Solving a New Keynesian model with Matlab

In [1]:

    

1+1


    

    




ans =

     2

In [3]:

    

mon1_txt = fileread(websave('Monomials_1.m', 'https://s3.amazonaws.com/clmm-resources/mfiles/Monomials_1.m'));
mon2_txt = fileread(websave('Monomials_2.m', 'https://s3.amazonaws.com/clmm-resources/mfiles/Monomials_2.m'));
poly_txt = fileread(websave('Ord_Polynomial_N.m', 'https://s3.amazonaws.com/clmm-resources/mfiles/Ord_Polynomial_N.m'));


    

In [5]:

    

mon1_txt


    

    




mon1_txt =

    '% Monomials_1.m is a routine that constructs integration nodes and weights 
     % under N-dimensional monomial (non-product) integration rule with 2N nodes;
     % see Judd, Maliar and Maliar, (2010), "A Cluster-Grid Projection Method:
     % Solving Problems with High Dimensionality", NBER Working Paper 15965
     % (henceforth, JMM, 2010).
     % -------------------------------------------------------------------------
     % Inputs:  "N" is the number of random variables; N>=1;
     %          "vcv" is the variance-covariance matrix; N-by-N
     
     % Outputs: "n_nodes" is the total number of integration nodes; 2*N;
     %          "epsi_nodes" are the integration nodes; n_nodes-by-N;
     %          "weight_nodes" are the integration weights; n_nodes-by-1
     % -------------------------------------------------------------------------
     % 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 [n_nodes,epsi_nodes,weight_nodes] = Monomials_1(N,vcv)
     
     n_nodes   = 2*N;       % Total number of integration nodes
     
     % 1. N-dimensional integration nodes for N uncorrelated random variables with
     % zero mean and unit variance
     % ---------------------------------------------------------------------------
     z1 = zeros(n_nodes,N); % A supplementary matrix for integration nodes;
                            % n_nodes-by-N
     
     for i = 1:N            % In each node, random variable i takes value either
                            % 1 or -1, and all other variables take value 0
         z1(2*(i-1)+1:2*i,i) = [1; -1];
     end                    % For example, for N = 2, z1 = [1 0; -1 0; 0 1; 0 -1]
     
     % z = z1*sqrt(N);      % Integration nodes
     
     % 2. N-dimensional integration nodes and weights for N correlated random
     % variables with zero mean and variance-covaraince matrix vcv
     % ----------------------------------------------------------------------
     sqrt_vcv = chol(vcv);  % Cholesky decomposition of the variance-covariance
                            % matrix
     
     R = sqrt(N)*sqrt_vcv;  % Variable R; see condition (20) in JMM (2010)
     
     epsi_nodes = z1*R;     % Integration nodes; see condition (20) in JMM (2010);
                            % n_nodes-by-N
     
     % 3. Integration weights
     %-----------------------
     weight_nodes = ones(n_nodes,1)/n_nodes;
                            % Integration weights are equal for all integration
                            % nodes; n_nodes-by-1; the weights are the same for
                            % the cases of correlated and uncorrelated random
                            % variables
     '

In [6]:

    

mon2_txt


    

    




mon2_txt =

    '% Monomials_2.m is a routine that constructs integration nodes and weights 
     % under N-dimensional monomial (non-product) integration rule with 2N^2+1
     % nodes; see Judd, Maliar and Maliar, (2010), "A Cluster-Grid Projection
     % Method: Solving Problems with High Dimensionality", NBER Working Paper
     % 15965 (henceforth, JMM, 2010).
     % -------------------------------------------------------------------------
     % Inputs:  "N" is the number of random variables; N>=1;
     %          "vcv" is the variance-covariance matrix; N-by-N;
     
     % Outputs: "n_nodes" is the total number of integration nodes; 2*N^2+1;
     %          "epsi_nodes" are the integration nodes; n_nodes-by-N;
     %          "weight_nodes" are the integration weights; n_nodes-by-1
     % -------------------------------------------------------------------------
     % 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 [n_nodes,epsi_nodes,weight_nodes] = Monomials_2(N,vcv)
     
     n_nodes = 2*N^2+1;    % Total number of integration nodes
     
     % 1. N-dimensional integration nodes for N uncorrelated random variables with
     % zero mean and unit variance
     % ---------------------------------------------------------------------------
     
     % 1.1 The origin point
     % --------------------
     z0 = zeros(1,N);       % A supplementary matrix for integration nodes: the
                            % origin point
     
     % 1.2 Deviations in one dimension
     % -------------------------------
     z1 = zeros(2*N,N);     % A supplementary matrix for integration nodes;
                            % n_nodes-by-N
     
     for i = 1:N            % In each node, random variable i takes value either
                            % 1 or -1, and all other variables take value 0
         z1(2*(i-1)+1:2*i,i) = [1; -1];
     end                    % For example, for N = 2, z1 = [1 0; -1 0; 0 1; 0 -1]
     
     % 1.3 Deviations in two dimensions
     % --------------------------------
     z2 = zeros(2*N*(N-1),N);   % A supplementary matrix for integration nodes;
                                % 2N(N-1)-by-N
     
     i=0;                       % In each node, a pair of random variables (p,q)
                                % takes either values (1,1) or (1,-1) or (-1,1) or
                                % (-1,-1), and all other variables take value 0
     for p = 1:N-1
         for q = p+1:N
             i=i+1;
             z2(4*(i-1)+1:4*i,p) = [1;-1;1;-1];
             z2(4*(i-1)+1:4*i,q) = [1;1;-1;-1];
         end
     end                 % For example, for N = 2, z2 = [1 1;1 -1;-1 1;-1 1]
     
     % z = [z0;z1*sqrt(N+2);z2*sqrt((N+2)/2)];   % Integration nodes
     
     % 2. N-dimensional integration nodes and weights for N correlated random
     % variables with zero mean and variance-covaraince matrix vcv
     % ----------------------------------------------------------------------
     sqrt_vcv = chol(vcv);            % Cholesky decomposition of the variance-
                                      % covariance matrix
     
     R = sqrt(N+2)*sqrt_vcv;          % Variable R; see condition (21) in JMM
                                      % (2010)
     
     S = sqrt((N+2)/2)* sqrt_vcv;     % Variable S; see condition (21) in JMM
                                      % (2010)
     
     epsi_nodes = [z0;z1*R;z2*S];
                                      % Integration nodes; see condition (21)
                                      % in JMM (2010); n_nodes-by-N
     
     % 3. Integration weights
     %-----------------------
     weight_nodes = [2/(N+2)*ones(size(z0,1),1);(4-N)/2/(N+2)^2*ones(size(z1,1),1);1/(N+2)^2*ones(size(z2,1),1)];
                                      % See condition (21) in JMM (2010);
                                      % n_nodes-by-1; the weights are the same
                                      % for the cases of correlated and
                                      % uncorrelated random variables
     '

