% This MATLAB software compares the computational cost of approximating 
% an expectation function of polynomial of the given degree by different 
% numerical integration methods, namely, Monte Carlo, Gauss-Hermite quadrature 
% rule and nonproduct monomial rules M1 and M2; 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).
%
% First version: July 3, 2013
% This version:  October 26, 2016

% ------------------------------------------------------------------------
% This software uses the following files: 
% ------------------------------------------------------------------------
% 1. "Cost_N.m"           compares costs of evaluating the expectation 
%                         function of polynomial of the given degree by 
%                         different numerical integration methods 
% 2. "Integrals_PREC_N.m" constructs the sets of precomputed integrals 
%                         defined by equation (53) in JMMT (2016) for ordinary 
%                         polynomials of the degrees from one to five, for  
%                         the N-country model
% 3. "Ord_Polynomial_N.m" constructs the sets of basis functions for ordinary 
%                         polynomials of the degrees from one to five, for
%                         the N-country model; borrowed from JMM (2011) 
% 4. "Monomials_1.m"      constructs integration nodes and weights for an N-
%                         dimensional monomial (non-product) integration rule 
%                         with 2N nodes; borrowed from Judd, Maliar and Maliar 
%                         (QE, 2011) 
% 5. "Monomials_2.m"      constructs integration nodes and weights for an N-
%                         dimensional monomial (non-product) integration rule 
%                         with 2N^2+1 nodes; borrowed from Judd, Maliar and
%                         Maliar (QE, 2011) 
% 6. "GH_Quadrature.m"    constructs integration nodes and weights for the 
%                         Gauss-Hermite rules with the number of nodes in 
%                         each dimension ranging from one to ten; borrowed 
%                         from Judd, Maliar and Maliar (QE, 2011)   
% 7."epsi30.mat"          contains a fixed 10,000-by-30 array of random 
%                         numbers drawn from a multivariate standard normal 
%                         distribution 
% -------------------------------------------------------------------------
% 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".
% -------------------------------------------------------------------------

clear all
clc
    
% 1. Parameters
% -------------
N = 1;                  % Number of stochastic shocks; note that the code 
                        % also works for N=1
sigma = 0.2;            % Standard deviation of shock to the log of the 
                        % productivity level
vcv = sigma^2*(eye(N)); % Variance-covariance matrix of uncorrelated shocks  
D_max = 5;              % Maximum degree of polynomial: the program compares 
                        % costs of evaluating expectations of ordinary polynomial
                        % of the degrees from one to D_max;
                        % (D_max can be from 1 to 5) 
                        
% 2. Steady state values of k and z
% ---------------------------------
ks = ones(1,N);         % Steady state of capital for N countries, equal to 1
zs = ones(1,N);         % Steady state of productivity level for N countries,
                        % equal to 1
                        
% 3. The matrix of the polynomial coefficients
% ---------------------------------------------  
for D = 1:D_max       % For the polynomial degrees from one to D_max...
    npol(D) = size(Ord_Polynomial_N([ks zs],D),2);      
end
                      % Construct polynomial bases and compute their
                      % number; the latter is needed for finding the number 
                      % of the polynomial coefficients in the policy functions)
BQ = zeros(npol(D_max),D_max); 
                      % Initialize the matrix of polynomial coefficients, b, 
                      % of an integrand function for the polynomial solutions 
                      % of the degrees from one to D_max; npol(D_max)-by-D_max                      
bq_mult = zeros(npol(D_max),D_max);
                      % Initialize the matrix of the precomputed integrals, 
                      % I, for the polynomial solutions of the degrees from
                      % one to D_max; npol(D_max)-by-D_max     
                                                           
% 4. Construct the sets of precomputed integrals
% ----------------------------------------------
                             
for D = 1:D_max;
    bq_mult(1:npol(D),D) = Integrals_PREC_N(N,N,vcv,D)';
                             % Compute the values of the integrals, I, of 
                             % future shocks, that multiply the polynomial
                             % coefficients b, i.e., I*b=b'; see equation 
                             % (53) in JMMT (2016) for a description of the 
                             % multidimensional integrals 
                           
                             
    % 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(1:npol(D),D) =  1 + eps;   % Polynomial coefficients of the integrand 
                                  % function; we add "eps" (the smallest 
                                  % number in MATLAB) to obtain
                                  % floating-point numbers
end                             
                                               
% 5. Evaluate the expectation function of polynomial of the given degree
% using precomputation
% ----------------------------------------------------------------------

time_prec = NaN(5,1);             % Initial computing times
EP = 0;                            % Initial value of the expectation function

