1+1 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')); mon1_txt mon2_txt poly_txt 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')); nkacc_txt nksim_txt % Start counting time % ------------------- tic; % 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) % 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); % 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 % 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 % 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 % 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) %-------------------------------------------------------------------------- % % 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)) % 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 % 10. Finish counting time % ------------------------ running_time = toc; % 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; % 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);