% Accuracy_Test_N.m is a routine for evaluating accuracy of solutions to the    
% multi-country model: it computes approximation errors in the optimality  
% conditions on a given set of points in the state space; see "Numerically 
% Stable and Accurate Stochastic Simulation Approaches for Solving Dynamic 
% Economic Models" by Kenneth L. Judd, Lilia Maliar and Serguei Maliar, 
% (2011), Quantitative Economics 2/2, 173–210 (henceforth, JMM, 2011)., and 
% Juillard and Villemot, (2011), "Multi-Country Real Business Cycle Models: 
% Accuracy Tests and Test Bench", Journal of Economic Dynamics and Control 
% 35, 178-185.
%
% This version: July 14, 2011. First version: August 27, 2009.
% -------------------------------------------------------------------------
% Inputs:    "k" and "a" are, respectively, current-period capital and 
%            productivity levels, in the given set of points on which the 
%            accuracy is tested; 
%            "bk" are the coefficients of the capital policy functions of N 
%            countries;
%            "IM" is the integration method for evaluating accuracy, 
%            IM=1,2,..,10=Gauss-Hermite quadrature rules with 1,2,...,10 
%            nodes in each dimension, respectively, 
%            IM=11=Monomial rule with 2N nodes,
%            IM=12=Monomial rule with 2N^2+1 nodes;
%            "alpha", "gam", "phi", "beta", "A", "tau", "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 approximation errors across all optimality conditions;
%            "Errors_max_EE", "Errors_max_MUC", "Errors_max_MUL", and 
%            "Errors_max_RC" are the maximum approximation errors disaggre- 
%            gated by optimality conditions 
% -------------------------------------------------------------------------
% Copyright © 2011 by Lilia Maliar and Serguei Maliar. 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] = Accuracy_Test_N(k,a,bk,D,IM,alpha,gam,delta,beta,A,tau,rho,vcv,discard)

tic                   % Start counting time needed to run the test        

[P,N] = size(a);      % Infer the number of points P on which the accuracy 
                      % is evaluated and the number of countries N
                      

% 2. Integration method for evaluating accuracy 
% ---------------------------------------------
if (IM>=1)&&(IM<=10)
    [n_nodes,epsi_nodes,weight_nodes] = GH_Quadrature(IM,N,vcv);
                             % Compute the number of integration nodes, 
                             % n_nodes, integration nodes, epsi_nodes, and 
                             % integration weights, weight_nodes, for Gauss-
                             % Hermite quadrature integration rule with IM 
                             % nodes in each dimension
elseif IM == 11
    [n_nodes,epsi_nodes,weight_nodes] = Monomials_1(N,vcv);
                             % Monomial integration rule with 2N nodes
elseif IM == 12
    [n_nodes,epsi_nodes,weight_nodes] = Monomials_2(N,vcv);
                             % Monomial integration rule with 2N^2+1 nodes
end    

% 3. Polynomial bases for the test
%---------------------------------
X = Ord_Polynomial_N([k a],D); % Form a complete ordinary polynomial of degree
                             % D on the given set of points 

