% MATLAB software that solves a new Keynesian model using an epsilon  
% distinguishable set algorithm (EDS) and cluster grid algorithm (CGA), as 
% described in the article "Merging simulation and projection approaches to  
% solve high-dimensional problems with an application to a new Keynesian 
% model" by Lilia Maliar and Serguei Maliar (2015), Quantitative Economics 
% 6/1, pages 1–47 (henceforth, MM, 2015). 
%
% This software is written to illustrate a general methodology proposed in 
% MM (2015) and to replicate the results in that paper. The running time was
% not the priority. The following modifications will speed up the code (the
% optimized code can be as fast as 5 seconds at only moderate accuracy loss):
% 
% 1. Use a shorter simulation length, for example, T = 10,000 (instead of 
% T=100,000) and do not remove the low density points. 
% 2. Use the simulated series produced by perturbation method for constructing
% the EDS / CGA grid and do not update the grid along iterations. 
% 3. Choose the number of grid points to be slightly larger than the number of 
% the polynomial coefficients, for example, use 50 grid points for 2nd degree 
% polynomial approximation with 45 coefficients. Do not iterate on "epsilon" 
% by using a bisection but use some prior experimentation to find the value 
% of "epsilon" that leads to the target number of grid points. 
% 4. Increase damping parameters "damp" from [0.1 0.1 0.1] to larger values
% that still leads to convergence (such values depend on a specific model).
% 
% This software is based on that of Lilia  Maliar and Serguei Maliar for 
% solving models using the generalized stochastic simulation algorithm (GSSA) 
% method, as described in the paper "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). The modifications made are concerned with 
% the construction of the grid.
%
% This version: March 19, 2015. First version: June 27, 2011.
% 
% ------------------------------------------------------------------------
% The software uses the following files: 
% ------------------------------------------------------------------------
% 1. "Main_NK_Degree2.m"         is a main file for computing the EDS and 
%                                CGA solutions to a new Keynesian model;
%                                it produces a polynomial approximation of 
%                                degree 2; it uses a Dynare routine  
%                                "NK_Dynare.mod" to find a perturbation
%                                solution
% 2. "NK_Degree3.m"              produces a polynomial approximation of 
%                                degree 3 using the polynomial approximation 
%                                of degree 2 generated by "Main_NK_Degree2.m" 
% 3. "NK_Dynare.mod"             is a Dynare routine that produces perturbation
%                                solutions to be used as an initial guess 
% 4. "NK_Simulation.m"           simulates a time series solution to the new 
%                                Keynesian model for a given sequence of 
%                                exogenous shocks
% 5. "NK_Accuracy.m"             computes residuals of the equilibrium
%                                conditions of the new Keynesian model on a   
%                                given set of points in the state space
% 6. "Density.m"                 estimates the density function from a 
%                                given set of points 
% 7. "Clusters.m"                constructs clusters from simulated series 
%                                and computes clusters' centers (to be used 
%                                as a grid) 
% 8. "EDS"                       constructs an epsilon distinguishable set 
%                                for a given set of data (to be used as a 
%                                grid)
% 9. "Ord_Polynomial_N.m"        constructs the sets of basis functions for 
%                                ordinary polynomials of the degrees from 
%                                one to five; borrowed from JMM (2011)
% 10. "Monomials_1.m"            constructs integration nodes and weights for 
%                                an N-dimensional monomial (non-product) 
%                                integration rule with 2N nodes; borrowed from 
%                                JMM (2011) 
% 11. "Monomials_2.m"            constructs integration nodes and weights for 
%                                an N-dimensional monomial (non-product) 
%                                integration rule with 2N^2+1 nodes; borrowed 
%                                from JMM (2011)
% 12. "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 JMM (2011)                     
% 13. "NK_shocks_100000.mat"     contains the series of the exogenous shocks 
%                                of length 100,000 for the solution procedure 
%                                and the series of the exogeneous shocks of 
%                                length 10,200 for the test 
% -------------------------------------------------------------------------
% Copyright © 2015 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".
% -------------------------------------------------------------------------

clc;
clear all;
                                 