for D = 1:D_max        
    tic           % Start counting time needed for evaluation of the expectation
    BQ_prime = BQ(1:npol(D),D).*bq_mult(1:npol(D),D);
                  % Compute the coefficients b'=b*I
    EP = Ord_Polynomial_N([ks,zs],D)*BQ_prime;
                  % Evaluate expectation of polynomial at a steady state;
                  % see equation (10) in JMMT (2016)
    time_prec(D,1) = toc;
                  % Time needed to evaluate expectation of polynomial 
                  % of degree D
end 

% 6. Evaluate the expectation function of a polynomial of a given degree
% using Monte Carlo, Gauss-Hermite, monomial integration rules
% ----------------------------------------------------------------------

% 6.1 Evaluate the expectation function of polynomial of a given degree
% using Monte Carlo integration with 100 nodes
% ----------------------------------------------------------------------

time_MC100 = NaN(5,1);        % Initial computing times
EP = 0;                       % Initial value of the expectation function
load('epsi30.mat','epsi30')   % Load previously saved sequence of shocks drawn 
                              % from a standard normal distribution

J = 100;                            % Number of the integration nodes  
epsi_nodes = sigma*epsi30(1:J,1:N); % Construct J integration nodes
z = exp(epsi_nodes);                % Productivity level at each integration 
                                    % node
                    
weight_nodes = ones(1,J)*1/J; 
                    % Construct J integration weights; Monte Carlo
                    % integration method assigns the same weight for each node

for D = 1:D_max        
    tic           % Start counting time needed for evaluation of the expectation
    for j = 1:J
        
        EP = EP + Ord_Polynomial_N([ks,z(j,:)],D)*BQ(1:npol(D),D)*weight_nodes(j);
                  % Evaluate expectation in the steady state
    end
    time_MC100(D,1) = toc;
                  % Time needed to evaluate expectation of polynomial 
                  % of degree D
end

FacMC100 = time_MC100./time_prec;
                  % Factor by which running time with Monte Carlo integration 
                  % with 100 nodes differs from precomputation

% 6.2 Evaluate the expectation function of polynomial of a given degree
% using Monte Carlo integration with 10,000 nodes
% ----------------------------------------------------------------------

time_MC10000 = NaN(5,1);      % Initial computing times

if N<=10
EP = 0;                        % Initial value of the expectation
load('epsi30.mat','epsi30')   % Load previously saved sequence of shocks drawn 
                              % from a standard normal distribution

J = 10000;                    % Number of integration nodes  
epsi_nodes = sigma*epsi30(1:J,1:N);  
                              % Construct J integration nodes
z = exp(epsi_nodes);          % Productivity level at each integration node
                    
weight_nodes = ones(1,J)*1/J; 
                              % Construct J integration weights; Monte Carlo
                              % integration method assigns the same weight 
                              % for each node

for D = 1:D_max        
    tic                       % Start counting time needed for evaluation 
                              % of the expectation
    for j=1:J
    EP = EP + Ord_Polynomial_N([ks,z(j,:)],D)*BQ(1:npol(D),D)*weight_nodes(j);
                              % Evaluate expectation at a steady state
    end
    time_MC10000(D,1)=toc;
                              % Time needed to evaluate expectation of 
                              % polynomial of degree D
end

end

FacMC10000 = time_MC10000./time_prec;
                  % Factor by which running time with Monte Carlo integration 
                  % with 10,000 nodes differs from precomputation
                  
% 6.3 Evaluate the expectation function of a polynomial of a given degree
% using Monomial (non-product) integration rule with 2N nodes
% ----------------------------------------------------------------------

time_M1 = NaN(5,1);           % Initial computing times
EP = 0;                       % Initial value of the expectation

[J,epsi_nodes,weight_nodes] = Monomials_1(N,vcv);
            % Compute the number of integration nodes, "n_nodes", integration 
            % nodes, "epsi_nodes", and integration weights, "weight_nodes",  
            % for monomial integration rule 
z = exp(epsi_nodes);  
            % Productivity level at each integration node

for D = 1:D_max        
    tic           % Start counting time needed for evaluation of the expectation
    for j = 1:J
    EP = EP + Ord_Polynomial_N([ks,z(j,:)],D)*BQ(1:npol(D),D)*weight_nodes(j);
                  % Evaluate expectation at a steady state
    end
    time_M1(D,1) = toc;
                  % Time needed to evaluate expectation of polynomial 
                  % of degree D
end

FacM1 = time_M1./time_prec;
                  % Factor by which running time with Monomial (non-product) 
                  % integration rule with 2N nodes differs from precomputation

% 6.4 Evaluate the expectation function of a polynomial of a given degree
% using Monomial (non-product) integration rule with 2N^2+1 nodes
% ----------------------------------------------------------------------