% 4. Given the solution for capital, compute consumption on the given set 
% of points  
%------------------------------------------------------------------------
for p = 1:P;                 % For each given point, ...     
    
       p                     % Display the point (with the purpose of  
                             % monitoring the progress) 
                                     
        % 4.1 Variables in point p 
        % ------------------------        
        k0 = k(p,1:N);       % N capital stocks of period t
        a0 = a(p,1:N);       % N productivity levels of period t
        X0 = X(p,:);         % Complete (second-degree) polynomial 
                             % bases at t
                                    
        % 4.2 Capital and consumption choices at t
        % ----------------------------------------
        k1 = X0*bk; 
        % Compute a row-vector of capital of period t+1 (chosen at t) using
        % the corresponding capital policy functions; 1-by-N
        
        C0 = (A*k0.^alpha.*a0 - 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; 1-by-1

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

        % 4.3 Capital and consumption choices at t+1
        %-------------------------------------------
        a1 = (ones(n_nodes,1)*a0).^rho.*exp(epsi_nodes);    
        % Compute the next-period productivity levels in each integration node
        % using condition (?) in the online appendix; n_nodes-by-N

        k1_dupl = ones(n_nodes,1)*k1; 
        % Duplicate k1 n_nodes times to create a matrix with n_nodes identical
        % rows; n_nodes-by-N 

        X1 = Ord_Polynomial_N([k1_dupl a1],D);
        % Form a complete polynomial of degree D (at t+1) in the given point 
        
        k2 = X1*bk; 
        % Compute capital of period t+2 (chosen at t+1) using the second-
        % degree capital policy functions; n_nodes-by-N 

        C1 = (A*k1_dupl.^alpha.*a1 - k2+k1_dupl*(1-delta))*ones(N,1);
        % Aggregate consumption is computed by summing up individual consumption, 
        % which in turn, is found from the individual budget constraints; 
        % n_nodes-by-1

        c1 = C1*ones(1,N)/N;    % Individual consumption is the same for 
                                % all countries; n_nodes-by-N

% 5. Approximation errors in point p
%-----------------------------------
        
        % 5.1 Lagrange multiplier associated with the aggregate resource
        % constraint
        %---------------------------------------------------------------
        for j = 1:N
            MUC0j(1,j) = tau*c0(1,j).^(-gam); 
            % Compute a country's marginal utility of consumption multiplied 
            % by its welfare weight
        end
        lambda0 = mean(MUC0j,2);
        % An optimality condition w.r.t. consumption of period t equates 
        % the Lagrange multiplier of the aggregate resource constraint of 
        % period t and each country's marginal utility of consumption 
        % multiplied by its welfare weight; to infer the Lagrange multiplier,  
        % we average across N countries; 1-by-1
        
        for j = 1:N
            MUC1j(1:n_nodes,j) = tau*c1(1:n_nodes,j).^(-gam);
            % Compute a country's marginal utility of consumption multiplied 
            % by its welfare weight
        end
        lambda1 = mean(MUC1j,2);
        % Similarly, the Lagrange multiplier of the aggregate resource 
        % constraint of period t+1 is equal to a country's marginal utility 
        % of consumption multiplied by its welfare weight; to infer the 
        % Lagrange multiplier, we average across N countries; 1-by-n_nodes
        
        % 5.2 Unit-free Euler-equation errors
        %------------------------------------
        for j = 1:N
            Errors(p,j) = 1-weight_nodes'*(beta*lambda1/lambda0.*(1-delta+alpha*A*k1(1,j)^(alpha-1)*a1(1:n_nodes,j)));
        % A unit-free Euler-equation approximation error of country j
        end
        
        % 5.2 Unit-free errors in the optimality conditions w.r.t. consumption
        %---------------------------------------------------------------------
        for j = 1:N;
            Errors(p,N+j) = 1-lambda0./(tau*c0(1,j)^(-gam)); 
        % A unit-free approximation error in the optimality condition w.r.t. 
        % consumption of country j (this condition equates marginal utility 
        % of consumption, multiplied by the welfare weight, and the 
        % Lagrange multiplier of the aggregate resource constraint)
        end
        
        % 5.3 Unit-free errors in the optimality conditions w.r.t. labor 
        %---------------------------------------------------------------
        Errors(p,2*N+1:3*N) = zeros(N,1);
        % These errors  are zero by construction 
        
        % 5.4 Unit-free approximation error in the aggregate resource constraint
        %-----------------------------------------------------------------------
        Errors(p,3*N+1) = 1-(c0(1,1:N) + k1(1,1:N)-k0(1,1:N)*(1-delta))*ones(N,1)/((A*k0(1,1:N).^alpha.*a0(1,1:N))*ones(N,1));
        % This error is a unit-free expression of the resource constraint  
        % (?) in the online appendix
        
        % 5.5 Approximation errors in the capital-accumulation equation
        %--------------------------------------------------------------
        Errors(p,3*N+2:4*N+1) = zeros(N,1);
        % These errors are always zero by construction
        
        % For this model, GSSA produces zero errors (with machine accuracy)
        % in all the optimality conditions except of the Euler equation. 
        % Here, the errors in all the optimality conditions are introduced 
        % to make our accuracy measures comparable to those in the February 
        % 2011 special issue of the Journal of Economic Dynamics and Control 
end

% 6. Mean and maximum approximation errors computed after discarding the  
% first "discard" observations
%-----------------------------------------------------------------------

        % 6.1 Approximation errors across all the optimality conditions
        %--------------------------------------------------------------
        Errors_mean = log10(mean(mean(abs(Errors(1+discard:end,:))))); 
        % Average absolute approximation error 
        
        Errors_max = log10(max(max(abs(Errors(1+discard:end,:)))));    
        % Maximum absolute approximation error
 
        % 6.2 Maximum approximation errors disaggregated by the optimality 
        % conditions
        %-----------------------------------------------------------------
        Errors_max_EE = log10(max(max(abs(Errors(1+discard:end,1:N)))));    
        % Across N Euler equations

        Errors_max_MUC = log10(max(max(abs(Errors(1+discard:end,N+1:2*N)))));    
        % Across N optimality conditions w.r.t. consumption (conditions on   
        % marginal utility of consumption, MUC)

        Errors_max_MUL = log10(max(max(abs(Errors(1+discard:end,2*N+1:3*N)))));    
        % Across N optimality conditions w.r.t. labor (conditions on marginal 
        % utility of labor, MUL)

        Errors_max_RC = log10(max(max(abs(Errors(1+discard:end,3*N+1)))));    
        % In the aggregate resource constraint 
 
% 7. Time needed to run the test
%-------------------------------
time_test = toc;         