% MC accuracy test based on dynare_symul
% @ Feb 2015 Gianni.lombardo giannilmbd@gmail.com
% This file used dynare++ to evaluate Euler equation residuals using MC
% integration.
% Contrary to the built in check in dynare++ here can control the
% simulation:e.g. using the same draws used in the generalized IRFs (thus
% get the error specific to that simulation). Here I'm using a truncated
% multivariate normal distribution to avoid extreme events (and
% explosiveness of series: rmvnrnd Copyright 2011-2013 Tim J. Benham,
% School of Mathematics and Physics,)
% To compute expectations I use Gaussian-quadrature (Function from CompEcon
% -- Miranda Fackler). This can create an inconsistency with the truncated
% multivariate normal used for the simulation.

% THERE ARE ONLY FEW PLACES TO FIX: E.G. FURTHER BELOW WHERE HAVE TO GIVE
% THE MODEL SPECIFIC EULER EQUATION.

%% Run dynare ++ NEED TO ADAPT TO YOUR MODEL (E.G. FILE NAME)
clear all
close all
fname='mickgiapp_QL';
param_mickgia;
% get number of processors

Nproc=java.lang.Runtime.getRuntime.availableProcessors;
tic
if 0 % zero if don't want to re-run
    
    % SET THE ORDER OF APPROXIMATION
    order=5;
    
    eval(['!dynare++ --seed 10001 --ss-tol 1.e-9 --per 0 --sim 0 --order ', num2str(order),' --no-irfs --threads ', num2str(Nproc),' --check pse --check-scale .01 --check-num 000 --check-evals 0000 ', fname,'.mod'])
end

disp('********************************************')
disp('time for dynare++ ')
disp('********************************************')
toc
eval(['load ',fname,'.mat'])
%% DRAW SHOCKS in advance (using truncated multivariate normal)
exo_nbr=size(dyn_vcov_exo,1);
scale=0.01;
% length of simulation
leng=60;
% REPETITIONS for simulations
repet=5;
% nodes for integration
repetN=100;
repetN=round(repetN^(1/3));
simsim=(repetN);
repetN=repetN^3;
[Hnodes,Hweights]= qnwnorm([simsim,simsim,simsim],[0,0,0],dyn_vcov_exo*scale^2);




shocks=zeros(leng+1,exo_nbr ,repet);

% See rmvnrnd for details on imput parameters: note that I'm using a
% symmetric square bound, where the lower bound can be given either as
% a negative identity matrix and positive parameter or vice-versa.
sqlimL=-.1;
sqlimU=.1;
%  Creat all shock sequences
parfor jj=1:repet
    
    %  shocks_tmp=mvnrnd(zeros(exo_nbr,1),dyn_vcov_exo,leng)*diag([0;0;1])*scale;
    shocks_tmp=rmvnrnd(zeros(exo_nbr,1),dyn_vcov_exo,leng+1,[-eye(exo_nbr); eye(exo_nbr)],[ones(exo_nbr,1)*sqlimU;ones(exo_nbr,1)*sqlimU]);%[rand(1000,exo_nbr)-0.5]*dyn_vcov_exo*1;
    % DEMEAN THE SERIES
    shocks(:,:,jj)=shocks_tmp-repmat(mean(shocks_tmp),size(shocks_tmp,1),1);
    
    %
    %     [Hnodes(:,:,jj),Hweights(:,jj)]= qnwnorm([simsim,simsim,simsim],[0,0,0],dyn_vcov_exo*scale);
    %     shocks(:,:,jj)=Hnodes(:,:,jj);
    
end

%% compute series
% initialize variables and gather information
zeroshock=zeros(exo_nbr,1);
varnames=dyn_vars;
Y=zeros(size(varnames,1),leng+1,repet);
Ye=zeros(leng+1,repet); % Matrix of expeted quantities
Yetmp=zeros(size(varnames,1),repetN);
%% timing assessment
tic;
Y(:,1,1)=dynare_simul(fname, shocks(1,:,1)');
g=toc;

tottime_1=g*(repet)/(60*60)/Nproc;
tottime_2=g*(repet*repetN/Nproc*leng)/(60*60);

fprintf('Approximate time to complete: %2.4f  hours \n', tottime_2+tottime_1)


%% Simulation
disp('PRESS ANY KEY TO CONTINUE')
pause
disp('Computing errors')
tic
% compute first period for all replications (can also be used to time replications)
% tic
parfor jj=1:repet
    Y(:,1,jj)=dynare_simul(fname, shocks(1,:,jj)');
end
h = waitbar(0,'Please wait...');
zz=1;
for kk=2:leng+1
    for uu=1:repet
        Y(:,kk,uu)=dynare_simul(fname, shocks(kk,:,uu)',Y(:,kk-1,uu));
        % integrate one period ahead Euler factor
        parfor jj=1:repetN
            Yetmp(:,jj)=dynare_simul(fname, Hnodes(jj,:)',Y(:,kk-1,uu));
            
        end
        
        Ye(kk,uu)=(((betta*(exp(Yetmp(dyn_i_C,:))).^(-siggma).*exp(Yetmp(dyn_i_R,:))./exp(Yetmp(dyn_i_pii,:))).^(-1/siggma)))*Hweights;
        
        res(:,kk,uu)=(Ye(kk,uu)/exp(Y(dyn_i_C,kk,uu)))-1;
         waitbar(zz/ (repet*leng))
zz=zz+1;
        
    end
end
close(h)

EER=log10(mean(mean(abs(res),3)));
MER=log10(max(max(abs(res),[],3)));
d=toc;
fprintf('Total time elapsed %3.5f hours\n', d/(60*60))
fprintf('Maximum absolute error our of %4i draws: %2.5f\n',repet*leng,MER)
fprintf('Average absolute error our of %4i draws: %2.5f\n',repet*leng,EER)