recsSimul Simulates a model from starting values given in s0 and for nper period
 
  SSIM = recsSimul(MODEL,INTERP,S0,NPER) simulates the model defined in the
  object MODEL, by using the interpolation structure defined in the structure
  INTERP. The simulation starts from the initial state S0 and lasts NPER (scalar)
  periods. S0 is a nrep-by-d matrix with nrep the number of scenarios to simulate,
  and d the number of state variables. recsSimul returns the nrep-by-d-by-nper
  array SSIM that contains the simulated state variables.
  MODEL is an object created by recsmodel.
  INTERP is a structure, which includes the following fields:
     ch      : coefficient matrix of the interpolation of the expectations
               function (optional, to be provided with method 'expfunapprox')
     cx      : coefficient matrix of the interpolation of the response variables
     cz      : coefficient matrix of the interpolation of the expectations variables
     fspace  : a definition structure for the interpolation family (created by
               the function fundef)
 
  SSIM = recsSimul(MODEL,INTERP,S0,NPER,SHOCKS) uses the nrep-by-q-by-(nper-1) array
  SHOCKS to simulate the model instead of drawing random numbers. In this case
  size(SHOCKS,3) supersedes NPER.
 
  SSIM = recsSimul(MODEL,INTERP,S0,NPER,SHOCKS,OPTIONS) simulates the model with
  the parameters defined by the structure OPTIONS. The fields of the structure are
     accuracy         : 1 to check accuracy on the asymptotic distribution (default: 0)
     display          : 1 (default) to display outputs
     eqsolver         : 'fsolve', 'lmmcp' (default), 'ncpsolve' or 'path'
     eqsolveroptions  : options structure to be passed to eqsolver
     extrapolate      : 1 if extrapolation is allowed outside the
                        interpolation space or 0 to forbid it (default: 1)
     functional       : 1 if the equilibrium equations are a functional equation
                        problem (default: 0)
     loop_over_s      : 0 (default) to solve all grid points at once, 1 to loop
                        over each grid points, or n to loop over n blocks of
                        grid points
     funapprox        : 'expapprox', 'expfunapprox', or 'resapprox' (default)
     simulmethod      : 'interpolation' (default) or 'solve'
     stat             : 1 to ouput summary statistics from the simulation
                        (default: 0)
     UseParallel      : 'always' (default) to use parallel calculation (require
                        Parallel Computing Toolbox)' or never'
 
  [SSIM,XSIM] = recsSimul(MODEL,INTERP,S0,NPER,...) returns the nrep-by-m-by-nper
  array XSIM that contains the simulated response variables.
 
  [SSIM,XSIM,ESIM] = recsSimul(MODEL,INTERP,S0,NPER,...) returns the
  nrep-by-q-by-nper array ESIM that contains the shocks.
 
  [SSIM,XSIM,ESIM,STAT] = recsSimul(MODEL,INTERP,S0,NPER,...) returns summary
  statistics as a structure that contains the number of observations (STAT.N),
  the moments (STAT.MOMENTS), the correlation between variables (STAT.COR), and
  the autocorrelation (STAT.ACOR). Asking recsSimul to return STAT forces the
  OPTIONS.STAT to 1.
 
  [SSIM,XSIM,ESIM,STAT,FSIM] = recsSimul(MODEL,INTERP,S0,NPER,...)  returns the
  nrep-by-m-by-nper array FSIM that contains the value of the equilibrium
  equations on the simulation.