time_M2 = NaN(5,1);           % Initial computing times
EP = 0;                       % Initial value of the expectation

[J,epsi_nodes,weight_nodes] = Monomials_2(N,vcv);
            % Compute the number of integration nodes, "n_nodes", integration 
            % nodes, "epsi_nodes", and integration weights, "weight_nodes",  
            % for monomial integration rule 
z = exp(epsi_nodes);  % Productivity level at each integration node

for D = 1:D_max        
    tic           % Start counting time needed for evaluation of the expectation
    for j = 1:J
    EP = EP + Ord_Polynomial_N([ks,z(j,:)],D)*BQ(1:npol(D),D)*weight_nodes(j);
                  % Evaluate expectation at a steady state
    end
    time_M2(D,1) = toc;
                  % Time needed to evaluate expectation of polynomial 
                  % of degree D
end 

FacM2 = time_M2./time_prec;
                  % Factor by which running time with Monomial (non-product) 
                  % integration rule with 2N^2+1 nodes differs from precomputation

% 6.5 Evaluate the expectation function of a polynomial of a given degree
% using Gauss-Hermite quadrature rule with 2^N nodes
% ----------------------------------------------------------------------

time_GH2 = NaN(5,1);        % Initial computing times

if N<=20
EP = 0;                     % Initial value of the expectation
IM = 2;                     % Number of nodes in each dimension
[J,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 
            % monomial integration rule 
z = exp(epsi_nodes);  % Productivity level at each integration node

for D = 1:D_max        
    tic           % Start counting time needed for evaluation of the expectation
    for j = 1:J
    EP = EP+Ord_Polynomial_N([ks,z(j,:)],D)*BQ(1:npol(D),D)*weight_nodes(j);
                  % Evaluate expectation at a steady state
    end
    time_GH2(D,1) = toc;
                  % Time needed to evaluate expectation of polynomial 
                  % of degree D
end 

end

FacGH2 = time_GH2./time_prec;
                  % Factor by which running time with Gauss-Hermite quadrature 
                  % rule with 2^N nodes differs from precomputation

% 6.6 Evaluate the expectation function of a polynomial of a given degree
% using Gauss-Hermite quadrature rule with 5^N nodes
% ----------------------------------------------------------------------

time_GH5 = NaN(5,1);        % Initial computing times

if N<=5 
EP = 0;                     % Initial value of the expectation
IM = 5;                     % Number of nodes in each dimension
[J,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 monomial integration rule 
z = exp(epsi_nodes);  % Productivity level at each integration node

for D = 1:D_max        
    tic           % Start counting time needed for evaluation of the expectation
    for j=1:J
    EP = EP+Ord_Polynomial_N([ks,z(j,:)],D)*BQ(1:npol(D),D)*weight_nodes(j);
                  % Evaluate expectation at a steady state
    end
    time_GH5(D,1) = toc;
                  % Time needed to evaluate expectation of polynomial 
                  % of degree D
end 

end

FacGH5 = time_GH5./time_prec;
                  % Factor by which running time with Gauss-Hermite quadrature 
                  % rule with 5^N nodes differs from precomputation

% 6.7 Evaluate the expectation function of a polynomial of a given degree
% using Gauss-Hermite quadrature rule with 10^N nodes
% ----------------------------------------------------------------------

time_GH10 = NaN(5,1);        % Initial computing times

if N<=3
EP = 0;                      % Initial value of the expectation
IM = 10;                     % Number of nodes in each dimension
[J,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 monomial integration rule 
z = exp(epsi_nodes);  % Productivity level at each integration node

for D = 1:D_max        
    tic           % Start counting time needed for evaluation of the expectation
    for j = 1:J
    EP = EP + Ord_Polynomial_N([ks,z(j,:)],D)*BQ(1:npol(D),D)*weight_nodes(j);
                  % Evaluate expectation at a steady state
    end
    time_GH10(D,1) = toc;
                  % Time needed to evaluate expectation of polynomial 
                  % of degree D
end 

end

FacGH10 = time_GH10./time_prec;
                  % Factor by which running time with Gauss-Hermite quadrature 
                  % rule with 10^N nodes differs from precomputation



fprintf(1,'NUMBER OF SHOCKS = %i\n',N);
disp(''); 
fprintf('COST OF NUMERICAL INTEGRATION METHODS \n')
disp('Precomputation (running time)')
disp([time_prec])
disp('THE DIFFERENCE IN RUNNING TIME VERSUS PRECOMPUTATION')
disp('     GH2       GH5       GH10       M1       M2')       
disp([FacGH2 FacGH5 FacGH10 FacM1 FacM2])
disp('    MC100' );   
disp([FacMC100])
disp('   MC10000');
disp([FacMC10000])
