% Residuals_PREC.m is a routine for evaluating the accuracy of the solutions
% to the neoclassical stochastic growth model; it computes  residuals in 
% the optimality conditions on a given set of points in the state space 
% by evaluating integrals exactly; see "How to Solve Dynamic Stochastic 
% Models Computing Expectations Just Once" by Kenneth L. Judd, Lilia Maliar, 
% Serguei Maliar and Inna Tsener, (2016), Quantitative Economics, forthcoming, 
% (henceforth, JMMT, 2016).
%
%-------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Inputs:    "k"     is current capital;
%            "z"     is current productivity
%            "A"     is the normalizing constant in production;
%            "alpha" is the share of capital in production;
%            "gam"   is the utility-function parameter;
%            "delta" is the depreciation rate;
%            "sigma" is the standard deviation of the productivity shock;
%            "rho"   is the autocorrelation coefficient in the process for
%                    shock;
%            "beta"  is the discount factor; 
%            "vf"    is the vector of polynomial coefficients in the 
%                    approximated value function;  
%            "D"     is the degree of polynomial approximation
%
% Output:    "Mean_Residuals" and "Max_Residuals" are, respectively, the 
%            mean and maximum absolute Euler equation residuals (in log10)
% -------------------------------------------------------------------------
% Copyright © 2016 by Lilia Maliar, Serguei Maliar and Inna Tsener. All 
% rights reserved. The code may be used, modified and redistributed under 
% the terms provided in the file "License_Agreement.txt".
% -------------------------------------------------------------------------

function [Mean_Residuals,Max_Residuals]  = Residuals_PREC(k,z,A,alpha,gam,delta,sigma,rho,beta,vf,D)

T     = size(k,1);   % Simulation length

% Precomputation of integrals in the accuracy evaluation procedure
%-----------------------------------------------------------------
vcv = sigma^2;                    % Variance covariance matrix
intV = Integrals_PREC_2d(D,vcv)'; % Precompute integrals for constructing 
                                  % expectation functions
                     
 for t=1:T;   
     
     k0 = k(t,1);                 % Capital of period t
     z0 = z(t,1);                 % Productivity level of period t 
             
     % Current period quantities in (k0,z0)
     % ------------------------------------     
     Vder0 = Polynomial_deriv_2d([k0 z0],D)*vf;  
                                  % Construct the derivative of value function 
     u0 = Vder0./(1-delta+A*alpha*z0*k0^(alpha-1));
                                  % Envelope condition
     c0 = u0.^(-1/gam);           % Consumption
     k1 = (1-delta)*k0+A*z0.*k0^alpha-c0;          
                                  % Next-period capital
     
     % Future period quantities in n_grid grid points (k1,z1)
     % -------------------------------------------------------------     
     vf1 = vf.*intV;              % A new vector of coefficients with 
                                  % precomputed integrals
     Vder1 = Polynomial_deriv_2d([k1 z0.^rho],D)*vf1;
                        % Derivative of value function evaluated in the 
                        % future state using the new vector of coefficients
     
     % Unit-free residuals in the Euler equation
     % -----------------------------------------
     Residuals(t,1) = (beta*Vder1)./c0.^(-gam)-1;
     
 end

 % Output of the accuracy evaluation test
 % --------------------------------------
 Mean_Residuals = log10(mean(mean(abs(Residuals(1:end))))); % Mean residuals 
 Max_Residuals = log10(max(max(abs(Residuals(1:end)))));    % Max  residuals
    
 
    