stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'pruned_state_space' into 'master' 

Speed up pruned_state_space_system.m by e.g. using persistent variables

See merge request Dynare/dynare!1749
time-shift
Sébastien Villemot 2020-07-16 09:39:01 +00:00
parent 77dbeab98a 8a00ee3dff
commit 499451d50a
4 changed files with 83 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
matlab/Q6_plication.m
View File

@ -63,7 +63,9 @@ for i1=1:p
end end
end end
end end
if nargout==2
DP6inv = (transpose(DP6)*DP6)\transpose(DP6); DP6inv = (transpose(DP6)*DP6)\transpose(DP6);
end
function m = mue(p,i1,i2,i3,i4,i5,i6) function m = mue(p,i1,i2,i3,i4,i5,i6)
% Auxiliary expression, see page 122 of Meijer (2005) % Auxiliary expression, see page 122 of Meijer (2005)

2
matlab/list_of_functions_to_be_cleared.m
View File

@ -1 +1 @@
list_of_functions = {'discretionary_policy_1', 'dsge_var_likelihood', 'dyn_first_order_solver', 'dyn_waitbar', 'ep_residuals', 'evaluate_likelihood', 'prior_draw_gsa', 'identification_analysis', 'computeDLIK', 'univariate_computeDLIK', 'metropolis_draw', 'flag_implicit_skip_nan', 'moment_function', 'mr_hessian', 'masterParallel', 'auxiliary_initialization', 'auxiliary_particle_filter', 'conditional_filter_proposal', 'conditional_particle_filter', 'gaussian_filter', 'gaussian_filter_bank', 'gaussian_mixture_filter', 'gaussian_mixture_filter_bank', 'Kalman_filter', 'online_auxiliary_filter', 'sequential_importance_particle_filter', 'solve_model_for_online_filter', 'perfect_foresight_simulation', 'prior_draw', 'priordens'}; list_of_functions = {'discretionary_policy_1', 'dsge_var_likelihood', 'dyn_first_order_solver', 'dyn_waitbar', 'ep_residuals', 'evaluate_likelihood', 'prior_draw_gsa', 'identification_analysis', 'computeDLIK', 'univariate_computeDLIK', 'metropolis_draw', 'flag_implicit_skip_nan', 'moment_function', 'mr_hessian', 'masterParallel', 'auxiliary_initialization', 'auxiliary_particle_filter', 'conditional_filter_proposal', 'conditional_particle_filter', 'gaussian_filter', 'gaussian_filter_bank', 'gaussian_mixture_filter', 'gaussian_mixture_filter_bank', 'Kalman_filter', 'online_auxiliary_filter', 'pruned_state_space_system', 'sequential_importance_particle_filter', 'solve_model_for_online_filter', 'perfect_foresight_simulation', 'prior_draw', 'priordens'};

56
matlab/pruned_state_space_system.m
View File

