% Euler equation residuals;
% and simulation of NP and Pruned solution
% take block from simult_.m
% to be run after dynare with pruning

accuracy=0; % if 0 simulate only

Equation='1-(beta*expected_part*exp(RR))/(exp(CC)^(-sigma))';
expected_part='exp(CCE).^(-sigma).*exp(one_def_)';

% 
% pos_cons=find_incell(M_.endo_names,'\<C\>');
% pos_R=find_incell(M_.endo_names,'\<R\>');
pos_one_def=find_incell(M_.endo_names,'\<y\>');
pos_sigma_a=find_incell(M_.param_names,'\<std_a\>');


iter=100;
sim=1;
scale=M_.params(pos_sigma_a);
dr=oo_.dr;
endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
ex_=zeros(iter,exo_nbr);randn(iter,exo_nbr)*diag([scale;0]);
ex_(1,1)=1;
% load ex_.mat
ex_=ex_(1:iter,:);
% ex_(:,2)=0;% hit first shock
% exe_=scale*randn(sim,exo_nbr);
% exe_(:,2)=0;% hit first shock
%  k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
k2 = dr.kstate(dr.kstate(:,2) <= M_.maximum_lag+1,[1 2]);
k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr;
order_var = dr.order_var;
% only with pruning
%         keyboard
y0=dr.ys;
ghx = dr.ghx;
ghu = dr.ghu;
ghxx = dr.ghxx;
ghxu = dr.ghxu;
ghuu = dr.ghuu;
ghs2 = dr.ghs2;
ghxxx = dr.ghxxx;
ghxxu = dr.ghxxu;
ghxuu = dr.ghxuu;
ghuuu = dr.ghuuu;
ghxss = dr.ghxss;
ghuss = dr.ghuss;
threads = options_.threads.kronecker.A_times_B_kronecker_C;
nspred = M_.nspred;
ipred = M_.nstatic+(1:nspred);
yhat1 = y0(order_var(k2))-dr.ys(order_var(k2));
yhat2 = zeros(nspred,1);
yhat3 = zeros(nspred,1);
yhat3NP = zeros(nspred,1);
y_=zeros(endo_nbr,iter);
yNP_=zeros(endo_nbr,iter);
ye_=zeros(endo_nbr,sim);
yeNP_=zeros(endo_nbr,sim);
yee_=zeros(endo_nbr,iter);
yexpected=zeros(1,iter);
yexpectedNP=zeros(1,iter);
% taken from Miranda-Fackler

if accuracy
    simsim=sqrt(sim);
[Hnodes,Hweights]= qnwnorm([simsim,simsim],[0,0],[scale,0;0,scale].^2);
shzer=(sum(ex_,1)~=0);
Hnodes=Hnodes*diag(shzer);
%zero shocks that are zero in simulation

Hweights=[0;Hweights];
end
%
tic
for i=2:iter+M_.maximum_lag
    %simulation block
    u = ex_(i-1,:)';
    [y_(:,i),yhat1,yhat2,yhat3]=evaluation_block_3ord(yhat1,yhat2,yhat3,u,i,order_var,ipred,ghx,ghu,ghxx,ghuu,ghxu,ghs2,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,threads);
    [yNP_(:,i),yhat3NP]=evaluation_block_3ord_NP(yhat3NP,u,i,order_var,ipred,ghx,ghu,ghxx,ghuu,ghxu,ghs2,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,threads);
    %expectation block
    yhat1e=yhat1;
    yhat2e=yhat2;
    yhat3e=yhat3;
    yhat3eNP=yhat3NP;
    if accuracy
        for j=2:sim+M_.maximum_lag
            ue=Hnodes(j-1,:)';
            [ye_(:,j),yhat1e,yhat2e,yhat3e]=evaluation_block_3ord(yhat1e,yhat2e,yhat3e,ue,j,order_var,ipred,ghx,ghu,ghxx,ghuu,ghxu,ghs2,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,threads);
            [yeNP_(:,j),yhat3eNP]=evaluation_block_3ord_NP(yhat3eNP,ue,i,order_var,ipred,ghx,ghu,ghxx,ghuu,ghxu,ghs2,ghxxx,ghuuu,ghxxu,ghxuu,ghxss,ghuss,threads);
            
        end
        CCE=ye_(pos_cons,:);
        one_def_=ye_(pos_one_def,:);
        %             yexpected(1,i)=mean(eval(expected_part)); % use if Montecarlo
        %             integration; if Hermitian integration use sum
        yexpected(1,i)=(eval(expected_part))*Hweights;
        CCE=yeNP_(pos_cons,:);
        one_def_=yeNP_(pos_one_def,:);
        %             yexpectedNP(1,i)=mean(eval(expected_part));
        yexpectedNP(1,i)=(eval(expected_part))*Hweights;
    end
end

% evaluation of Euler equation
if accuracy
    res=zeros(length(y_),1);
    resNP=zeros(length(y_),1);
    for jj=1:length(y_)
        CC=y_(pos_cons,jj);
        RR=y_(pos_R,jj);
        expected_part=yexpected(1,jj);
        
        res(jj,1)=eval(Equation);
        
        CC=yNP_(pos_cons,jj);
        RR=yNP_(pos_R,jj);
        expected_part=yexpectedNP(1,jj);
        
        resNP(jj,1)=eval(Equation);
        
    end
    
    
    figure;plot(log10(res(2:end)))
    title('Euler Residual (log10): Pruning');
    
    figure;plot(log10(resNP(2:end)))
    title('Euler Residual (log10): NO Pruning');
end
toc
figure;hist(exp(y_(pos_one_def,:)),10)
title('distribution of 1-defaults (level): SE')
figure;hist(exp(yNP_(pos_one_def,:)),10)
title('distribution of 1-defaults (level): NP')