In [7]:

    

poly_txt


    

    




poly_txt =

    '% Ord_Polynomial_N.m is a routine that constructs the basis functions of 
     % complete ordinary polynomial of the degrees from one to five for the
     % multi-dimensional case; 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).
     %
     % This version: July 14, 2011. First version: August 27, 2009.
     % -------------------------------------------------------------------------
     % Inputs:  "z" is the data points on which the polynomial basis functions
     %               must be constructed; n_rows-by-dimen;
     %          "D" is the degree of the polynomial whose basis functions must
     %               be constructed; (can be 1,2,3,4 or 5)
     %
     % Output:  "basis_fs" is the matrix of basis functions of a complete
     %               polynomial of the given degree
     % -------------------------------------------------------------------------
     % 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.pdf".
     % -------------------------------------------------------------------------
     
     function basis_fs = Ord_Polynomial_N(z,D)
     
     
     % A polynomial is given by the sum of polynomial basis functions, phi(i),
     % multiplied by the coefficients; see condition (13) in JMM (2011). By
     % convention, the first basis function is one.
     
     [n_rows,dimen] = size(z); % Compute the number of rows, n_rows, and the
                               % number of variables (columns), dimen, in the
                               % data z on which the polynomial basis functions
                               % must be constructed
     
         % 1. The matrix of the basis functions of the first-degree polynomial
         % (the default option)
         % -------------------------------------------------------------------
         basis_fs = [ones(n_rows,1) z];  % The matrix includes a column of ones
                                         % (the first basis function is one for
                                         % n_rows points) and linear basis
                                         % functions
         i = dimen+1; % Index number of a polynomial basis function; the first
                      % basis function (equal to one) and linear basis functions
                      % are indexed from 1 to dimen+1, and subsequent polynomial
                      % basis functions will be indexed from dimen+2 and on
     
         % 2. The matrix of the basis functions of the second-degree polynomial
         % --------------------------------------------------------------------
         if D == 2
     
     % Version one (not vectorized):
             for j1 = 1:dimen
                 for j2 = j1:dimen
                     i = i+1;
                     basis_fs(:,i) = z(:,j1).*z(:,j2);
                 end
             end
     
     % Version 2 (vectorized): Note that this version works only for a second-degree
     % polynomial in which all state variables take non-zero values
     %        for r = 1:n_rows
     %            basis_fs(r,2+dimen:1+dimen+dimen*(dimen+1)/2) = [nonzeros(tril(z(r,:)'*z(r,:)))'];
                 % Compute linear and quadratic polynomial basis functions for
                 % each row r; "tril" extracts a lower triangular part of z'z
                 % (note that matrix z'z is symmetric so that an upper triangular
                 % part is redundant); "nonzeros" forms a column vector by
                 % stacking the columns of the original matrix one after another
                 % and by eliminating zero terms
      %       end
     
         % 3. The matrix of the basis functions of the third-degree polynomial
         % -------------------------------------------------------------------
         elseif D == 3
     
             for j1 = 1:dimen
                 for j2 = j1:dimen
                     i = i+1;
                     basis_fs(:,i) = z(:,j1).*z(:,j2);
                     for j3 = j2:dimen
                         i = i+1;
                         basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3);
                     end
                 end
             end
     
         % 4. The matrix of the basis functions of the fourth-degree polynomial
         % -------------------------------------------------------------------
         elseif D == 4
     
             for j1 = 1:dimen
                 for j2 = j1:dimen
                     i = i+1;
                     basis_fs(:,i) = z(:,j1).*z(:,j2);
                     for j3 = j2:dimen
                         i = i+1;
                         basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3);
                         for j4 = j3:dimen
                             i = i+1;
                             basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3).*z(:,j4);
                         end
                     end
                 end
             end
     
         % 5. The matrix of the basis functions of the fifth-degree polynomial
         % -------------------------------------------------------------------
     
         elseif D == 5
     
             for j1 = 1:dimen
                 for j2 = j1:dimen
                     i = i+1;
                     basis_fs(:,i) = z(:,j1).*z(:,j2);
                     for j3 = j2:dimen
                         i = i+1;
                         basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3);
                         for j4 = j3:dimen
                             i = i+1;
                             basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3).*z(:,j4);
                             for j5 = j4:dimen
                                 i = i+1;
                                 basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3).*z(:,j4).*z(:,j5);
                             end
                         end
                     end
                 end
             end
     
         end                         % end of "if"/"elseif"
     '

