% conditional simulation 
% need to run dynare first
% Could use Hermitian quadrature (see simul_higher_order_dyn.m)

rng(161616)
repet=100;
leng=100;
burn=round(0.2*leng);
nseq=1;
pos_sigma=find_incell(M_.param_names,'\<sigma_a\>');
% note that innovation have been normalized to have unit variance
intensity=1;
seq=repmat(intensity,nseq,1);%*abs(randn(5,1));
Y=zeros(size(M_.endo_names,1),leng,repet);
shocks=zeros(leng,M_.exo_nbr,repet);
pos_sh=1;
tic
for jj=1:repet
        shocks(:,:,jj)=randn(leng,M_.exo_nbr);
        [r]=dynare_simul_matlab(oo_,M_,options_,shocks(:,:,jj));
        shocks(burn+1:burn+nseq,pos_sh,jj)=shocks(burn+1:burn+nseq,pos_sh,jj)+seq;
         [r_seq]=dynare_simul_matlab(oo_,M_,options_,shocks(:,:,jj));
        Y(:,:,jj)=r_seq-r;
end
toc
YY=mean(Y,3);
pos_y=find_incell(M_.endo_names,'\<Y\>');
var_inter=YY(pos_y,burn-1:end);
figure;plot(var_inter);
% compare to first order
options_.irf=10;
irf_dynare;
var_inter(3:13)-oo_.irfs.Y_e_a*intensity