% accuracy.m
%
% this program checks the loglinear solution for accuracy
%

%indicate which column of decision corresponds with which variable
%note that for each application you have to check that these column
%indications are correct. Here the labels are the same for each program
iiicap=1;
iiicons=3;

dynare modelcloglinear noclearall
% now we get a policy rule for the LOG of capital and the LOG of consumption
% note that productivity is always in logs
load dynarerocks
% this loads the matrix "decision" with coefficients into memory

%==========================================================================
%
% PRELIMINARIES - DEFINITION OF POLICY FUNCTION
%
% To make the program more readable I wrote a function polfunction(.) with
% which you can calculate the policy choices as a function of the inputs
% The format is as follows
%
% cap = polfunction(caplag,prodlag,shocks,iiicap)
% con = polfunction(caplag,prodlag,shocks,iiicons)
%
global decision capstst
%
% the policy function uses these two inputs 
%==========================================================================

%==========================================================================
%
%  DenHaan-Marcet Accuracy test
%
%==========================================================================

% step 1: Simulate the data

T=1000;
cap    = zeros(T,1);
con    = zeros(T,1);
prod   = zeros(T,1);	
randn('state',2008);
sig = 0.01; % this value of sigma has to equal the one given in the Dynare programs (for 2nd-order)
shocks = sig*randn(T,1); 
prod(1)= 0;

capstst = decision(1,iiicap); 	
cap(1)  = capstst; % initial value is steady state level
con(1)  = exp(alpha*capstst) - delta*exp(capstst);
con(1)  = log(con(1));
for t = 2:T	
prod(t)=rho*prod(t-1)+shocks(t);
cap(t) =polfunction(cap(t-1),prod(t-1),shocks(t),iiicap);
con(t) =polfunction(cap(t-1),prod(t-1),shocks(t),iiicons);
end	

% step 2: define error term

error = -exp(-nu*con(1:T-1)) + beta*exp(-nu*con(2:T)).*(alpha*exp(prod(2:T)+(alpha-1)*cap(2:T))+1-delta);
% this is just the realization of the Euler equation itself. You could
% multiply this with any term in the information set