In [8]:

    

nkacc_txt = fileread(websave('NK_accuracy.m', 'https://s3.amazonaws.com/clmm-resources/mfiles/NK_accuracy.m'));
nksim_txt = fileread(websave('NK_simulation.m', 'https://s3.amazonaws.com/clmm-resources/mfiles/NK_simulation.m'));


    

In [9]:

    

nkacc_txt


    

    




nkacc_txt =

    '% NK_accuracy.m is a routine for evaluating accuracy of the solutions    
     % to the new Keynesian model considerd 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,
     % Quantitative Economics 6, 1-47 (2015) (henceforth, MM, 2015).
     % -------------------------------------------------------------------------
     % Inputs:    "nua", "nuL", "nuR", "nuG", "nuB", "nuu" are the time series
     %            of shocks generated with the innovations for the test;
     %            "delta", "L", "Y", "Yn", "pie", "S", "F", "C" are the time
     %            series solution generated with the innovations for test;
     %            "rho_nua", "rho_nuL", "rho_nuR", "rho_nuu", "rho_nuB", "rho_nuG"
     %            are the parameters of the laws of motion for shocks;
     %            "mu", "gam", "epsil", "vartheta", "beta", "A", "tau", "rho",
     %            "vcv", "beta", "phi_y", "phi_pie", "theta", "piestar" and "Gbar"
     %            are the parameters of the model;
     %            "vk_2d" is the matrix of coefficients of the GSSA solution;
     %            "discard" is the number of data points to discard;
     %            "Degree" is the degree of polynomial approximation
     %
     % Output:    "Residuals_mean" and "Residuals_max" are, respectively, the mean
     %            and maximum absolute residuals across all points and all equi-
     %            librium conditions; "Residuals_max_E" is the maximum absolute
     %            residuals across all points, disaggregated  by optimality conditions;
     %            "Residuals" are absolute residuals disaggregated by the equi-
     %            librium conditions
     % -------------------------------------------------------------------------
     % Copyright � 2013-2016 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 [Residuals_mean, Residuals_max, Residuals_max_E, Residuals] = NK_accuracy(nua,nuL,nuR,nuG,nuB,nuu,R,delta,L,Y,Yn,pie,S,F,C,rho_nua,rho_nuL,rho_nuR,rho_nuu,rho_nuB,rho_nuG,gam,vartheta,epsil,beta,phi_y,phi_pie,mu,theta,piestar,vcv,discard,vk_2d,Gbar,zlb,Degree)
     
     tic                     % Start counting time for running the test
     
     [T] = size(nua,1);      % Infer the number of points on which accuracy is
                             % evaluated
     Residuals = zeros(T,6); % Allocate memory to the matrix of residuals; T-by-6
     
     % Integration method for evaluating accuracy
     % ------------------------------------------
     [n_nodes,epsi_nodes,weight_nodes] = Monomials_2(6,vcv);
                                  % Monomial integration rule with 2N^2+1 nodes
     
     % Compute variables on the given set of points
     %---------------------------------------------
     
     for t = 1:T;                 % For each given point,
             t;
     
        % Take the corresponding value for shocks at t
        %---------------------------------------------
        nuR0 = nuR(t,1); % nuR(t)
        nua0 = nua(t,1); % nua(t)
        nuL0 = nuL(t,1); % nuL(t)
        nuu0 = nuu(t,1); % nuu(t)
        nuB0 = nuB(t,1); % nuB(t)
        nuG0 = nuG(t,1); % nuG(t)
     
        % Compute shocks at t+1 in all future nodes using their laws of motion
        %---------------------------------------------------------------------
        % Note that we do not premultiply by standard deviations as epsi_nodes
        % already include them
        nuR1(1:n_nodes,1) = (ones(n_nodes,1)*nuR0)*rho_nuR + epsi_nodes(:,1);
        % nuR(t+1); n_nodes-by-1
        nua1(1:n_nodes,1) = (ones(n_nodes,1)*nua0)*rho_nua + epsi_nodes(:,2);
        % nua(t+1); n_nodes-by-1
        nuL1(1:n_nodes,1) = (ones(n_nodes,1)*nuL0)*rho_nuL + epsi_nodes(:,3);
        % nuL(t+1); n_nodes-by-1
        nuu1(1:n_nodes,1) = (ones(n_nodes,1)*nuu0)*rho_nuu + epsi_nodes(:,4);
        % nuu(t+1); n_nodes-by-1
        nuB1(1:n_nodes,1) = (ones(n_nodes,1)*nuB0)*rho_nuB + epsi_nodes(:,5);
        % nuB(t+1); n_nodes-by-1
        nuG1(1:n_nodes,1) = (ones(n_nodes,1)*nuG0)*rho_nuG + epsi_nodes(:,6);
        % nuG(t+1); n_nodes-by-1
     
        R0  = R(t,1);            % R(t-1)
        delta0  = delta(t,1);    % delta(t-1)
        R1  = R(t+1,1);          % R(t)
        delta1  = delta(t+1,1);  % delta(t)
     
        L0 = L(t,1);             % L(t)
        Y0 = Y(t,1);             % Y(t)
        Yn0 = Yn(t,1);           % Yn(t)
        pie0 = pie(t,1);         % pie(t)
        S0 = S(t,1);             % S(t)
        F0 = F(t,1);             % F(t)
        C0 = C(t,1);             % C(t)
     
        % Future choices at t+1
        %----------------------
        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_nodes-by-1
     
        X1 = Ord_Polynomial_N([log(R1_dupl) log(delta1_dupl) nuR1 nua1 nuL1 nuu1 nuB1 nuG1],Degree);
        % Form a complete polynomial of degree "Degree" (at t+1) on future state
        % variables; n_nodes-by-npol_2d
     
        S1 = X1*vk_2d(:,1);             % Compute S(t+1) in all nodes using vk_2d
        F1 = X1*vk_2d(:,2);             % Compute F(t+1) in all nodes using vk_2d
        C1 = (X1*vk_2d(:,3)).^(-1/gam); % Compute C(t+1) in all nodes using vk_2d
        pie1 = ((1-(1-theta)*(S1./F1).^(1-epsil))/theta).^(1/(epsil-1));
                                        % Compute pie(t+1) using condition (35)
                                        % in MM (2015)
     
         % Compute residuals for each of the 9 equilibrium conditions
         %-----------------------------------------------------------
         Residuals(t,1) = 1-weight_nodes'*(exp(nuu0)*exp(nuL0)*L0^vartheta*Y0/exp(nua0) + beta*theta*pie1.^epsil.*S1)/S0;
         Residuals(t,2) = 1-weight_nodes'*(exp(nuu0)*C0^(-gam)*Y0 + beta*theta*pie1.^(epsil-1).*F1)/F0;
         Residuals(t,3) = 1-weight_nodes'*(beta*exp(nuB0)/exp(nuu0)*R1*exp(nuu1).*C1.^(-gam)./pie1)/C0^(-gam);
         Residuals(t,4) = 1-((1-theta*pie0^(epsil-1))/(1-theta))^(1/(1-epsil))*F0/S0;
         Residuals(t,5) = 1-((1-theta)*((1-theta*pie0^(epsil-1))/(1-theta))^(epsil/(epsil-1)) + theta*pie0^epsil/delta0)^(-1)/delta1;
         Residuals(t,6) = 1-exp(nua0)*L0*delta1/Y0;
         Residuals(t,7) = 1-(1-Gbar/exp(nuG0))*Y0/C0;
         Residuals(t,8) = 1-(exp(nua0)^(1+vartheta)*(1-Gbar/exp(nuG0))^(-gam)/exp(nuL0))^(1/(vartheta+gam))/Yn0;
         Residuals(t,9) = 1-piestar/beta*(R0*beta/piestar)^mu*((pie0/piestar)^phi_pie * (Y0/Yn0)^phi_y)^(1-mu)*exp(nuR0)/R1;   % Taylor rule
         if (zlb==1); Residuals(t,9) = Residuals(t,9)*(R1>1);end
         % If the ZLB is imposed and R>1, the residuals in the Taylor rule (the
         % 9th equation) are zero
     
     end
     
     % Residuals across all the equilibrium conditions and all test points
     %--------------------------------------------------------------------
        Residuals_mean = log10(mean(mean(abs(Residuals(1+discard:end,:)))));
        % Mean absolute residuals computed after discarding the first
        % "discard" observations
     
        Residuals_max = log10(max(max(abs(Residuals(1+discard:end,:)))));
        % Maximum absolute residuals computed after discarding the first
        % "discard" observations
     
     % Residuals disaggregated by the eqiulibrium conditions
     %------------------------------------------------------
      Residuals_max_E = log10(max(abs(Residuals(1+discard:end,:))))';
        % Maximum absolute residuals across all test points for each of the 9
        % equilibrium conditions computed after discarding the first "discard"
        % observations;
     
     time_test = toc;     % Time needed to run the test
     '

