% Residuals_PREC_N.m is a routine for evaluating the accuracy of solutions 
% to the multi-country 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" and "z" are, respectively, current-period capital and 
%                 productivity levels, in the given set of points on 
%                 which the accuracy is tested; 
%            "bk" is the coefficients of the capital policy functions of N 
%                 countries;
%            "D"  is the degree of polynomial approximation;
%            "N"  is the number of countries;
%            "alpha", "gam","delta", "beta", "A", "rho" and "vcv" are the 
%                 parameters of the model;
%            "discard" is the number of data points to discard 
%
% Outputs:   "Errors_mean" and "Errors_max" 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 [Errors_mean, Errors_max, time_test] = Residuals_PREC_N(k,z,bk,N,D,alpha,gam,delta,beta,A,rho,vcv,discard)

tic                   % Start counting time needed to run the test        

[P,N1] = size(k);     % Infer the number of points P on which the accuracy 
                      % is evaluated and the number of endogenous state
                      % variables N1

N2 = size(z,2);       % Infer the number of exogenous state variables N2
   
integrals = Integrals_PREC_N(N1,N2,vcv,D)';
                      % Compute the values of the integrals of future 
                      % shocks defined by equation (50) for a
                      % country j
                             
% Note that "Integrals_PREC_N.m" is a function that must be used together 
% with "Ord_Polynomial_N.m". The input of this polynomial routine must 
% be in the form s = [k, z], where k and z are n1 and n2-dimensional vectors, 
% respectively. This ordering of state variables ensures that precomputed  
% integrals match the corresponding polynomial terms in the approximating 
% function of degree D

bq_mult = integrals*ones(1,N);
                    % Matrix of the precomputed integrals for the polynomial
                    % solutions 
                          
% 1. Polynomial bases for the test
%---------------------------------
X = Ord_Polynomial_N([k z],D); % Form a complete ordinary polynomial of degree
                               % D on the given set of points 

% 2. Given the solution for capital, compute consumption on the given set 
% of points  
%------------------------------------------------------------------------
        
        % 2.1 Compute consumption series of all countries 
        % -----------------------------------------------        
        k0 = k(1:P,:);       % N capital stocks of period t
        z0 = z(1:P,:);       % N productivity levels of period t
        X0 = X(1:P,:);       % The basis functions of a polynomial of 
                             % degree D
                                    
        % 2.2 Capital and consumption choices at t
        % ----------------------------------------
        k1 = X0*bk; 
        % Compute capital of period t+1 (chosen at t) using the corresponding 
        % capital policy functions; P-by-N
        
        C0 = (A*k0.^alpha.*z0 - k1+k0*(1-delta))*ones(N,1);
        % C is computed by summing up individual consumption, which in turn, is 
        % found from the individual budget constraints; P-by-1

        c0 = C0*ones(1,N)/N;  % Individual consumption is the same for all
                              % countries; P-by-N 

        % 2.3 Approximate N integrand functions
        %--------------------------------------
        
        Q0 = c0.^(-gam).*(1-delta+A*alpha*z0.*k0.^(alpha-1));                         
                                    % Construct derivative of value function
        bq = X0\Q0;                 % Compute the coefficients
        
        % 2.4 Approximate expectations of N integrand functions 
        %------------------------------------------------------
        
        bq1 = bq.*bq_mult;          % A new vector of coefficients with 
                                    % precomputed integrals
        EQ1 = Ord_Polynomial_N([k1 z0.^rho],D)*bq1;
                                    % Approximate N expectations

        % 2.5 Mean and maximum absolute Euler equation residuals (in log10); 
        % errrors are computed after discarding the first "discard" observations
        % -----------------------------------------------------------------
        
        Errors(1:P,1:N) = beta*EQ1./c0.^(-gam)-1;
                                    % Compute Euler equation errors; P-by-N
        Errors_mean = log10(mean(mean(abs(Errors(1+discard:end,1:N))))); 
                                    % Average absolute approximation error 
        Errors_max = log10(max(max(abs(Errors(1+discard:end,1:N))))); 
                                    % Maximum absolute approximation error
  
% 3. Time needed to run the test
%-------------------------------
time_test = toc;  