% step 3: test whether the mean of the error term is zero

    variance = ((error-mean(error))'*(error-mean(error)))/size(error,1);
    dhmtest = size(error,1)*mean(error)^2/variance;

    fprintf(' \n')
    fprintf(' ================================= \n')
    fprintf(' \n')
    fprintf(' Beginning of Accuracy Tests \n')
    fprintf(' \n')
    fprintf(' ================================= \n')
    fprintf(' \n')    
    fprintf(' DHM statistic \n');
    fprintf(' \n')
    fprintf(' lower 5perc & upper 5perc & DHM statistic %2.5f %2.5f %2.5f \n',chi2inv(0.05,1),chi2inv(0.95,1),dhmtest);
    fprintf(' \n')


%==========================================================================
%
%  EULER EQUATION ERRORS WHEN INNOVATIONS HAVE DISCRETE SUPPORT
%
%==========================================================================

%==========================================================================
% simple example
%==========================================================================

%
% First, I show how to do the accuracy test when the initial conditions 
% correspond to steady state values
% Of course, a good accuracy test is done at many points in the state
% space, which is done below


capstst = decision(1,iiicap); 	
caplag  = capstst; 
prodlag = 0;
shocks  = +sig; % I arbitrarily pick one of the two realizations

% consumption as implied by the numerical solution
con = polfunction(caplag,prodlag,shocks,iiicons);

%
% now we are going to calculate the implied consumption value
% if the numerical solution is accurate then we should get a very similar
% value
%

prod =rho*prodlag + shocks;
cap = polfunction(caplag,prodlag,shocks,iiicap);
% these two values are independent of what will happen next period

% start calculating the conditional expectation

% calculate realization when next period's realization is high

shocksnext = sig;
prodnext = rho*prod+shocksnext;
connext = polfunction(cap,prod,shocksnext,iiicons);
realizationhigh = beta*exp(-nu*connext)*(alpha*exp(prodnext+(alpha-1)*cap)+1-delta);

% calculate realization when next period's realization is low

shocksnext = -sig;
prodnext = rho*prod+shocksnext;
connext = polfunction(cap,prod,shocksnext,iiicons);
realizationlow  = beta*exp(-nu*connext)*(alpha*exp(prodnext+(alpha-1)*cap)+1-delta);

% weigh two possibilities with probabilities

condexp = 0.5*realizationhigh + 0.5*realizationlow;

conimplied = condexp^(-1/nu); 
conimplied = log(conimplied); % so we are comparing log of consumpion

percerror = 100*(con - conimplied); % since we have logs perc. error is simply the difference

    fprintf(' Euler Errors when shocks have discrete support \n');
    fprintf(' \n')
    fprintf(' error at one particular point %2.5f \n', percerror);

%==========================================================================
% THE REAL ACCURACY TEST WHEN INNOVATIONS ARE DISCRETE
% that is, do the accuracy test at many points in the state space
%==========================================================================

%
% Now do the real accuracy test at many points in the state space
%

% construct grid

% number of grid points
numgridk = 101;
numgridp = 101;
numgrids =   2;
error = zeros(numgridk,numgridp,numgrids);

% bounds of the grid (reasonable bounds can e.g. be obtained by looking at
% a simulation
klow = 1.7;
khigh = 2.05;
plow = -0.1;
phigh = 0.1;

gridk = (klow: (khigh-klow)/(numgridk-1) : khigh)';
gridp = (plow: (phigh-plow)/(numgridp-1) : phigh)';
grids = [-sig;sig];

% The Dynare solution uses prod(t-1) and shocks(t) as separate inputs
% typically, however, prod(t) captures the information of both. You can 
% check this in our example by noting that
% decision(3,iii)*prod(t-1)+decision(4,iii)*shocks(t)
% =
% decision(4,iii)*prod(t)
%
% The reason for this equality is that (for our and many other models)
% decision(3,iii)=rho*decision(4,iii) 
%
% The consequence of this is that the accuracy test could have been done
% using only lagged capital and the current value of productivity but to
% stick close to the Dynare specification I keep lagged productivity and
% the innovation as separate.

for ik = 1:numgridk
    for ip = 1:numgridp
        for is = 1:numgrids

        % set initial values
        caplag  = gridk(ik); 
        prodlag = gridp(ip);
        shocks  = grids(is);
     
        % now do the exact same thing as was done above

        con = polfunction(caplag,prodlag,shocks,iiicons);

        prod =rho*prodlag + shocks;
        cap = polfunction(caplag,prodlag,shocks,iiicap);

        shocksnext = sig;
        prodnext = rho*prod+shocksnext;
        connext = polfunction(cap,prod,shocksnext,iiicons);
        realizationhigh = beta*exp(-nu*connext)*(alpha*exp(prodnext+(alpha-1)*cap)+1-delta);
        shocksnext = -sig;
        prodnext = rho*prod+shocksnext;
        connext = polfunction(cap,prod,shocksnext,iiicons);
        realizationlow  = beta*exp(-nu*connext)*(alpha*exp(prodnext+(alpha-1)*cap)+1-delta);
               
        condexp = 0.5*realizationhigh + 0.5*realizationlow;
        conimplied = condexp^(-1/nu); 
        conimplied = log(conimplied); % so we are comparing log of consumpion

        error(ik,ip,is) = 100*abs((con - conimplied)); % since we have logs perc. error is simply the difference

        end
    end
end

    fprintf(' mean of perc. errors %2.5f \n', mean(mean(mean(error))));
    fprintf(' max  of perc. errors %2.5f \n', max(max(max(error))));

%==========================================================================
%
%  EULER EQUATION ERRORS WHEN SHOCKS ARE NORMALLY DISTRIBUTED
%
%==========================================================================

numnodes = 10; % number of Hermite nodes used in integration
[Hnodes , Hweights] = hernodes(numnodes);


% Again, I first show how to do the accuracy test when the initial conditions 
% correspond to steady state values

caplag  = capstst; 
prodlag = 0;
shocks  = 0;

% consumption as implied by the numerical solution
con = polfunction(caplag,prodlag,shocks,iiicons);

%
% now we are going to calculate the implied consumption value
% if the numerical solution is accurate then we should get a very similar
% value
%

prod =rho*prodlag + shocks;
cap = polfunction(caplag,prodlag,shocks,iiicap);
% these two values are independent of what will happen next period

% start calculating the conditional expectation

% calculate for each hermite node the term inside the integrand (taking
% into account the change in variables to get correct expression for
% Gaussian Hermite quadrature

condexp = 0; %initialize value
for i = 1:numnodes
    shocksnext  = sqrt(2)*sig*Hnodes(i);
    prodnext    = rho*prod+shocksnext;
    connext     = polfunction(cap,prod,shocksnext,iiicons);
    realization = beta*exp(-nu*connext)*(alpha*exp(prodnext+(alpha-1)*cap)+1-delta);
    condexp     = condexp + Hweights(i)*realization/sqrt(pi);
end

% the 1/sqrt(pi) term above is what remains of the constant of the normal density
% after the multiplication with the Jacobian due to the change in variables

conimplied = condexp^(-1/nu); 
conimplied = log(conimplied); % so we are comparing log of consumpion
percerror = 100*(con - conimplied); % since we have logs perc. error is simply the difference

    fprintf(' \n')
    fprintf(' \n')
    fprintf(' Euler Errors when shocks are Normally distributed \n');
    fprintf(' \n')
    fprintf(' error at one particular point %2.5f \n', percerror);

%==========================================================================
% THE REAL ACCURACY TEST WHEN INNOVATIONS ARE NORMAL
% that is, do the accuracy test at many points in the state space
%==========================================================================

% construct grid

% number of grid points
numgridk =  51;
numgridp =  51;
numgrids =  51;
error = zeros(numgridk,numgridp,numgrids);

% bounds of the grid (reasonable bounds can e.g. be obtained by looking at
% a simulation
klow = 1.7;
khigh = 2.05;
plow = -0.1;
phigh = 0.1;
slow  = -3*sig;
shigh =  3*sig;

gridk = (klow: (khigh-klow)/(numgridk-1) : khigh)';
gridp = (plow: (phigh-plow)/(numgridp-1) : phigh)';
grids = (slow: (shigh-slow)/(numgrids-1) : shigh)';


% The Dynare solution uses prod(t-1) and shocks(t) as separate inputs
% typically, however, prod(t) captures the information of both. You can 
% check this in our example by noting that
% decision(3,iii)*prod(t-1)+decision(4,iii)*shocks(t)
% =
% decision(4,iii)*prod(t)
%
% The reason for this equality is that (for our and many other models)
% decision(3,iii)=rho*decision(4,iii) 
%
% The consequence of this is that the accuracy test could have been done
% using only lagged capital and the current value of productivity but to
% stick close to the Dynare specification I keep lagged productivity and
% the innovation as separate.

for ik = 1:numgridk
    for ip = 1:numgridp
        for is = 1:numgrids

        % set initial values
        caplag  = gridk(ik); 
        prodlag = gridp(ip);
        shocks  = grids(is);
     
        % now do the exact same thing as was done above

        con = polfunction(caplag,prodlag,shocks,iiicons);

        prod =rho*prodlag + shocks;
        cap = polfunction(caplag,prodlag,shocks,iiicap);

        condexp = 0; 
        for i = 1:numnodes
            shocksnext  = sqrt(2)*sig*Hnodes(i);
            prodnext    = rho*prod+shocksnext;
            connext     = polfunction(cap,prod,shocksnext,iiicons);
            realization = beta*exp(-nu*connext)*(alpha*exp(prodnext+(alpha-1)*cap)+1-delta);
            condexp     = condexp + Hweights(i)*realization/sqrt(pi);
        end

        conimplied = condexp^(-1/nu); 
        conimplied = log(conimplied); % so we are comparing log of consumpion
        error(ik,ip,is) = 100*abs((con - conimplied)); % since we have logs perc. error is simply the difference
 
    
        end
    end
end

    fprintf(' mean of perc. errors %2.5f \n', mean(mean(mean(error))));
    fprintf(' max  of perc. errors %2.5f \n', max(max(max(error))));