@ -10,12 +10,13 @@ function pruned_state_space = pruned_state_space_system(M, options, dr, indy, nl
% Econometrics and Statistics, Volume 6, Pages 44-56. % Econometrics and Statistics, Volume 6, Pages 44-56.
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
% M: [structure] storing the model information % M: [structure] storing the model information
% options: [structure] storing the options % options: [structure] storing the options
% dr: [structure] storing the results from perturbation approximation % dr: [structure] storing the results from perturbation approximation
% indy: [vector] index of control variables in DR order % indy: [vector] index of control variables in DR order
% nlags: [integer] number of lags in autocovariances and autocorrelations % nlags: [integer] number of lags in autocovariances and autocorrelations
% useautocorr: [boolean] true: compute autocorrelations % useautocorr: [boolean] true: compute autocorrelations
% compute_derivs: [boolean] true: compute derivatives
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% OUTPUTS % OUTPUTS
% pruned_state_space: [structure] with the following fields: % pruned_state_space: [structure] with the following fields:
@ -240,6 +241,7 @@ function pruned_state_space = pruned_state_space_system(M, options, dr, indy, nl
% See code below how z and inov are defined at first, second, and third order, % See code below how z and inov are defined at first, second, and third order,
% and how to set up A, B, C, D and compute unconditional first and second moments of inov, z and y % and how to set up A, B, C, D and compute unconditional first and second moments of inov, z and y
persistent QPu COMBOS4 Q6Pu COMBOS6 K_u_xx K_u_ux K_xx_x
%% Auxiliary indices and objects %% Auxiliary indices and objects
order = options.order; order = options.order;
@ -418,8 +420,10 @@ if order > 1
%Compute unique fourth order product moments of u, i.e. unique(E[kron(kron(kron(u,u),u),u)],'stable') %Compute unique fourth order product moments of u, i.e. unique(E[kron(kron(kron(u,u),u),u)],'stable')
u_nbr4 = u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4; u_nbr4 = u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4;
QPu = quadruplication(u_nbr); if isempty(QPu)
COMBOS4 = flipud(allVL1(u_nbr, 4)); %all possible (unique) combinations of powers that sum up to four QPu = quadruplication(u_nbr);
COMBOS4 = flipud(allVL1(u_nbr, 4)); %all possible (unique) combinations of powers that sum up to four
end
E_u_u_u_u = zeros(u_nbr4,1); %only unique entries E_u_u_u_u = zeros(u_nbr4,1); %only unique entries
if compute_derivs && (stderrparam_nbr+corrparam_nbr>0) if compute_derivs && (stderrparam_nbr+corrparam_nbr>0)
dE_u_u_u_u = zeros(u_nbr4,stderrparam_nbr+corrparam_nbr); dE_u_u_u_u = zeros(u_nbr4,stderrparam_nbr+corrparam_nbr);
@ -588,8 +592,11 @@ if order > 1
if order > 2 if order > 2
% Some common and useful objects for order > 2 % Some common and useful objects for order > 2
K_u_xx = commutation(u_nbr,x_nbr^2,1); if isempty(K_u_xx)
K_u_ux = commutation(u_nbr,u_nbr*x_nbr,1); K_u_xx = commutation(u_nbr,x_nbr^2,1);
K_u_ux = commutation(u_nbr,u_nbr*x_nbr,1);
K_xx_x = commutation(x_nbr^2,x_nbr);
end
hx_hss2 = kron(hx,1/2*hss); hx_hss2 = kron(hx,1/2*hss);
hu_hss2 = kron(hu,1/2*hss); hu_hss2 = kron(hu,1/2*hss);
hx_hxx2 = kron(hx,1/2*hxx); hx_hxx2 = kron(hx,1/2*hxx);
@ -658,9 +665,11 @@ if order > 1
end end
% Compute unique sixth-order product moments of u, i.e. unique(E[kron(kron(kron(kron(kron(u,u),u),u),u),u)],'stable') % Compute unique sixth-order product moments of u, i.e. unique(E[kron(kron(kron(kron(kron(u,u),u),u),u),u)],'stable')
u_nbr6 = u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4*(u_nbr+4)/5*(u_nbr+5)/6; u_nbr6 = u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4*(u_nbr+4)/5*(u_nbr+5)/6;
Q6Pu = Q6_plication(u_nbr); if isempty(Q6Pu)
COMBOS6 = flipud(allVL1(u_nbr, 6)); %all possible (unique) combinations of powers that sum up to six Q6Pu = Q6_plication(u_nbr);
COMBOS6 = flipud(allVL1(u_nbr, 6)); %all possible (unique) combinations of powers that sum up to six
end
E_u_u_u_u_u_u = zeros(u_nbr6,1); %only unique entries E_u_u_u_u_u_u = zeros(u_nbr6,1); %only unique entries
if compute_derivs && (stderrparam_nbr+corrparam_nbr>0) if compute_derivs && (stderrparam_nbr+corrparam_nbr>0)
dE_u_u_u_u_u_u = zeros(u_nbr6,stderrparam_nbr+corrparam_nbr); dE_u_u_u_u_u_u = zeros(u_nbr6,stderrparam_nbr+corrparam_nbr);
@ -798,7 +807,7 @@ if order > 1
E_inovzlag1 = zeros(inov_nbr,z_nbr); % Attention: E[inov*z(-1)'] is not equal to zero for a third-order approximation due to kron(kron(xf(-1),u),u) E_inovzlag1 = zeros(inov_nbr,z_nbr); % Attention: E[inov*z(-1)'] is not equal to zero for a third-order approximation due to kron(kron(xf(-1),u),u)
E_inovzlag1(id_inov6_xf_u_u , id_z1_xf ) = kron(E_xfxf,E_uu(:)); E_inovzlag1(id_inov6_xf_u_u , id_z1_xf ) = kron(E_xfxf,E_uu(:));
E_inovzlag1(id_inov6_xf_u_u , id_z4_xrd ) = kron(E_xrdxf',E_uu(:)); E_inovzlag1(id_inov6_xf_u_u , id_z4_xrd ) = kron(E_xrdxf',E_uu(:));
E_inovzlag1(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu)) ; E_inovzlag1(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(reshape(K_xx_x*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu)) ;
E_inovzlag1(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:)); E_inovzlag1(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:));
Binovzlag1A= B*E_inovzlag1*transpose(A); Binovzlag1A= B*E_inovzlag1*transpose(A);
@ -980,7 +989,7 @@ if order > 1
dE_inovzlag1(id_inov6_xf_u_u , id_z1_xf , jp3) = kron(dE_xfxf_jp3,E_uu(:)) + kron(E_xfxf,dE_uu_jp3(:)); dE_inovzlag1(id_inov6_xf_u_u , id_z1_xf , jp3) = kron(dE_xfxf_jp3,E_uu(:)) + kron(E_xfxf,dE_uu_jp3(:));
dE_inovzlag1(id_inov6_xf_u_u , id_z4_xrd , jp3) = kron(dE_xrdxf_jp3',E_uu(:)) + kron(E_xrdxf',dE_uu_jp3(:)); dE_inovzlag1(id_inov6_xf_u_u , id_z4_xrd , jp3) = kron(dE_xrdxf_jp3',E_uu(:)) + kron(E_xrdxf',dE_uu_jp3(:));
dE_inovzlag1(id_inov6_xf_u_u , id_z5_xf_xs , jp3) = kron(reshape(commutation(x_nbr^2,x_nbr)*vec(dE_xsxf_xf_jp3),x_nbr,x_nbr^2),vec(E_uu)) + kron(reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(dE_uu_jp3)) ; dE_inovzlag1(id_inov6_xf_u_u , id_z5_xf_xs , jp3) = kron(reshape(K_xx_x*vec(dE_xsxf_xf_jp3),x_nbr,x_nbr^2),vec(E_uu)) + kron(reshape(K_xx_x*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(dE_uu_jp3)) ;
dE_inovzlag1(id_inov6_xf_u_u , id_z6_xf_xf_xf , jp3) = kron(reshape(dE_xf_xfxf_xf_jp3,x_nbr,x_nbr^3),E_uu(:)) + kron(reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),dE_uu_jp3(:)); dE_inovzlag1(id_inov6_xf_u_u , id_z6_xf_xf_xf , jp3) = kron(reshape(dE_xf_xfxf_xf_jp3,x_nbr,x_nbr^3),E_uu(:)) + kron(reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),dE_uu_jp3(:));
dBinovzlag1A_jp3 = dB(:,:,jp3)*E_inovzlag1*transpose(A) + B*dE_inovzlag1(:,:,jp3)*transpose(A) + B*E_inovzlag1*transpose(dA(:,:,jp3)); dBinovzlag1A_jp3 = dB(:,:,jp3)*E_inovzlag1*transpose(A) + B*dE_inovzlag1(:,:,jp3)*transpose(A) + B*E_inovzlag1*transpose(dA(:,:,jp3));
@ -1029,8 +1038,8 @@ else
+ C(stationary_vars,:)*transpose(E_inovzlag1)*D(stationary_vars,:)'... + C(stationary_vars,:)*transpose(E_inovzlag1)*D(stationary_vars,:)'...
+ D(stationary_vars,:)*Varinov*D(stationary_vars,:)'; + D(stationary_vars,:)*Varinov*D(stationary_vars,:)';
end end
indzeros = find(abs(Var_y) < 1e-12); %find values that are numerical zero
Var_y(indzeros) = 0; Var_y(abs(Var_y) < 1e-12) = 0; %find values that are numerical zero
if useautocorr if useautocorr
sdy = sqrt(diag(Var_y)); %theoretical standard deviation sdy = sqrt(diag(Var_y)); %theoretical standard deviation
sdy = sdy(stationary_vars); sdy = sdy(stationary_vars);
@ -1056,8 +1065,7 @@ if compute_derivs
+ dC(stationary_vars,:,jpV)*transpose(E_inovzlag1)*D(stationary_vars,:)' + C(stationary_vars,:)*transpose(dE_inovzlag1(:,:,jpV))*D(stationary_vars,:)' + C(stationary_vars,:)*transpose(E_inovzlag1)*dD(stationary_vars,:,jpV)'... + dC(stationary_vars,:,jpV)*transpose(E_inovzlag1)*D(stationary_vars,:)' + C(stationary_vars,:)*transpose(dE_inovzlag1(:,:,jpV))*D(stationary_vars,:)' + C(stationary_vars,:)*transpose(E_inovzlag1)*dD(stationary_vars,:,jpV)'...
+ dD(stationary_vars,:,jpV)*Varinov*D(stationary_vars,:)' + D(stationary_vars,:)*dVarinov(:,:,jpV)*D(stationary_vars,:)' + D(stationary_vars,:)*Varinov*dD(stationary_vars,:,jpV)'; + dD(stationary_vars,:,jpV)*Varinov*D(stationary_vars,:)' + D(stationary_vars,:)*dVarinov(:,:,jpV)*D(stationary_vars,:)' + D(stationary_vars,:)*Varinov*dD(stationary_vars,:,jpV)';
end end
indzeros = find(abs(dVar_y_tmp) < 1e-12); %find values that are numerical zero dVar_y_tmp(abs(dVar_y_tmp) < 1e-12) = 0; %find values that are numerical zero
dVar_y_tmp(indzeros) = 0;
dVar_y(stationary_vars,stationary_vars,jpV) = dVar_y_tmp; dVar_y(stationary_vars,stationary_vars,jpV) = dVar_y_tmp;
if useautocorr if useautocorr
dsy = 1/2./sdy.*diag(dVar_y(:,:,jpV)); dsy = 1/2./sdy.*diag(dVar_y(:,:,jpV));
@ -1089,7 +1097,7 @@ for i = 1:nlags
E_inovzlagi = zeros(inov_nbr,z_nbr); E_inovzlagi = zeros(inov_nbr,z_nbr);
E_inovzlagi(id_inov6_xf_u_u , id_z1_xf ) = kron(hxi*E_xfxf,E_uu(:)); E_inovzlagi(id_inov6_xf_u_u , id_z1_xf ) = kron(hxi*E_xfxf,E_uu(:));
E_inovzlagi(id_inov6_xf_u_u , id_z4_xrd ) = kron(hxi*E_xrdxf',E_uu(:)); E_inovzlagi(id_inov6_xf_u_u , id_z4_xrd ) = kron(hxi*E_xrdxf',E_uu(:));
E_inovzlagi(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(hxi*reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu)); E_inovzlagi(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(hxi*reshape(K_xx_x*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu));
E_inovzlagi(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(hxi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:)); E_inovzlagi(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(hxi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:));
Var_yi(stationary_vars,stationary_vars,i) = C(stationary_vars,:)*Var_zi*C(stationary_vars,:)' + C(stationary_vars,:)*Ai*tmp + D(stationary_vars,:)*E_inovzlagi*C(stationary_vars,:)'; Var_yi(stationary_vars,stationary_vars,i) = C(stationary_vars,:)*Var_zi*C(stationary_vars,:)' + C(stationary_vars,:)*Ai*tmp + D(stationary_vars,:)*E_inovzlagi*C(stationary_vars,:)';
end end
@ -1125,12 +1133,12 @@ if compute_derivs
E_inovzlagi = zeros(inov_nbr,z_nbr); E_inovzlagi = zeros(inov_nbr,z_nbr);
E_inovzlagi(id_inov6_xf_u_u , id_z1_xf ) = kron(hxi*E_xfxf,E_uu(:)); E_inovzlagi(id_inov6_xf_u_u , id_z1_xf ) = kron(hxi*E_xfxf,E_uu(:));
E_inovzlagi(id_inov6_xf_u_u , id_z4_xrd ) = kron(hxi*E_xrdxf',E_uu(:)); E_inovzlagi(id_inov6_xf_u_u , id_z4_xrd ) = kron(hxi*E_xrdxf',E_uu(:));
E_inovzlagi(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(hxi*reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu)); E_inovzlagi(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(hxi*reshape(K_xx_x*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu));
E_inovzlagi(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(hxi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:)); E_inovzlagi(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(hxi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:));
dE_inovzlagi_jpVi = zeros(inov_nbr,z_nbr); dE_inovzlagi_jpVi = zeros(inov_nbr,z_nbr);
dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z1_xf ) = kron(dhxi_jpVi*E_xfxf,E_uu(:)) + kron(hxi*dE_xfxf(:,:,jpVi),E_uu(:)) + kron(hxi*E_xfxf,vec(dE_uu(:,:,jpVi))); dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z1_xf ) = kron(dhxi_jpVi*E_xfxf,E_uu(:)) + kron(hxi*dE_xfxf(:,:,jpVi),E_uu(:)) + kron(hxi*E_xfxf,vec(dE_uu(:,:,jpVi)));
dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z4_xrd ) = kron(dhxi_jpVi*E_xrdxf',E_uu(:)) + kron(hxi*dE_xrdxf(:,:,jpVi)',E_uu(:)) + kron(hxi*E_xrdxf',vec(dE_uu(:,:,jpVi))); dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z4_xrd ) = kron(dhxi_jpVi*E_xrdxf',E_uu(:)) + kron(hxi*dE_xrdxf(:,:,jpVi)',E_uu(:)) + kron(hxi*E_xrdxf',vec(dE_uu(:,:,jpVi)));
dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(dhxi_jpVi*reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu)) + kron(hxi*reshape(commutation(x_nbr^2,x_nbr)*vec(dE_xsxf_xf(:,:,jpVi)),x_nbr,x_nbr^2),vec(E_uu)) + kron(hxi*reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(dE_uu(:,:,jpVi))); dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(dhxi_jpVi*reshape(K_xx_x*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu)) + kron(hxi*reshape(K_xx_x*vec(dE_xsxf_xf(:,:,jpVi)),x_nbr,x_nbr^2),vec(E_uu)) + kron(hxi*reshape(K_xx_x*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(dE_uu(:,:,jpVi)));
dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(dhxi_jpVi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:)) + kron(hxi*reshape(dE_xf_xfxf_xf(:,:,jpVi),x_nbr,x_nbr^3),E_uu(:)) + kron(hxi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),vec(dE_uu(:,:,jpVi))); dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(dhxi_jpVi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:)) + kron(hxi*reshape(dE_xf_xfxf_xf(:,:,jpVi),x_nbr,x_nbr^3),E_uu(:)) + kron(hxi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),vec(dE_uu(:,:,jpVi)));
dVar_yi(stationary_vars,stationary_vars,i,jpVi) = dC(stationary_vars,:,jpVi)*Var_zi*C(stationary_vars,:)' + C(stationary_vars,:)*dVar_zi_jpVi*C(stationary_vars,:)' + C(stationary_vars,:)*Var_zi*dC(stationary_vars,:,jpVi)'... dVar_yi(stationary_vars,stationary_vars,i,jpVi) = dC(stationary_vars,:,jpVi)*Var_zi*C(stationary_vars,:)' + C(stationary_vars,:)*dVar_zi_jpVi*C(stationary_vars,:)' + C(stationary_vars,:)*Var_zi*dC(stationary_vars,:,jpVi)'...
+ dC(stationary_vars,:,jpVi)*Ai*tmp + C(stationary_vars,:)*dAi_jpVi*tmp + C(stationary_vars,:)*Ai*dtmp_jpVi... + dC(stationary_vars,:,jpVi)*Ai*tmp + C(stationary_vars,:)*dAi_jpVi*tmp + C(stationary_vars,:)*Ai*dtmp_jpVi...
@ -1179,7 +1187,7 @@ if compute_derivs
end end
end end
end end
non_stationary_vars = setdiff(1:y_nbr,stationary_vars); non_stationary_vars = ~ismember((1:y_nbr)',stationary_vars);
E_y(non_stationary_vars) = NaN; E_y(non_stationary_vars) = NaN;
if compute_derivs if compute_derivs
dE_y(non_stationary_vars,:) = NaN; dE_y(non_stationary_vars,:) = NaN;
@ -1195,7 +1203,7 @@ pruned_state_space.D = D;
pruned_state_space.c = c; pruned_state_space.c = c;
pruned_state_space.d = d; pruned_state_space.d = d;
pruned_state_space.Varinov = Varinov; pruned_state_space.Varinov = Varinov;
pruned_state_space.Var_z = Var_z; %remove in future [@wmutschl] % pruned_state_space.Var_z = Var_z; %
pruned_state_space.Var_y = Var_y; pruned_state_space.Var_y = Var_y;
pruned_state_space.Var_yi = Var_yi; pruned_state_space.Var_yi = Var_yi;
if useautocorr if useautocorr

96
tests/pruning/AnSchorfheide_pruned_state_space.mod
View File

@ -143,51 +143,51 @@ for iorder = 1:3
error('Something wrong with pruned_state_space.m compared to Andreasen et al 2018 Toolbox v2 at order %d.',iorder); error('Something wrong with pruned_state_space.m compared to Andreasen et al 2018 Toolbox v2 at order %d.',iorder);
end end
end end
skipline(); % skipline();
fprintf('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'); % fprintf('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');
fprintf('EXAMPLE:\n') % fprintf('EXAMPLE:\n')
fprintf(' 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'); % fprintf(' 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');
fprintf(' 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'); % fprintf(' 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');
fprintf(' which can be computed using the prodmom.m function by providing E[xf*xf''] as covariance matrix.\n'); % fprintf(' which can be computed using the prodmom.m function by providing E[xf*xf''] as covariance matrix.\n');
fprintf(' In order to replicate this you have to change UnconditionalMoments_3rd_Lyap.m to also output Var_z.\n') % fprintf(' In order to replicate this you have to change UnconditionalMoments_3rd_Lyap.m to also output Var_z.\n')
%
dynare_nx = M_.nspred; % dynare_nx = M_.nspred;
dynare_E_xf2 = pruned_state_space.order_3.Var_z(1:dynare_nx,1:dynare_nx); % dynare_E_xf2 = pruned_state_space.order_3.Var_z(1:dynare_nx,1:dynare_nx);
dynare_E_xf6 = pruned_state_space.order_3.Var_z((end-dynare_nx^3+1):end,(end-dynare_nx^3+1):end); % dynare_E_xf6 = pruned_state_space.order_3.Var_z((end-dynare_nx^3+1):end,(end-dynare_nx^3+1):end);
dynare_E_xf6 = dynare_E_xf6(:); % dynare_E_xf6 = dynare_E_xf6(:);
%
Andreasen_nx = M_.nspred+M_.exo_nbr; % Andreasen_nx = M_.nspred+M_.exo_nbr;
Andreasen_E_xf2 = outAndreasenetal.order_3.Var_z(1:Andreasen_nx,1:Andreasen_nx); % Andreasen_E_xf2 = outAndreasenetal.order_3.Var_z(1:Andreasen_nx,1:Andreasen_nx);
Andreasen_E_xf6 = outAndreasenetal.order_3.Var_z((end-Andreasen_nx^3+1):end,(end-Andreasen_nx^3+1):end); % Andreasen_E_xf6 = outAndreasenetal.order_3.Var_z((end-Andreasen_nx^3+1):end,(end-Andreasen_nx^3+1):end);
Andreasen_E_xf6 = Andreasen_E_xf6(:); % Andreasen_E_xf6 = Andreasen_E_xf6(:);
%
fprintf('Second-order product moments of xf and u are the same:\n') % fprintf('Second-order product moments of xf and u are the same:\n')
norm_E_xf2 = norm(dynare_E_xf2-Andreasen_E_xf2(1:M_.nspred,1:M_.nspred),Inf) % norm_E_xf2 = norm(dynare_E_xf2-Andreasen_E_xf2(1:M_.nspred,1:M_.nspred),Inf)
norm_E_uu = norm(M_.Sigma_e-Andreasen_E_xf2(M_.nspred+(1:M_.exo_nbr),M_.nspred+(1:M_.exo_nbr)),Inf) % norm_E_uu = norm(M_.Sigma_e-Andreasen_E_xf2(M_.nspred+(1:M_.exo_nbr),M_.nspred+(1:M_.exo_nbr)),Inf)
%
% Compute unique sixth-order product moments of xf, i.e. unique(E[kron(kron(kron(kron(kron(xf,xf),xf),xf),xf),xf)],'stable') % % Compute unique sixth-order product moments of xf, i.e. unique(E[kron(kron(kron(kron(kron(xf,xf),xf),xf),xf),xf)],'stable')
dynare_nx6 = dynare_nx*(dynare_nx+1)/2*(dynare_nx+2)/3*(dynare_nx+3)/4*(dynare_nx+4)/5*(dynare_nx+5)/6; % dynare_nx6 = dynare_nx*(dynare_nx+1)/2*(dynare_nx+2)/3*(dynare_nx+3)/4*(dynare_nx+4)/5*(dynare_nx+5)/6;
dynare_Q6Px = Q6_plication(dynare_nx); % dynare_Q6Px = Q6_plication(dynare_nx);
dynare_COMBOS6 = flipud(allVL1(dynare_nx, 6)); %all possible (unique) combinations of powers that sum up to six % dynare_COMBOS6 = flipud(allVL1(dynare_nx, 6)); %all possible (unique) combinations of powers that sum up to six
dynare_true_E_xf6 = zeros(dynare_nx6,1); %only unique entries % dynare_true_E_xf6 = zeros(dynare_nx6,1); %only unique entries
for j6 = 1:size(dynare_COMBOS6,1) % for j6 = 1:size(dynare_COMBOS6,1)
dynare_true_E_xf6(j6) = prodmom(dynare_E_xf2, 1:dynare_nx, dynare_COMBOS6(j6,:)); % dynare_true_E_xf6(j6) = prodmom(dynare_E_xf2, 1:dynare_nx, dynare_COMBOS6(j6,:));
end % end
dynare_true_E_xf6 = dynare_Q6Px*dynare_true_E_xf6; %add duplicate entries % dynare_true_E_xf6 = dynare_Q6Px*dynare_true_E_xf6; %add duplicate entries
norm_dynare_E_xf6 = norm(dynare_true_E_xf6 - dynare_E_xf6, Inf); % norm_dynare_E_xf6 = norm(dynare_true_E_xf6 - dynare_E_xf6, Inf);
%
Andreasen_nx6 = Andreasen_nx*(Andreasen_nx+1)/2*(Andreasen_nx+2)/3*(Andreasen_nx+3)/4*(Andreasen_nx+4)/5*(Andreasen_nx+5)/6; % Andreasen_nx6 = Andreasen_nx*(Andreasen_nx+1)/2*(Andreasen_nx+2)/3*(Andreasen_nx+3)/4*(Andreasen_nx+4)/5*(Andreasen_nx+5)/6;
Andreasen_Q6Px = Q6_plication(Andreasen_nx); % Andreasen_Q6Px = Q6_plication(Andreasen_nx);
Andreasen_COMBOS6 = flipud(allVL1(Andreasen_nx, 6)); %all possible (unique) combinations of powers that sum up to six % Andreasen_COMBOS6 = flipud(allVL1(Andreasen_nx, 6)); %all possible (unique) combinations of powers that sum up to six
Andreasen_true_E_xf6 = zeros(Andreasen_nx6,1); %only unique entries % Andreasen_true_E_xf6 = zeros(Andreasen_nx6,1); %only unique entries
for j6 = 1:size(Andreasen_COMBOS6,1) % for j6 = 1:size(Andreasen_COMBOS6,1)
Andreasen_true_E_xf6(j6) = prodmom(Andreasen_E_xf2, 1:Andreasen_nx, Andreasen_COMBOS6(j6,:)); % Andreasen_true_E_xf6(j6) = prodmom(Andreasen_E_xf2, 1:Andreasen_nx, Andreasen_COMBOS6(j6,:));
end % end
Andreasen_true_E_xf6 = Andreasen_Q6Px*Andreasen_true_E_xf6; %add duplicate entries % Andreasen_true_E_xf6 = Andreasen_Q6Px*Andreasen_true_E_xf6; %add duplicate entries
norm_Andreasen_E_xf6 = norm(Andreasen_true_E_xf6 - Andreasen_E_xf6, Inf); % norm_Andreasen_E_xf6 = norm(Andreasen_true_E_xf6 - Andreasen_E_xf6, Inf);
%
fprintf('Sixth-order product moments of xf and u are not the same!\n'); % fprintf('Sixth-order product moments of xf and u are not the same!\n');
fprintf(' Dynare maximum absolute deviations of sixth-order product moments of xf: %d\n',norm_dynare_E_xf6) % fprintf(' Dynare maximum absolute deviations of sixth-order product moments of xf: %d\n',norm_dynare_E_xf6)
fprintf(' Andreasen et al maximum absolute deviations of sixth-order product moments of xf: %d\n',norm_Andreasen_E_xf6) % fprintf(' Andreasen et al maximum absolute deviations of sixth-order product moments of xf: %d\n',norm_Andreasen_E_xf6)
skipline(); % skipline();
fprintf('Note that the standard deviations are set quite high to make the numerical differences more apparent.\n'); % fprintf('Note that the standard deviations are set quite high to make the numerical differences more apparent.\n');