Solving a New Keynesian model with Matlab

This notebook is part of a computational appendix that accompanies the paper.

MATLAB, Python, Julia: What to Choose in Economics?

Coleman, Lyon, Maliar, and Maliar (2017)

In order to run the codes in this notebook you will need to install and configure Matlab and the a jupyter Matlab kernel. We assume you already have Matlab installed and show how to install the jupyter matlab kernel in a few steps:

  1. Install Python version 3.5. We recommend following the instructions on quantecon and using the Anaconda Python distribution. Make sure you don't get a python version above version 3.5.
  2. Install the MATLAB Engine API for Python. You can follow the official instructions from MathWorks at the link in the previous sentence.
  3. Install the imatlab kernel by doing the following from the command prompt (or terminal prompt for OSX/Linux userse):
    python -m pip install matlab_kernel
    python -m matlab_kernel install
    

Numerical Tools

In order to implement our routines, we need a few functions defining numerical tools.

Do to limitations in the Matlab Juptyer kernel, we cannot define functions in the notebook itself. So, we will download the file we need from online, save them to the same directory as this notebook, and then add this directory to the Matlab path.

So that we know what these files are doing, we will print the text of these files here in the notebook.

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"
     '

Accuracy and Simulation

Before proceeding to the main code, we will also need to have routines for simulating and checking the accuracy of our solution to the model.

We also download these from online as follows:

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
     '

Main routine

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