In [10]:

    

nksim_txt


    

    




nksim_txt =

    '% NK_simulation.m is a routine for simulating the solution to the new 
     % Keynesian model considerd 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, Quantitative Economics
     % 6, 1-47 (2015) (henceforth, MM, 2015).
     % -------------------------------------------------------------------------
     
     % -------------------------------------------------------------------------
     % Inputs:    "vk" is the matrix of coefficients of the GSSA solution;
     %            "nua", "nuL", "nuR", "nuG", "nuB", "nuu" are the time series
     %            of shocks generated with the innovations for the test;
     %            "R_init" and "delta_init" are initial values for R and delta
     %            "gam", "vartheta", "epsil", "betta", "phi_y", "phi_pie", "mu",
     %            "theta", "piestar", "Gbar"  are the parameters of the model;
     %            "zlb" is a dummy parameter which is equal to 0 when ZLB is not
     %            imposed and is equal to 1 when it is imposed;
     %            "Degree" is the degree of polynomial approximation
     %
     % Output:    "S", "F", "delta", "C", "Y", "Yn", "L", "R" and "pie2" are the
     %             simulated series
     % -------------------------------------------------------------------------
     % Copyright � 2013-2016 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 [S F delta C Y Yn L R pie w] = NK_simulation(vk,nuR,nua,nuL,nuu,nuB,nuG,R_init,delta_init,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,piestar,Gbar,zlb,Degree)
     
     
     [T] = size(nua,1);       % Infer the number of points on which the accuracy
                              % is evaluated
     
     delta = ones(T+1,1);     % Allocate memory for the time series of delta(t)
     R = ones(T+1,1);         % Allocate memory for the time series of R(t)
     S = ones(T,1);           % Allocate memory for the time series of S(t)
     F = ones(T,1);           % Allocate memory for the time series of F(t)
     C = ones(T,1);           % Allocate memory for the time series of C(t)
     pie = ones(T,1);         % Allocate memory for the time series of pie(t)
     Y = ones(T,1);           % Allocate memory for the time series of Y(t)
     L = ones(T,1);           % Allocate memory for the time series of L(t)
     Yn = ones(T,1);          % Allocate memory for the time series of Yn(t)
     w  = ones(T,1);
     state = ones(T,8);       % Allocate memory for the time series of the state
                              % variables; T-by-8
     
     delta(1,1) = delta_init; % Initial condition for delta(t-1), i.e., delta(-1)
     R(1,1) = R_init;         % Initial condition for R(t-1)
     
     for t = 1:T
         pol_bases = Ord_Polynomial_N([log(R(t,1)) log(delta(t,1)) nuR(t,1) nua(t,1) nuL(t,1) nuu(t,1) nuB(t,1) nuG(t,1)],Degree);
         % Construct the matrix of explanatory variables "pol_bases" on the series
         % of state variables; columns of "pol_bases" are given by the basis
         % functions of the polynomial of degree "Degree"
         S(t,1) = pol_bases*vk(:,1);             % Compute S(t) using vk
         F(t,1) = pol_bases*vk(:,2);             % Compute F(t) using vk
         C(t,1) = (pol_bases*vk(:,3)).^(-1/gam); % Compute C(t) using vk
         pie(t,:) = ((1-(1-theta)*(S(t,:)/F(t,:))^(1-epsil))/theta)^(1/(epsil-1));
                              % Compute pie(t) from condition (35) in MM (2015)
         delta(t+1,:) = ((1-theta)*((1-theta*pie(t,1)^(epsil-1))/(1-theta))^(epsil/(epsil-1))+theta*pie(t,1)^epsil/delta(t,1))^-1;
                               % Compute delta(t) from condition (36) in MM (2015)
         Y(t,:) = C(t,1)/(1-Gbar/exp(nuG(t,1)));
                              % Compute Y(t) from condition (38) in MM (2015)
         L(t,1) = Y(t,1)/delta(t+1,1)/exp(nua(t,1));
                              % Compute L(t) from condition (37) in MM (2015)
         Yn(t,:) = (exp(nua(t,1))^(1+vartheta)*(1-Gbar/exp(nuG(t,1)))^(-gam)/exp(nuL(t,1)))^(1/(vartheta+gam));
                              %  Compute Yn(t) from condition (31) in MM (2015)
         R(t+1,:) = piestar/betta*(R(t,:)*betta/piestar)^mu*((pie(t,:)/piestar)^phi_pie*(Y(t,:)/Yn(t,:))^phi_y)^(1-mu)*exp(nuR(t,1));
              % Compute R(t) from conditions (27), (39) in MM (2015)
         w(t,1) = exp(nuL(t,1))*(L(t,1)^vartheta)*(C(t,1)^gam);
              % Compute real wage
     
         if zlb == 1
             R(t+1,:) = max(R(t+1,:),1);
             % If ZLB is imposed, set R(t)=1 if ZLB binds
         end;
     
     end
     '

