1% This mod file compares the functionality of Dynare's pruned_state_space.m with the 2% external Dynare pruning toolbox of Andreasen, Fernández-Villaverde and Rubio-Ramírez (2018): 3% "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", 4% Review of Economic Studies, Volume 85, Issue 1, Pages 1–49. 5% The model under study is taken from An and Schorfheide (2007): "Bayesian Analysis of DSGE Models", 6% Econometric Reviews, Volume 26, Issue 2-4, Pages 113-172. 7% Note that we use version 2 of the toolbox, i.e. the one which is called 8% "Third-order GMM estimate package for DSGE models (version 2)" and can be 9% downloaded from https://sites.google.com/site/mandreasendk/home-1 10% 11% Created by @wmutschl (Willi Mutschler, willi@mutschler.eu) 12% 13% ========================================================================= 14% Copyright (C) 2020 Dynare Team 15% 16% This file is part of Dynare. 17% 18% Dynare is free software: you can redistribute it and/or modify 19% it under the terms of the GNU General Public License as published by 20% the Free Software Foundation, either version 3 of the License, or 21% (at your option) any later version. 22% 23% Dynare is distributed in the hope that it will be useful, 24% but WITHOUT ANY WARRANTY; without even the implied warranty of 25% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 26% GNU General Public License for more details. 27% 28% You should have received a copy of the GNU General Public License 29% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 30% ========================================================================= 31 32% set this to 1 if you want to recompute using the Andreasen et al toolbox 33% otherwise the results are loaded from Andreasen_et_al_2018_Dynare44Pruning_v2.mat 34@#define Andreasen_et_al_toolbox = 0 35 36var YGR INFL INT 37 c p R g y z; %if ordering of var is changed comparison code below needs to be adapted 38varexo e_r e_g e_z; 39parameters tau nu kap cyst psi1 psi2 rhor rhog rhoz rrst pist gamst; 40 41tau = 2; 42nu = 0.1; 43kap = 0.33; 44cyst = 0.85; 45psi1 = 1.5; 46psi2 = 0.125; 47rhor = 0.75; 48rhog = 0.95; 49rhoz = 0.9; 50rrst = 1; 51pist = 3.2; 52gamst = 0.55; 53sig_r = .2; 54sig_g = .6; 55sig_z = .3; 56 57model; 58#pist2 = exp(pist/400); 59#rrst2 = exp(rrst/400); 60#bet = 1/rrst2; 61#phi = tau*(1-nu)/nu/kap/pist2^2; 62#gst = 1/cyst; 63#cst = (1-nu)^(1/tau); 64#yst = cst*gst; 65#dy = y-y(-1); 661 = exp(-tau*c(+1)+tau*c+R-z(+1)-p(+1)); 67(1-nu)/nu/phi/(pist2^2)*(exp(tau*c)-1) = (exp(p)-1)*((1-1/2/nu)*exp(p)+1/2/nu) - bet*(exp(p(+1))-1)*exp(-tau*c(+1)+tau*c+y(+1)-y+p(+1)); 68exp(c-y) = exp(-g) - phi*pist2^2*gst/2*(exp(p)-1)^2; 69R = rhor*R(-1) + (1-rhor)*psi1*p + (1-rhor)*psi2*(y-g) + e_r; 70g = rhog*g(-1) + e_g; 71z = rhoz*z(-1) + e_z; 72YGR = gamst+100*(dy+z); 73INFL = pist+400*p; 74INT = pist+rrst+4*gamst+400*R; 75end; 76 77shocks; 78var e_r = sig_r^2; 79var e_g = sig_g^2; 80var e_z = sig_z^2; 81end; 82 83steady_state_model; 84y = 0; 85R = 0; 86g = 0; 87z = 0; 88c = 0; 89p = 0; 90YGR = gamst; 91INFL = pist; 92INT = pist + rrst + 4*gamst; 93end; 94 95steady; check; model_diagnostics; 96 97@#for orderApp in [1, 2, 3] 98 stoch_simul(order=@{orderApp},pruning,irf=0,periods=0); 99 pruned_state_space.order_@{orderApp} = pruned_state_space_system(M_, options_, oo_.dr, [], options_.ar, 1, 0); 100 @#if Andreasen_et_al_toolbox 101 addpath('Dynare44Pruning_v2/simAndMoments3order'); %provide path to toolbox 102 optPruning.orderApp = @{orderApp}; 103 outAndreasenetal.order_@{orderApp} = RunDynarePruning(optPruning,oo_,M_,[oo_.dr.ghx oo_.dr.ghu]); 104 rmpath('Dynare44Pruning_v2/simAndMoments3order'); 105 close all; 106 @#endif 107@#endfor 108 109@#if Andreasen_et_al_toolbox 110 delete Andreasen_et_al_2018_Dynare44Pruning_v2.mat; 111 pause(3); 112 save('Andreasen_et_al_2018_Dynare44Pruning_v2.mat', 'outAndreasenetal') 113 pause(3); 114@#endif 115 116load('Andreasen_et_al_2018_Dynare44Pruning_v2.mat') 117 118% Make comparisons only at orders 1 and 2 119for iorder = 1:3 120 fprintf('ORDER %d:\n',iorder); 121 pruned = pruned_state_space.(sprintf('order_%d',iorder)); 122 outAndreasen = outAndreasenetal.(sprintf('order_%d',iorder)); 123 %make sure variable ordering is correct 124 if ~isequal(M_.endo_names,[outAndreasen.label_y; outAndreasen.label_v(1:M_.nspred)]) 125 error('variable ordering is not the same, change declaration order'); 126 end 127 norm_E_yx = norm(pruned.E_y(oo_.dr.inv_order_var) - [outAndreasen.Mean_y; outAndreasen.Mean_v(1:M_.nspred)] , Inf); 128 fprintf('max(sum(abs(E[y;x]''))): %d\n',norm_E_yx); 129 norm_Var_y = norm(pruned.Var_y(oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred)),oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred))) - outAndreasen.Var_y , Inf); 130 fprintf('max(sum(abs(Var[y]''))):: %d\n',norm_Var_y); 131 norm_Var_x = norm(pruned.Var_y(M_.nstatic+(1:M_.nspred),M_.nstatic+(1:M_.nspred)) - outAndreasen.Var_v(1:M_.nspred,1:M_.nspred) , Inf); 132 fprintf('max(sum(abs(Var[x]''))): %d\n',norm_Var_x); 133 norm_Corr_yi1 = norm(pruned.Corr_yi(oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred)),oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred)),1) - outAndreasen.Corr_y(:,:,1) , Inf); 134 fprintf('max(sum(abs(Corr[y,y(-1)]''))): %d\n',norm_Corr_yi1); 135 norm_Corr_yi2 = norm(pruned.Corr_yi(oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred)),oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred)),2) - outAndreasen.Corr_y(:,:,2) , Inf); 136 fprintf('max(sum(abs(Corr[y,y(-2)]''))): %d\n',norm_Corr_yi2); 137 norm_Corr_xi1 = norm(pruned.Corr_yi(M_.nstatic+(1:M_.nspred),M_.nstatic+(1:M_.nspred),1) - outAndreasen.Corr_v(1:M_.nspred,1:M_.nspred,1) , Inf); 138 fprintf('max(sum(abs(Corr[x,x(-1)]''))): %d\n',norm_Corr_xi1); 139 norm_Corr_xi2 = norm(pruned.Corr_yi(M_.nstatic+(1:M_.nspred),M_.nstatic+(1:M_.nspred),2) - outAndreasen.Corr_v(1:M_.nspred,1:M_.nspred,2) , Inf); 140 fprintf('max(sum(abs(Corr[x,x(-2)]''))): %d\n',norm_Corr_xi2); 141 142 if iorder < 3 && any([norm_E_yx norm_Var_y norm_Var_x norm_Corr_yi1 norm_Corr_yi2 norm_Corr_xi1 norm_Corr_xi2] > 1e-5) 143 error('Something wrong with pruned_state_space.m compared to Andreasen et al 2018 Toolbox v2 at order %d.',iorder); 144 end 145end 146skipline(); 147fprintf('Note that at third order, there is an error in the computation of Var_z in Andreasen et al (2018)''s toolbox, @wmutschl is in contact to clarify this.\n'); 148fprintf('EXAMPLE:\n') 149fprintf(' Consider Var[kron(kron(xf,xf),xf)] = E[kron(kron(kron(kron(kron(xf,xf),xf),xf),xf),xf)] - E[kron(kron(xf,xf),xf)]*E[kron(kron(xf,xf),xf)].''\n'); 150fprintf(' Now note that xf=hx*xf(-1)+hu*u is Gaussian, that is E[kron(kron(xf,xf),xf)]=0, and Var[kron(kron(xf,xf),xf)] are the sixth-order product moments\n'); 151fprintf(' which can be computed using the prodmom.m function by providing E[xf*xf''] as covariance matrix.\n'); 152fprintf(' In order to replicate this you have to change UnconditionalMoments_3rd_Lyap.m to also output Var_z.\n') 153 154dynare_nx = M_.nspred; 155dynare_E_xf2 = pruned_state_space.order_3.Var_z(1:dynare_nx,1:dynare_nx); 156dynare_E_xf6 = pruned_state_space.order_3.Var_z((end-dynare_nx^3+1):end,(end-dynare_nx^3+1):end); 157dynare_E_xf6 = dynare_E_xf6(:); 158 159Andreasen_nx = M_.nspred+M_.exo_nbr; 160Andreasen_E_xf2 = outAndreasenetal.order_3.Var_z(1:Andreasen_nx,1:Andreasen_nx); 161Andreasen_E_xf6 = outAndreasenetal.order_3.Var_z((end-Andreasen_nx^3+1):end,(end-Andreasen_nx^3+1):end); 162Andreasen_E_xf6 = Andreasen_E_xf6(:); 163 164fprintf('Second-order product moments of xf and u are the same:\n') 165norm_E_xf2 = norm(dynare_E_xf2-Andreasen_E_xf2(1:M_.nspred,1:M_.nspred),Inf) 166norm_E_uu = norm(M_.Sigma_e-Andreasen_E_xf2(M_.nspred+(1:M_.exo_nbr),M_.nspred+(1:M_.exo_nbr)),Inf) 167 168% Compute unique sixth-order product moments of xf, i.e. unique(E[kron(kron(kron(kron(kron(xf,xf),xf),xf),xf),xf)],'stable') 169dynare_nx6 = dynare_nx*(dynare_nx+1)/2*(dynare_nx+2)/3*(dynare_nx+3)/4*(dynare_nx+4)/5*(dynare_nx+5)/6; 170dynare_Q6Px = Q6_plication(dynare_nx); 171dynare_COMBOS6 = flipud(allVL1(dynare_nx, 6)); %all possible (unique) combinations of powers that sum up to six 172dynare_true_E_xf6 = zeros(dynare_nx6,1); %only unique entries 173for j6 = 1:size(dynare_COMBOS6,1) 174 dynare_true_E_xf6(j6) = prodmom(dynare_E_xf2, 1:dynare_nx, dynare_COMBOS6(j6,:)); 175end 176dynare_true_E_xf6 = dynare_Q6Px*dynare_true_E_xf6; %add duplicate entries 177norm_dynare_E_xf6 = norm(dynare_true_E_xf6 - dynare_E_xf6, Inf); 178 179Andreasen_nx6 = Andreasen_nx*(Andreasen_nx+1)/2*(Andreasen_nx+2)/3*(Andreasen_nx+3)/4*(Andreasen_nx+4)/5*(Andreasen_nx+5)/6; 180Andreasen_Q6Px = Q6_plication(Andreasen_nx); 181Andreasen_COMBOS6 = flipud(allVL1(Andreasen_nx, 6)); %all possible (unique) combinations of powers that sum up to six 182Andreasen_true_E_xf6 = zeros(Andreasen_nx6,1); %only unique entries 183for j6 = 1:size(Andreasen_COMBOS6,1) 184 Andreasen_true_E_xf6(j6) = prodmom(Andreasen_E_xf2, 1:Andreasen_nx, Andreasen_COMBOS6(j6,:)); 185end 186Andreasen_true_E_xf6 = Andreasen_Q6Px*Andreasen_true_E_xf6; %add duplicate entries 187norm_Andreasen_E_xf6 = norm(Andreasen_true_E_xf6 - Andreasen_E_xf6, Inf); 188 189fprintf('Sixth-order product moments of xf and u are not the same!\n'); 190fprintf(' Dynare maximum absolute deviations of sixth-order product moments of xf: %d\n',norm_dynare_E_xf6) 191fprintf(' Andreasen et al maximum absolute deviations of sixth-order product moments of xf: %d\n',norm_Andreasen_E_xf6) 192skipline(); 193fprintf('Note that the standard deviations are set quite high to make the numerical differences more apparent.\n'); 194