% DPSIMUL  Monte Carlo simulation of discrete time controlled Markov process
% USAGE
%   [ssim,xsim] = dpsimul(model,ss,nper,sres,xres)
% INPUTS
%   model   : model structure variable
%   ss      : k by d vector of initial states
%   nper    : number of simulated time periods
%   sres    : coordinates of the evaluation grid (from dpsolve)
%   xres    : optimal control function values 
%               at grid defined by sres
% OUTPUTS
%   ssim    : k by d by nper+1 vector of simulated states       
%   xsim    : k by d by nper+1 vector of simulated actions
% For finite horizon problems, xsim is k by d by nper
% If d=1, ssim and xsim are k by nper+1
%
% USES: DISCRAND

% Copyright (c) 1997-2002, Paul L. Fackler & Mario J. Miranda
% paul_fackler@ncsu.edu, miranda.4@osu.edu

function [ssim,xsim] = mydpsimul(model,ss,nper,sres,xres,noshocks,impulse_shock)

if noshocks == true;fprintf('Simulating without shocks! (converge to stochastic steady state)\n'); end

% check that impulse (shock) has correct size
if nargin > 6
    if size(impulse_shock,2) ~= size(model.e,2)
        impulse_shock = impulse_shock';
    end
    assert(size(impulse_shock,2)== size(model.e,2),'Impulse shock does not have the right dimension')
else
    impulse_shock = [];
end
n_shocks = size(impulse_shock,1);


% DUMMY SHOCKS/WEIGHTS IF MODEL DETERMINISTIC  
  if isfield(model,'horizon'), T=model.horizon; else, T=inf; end;
  if isfield(model,'e'), e=model.e; else, e=0; end;
  if isfield(model,'w'), w=model.w; else, w=1; end;
 
  func=model.func;
  params=model.params;
  
  nper = min(nper,T);
  nrep = size(ss,1);
  ds   = size(ss,2);
  st   = gridmake(sres);
  if T==inf
    dx = ds*length(xres(:))/length(st(:));
  else
    dx = ds*length(xres(:))/(T*length(st(:)));
  end
  ssim = zeros(nrep,ds,nper+1);  
  xsim = zeros(nrep,dx,nper+1);  
 
    nx = numel(xres)/dx;
    xres = reshape(xres,nx,dx);
    xmean = zeros(1,size(xres,2));
    bar = waitbar(0,'Simulating...');
    for t=1:nper+1
      xx = minterp(sres,xres,ss);
      ssim(:,:,t) = ss;      
      xsim(:,:,t) = xx; 
      if noshocks == false
          ee = e(discrand(nrep,w),:);
      else
          ee = zeros(1,size(e,2));
      end
      if t <= n_shocks
          ee = impulse_shock(t,:);
      end
      ss = feval(func,'g',ss,xx,ee,params{:}); 
    waitbar(t / nper)
    end   
    close(bar)
    xsim=xsim(:,:,1:t); 


 if ds==1; ssim=squeeze(ssim); end 
 if dx==1; xsim=squeeze(xsim); end