In [24]:

    

% Start counting time
% -------------------
tic;


    

In [25]:

    

% 1. Parameter values
% -------------------
zlb        = 0;           % Impose ZLB on nominal interest rate
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.0;         % 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.0;         % See process (28) in MM (2015)
rho_nuu    = 0.92;        % See process (15) in MM (2015)
rho_nuB    = 0.0;         % 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.4054;      % 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.0010;      % See process (17) in MM (2015)
sigma_nuG  = 0.0038;      % See process (26) in MM (2015)


    

In [26]:

    

% 2. Steady state of the model (see page 19 in Supplement to MM (2015))
% -------------------------------------------------------------------
Yn_ss     = exp(Gbar)^(gam/(vartheta+gam));
Y_ss      = Yn_ss;
pie_ss    = 1;
delta_ss  = 1;
L_ss      = Y_ss/delta_ss;
C_ss      = (1-Gbar)*Y_ss;
F_ss      = C_ss^(-gam)*Y_ss/(1-betta*theta*pie_ss^(epsil-1));
S_ss      = L_ss^vartheta*Y_ss/(1-betta*theta*pie_ss^epsil);
R_ss      = pie_ss/betta;
w_ss      = (L_ss^vartheta)*(C_ss^gam);


    

In [27]:

    

% 3. Construct a grid for computing a solution
% --------------------------------------------
m = 200;        % Choose the number of grid points
grid_type = 1;  % Choose a grid type; grid_type = 1 corresponds to a uniformely
                % distribued random grid, and  grid_type = 2 corresponds to
                % a quasi Monte Carlo grid

