% This MATLAB software compares the accuracy of approximating 
% an integral I=E(exp(l*e)) (here, l=1,...,5 and e is a normally distributed 
% multivariate random variable) by different integration methods; see 
% "How to Solve Dynamic tochastic 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. "Accuracy_1.m"       compares the accuracy of approximating a unidimen-
%                         sional integral I=E(exp(l*e)) by different
%                         integration methods
% 2. "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) 
% 3. "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) 
% 4. "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)   
% 5."epsi.mat"            contains a fixed sequence of random numbers of 
%                         length 10,000 drawn from a 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".
% -------------------------------------------------------------------------

clc;
clear all;

% We approximate an integral I=E(exp(l*e)), where l=1,2,3,4,5 and e is a 
% normally distributed multivariate random shock to the productivity level 
% with mean 0 and variance-covariance matrix vcv

% 1. Parameters
% -------------
nshocks = 1;         % Number of stochastic shocks is 1
sigma = 0.2;         % Standard deviation of shock to the log of the 
                     % productivity level
vcv = sigma^2;       % Variance of shock

% 2. Compute exact value of integrals using equaton (22) of JMMT (2016)
% ---------------------------------------------------------------------
v = ones(5,1);       % Initial values of the integrals I = E(exp(l*e)), 
                     % l=1,2,3,4,5

for i=1:5
    v(i,1) = exp(i^2*vcv/2);  
                     % Compute exact values of I = E(exp(l*e)) 
end

% 3. Approximate the value of integrals using Monte Carlo, Gauss-Hermite, 
% M1 and M2 integration rules
% ----------------------------------------------------------------------

% 3.1 Monte-Carlo integration with 100 nodes
% ------------------------------------------
load('epsi.mat')    % Load previously saved sequence of shocks drawn from a 
                    % standard normal distribution

J = 100;            % Number of integration nodes  
epsi_nodes = sigma*epsi(1:J,1);
                    % Construct J integration nodes
weight_nodes = ones(1,J)*1/J; 
                    % Construct J integration weights; Monte Carlo
                    % integration method assigns the same weight for each node
vMC100 = zeros(5,1);% Initial approximated values of the integrals 

for i=1:5
    vMC100(i,1) = weight_nodes*exp(i*epsi_nodes);
                    % Compute approximated values of I = E(exp(l*e))
end

% Compute errors in approximation of I = E(exp(l*e)) by Monte Carlo 
% integration method with 100 nodes (measured in percent)

difMC100=(vMC100./v-1)*100;   % Errors in approximation 


% 3.2 Monte-Carlo integration with 10,000 nodes
% ---------------------------------------------
J = 10000;          % Number of integration nodes  
epsi_nodes = sigma*epsi(1:J,1);
                    % Construct J integration nodes
weight_nodes = ones(1,J)*1/J; 
                    % Construct J integration weights; Monte Carlo
                    % integration method assigns the same weight for each node
vMC10000 = zeros(5,1);% Initial approximated values of the integrals 

for i=1:5
    vMC10000(i,1) = weight_nodes*exp(i*epsi_nodes);
                    % Compute approximated values of I = E(exp(l*e))
end

% Compute errors in approximation of I = E(exp(l*e)) by Monte Carlo 
% integration method with 10,000 nodes (measured in percent)

difMC10000=(vMC10000./v-1)*100;   % Errors in approximation 


% 3.3 Monomial (non-product) integration rule with 2 nodes
% ---------------------------------------------------------

[n_nodes,epsi_nodes,weight_nodes] = Monomials_1(nshocks,vcv);
            % Compute the number of integration nodes, "n_nodes", integration 
            % nodes, "epsi_nodes", and integration weights, "weight_nodes", for 
            % monomial integration rule 

vM1 = zeros(5,1);   % Initial approximated values of the integrals 

for i=1:5
    vM1(i,1) = weight_nodes'*exp(i*epsi_nodes);
                    % Compute approximated values of I = E(exp(l*e))
end