% 1. Parameter values 
% -------------------
gam        = 1;           % Utility-function parameter
betta      = 0.99;        % Discount factor
vartheta   = 2.09;        % Utility-function parameter
epsil      = 4.45;        % Parameter in the Dixit-Stiglitz aggregator
phi_y      = 0.07;        % Parameter of the Taylor rule
phi_pie    = 2.21;        % Parameter of the Taylor rule
mu         = 0.82;        % Parameter of the Taylor rule
theta      = 0.83;        % Share of non-reoptimizing firms (Calvo's pricing)
piestar    = 1;           % Target (gross) inflation rate 
Gbar       = 0.23;        % Steady-state share of government spending in output

%  Autocorrelation coefficients in the processes for shocks
%----------------------------------------------------------
rho_nua    = 0.95;        % See process (22) in MM (2015)
rho_nuL    = 0.25;        % See process (16) in MM (2015)
rho_nuR    = 0.15;        % See process (28) in MM (2015)
rho_nuu    = 0.92;        % See process (15) in MM (2015)
rho_nuB    = 0.22;        % See process (17) in MM (2015)
rho_nuG    = 0.95;        % See process (26) in MM (2015)

%  Standard deviations of the innovations in the processes for shocks
%--------------------------------------------------------------------
sigma_nua  = 0.0045;      % See process (22) in MM (2015) 
sigma_nuL  = 0.1821;      % See process (16) in MM (2015)
sigma_nuR  = 0.0028;      % See process (28) in MM (2015)
sigma_nuu  = 0.0054;      % See process (15) in MM (2015)
sigma_nuB  = 0.0023;      % See process (17) in MM (2015)
sigma_nuG  = 0.0038;      % See process (26) in MM (2015)

% Save the parameters into "parametrfile" that will be used by "NK_Dynare.mod"
%----------------------------------------------------------------------------
save parameterfile  betta piestar Gbar gam vartheta epsil phi_y phi_pie mu theta ...
	rho_nua rho_nuu rho_nuG rho_nuL rho_nuR rho_nuB ...
	sigma_nua sigma_nuu sigma_nuG sigma_nuL sigma_nuR sigma_nuB;

% 2. Compute a Dynare solution using "NK_Dynare.mod"
% --------------------------------------------------
path(path,'4.2.1\matlab')  % Add a Dynare path
dynare NK_Dynare;          % Compute the first- and second-order perturbation 
                           % solutions 

% 3. Recover the matrices that are components of the decision rules produced 
% by Dynare 
%---------------------------------------------------------------------------
Order = oo_.dr.order_var;  % Check that the order in which the vaiables 
                           % appear in decision rules (DR-order) coincides 
                           % with the order in which the variables are declared

% 3.1 Recover the matrices relevant for the first-order perturbation solution 
% (PER1)
%---------------------------------------------------------------------------
% The first-order approximation has the form y(t)=SS+A*(y(t-1)-SS)+B*u(t),
% where y(t) and u(t) are endogenous and exogenous (i.e., shocks) variables,
% resp.                           
SS    = oo_.dr.ys;         % The steady state values of all variables appearing
                           % in the order in which the endogenous variables 
                           % are declared in the Dynare program (there are 15 
                           % variables); 15-by-1

A     = oo_.dr.ghx;        % A matrix of coefficients whose rows correspond
                           % to all the endogenous variables in the DR-order,
                           % and whose columns correspond to the state
                           % variables in the DR-order; 15-by-8
B     = oo_.dr.ghu;        % A matrix of coefficients whose rows correspond
                           % to all the endogenous variables in the DR-order,
                           % and whose columns correspond to the exogenous
                           % variables (shocks) in the DR-order; 15-by-6

% 3.2 Recover the matrices relevant for the second-order perturbation 
% solution (PER2)
%---------------------------------------------------------------------                          
% The second-order approximation has the form: 
% y(t)=SS+0.5*del2+A(y(t-1)-SS)+B*u(t)+0.5*C*(y(t-1)-SS)*kron*(y(t-1)-SS)+0.5*D*
% *(u(t)*kron*u(t)+E*(y(t-1)-SS)*kron*u(t))

del2  = oo_.dr.ghs2;       % A column vector whose rows correspond to all 
                           % endogenous variables in the DR-order; 15-by-1
C     = oo_.dr.ghxx;       % A matrix of coefficients whose row correspond
                           % to all the endogenous variables in the
                           % DR-order, and whose columns correspond to the
                           % Kronecker product of the vector of state
                           % variables (there are 6 of them) in the DR-order; 
                           % 15-by-64
D     = oo_.dr.ghuu;       % A matrix of coefficients whose row correspond
                           % to all the endogenous variables in the
                           % DR-order, and whose columns correspond to the
                           % Kronecker product of the exogenous variables
                           % (there are 6 of them) in the declaration order; 
                           % 16-by-36
E     = oo_.dr.ghxu;       % A matrix of coefficients whose row correspond
                           % to all the endogenous variables in the
                           % DR-order, and whose columns correspond to the
                           % Kronecker product of the vector of state
                           % variables (there are 6 of them) in the DR-order 
                           % by the vector of exogenous variables (there are
                           % 6 of them) in the declaration order; 15-by-48


% 4. Load shocks for simulating the exogenous state variables 
%------------------------------------------------------------
load NK_shocks_100000;     % The matrix of shocks used for simulation of the 
                           % model is called "epsi", and that used for
                           % the accuracy test is called "epsi_test" 


% 5. Simulate the time-series perturbation solutions
%---------------------------------------------------                           
% We simulate the Dynare solution "manually" instead of using a Dynare
% simulation routine because we need to impose the zero lower bound (ZLB)  
% on nominal interest rate in the simulated solution 

T = 10000;               % The simulation length 

% 5.1 Simulate a perturbation solution for T observations without ZLB imposed
%----------------------------------------------------------------------------
zlb = 0;                 % Set a dummy variable "zlb" to 0 (no ZLB imposed) 
[PER1 PER2]  = NK_PER12_simulation(SS,del2,A,B,C,D,E,epsi(1:T,:),sigma_nua,sigma_nuL,sigma_nuR,sigma_nuu,sigma_nuB,sigma_nuG,T,zlb);
                         % Simulate the model using the shocks epsi(1:T,:)

% 5.2 Simulate a perturbation solution for T observations with ZLB imposed
%-------------------------------------------------------------------------
zlb = 1;                 % Set a dummy variable "zlb" to 1 (ZLB is imposed)
[PER1_zlb PER2_zlb]  = NK_PER12_simulation(SS,del2,A,B,C,D,E,epsi(1:T,:),sigma_nua,sigma_nuL,sigma_nuR,sigma_nuu,sigma_nuB,sigma_nuG,T,zlb);
                         % Simulate the model using the shocks epsi(1:T,:)

% _____________________________________________________________
%
% The second-degree polynomial solution under the EDS algorithm
% _____________________________________________________________

time0 = cputime;         % Start counting time for computing the EDS solution

% 6. The number of state variables, simulation length, zlb dummy,
% polynomial degree
% ---------------------------------------------------------------
N     = 8;       % Indicate the number of state variables
T     = 10000;   % Choose the simulation length for the solution procedure,
                 % T<=10,000   

zlb   = 1;       % Set the dummy variable "zlb" to:
                 %                     =0 (no ZLB imposed), 
                 %                     =1 (zlb imposed) 
                
Degree = 2;      % Consider the second-degree polynomial approximation

% 7. Initial condition from the perturbation solutions
% ----------------------------------------------------
% The order of variables in the perturbation solution is
% L Y Yn nua R delta nuL nuR nuG nuB nuu pie S F  C

% The correspondence between Dynare and EDS variables is tricky (we should 
% pay attention to the time indices)

 % The series for shocks
 nua0 = PER1(4,2:T)'; % The variable "nua" is the 4th in PER1
 nuL0 = PER1(7,2:T)';
 nuR0 = PER1(8,2:T)';
 nuG0 = PER1(9,2:T)';
 nuB0 = PER1(10,2:T)';
 nuu0 = PER1(11,2:T)';

 % The series for the current-period variables 
 if zlb==1            % ZLB is imposed
    R_nl = PER2_zlb(5,1:T-1)';
    delta_nl = PER2_zlb(6,1:T-1)';
    S_per2    = PER2_zlb(13,2:T+1)';
    F_per2    = PER2_zlb(14,2:T+1)';
    C_per2    = PER2_zlb(15,2:T)';
 else                 % ZLB is not imposed
    R_nl = PER2(5,1:T-1)';
    delta_nl = PER2(6,1:T-1)';
    S_per2 = PER2(13,2:T+1)';
    F_per2 = PER2(14,2:T+1)';
    C_per2 = PER2(15,2:T)';
 end

% 8. Compute an initial condition for the polynomial coefficients
%----------------------------------------------------------------
Data = [log(R_nl(1:9999,1)) log(delta_nl(1:9999,1)) nuR0 nua0 nuL0 nuu0 nuB0 nuG0];
% Form a matrix of time series of the state variables; (T-1)-by-N

X0 = Ord_Polynomial_N(Data(1:T-1,:),Degree);
% Construct the matrix of explanatory variables X0 on the series of state 
% variables from the previously computed time-series perturbation solution; 
% columns of X0 are given by the basis functions of the polynomial of degree 
% "Degree" 

npol_2d = size(X0,2);       % Number of coefficients in polynomial of degree 
                            % "Degree" (recall that we chose Degree=2)
                            
VK = zeros(npol_2d,3,2);    % Initialize the matrix of polynomial coefficients; 
                            % npol_2d-by-3-by-2, where 3 is the number of 
                            % decision rules we compute, i.e., for S, F,  
                            % (C^(-gam) (i.e., marginal utility), and 2 
                            % stands for 2 degrees of approximations we  
                            % consider, of degree 1 and degree 2; npol_2d-by-3-by-2 
                            
VK(:,1:3,1)   = X0\[S_per2(1:T-1,1) F_per2(1:T-1,1) C_per2(1:T-1,1).^(-gam)];
                            % Use a backslach operator to regress [S F C^(-gam)]
                            % on X0 to obtain the coefficients of the
                            % decision rules

% 9. Generate long series of shocks 
%----------------------------------
Tepsi = 100000;             % Choose a long simulation length; the shocks 
                            % generated with this simulation length will be 
                            % used to simulate the solution with the aim of 
                            % creating an EDS grid

% Initialize the series for our 6 shocks; as an initial value of a shock, use 
% the same  value as in PER1 solution;
nua = zeros(Tepsi,1) + PER1(4,1);  
nuL = zeros(Tepsi,1) + PER1(7,1);
nuR = zeros(Tepsi,1) + PER1(8,1);
nuG = zeros(Tepsi,1) + PER1(9,1);
nuB = zeros(Tepsi,1) + PER1(10,1);
nuu = zeros(Tepsi,1) + PER1(11,1);

% Use the processes for shocks to generate the values of shocks in all
% periods
for t = 1:Tepsi-1;
        nua(t+1,1) = nua(t,1)*rho_nua + epsi(t+1,2)*sigma_nua; 
                             % See process (22) in MM (2015)
        nuL(t+1,1) = nuL(t,1)*rho_nuL + epsi(t+1,3)*sigma_nuL; 
                             % See process (16) in MM (2015)
        nuR(t+1,1) = nuR(t,1)*rho_nuR + epsi(t+1,1)*sigma_nuR; 
                             % See process (28) in MM (2015)
        nuG(t+1,1) = nuG(t,1)*rho_nuG + epsi(t+1,6)*sigma_nuG;
                             % See process (15) in MM (2015)
        nuB(t+1,1) = nuB(t,1)*rho_nuB + epsi(t+1,5)*sigma_nuB; 
                             % See process (17) in MM (2015)
        nuu(t+1,1) = nuu(t,1)*rho_nuu + epsi(t+1,4)*sigma_nuu; 
                             % See process (26) in MM (2015)
end

% 10. Initialize the second-degree capital policy functions for EDS
% -----------------------------------------------------------------
vk_2d = VK(:,:,1);         % Start from the second-degree polynomial solution
                           % obtained in paragraph 8

% 11. Select the integration method 
% ---------------------------------

vcv = diag([sigma_nuR^2 sigma_nua^2 sigma_nuL^2 sigma_nuu^2 sigma_nuB^2 sigma_nuG^2]);
                             % Variance covariance matrix

[n_nodes,epsi_nodes,weight_nodes] = Monomials_1(6,vcv);
                             % Monomial integration rule with 2N nodes
                             
%[n_nodes,epsi_nodes,weight_nodes] = Monomials_2(N,vcv);
                             % Monomial integration rule with 2N^2+1 nodes

%[n_nodes,epsi_nodes,weight_nodes] = Monomials_3(N,vcv);
                             % Monomial integration rule with 2^N nodes;  
                             % this rule coincides with Gauss-Hermite 
                             % quadrature (product) integration rule with 
                             % 2 nodes in each dimension 

%n_nodes = 1; epsi_nodes = zeros(1,N); weight_nodes = 1;                             
                             % Gauss-Hermite quadrature integration rule 
                             % with one node 

% 12. Initial values of the parameters of the EDS algorithm 
% ----------------------------------------------------------
epsilon = 0;                 % Initial value of the parameter that measures 
                             % the minimum distance between points in the 
                             % EDS grid
Di_G_min2 = 1e10;            % The minimum distance between a new  grid
                             % point and all old grid points (needed to
                             % check the convergence of the grid); see 
                             % discussion on page 22 of MM (2015)


% 13. The main iterative cycle of the EDS algorithm
% --------------------------------------------------
j = 0;                       % Index of an iteration on the grid construction; 
                             % iterations are performed util the grid converges; 
                             % initially, it is 0

while max(Di_G_min2)>(2*epsilon);    
                             % The criterion for convergence of the grid is
                             % that the distance from each point of the old 
                             % grid to some point of the new grid is smaller 
                             % than 2*epsilon
                                               
j = j+1;    
 
% 13.1 Construct an EDS grid
% --------------------------- 
% A grid of representative points will be used to compute a solution to 
% the model. An EDS grid is constructed from a given set of simulated points 
% using "EDS.m" (it constructs an epsilon distinguishable subset of the simulated 
% data points such that any two points in this subset are located at the 
% distance at least epsilon). Another option is to construct a cluster grid  
% using "Clusters.m"    

  % 13.1.1 Create by simulation a matrix of data on which the EDS grid will 
  % be constructed
  % ------------------------------------------------------------------------
  [S_nl F_nl delta_nl C_nl Y_nl Yn_nl L_nl R_nl pie_nl] = NK_EDS_simulation(vk_2d,nuR,nua,nuL,nuu,nuB,nuG,R_nl(1,1),delta_nl(1,1),gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,piestar,Gbar,zlb,Degree);
                            % Simulate the solution using the vector of
                            % coefficients vk_2d (see paragraph 10)
  kth = round(Tepsi/10000); % Consider each kth simulated point (for Tepsi=
                            % 100000, this corresponds to each 10th point)


   for t = 1:10000    
        Data(t,:) = [log(R_nl(t*kth,1)) log(delta_nl(t*kth,1)) nuR(t*kth,1) nua(t*kth,1) nuL(t*kth,1) nuu(t*kth,1) nuB(t*kth,1) nuG(t*kth,1)]; 
                            % From the simulated series, construct a matrix 
                            % of data that will be used to obtain the EDS
                            % grid
   end


  
  % 13.1.2 Estimate the density function and remove low-density points 
  % ------------------------------------------------------------------   
  [density,PCn,Di_min] = Density(Data,Data);
                       % Estimate the density function "density" in all 
                       % simulated points; "PCn" is the normalized principal 
                       % components of the matrix of data for which the 
                       % density is estimated; "Di_min" is the distance 
                       % from each point in "Data" to the closest neighbors                    
                        
  Data_sort = sortrows([density Di_min Data],1);
                       % Sort the simulated points by the density function 
                        
  cutoff = round(Tepsi/kth*0.01);
                       % Cutoff level is 1% of points with the lowest density

  Data1 = Data_sort(1+cutoff:end,3:end);
                       % Remove the low density points with the given cutoff 
                       % level

  % 13.1.3 Construct EDS with a target number of M grid points using bisection
  % --------------------------------------------------------------------------  
    
    % Construct EDS sets with the parameters "epsilon1" and "epsilon2" that 
    % correspond to the limiting hyperspheres that contain none and all data 
    % points, respectively; these are the initial values that are necessary 
    % for bisection; see Algorithm M_bar of MM (2015)
    % ---------------------------------------------------------------------

    M = 500;           % Choose the number of grid points (must be larger 
                       % than the number of polynomial terms
                       % (coefficients); for the EDS, this is a target number
                       % of grid points (for the cluster grid, it is an exact
                       % number)

    N = size(Data,2);  % Number of state variables;                  
                       

    PCn_sort = sortrows([density Di_min PCn],1); 
                      % Sort the principal components "PCn" by the density 
                      % function
    r1 = min(sqrt(sum(PCn_sort(1+cutoff:end,3:end).^2,2)));
                      % The distance from the center to the closest point;  
                      % a ball with the radius r1 contains no data points
    r2 = max(sqrt(sum(PCn_sort(1+cutoff:end,3:end).^2,2)));
                      % The distance from the center to the furthest point; 
                      % a ball with the radius r2 contains all data points
    epsilon1 = r1/2/M^(1/N);          
                      % The lower bound on epsilon used in bisection
    epsilon2 = r2/(M^(1/N)-1);  
                      % The upper bound on epsilon used in bisection
    
    % Construct an EDS set for a given "epsilon" that contains a target 
    % (i.e., approximate) number of points equal to M 
    % ---------------------------------------------------------------------
    size_Grid_old = 0;% The number of points in the old grid; initially,
                      % it is 0
    
    Grid = zeros(M,N);% Initialize the EDS grid; initially, consists of M 
                      % points   

    while abs(size_Grid_old - size(Grid,1))>0;  
                      % Compute the difference in the number of points in 
                      % the new and old grids
                                              
        size_Grid_old = size(Grid,1);         
                      % Store the number of points in the old grid
        epsilon = (epsilon1 + epsilon2)/2;          
                      % Compute epsilon for constructing the new grid
        [Grid] = EDS(Data1,epsilon);          
                      % Construct the new grid
     if size(Grid,1)>M;% If the new grid has more points than needed, ...
        epsilon1 = epsilon;                     
                      % Then, increase the lower bound
     else             % Otherwise, ...  
        epsilon2 = epsilon;                     
                      % Decrease the upper bound
     end
    end  

    
 % 13.2 Compute the number of grid points
 %----------------------------------
 n_G = size(Grid,1);    % The number of points in the constructed EDS grid
                        % that corresponds to the current iteration; the 
                        % obtained grid is of size n_G-by-N
 N_G(j) = size(Grid,1); % The number of points in the grid obtained on the 
                        % jth iteration (needed to check the convergence of 
                        % the grid)

 % 13.3 Allocate memory to the integrals in the right side of the 3 Euler 
 % equations that we parameterize; see conditions (32), (33), (34) in MM (2015)
 %-------------------------------------------------------------------------
 e = zeros(n_G,3);     

 % 13.4 Form complete polynomial of a degree "Degree" on the grid points for  
 % future regressions
 %-------------------------------------------------------------------------
 X0_G = Ord_Polynomial_N(Grid,Degree);   
                              % For Degree=2, the matrix "X0_G" consists of 
                              % a column of ones, the linear and quadratic  
                              % polynomial bases on the grid points

 % 13.5 Allocate memory to the next-period values of S, F and C on the grid 
 %------------------------------------------------------------------------
   % Allocate memory to S, F, C obtained in the previous iteration (to check
   % convergence)
   S0_old_G = ones(n_G,1);
   F0_old_G = ones(n_G,1);
   C0_old_G = ones(n_G,1);
   
   % Allocate memory to S, F, C obtained in the current iteration (to check
   % convergence) 
   S0_new_G = ones(n_G,1);     
   F0_new_G = ones(n_G,1);     
   C0_new_G = ones(n_G,1);     
                             
 % 13.6. The parameters of the EDS algorithm 
 % -----------------------------------------
 damp     = [0.1 0.1 0.1];% Damping parameter for (fixed-point) iteration on 
                          % the coefficients of the 3 policy functions (for
                          % S, F and C^(-gam))
 dif_EDS	  = 1e+10;    % Convergence criterion (initially is not satisfied)

 % 13.7 The main iterative cycle of the EDS algorithm
 % --------------------------------------------------              
 while dif_EDS > 1e-7;% The convergence criterion (which is unit free 
                      % because dif_EDS is unit free)
                       
    for i=1:n_G;      % For each grid point, ... 

        % 13.7.1 Variables in a grid point i 
        % ----------------------------------   
        
         % Endogenous state variables on the grid
         %---------------------------------------
         R0  = exp(Grid(i,1));
         delta0  = exp(Grid(i,2));
        
         % Exogenous state variables on the grid
         %--------------------------------------
         nuR0 = Grid(i,3);          % The variable "nuR" appears in the 3d
                                    % column of "Grid"
         nua0 = Grid(i,4);            
         nuL0 = Grid(i,5);           
         nuu0 = Grid(i,6);            
         nuB0 = Grid(i,7);           
         nuG0 = Grid(i,8);            
        
         % Compute future shocks in all grid points and all integration nodes
         % ------------------------------------------------------------------                           
         nuR1(1:n_nodes,1) = (ones(n_nodes,1)*nuR0)*rho_nuR+epsi_nodes(:,1); 
         nua1(1:n_nodes,1) = (ones(n_nodes,1)*nua0)*rho_nua+epsi_nodes(:,2); 
         nuL1(1:n_nodes,1) = (ones(n_nodes,1)*nuL0)*rho_nuL+epsi_nodes(:,3); 
         nuu1(1:n_nodes,1) = (ones(n_nodes,1)*nuu0)*rho_nuu+epsi_nodes(:,4); 
         nuB1(1:n_nodes,1) = (ones(n_nodes,1)*nuB0)*rho_nuB+epsi_nodes(:,5); 
         nuG1(1:n_nodes,1) = (ones(n_nodes,1)*nuG0)*rho_nuG+epsi_nodes(:,6);
         % The size of each od these vectors is 1-by-n_nodes


        % Give a name to an ith row of previously created complete polynomial 
        % "X0_G" (the one that corresponds to the grid point i)
        %------------------------------------------------------------------
        X0 = X0_G(i,:);                 % This is done for convenience 
                                    
        % 13.7.2 Current-period choices in a grid point i
        % -----------------------------------------------
        S0 = X0*vk_2d(:,1);              % Compute S(i) using vk_2d
        F0 = X0*vk_2d(:,2);              % Compute F(i) using vk_2d
        C0 = (X0*vk_2d(:,3)).^(-1/gam);  % Compute C(i) using vk_2d 
        pie0 = ((1-(1-theta)*(S0/F0)^(1-epsil))/theta)^(1/(epsil-1));
                                         % Compute pie(i) from condition 
                                         % (35) in MM (2015)
        delta1 = ((1-theta)*((1-theta*pie0^(epsil-1))/(1-theta))^(epsil/(epsil-1))+theta*pie0^epsil/delta0)^(-1);
                                         % Compute delta(i) from condition 
                                         % (36) in MM (2015)
        Y0 = C0/(1-Gbar/exp(nuG0));      % Compute Y(i) from condition (38)
                                         % in MM (2015)
        L0 = Y0/exp(nua0)/delta1;        % Compute L(i) from condition (37) 
                                         % in MM (2015)
        Yn0 = (exp(nua0)^(1+vartheta)*(1-Gbar/exp(nuG0))^(-gam)/exp(nuL0))^(1/(vartheta+gam));
                                         %  Compute Yn(i) from condition (31) 
                                         % in MM (2015)
        R1 = piestar/betta*(R0*betta/piestar)^mu*((pie0/piestar)^phi_pie * (Y0/Yn0)^phi_y)^(1-mu)*exp(nuR0);   
                                         % Compute R(i) from conditions (27),
                                         % (39) in MM (2015)
        if zlb==1;R1=R1.*(R1>=1)+(R1<1); end
                                         % If ZLB is imposed, set R(i)=1 
                                         % if ZLB binds
        
        % 13.7.3 Next-period choices in grid point i
        %-------------------------------------------
        delta1_dupl = ones(n_nodes,1)*delta1; 
        R1_dupl = ones(n_nodes,1)*R1;
        % Duplicate "delta1" and "R1" n_nodes times to create a matrix with
        % n_nodes identical rows; n_G-by-n_nodes
                

        X1 = Ord_Polynomial_N([log(R1_dupl) log(delta1_dupl) nuR1 nua1 nuL1 nuu1 nuB1 nuG1],Degree);
        % Form complete polynomial of degree "Degree" on next-period state 
        % variables; n_nodes-by-npol_2d 
       
        S1 = X1*vk_2d(:,1);             % Compute next-period S using the 
                                        % corresponding vk_2d
        F1 = X1*vk_2d(:,2);             % Compute next-period F using the 
                                        % corresponding vk_2d
        C1 = (X1*vk_2d(:,3)).^(-1/gam); % Compute next-period C using the 
                                        % corresponding vk_2d 
        pie1 = ((1-(1-theta)*(S1./F1).^(1-epsil))/theta).^(1/(epsil-1));
                                        % Compute next-period pie using 
                                        % condition (35) in MM (2015)

        % 13.7.4. Evaluate conditional expectations in the Euler equations
        %------------------------------------------------------------------
        e(i,1) = (exp(nuu0)*exp(nuL0)*L0^vartheta*Y0/exp(nua0) + betta*theta*pie1.^epsil.*S1)'*weight_nodes; 
        e(i,2) = (exp(nuu0)*C0^(-gam)*Y0 + betta*theta*pie1.^(epsil-1).*F1)'*weight_nodes;
        e(i,3) = (betta*exp(nuB0)/exp(nuu0)*R1*exp(nuu1).*C1.^(-gam)./pie1)'*weight_nodes;

        % 13.7.5 Variables of the current iteration 
        %------------------------------------------
        S0_new_G(i,1) = S0(1,1);     
        F0_new_G(i,1) = F0(1,1);
        C0_new_G(i,1) = C0(1,1);

    end  
    
 % 13.8 Compute and update the coefficients of the policy functions 
 % ----------------------------------------------------------------
 vk_hat_2d = X0_G\e; % Compute the new coefficients of the 3 policy functions                          
 vk_2d(:,1) = damp(1)*vk_hat_2d(:,1) + (1-damp(1))*vk_2d(:,1); 
                     % Update the coefficients using damping
 vk_2d(:,2) = damp(2)*vk_hat_2d(:,2) + (1-damp(2))*vk_2d(:,2); 
                     % Update the coefficients using damping
 vk_2d(:,3) = damp(3)*vk_hat_2d(:,3) + (1-damp(3))*vk_2d(:,3); 
                     % Update the coefficients using damping
  
 % 13.9 Evaluate the percentage (unit-free) difference between the values on  
 % the grid from the previous and current iterations
 % -------------------------------------------------------------------------
 dif_EDS = mean(mean(abs(1-S0_new_G./S0_old_G)))/damp(1)+mean(mean(abs(1-F0_new_G./F0_old_G)))/damp(2)+mean(mean(abs(1-C0_new_G./C0_old_G)))/damp(3)
                   % The convergence criterion is adjusted to the damping 
                   % parameters   
                                          
 % 13.10 Store the values on the grid to be used on the subsequent iteration 
 %-------------------------------------------------------------------------
 S0_old_G = S0_new_G; 
 F0_old_G = F0_new_G; 
 C0_old_G = C0_new_G; 
end

 % 13.11 Check the criterion for he grid convergence, i.e., the distance 
 % from each point of the new grid to some point point of the old grid is 
 % smaller than 2*epsilon
 %-------------------------------------------------------------------------
 [densityG,PCnG,Di_min] = Density(Data1,Grid);  
                          % Estimate the density "densityG" in the grid 
                          % points, given by the matrix "Grid"; "Data1" is 
                          % the matrix of all data from which we have 
                          % previously constructed the grid; "PCnG" is 
                          % the normalized principal components of the matrix 
                          % "Grid"; "Di_min" is the distance from each point 
                          % in "Data1" to the closest neighbors
    if j>1;                         
        Di_G_min2 = zeros(N_G(j),1);
                          % Initialize the minimum distance between the new 
                          % and old grid points (needed to check the 
                          % convergence of the grid)
        for ii = 1:N_G(j) % For each grid point, ...                
            Di_2 = (ones(N_G(j-1),1)*PCnG(ii,:)-PCnG_old).^2*ones(8,1);
                          % Compute the squared distance between a new 
                          % grid point ii and all old grid points (there are 
                          % N_G(j-1) of them); N_G(j-1)-by-1
            Di_2(ii,1) = inf;   
                          % Set the squared distance between the point  
                          % ii and itself to infinity 
            Di_G_min2(ii,1) = sqrt(min(Di_2)); 
                          % Find the distance from the point ii to its
                          % closest neighbor
        end        
    end        
 PCnG_old = PCnG;        % Store the normalized principal components of the 
                         % constructed grid (to check its convergence) 
            
end
    

VK(:,:,2) = vk_2d;           % Fill in the matrix of coefficients obtained 
                             % in the main iterative cycle into VK


time_CGA = cputime - time0;  % Compute time used to solve the model


% 14. Comparison between perturbation and EDS and accuracy tests
%---------------------------------------------------------------
T_test = 10200; % The simulation length for the accuracy test

% 14.1 Simulate a perturbation solution without ZLB imposed and compute the
% residuals
%--------------------------------------------------------------------------
zlb = 0;        % ZLB is not imposed

[PER1_test PER2_test]  = NK_PER12_simulation(SS,del2,A,B,C,D,E,epsi_test,sigma_nua,sigma_nuL,sigma_nuR,sigma_nuu,sigma_nuB,sigma_nuG,T_test,zlb);
                % Simulate PER1 and PER2 solutions

order = 1;      % Consider PER1
discard = 200;  % The number of observations to discard
[Residuals_mean(1) Residuals_max(1) Residuals_max_E(1:9,1)] = NK_PER12_accuracy(SS,del2,A,B,C,D,E,PER1_test,order,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,Gbar,piestar,vcv,discard,zlb);
                % Compute the mean and maximum residuals, "Residuals_mean" 
                % and "Residuals_max", across all points and all equilibrium
                %  conditions, as well as the maximum absolute  residuals, 
                % "Residuals_max_E", across all points, disaggregated  
                % by optimality conditions

order = 2;      % Consider PER2
discard = 200;  % The number of observations to discard
[Residuals_mean(2) Residuals_max(2) Residuals_max_E(1:9,2)] = NK_PER12_accuracy(SS,del2,A,B,C,D,E,PER2_test,order,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,Gbar,piestar,vcv,discard,zlb);
                % Compute the mean and maximum residuals, "Residuals_mean" 
                % and "Residuals_max", across all points and all equilibrium
                %  conditions, as well as the maximum absolute  residuals, 
                % "Residuals_max_E", across all points, disaggregated  
                % by optimality conditions

% 14.2 Simulate a perturbation solution with ZLB imposed and compute the
% residuals
%------------------------------------------------------------------------
zlb = 1;        % ZLB is imposed    
[PER1_test_zlb PER2_test_zlb]  = NK_PER12_simulation(SS,del2,A,B,C,D,E,epsi_test,sigma_nua,sigma_nuL,sigma_nuR,sigma_nuu,sigma_nuB,sigma_nuG,T_test,zlb);
                % Simulate PER1 and PER2 solutions
                
order = 1;      % Consider PER1
discard = 200;  % The number of observations to discard
[Residuals_mean(3) Residuals_max(3) Residuals_max_E(1:9,3)] = NK_PER12_accuracy(SS,del2,A,B,C,D,E,PER1_test_zlb,order,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,Gbar,piestar,vcv,discard,zlb);
                % Compute the mean and maximum residuals, "Residuals_mean" 
                % and "Residuals_max", across all points and all equilibrium
                %  conditions, as well as the maximum absolute  residuals, 
                % "Residuals_max_E", across all points, disaggregated  
                % by optimality conditions

order = 2;      % Consider PER2
discard = 200;  % The number of observations to discard
[Residuals_mean(4) Residuals_max(4) Residuals_max_E(1:9,4)] = NK_PER12_accuracy(SS,del2,A,B,C,D,E,PER2_test_zlb,order,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,Gbar,piestar,vcv,discard,zlb);
                % Compute the mean and maximum residuals, "Residuals_mean" 
                % and "Residuals_max", across all points and all equilibrium
                %  conditions, as well as the maximum absolute  residuals, 
                % "Residuals_max_E", across all points, disaggregated  
                % by optimality conditions
                

% 14.3 Simulate the EDS solution 
%-------------------------------
 
% Generate the series for shocks using PER1
% -----------------------------------------
nua_test = PER1_test(4,2:end)'; 
nuL_test = PER1_test(7,2:end)';
nuR_test = PER1_test(8,2:end)';
nuG_test = PER1_test(9,2:end)';
nuB_test = PER1_test(10,2:end)';
nuu_test = PER1_test(11,2:end)';

% Initial values for the endogenous state variables
%--------------------------------------------------
R_init  = PER1_test(5,1);    % Nominal interest rate in the initial period
delta_init  = PER1_test(6,1);% Price dispersion in the initial period


[S_test F_test delta_test C_test Y_test Yn_test L_test R_test pie_test] = NK_EDS_simulation(vk_2d,nuR_test,nua_test,nuL_test,nuu_test,nuB_test,nuG_test,R_init,delta_init,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,piestar,Gbar,zlb,Degree);%end
                % Simulate the EDS solution

discard = 200;  % The number of observations to discard
[Residuals_mean(5) Residuals_max(5) Residuals_max_E(1:9,5)] = NK_EDS_accuracy(nua_test,nuL_test,nuR_test,nuG_test,nuB_test,nuu_test,R_test,delta_test,L_test,Y_test,Yn_test,pie_test,S_test,F_test,C_test,rho_nua,rho_nuL,rho_nuR,rho_nuu,rho_nuB,rho_nuG,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,piestar,vcv,discard,vk_2d,Gbar,zlb,Degree);
                % Compute the mean and maximum residuals, "Residuals_mean" 
                % and "Residuals_max", across all points and all equilibrium
                %  conditions, as well as the maximum absolute  residuals, 
                % "Residuals_max_E", across all points, disaggregated  
                % by optimality conditions

Residuals_mean  % Display the mean residuals; 1-by-5
Residuals_max   % Display the maximum residuals; 1-by-5
Residuals_max_E % Display the maximum residuals by equilibrium conditions; 9-by-5
    
save pie1_D2;   % Save the results
CPU = toc       % Save the elapsed time in CPU