% Random grid
% -----------
if grid_type == 1 % Choose the grid type
nuR0    = (-2*sigma_nuR+4*sigma_nuR*rand(m,1))/sqrt(1-rho_nuR^2);
nua0    = (-2*sigma_nua+4*sigma_nua*rand(m,1))/sqrt(1-rho_nua^2);
nuL0    = (-2*sigma_nuL+4*sigma_nuL*rand(m,1))/sqrt(1-rho_nuL^2);
nuu0    = (-2*sigma_nuu+4*sigma_nuu*rand(m,1))/sqrt(1-rho_nuu^2);
nuB0    = (-2*sigma_nuB+4*sigma_nuB*rand(m,1))/sqrt(1-rho_nuB^2);
nuG0    = (-2*sigma_nuG+4*sigma_nuG*rand(m,1))/sqrt(1-rho_nuG^2);
    % Values of exogenous state variables are distributed uniformly
    % in the interval +/- std/sqrt(1-rho_nu^2)

R0      = 1+0.05*rand(m,1);
delta0  = 0.95+0.05*rand(m,1);
    % Values of endogenous state variables are distributed uniformly
    % in the intervals [1 1.05] and [0.95 1], respectively
end

% Quasi Monte-Carlo grid
%-----------------------
if grid_type == 2  % Choose the grid type
dimensionality = 8; % The number of state variables (exogenous and endogenous)
Def_sobol = sobolset(dimensionality);
                    % Constructs a Sobol sequence point set in
                    % "dimensionality" dimensions

Sob = net(Def_sobol,m);
                    % Get the first m points

nuR0    = (-2*sigma_nuR+4*(max(Sob(:,1))-Sob(:,1))/(max(Sob(:,1))-min(Sob(:,1)))*sigma_nuR)/sqrt(1-rho_nuR^2);
nua0    = (-2*sigma_nua+4*(max(Sob(:,2))-Sob(:,2))/(max(Sob(:,2))-min(Sob(:,2)))*sigma_nua)/sqrt(1-rho_nua^2);
nuL0    = (-2*sigma_nuL+4*(max(Sob(:,3))-Sob(:,3))/(max(Sob(:,3))-min(Sob(:,3)))*sigma_nuL)/sqrt(1-rho_nuL^2);
nuu0    = (-2*sigma_nuu+4*(max(Sob(:,4))-Sob(:,4))/(max(Sob(:,4))-min(Sob(:,4)))*sigma_nuu)/sqrt(1-rho_nuu^2);
nuB0    = (-2*sigma_nuB+4*(max(Sob(:,5))-Sob(:,5))/(max(Sob(:,5))-min(Sob(:,5)))*sigma_nuB)/sqrt(1-rho_nuB^2);
nuG0    = (-2*sigma_nuG+4*(max(Sob(:,6))-Sob(:,6))/(max(Sob(:,6))-min(Sob(:,6)))*sigma_nuG)/sqrt(1-rho_nuG^2);
    % Values of exogenous state variables are in the interval +/- std/sqrt(1-rho^2)

R0      = 1+0.05*(max(Sob(:,7))-Sob(:,7))/(max(Sob(:,7))-min(Sob(:,7)));
delta0  = 0.95+0.05*(max(Sob(:,8))-Sob(:,8))/(max(Sob(:,8))-min(Sob(:,8)));
    % Values of endogenous state variables are in the intervals [1 1.05] and
    % [0.95 1], respectively
end

if zlb == 1; R0=max(R0,1); end
             % If ZLB is imposed, set R(t)=1 if ZLB binds

Grid    = [log(R0(1:m,1)) log(delta0(1:m,1)) nuR0 nua0 nuL0 nuu0 nuB0 nuG0];
    % Construct the matrix of grid points; m-by-dimensionality


    

In [28]:

    

% 4. Constructing polynomial on the grid
% --------------------------------------
Degree  = 2;         % Degree of polynomial approximation
X0_Gs{1} = Ord_Polynomial_N(Grid, 1);
X0_Gs{Degree} = Ord_Polynomial_N(Grid, Degree);
                     % Construct the matrix of explanatory variables X0_Gs
                     % on the grid of state variables; the columns of X0_Gs
                     % are given by the basis functions of polynomial of
                     % degree "Degree"
npol = size(X0_Gs{Degree},2); % Number of coefficients in polynomial of degree
                              % "Degree"; it must be smaller than the number of grid
                              % points


    

In [29]:

    

% 5. Integration in the solution procedure
% ----------------------------------------
N = 6;               % Total number of exogenous shocks
vcv = diag([sigma_nuR^2 sigma_nua^2 sigma_nuL^2 sigma_nuu^2 sigma_nuB^2 sigma_nuG^2]);
                     % Variance covariance matrix