% Compute errors in approximation of I = E(exp(l*e)) by monomial 
% integration rule with 2 nodes (measured in percent)

difM1=(vM1./v-1)*100;   % Errors in approximation 

% 3.4 Monomial (non-product) integration rule with 3 nodes
% ---------------------------------------------------------

[n_nodes,epsi_nodes,weight_nodes] = Monomials_2(nshocks,vcv);
            % Compute the number of integration nodes, "n_nodes", integration 
            % nodes, "epsi_nodes", and integration weights, "weight_nodes", for 
            % monomial integration rule 

vM2 = zeros(5,1);   % Initial approximated values of the integrals 

for i=1:5
    vM2(i,1) = weight_nodes'*exp(i*epsi_nodes);
                    % Compute approximated values of I = E(exp(l*e))
end

% Compute errors in approximation of I = E(exp(l*e)) by monomial 
% integration rule with 3 nodes (measured in percent)

difM2=(vM2./v-1)*100;   % Errors in approximation 

% 3.5 Gauss-Hermite quadrature rule with 2 nodes
% ---------------------------------------------------------
IM = 2;     % Number of nodes
[n_nodes,epsi_nodes,weight_nodes] = GH_Quadrature(IM,nshocks,vcv);
            % Compute the number of integration nodes, "n_nodes", integration 
            % nodes, "epsi_nodes", and integration weights, "weight_nodes", for 
            % monomial integration rule 

vGH2 = zeros(5,1);   % Initial approximated values of the integrals 

for i=1:5
    vGH2(i,1) = weight_nodes'*exp(i*epsi_nodes);
                    % Compute approximated values of I = E(exp(l*e))
end

% Compute errors in approximation of I = E(exp(l*e)) by Gauss-Hermite 
% quadrature rule with 2 nodes (measured in percent)

difGH2=(vGH2./v-1)*100;   % Errors in approximation 

% 3.6 Gauss-Hermite quadrature rule with 5 nodes
% ---------------------------------------------------------
IM = 5;     % Number of nodes
[n_nodes,epsi_nodes,weight_nodes] = GH_Quadrature(IM,nshocks,vcv);
            % Compute the number of integration nodes, "n_nodes", integration 
            % nodes, "epsi_nodes", and integration weights, "weight_nodes",
            % for monomial integration rule 

vGH5 = zeros(5,1);   % Initial approximated values of the integrals 

for i = 1:5
    vGH5(i,1) = weight_nodes'*exp(i*epsi_nodes);
                    % Compute approximated values of I = E(exp(l*e))
end

% Compute errors in approximation of I = E(exp(l*e)) by Gauss-Hermite 
% quadrature rule with 5 nodes (measured in percent)

difGH5 = (vGH5./v-1)*100;   % Errors in approximation 

% 3.7 Gauss-Hermite quadrature rule with 10 nodes
% ---------------------------------------------------------
IM = 10;    % Number of nodes
[n_nodes,epsi_nodes,weight_nodes] = GH_Quadrature(IM,nshocks,vcv);
            % Compute the number of integration nodes, "n_nodes", integration 
            % nodes, "epsi_nodes", and integration weights, "weight_nodes", 
            % for monomial integration rule 

vGH10 = zeros(5,1);   % Initial approximated values of the integrals 

for i = 1:5
    vGH10(i,1) = weight_nodes'*exp(i*epsi_nodes);
                    % Compute approximated values of I = E(exp(l*e))
end

% Compute errors in approximation of I = E(exp(l*e)) by Gauss-Hermite 
% quadrature rule with 10 nodes (measured in percent)

difGH10 = (vGH10./v-1)*100;   % Errors in approximation 

% 4. Display accuracy results
% ---------------------------
fprintf('ACCURACY OF UNIVARIATE NUMERICAL INTEGRATION METHODS \n')
disp('Exact value    GH2       GH5     GH10      MC100    MC10000'); 
disp([v difGH2 difGH5 difGH10 difMC100 difMC10000 ])