% Compute the number of integration nodes, their values and weights
%------------------------------------------------------------------
[n_nodes,epsi_nodes,weight_nodes] = Monomials_1(N,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

nuR1(:,1:n_nodes) = (nuR0*ones(1,n_nodes)).*rho_nuR + ones(m,1)*epsi_nodes(:,1)';
nua1(:,1:n_nodes) = (nua0*ones(1,n_nodes)).*rho_nua + ones(m,1)*epsi_nodes(:,2)';
nuL1(:,1:n_nodes) = (nuL0*ones(1,n_nodes)).*rho_nuL + ones(m,1)*epsi_nodes(:,3)';
nuu1(:,1:n_nodes) = (nuu0*ones(1,n_nodes)).*rho_nuu + ones(m,1)*epsi_nodes(:,4)';
nuB1(:,1:n_nodes) = (nuB0*ones(1,n_nodes)).*rho_nuB + ones(m,1)*epsi_nodes(:,5)';
nuG1(:,1:n_nodes) = (nuG0*ones(1,n_nodes)).*rho_nuG + ones(m,1)*epsi_nodes(:,6)';
            % Compute future shocks in all grid points and all integration
            % nodes; the size of each of these matrices is m-by-n_nodes


    

In [30]:

    

% 6. Allocate memory
%--------------------
% TODO: pick up here!!!!
e = zeros(m,3);
% Allocate memory to integrals in the right side of 3 Euler equations

S0_old_G = ones(m,1);
F0_old_G = ones(m,1);
C0_old_G = ones(m,1);
% Allocate memory to S, F and C from the previous iteration (to check
% convergence)

S0_new_G = ones(m,1);
F0_new_G = ones(m,1);
C0_new_G = ones(m,1);
% Allocate memory to S, F, C from the current iteration (to check
% convergence)


    

In [31]:

    

%--------------------------------------------------------------------------
%
% The main iterative cycle
%
% -------------------------------------------------------------------------

% 7. Algorithm parameters
%------------------------
damp     = 0.1;           % Damping parameter for (fixed-point) iteration on
                          % the coefficients of 3 decision functions (for
                          % S, F and C^(-gam))


    

In [32]:

    

% 8. The loop over the polynomial coefficients
% --------------------------------------------
for deg = [1, Degree]
    diff = 1e+10;
    it = 0;
    X0_G = X0_Gs{deg};

    % 9. Initial guess for coefficients of the decision functions for the
    % variables S and F and marginal utility MU
    % -------------------------------------------------------------------
    if deg == 1
         vk       =  ones(size(Grid, 2)+1, 3)*1e-5;  % Initialize first all the coefficients
                                                     % at 1e-5
         vk(1,:)  = [S_ss F_ss C_ss.^(-gam)];        % Set the initial values of the constant
                                                     % terms in the decision rules for S,
                                                     % F and MU to values that give the
                                                     % deterministic steady state
    else
        % For degree > 1, initial guess for coefficients is given by
        % regressing final state matrix from degree 1 solution (e) on the
        % complete polynomial basis matrix
        vk = X0_G\e;
    end

    while diff > 1e-7        % The convergence criterion (which is unit free
                              % because diff is unit free)

        % Current choices (at t)
        % ------------------------------
        S0 = X0_G*vk(:,1);              % Compute S(t) using vk
        F0 = X0_G*vk(:,2);              % Compute F(t) using vk
        C0 = (X0_G*vk(:,3)).^(-1/gam);  % Compute C(t) using vk

        pie0 = ((1-(1-theta)*(S0./F0).^(1-epsil))/theta).^(1/(epsil-1));
                   % Compute pie(t) 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(t) from condition (36) in MM (2015)
        Y0 = C0./(1-Gbar./exp(nuG0));
                   % Compute Y(t) from condition (38) in MM (2015)
        L0 = Y0./exp(nua0)./delta1;
                   % Compute L(t) from condition (37) in MM (2015)
        Yn0 = (exp(nua0).^(1+vartheta).*(1-Gbar./exp(nuG0)).^(-gam)./exp(nuL0)).^(1/(vartheta+gam));
                   %  Compute Yn(t) 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);    % Taylor rule
                   % Compute R(t) from conditions (27), (39) in MM (2015)
        if zlb == 1; R1=max(R1,1); end
                   % If ZLB is imposed, set R(t)=1 if ZLB binds

        % Future choices (at t+1)
        %--------------------------------
        for u = 1:n_nodes

            X1 = Ord_Polynomial_N([log(R1) log(delta1) nuR1(:,u) nua1(:,u) nuL1(:,u) nuu1(:,u) nuB1(:,u) nuG1(:,u)],deg);
            % Form complete polynomial of degree "Degree" (at t+1) on future state
            % variables; n_nodes-by-npol

            S1(:,u) = X1*vk(:,1);             % Compute S(t+1) in all nodes using vk
            F1(:,u) = X1*vk(:,2);             % Compute F(t+1) in all nodes using vk
            C1(:,u) = (X1*vk(:,3)).^(-1/gam); % Compute C(t+1) in all nodes using vk

        end

        pie1 = ((1-(1-theta)*(S1./F1).^(1-epsil))/theta).^(1/(epsil-1));
                                          % Compute next-period pie using condition
                                          % (35) in MM (2015)


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


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

       % Compute and update the coefficients of the decision functions
       % -------------------------------------------------------------
       vk_hat_2d = X0_G\e;   % Compute the new coefficients of the decision
                             % functions using a backslash operator

       vk = damp*vk_hat_2d + (1-damp)*vk;
                             % Update the coefficients using damping

       % Evaluate the percentage (unit-free) difference between the values
       % on the grid from the previous and current iterations
       % -----------------------------------------------------------------
       diff = mean(mean(abs(1-S0_new_G./S0_old_G)))+mean(mean(abs(1-F0_new_G./F0_old_G)))+mean(mean(abs(1-C0_new_G./C0_old_G)));
                       % The convergence criterion is adjusted to the damping
                       % parameters

       % Store the obtained values for S(t), F(t), C(t) on the grid to
       % be used on the subsequent iteration in Section 10.2.6
       %-----------------------------------------------------------------------
       S0_old_G = S0_new_G;
       F0_old_G = F0_new_G;
       C0_old_G = C0_new_G;       
    end
end


    

In [33]:

    

% 10. Finish counting time
% ------------------------
running_time  = toc;


    

In [34]:

    

% 11. Simualating a time-series solution
%---------------------------------------
T = 10201; % The length of stochastic simulation


% Initialize the values of 6 exogenous shocks and draw innovations
%-----------------------------------------------------------------
nuR = zeros(T,1); eps_nuR = randn(T,1)*sigma_nuR;
nua = zeros(T,1); eps_nua = randn(T,1)*sigma_nua;
nuL = zeros(T,1); eps_nuL = randn(T,1)*sigma_nuL;
nuu = zeros(T,1); eps_nuu = randn(T,1)*sigma_nuu;
nuB = zeros(T,1); eps_nuB = randn(T,1)*sigma_nuB;
nuG = zeros(T,1); eps_nuG = randn(T,1)*sigma_nuG;

% Generate the series for shocks
%-------------------------------
for t = 1:T-1
    nuR(t+1,1) = rho_nuR*nuR(t,1) + eps_nuR(t);
    nua(t+1,1) = rho_nua*nua(t,1) + eps_nua(t);
    nuL(t+1,1) = rho_nuL*nuL(t,1) + eps_nuL(t);
    nuu(t+1,1) = rho_nuu*nuu(t,1) + eps_nuu(t);
    nuB(t+1,1) = rho_nuB*nuB(t,1) + eps_nuB(t);
    nuG(t+1,1) = rho_nuG*nuG(t,1) + eps_nuG(t)  ;
end

% Initial values of two endogenous state variables
%-------------------------------------------------
R_initial      = 1; % Nominal interest rate R
delta_initial  = 1; % Price dispersion "delta"

tic;
% Simulate the model
%-------------------
[S, F, delta, C, Y, Yn, L, R, pie, w] = NK_simulation(vk,nuR,nua,nuL,nuu,nuB,nuG,R_initial,delta_initial,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,piestar,Gbar,zlb,Degree);

simulation_time = toc;


    

In [35]:

    

% 12. Compute unit free Euler equation residuals on simulated points
%-------------------------------------------------------------------
discard = 200; % The number of observations to discard
time0 = tic;

[Residuals_mean(1), Residuals_max(1), Residuals_max_E(1:9,1), Residuals] = NK_accuracy(nua,nuL,nuR,nuG,nuB,nuu,R,delta,L,Y,Yn,pie,S,F,C,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,Gbar,zlb,Degree);

accuracy_time = toc;
% Compute the residuals; Residuals_mean, Residuals_max and Residuals_max_E(1:9,5)
% are filled into the 5h columns of the corresponding vectors/matrix, while
% the first 4 columns correspond to PER1, PER2 with and without ZLB

% ------------------------------------------------------------------------
% disp(' '); disp('           OUTPUT:'); disp(' ');
% disp('RUNNING TIME (in seconds):'); disp('');
% display(running_time);
% disp('SIMULATION TIME (in seconds):'); disp('');
% display(simulation_time);
% disp('ACCURACY TIME (in seconds):'); disp('');
% display(accuracy_time);
% disp('APPROXIMATION ERRORS (log10):'); disp('');
% disp('a) mean error in the model equations');
% disp(Residuals_mean)
% disp('b) max error in the model equations');
% disp(Residuals_max)
% disp('b) max error in by equation');
% disp(Residuals_max_E(1:9,1))

l1 = log10(sum(max(abs(Residuals), [], 1)));
tot_time = running_time + simulation_time + accuracy_time;
fprintf('Solver time: %4.6f seconds\n', running_time);
fprintf('Simulation time: %4.6f seconds\n', simulation_time);
fprintf('Accuracy time: %4.6f seconds\n', accuracy_time);
fprintf('Total time: %4.6f seconds\n', running_time + simulation_time + accuracy_time)
fprintf('pistar: %3.4f    sigma_L: %3.4f\n', piestar, sigma_nuL);
fprintf('zlb : %d    grid: %d\n', zlb, grid_type);
fprintf('tex line: %0.2f & %0.2f & %0.2f\n', l1, Residuals_max, tot_time);


    

    




Solver time: 3.988508 seconds
Simulation time: 0.194495 seconds
Accuracy time: 0.922081 seconds
Total time: 5.105085 seconds
pistar: 1.0000    sigma_L: 0.4054
zlb : 0    grid: 1
tex line: -1.08 & -1.27 & 5.11

In [ ]: