Commit cdb7b018 authored by Johannes Pfeifer 's avatar Johannes Pfeifer
Browse files

Implement Andreasen et al. 2013 pruning at order

Implements and documents the Andreasen et al. 2013 pruning at order=3
and sets it as the default at this order. Michel's pruning based on the
approximation of the forecast function has been assigned the option
pruning_forecast_approximation. The preprocessor-interface still needs
to be added for this option. Moreover, more documentation/a reference
what this option does is needed. At a later point, we might change the
default to Michel's approach.
parent efaa3c71
......@@ -3234,11 +3234,14 @@ period(s). The periods must be strictly positive. Conditional variances are give
decomposition provides the decomposition of the effects of shocks upon
impact. The results are stored in
@code{oo_.conditional_variance_decomposition}
(@pxref{oo_.conditional_variance_decomposition}). The variance decomposition is only conducted, if theoretical moments are requested, i.e. using the @code{periods=0}-option. In case of @code{order=2}, Dynare provides a second-order accurate approximation to the true second moments based on the linear terms of the second-order solution (see @cite{Kim, Kim, Schaumburg and Sims (2008)}).
(@pxref{oo_.conditional_variance_decomposition}). The variance decomposition is only conducted, if theoretical moments are requested, i.e. using the @code{periods=0}-option. In case of @code{order=2}, Dynare provides a second-order accurate approximation to the true second moments based on the linear terms of the second-order solution (see @cite{Kim, Kim, Schaumburg and Sims (2008)}).
@item pruning
Discard higher order terms when iteratively computing simulations of
the solution, as in @cite{Kim, Kim, Schaumburg and Sims (2008)}.
Discard higher order terms when iteratively computing simulations of
the solution. At second order, Dynare uses the algorithm of @cite{Kim, Kim, Schaumburg and Sims (2008)}, while at third order its generalization by @cite{Andreasen, Fernández-Villaverde and Rubio-Ramírez (2013)} is used.
@item pruning_forecast_approximation
This type of pruning generalizes the @cite{Kim, Kim, Schaumburg and Sims (2008)} approach to higher orders based on the k-order approximation of the forecasting function.
@item partial_information
@anchor{partial_information}
......@@ -7725,6 +7728,9 @@ internals --test particle/local_state_iteration
Aguiar, Mark and Gopinath, Gita (2004): ``Emerging Market Business
Cycles: The Cycle is the Trend,'' @i{NBER Working Paper}, 10734
@item
Andreasen, Martin M., Fernández-Villaverde, Jesús and Juan Rubio-Ramírez (2013): ``The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications,'' @i{NBER Working Paper}, 18983
@item
Backus, David K., Patrick J. Kehoe, and Finn E. Kydland (1992):
``International Real Business Cycles,'' @i{Journal of Political
......
......@@ -259,6 +259,8 @@ options_.minimal_solving_periods = 1;
% Solution
options_.order = 2;
options_.pruning = 0;
options_.pruning = 0;
options_.pruning_forecast_approximation = 0;
options_.solve_algo = 2;
options_.linear = 0;
options_.replic = 50;
......
......@@ -39,7 +39,7 @@ switch(order)
dr.g_1 = g_1;
dr.g_2 = g_2;
case 3
if options.pruning
if options.pruning || options_.pruning_forecast_approximation
[err, g_0, g_1, g_2, g_3, derivs] = k_order_perturbation(dr, ...
M,options);
dr.ghx = derivs.gy;
......@@ -67,7 +67,7 @@ switch(order)
error('order > 3 isn''t implemented')
end
if options.pruning
if options.pruning || options_.pruning_forecast_approximation
return
end
......
......@@ -53,7 +53,7 @@ if ~options_.k_order_solver
end
end
if options_.k_order_solver && ~options_.pruning % Call dynare++ routines.
if options_.k_order_solver && ~options_.pruning && ~options_.pruning_forecast_approximation% Call dynare++ routines.
ex_ = [zeros(1,exo_nbr); ex_];
switch options_.order
case 1
......@@ -107,7 +107,7 @@ else
y_ = bsxfun(@plus,y_,dr.ys);
case 2
constant = dr.ys(order_var)+.5*dr.ghs2;
if options_.pruning
if options_.pruning || options_.pruning_forecast_approximation
y__ = y0;
for i = 2:iter+M_.maximum_lag
yhat1 = y__(order_var(k2))-dr.ys(order_var(k2));
......@@ -157,38 +157,87 @@ else
threads = options_.threads.kronecker.A_times_B_kronecker_C;
nspred = M_.nspred;
ipred = M_.nstatic+(1:nspred);
yhat1 = y0(order_var(k2))-dr.ys(order_var(k2));
yhat2 = zeros(nspred,1);
yhat3 = zeros(nspred,1);
for i=2:iter+M_.maximum_lag
u = ex_(i-1,:)';
[gyy, err] = A_times_B_kronecker_C(ghxx,yhat1,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[guu, err] = A_times_B_kronecker_C(ghuu,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[gyu, err] = A_times_B_kronecker_C(ghxu,yhat1,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
y2a = kron(yhat1,yhat1);
[gyyy, err] = A_times_B_kronecker_C(ghxxx,y2a,yhat1,threads);
mexErrCheck('A_times_B_kronecker_C', err);
u2a = kron(u,u);
[guuu, err] = A_times_B_kronecker_C(ghuuu,u2a,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
yu = kron(yhat1,u);
[gyyu, err] = A_times_B_kronecker_C(ghxxu,yhat1,yu,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[gyuu, err] = A_times_B_kronecker_C(ghxuu,yu,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[gyy12, err] = A_times_B_kronecker_C(ghxx,yhat1,yhat2,threads);
mexErrCheck('A_times_B_kronecker_C', err);
yhat3 = ghx*yhat3 + gyyy + guuu + 3*(gyyu + gyuu ...
+ gyy12 + ghxss*yhat1 + ghuss*u);
yhat2 = ghx*yhat2 + gyy + guu + 2*gyu + ghs2;
yhat1 = ghx*yhat1 + ghu*u;
y_(order_var,i) = constant + yhat1 + (1/2)*yhat2 + (1/6)*yhat3;
yhat1 = yhat1(ipred);
yhat2 = yhat2(ipred);
yhat3 = yhat3(ipred);
if options_.pruning_forecast_approximation
yhat1 = y0(order_var(k2))-dr.ys(order_var(k2));
yhat2 = zeros(nspred,1);
yhat3 = zeros(nspred,1);
for i=2:iter+M_.maximum_lag
u = ex_(i-1,:)';
[gyy, err] = A_times_B_kronecker_C(ghxx,yhat1,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[guu, err] = A_times_B_kronecker_C(ghuu,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[gyu, err] = A_times_B_kronecker_C(ghxu,yhat1,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
y2a = kron(yhat1,yhat1);
[gyyy, err] = A_times_B_kronecker_C(ghxxx,y2a,yhat1,threads);
mexErrCheck('A_times_B_kronecker_C', err);
u2a = kron(u,u);
[guuu, err] = A_times_B_kronecker_C(ghuuu,u2a,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
yu = kron(yhat1,u);
[gyyu, err] = A_times_B_kronecker_C(ghxxu,yhat1,yu,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[gyuu, err] = A_times_B_kronecker_C(ghxuu,yu,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[gyy12, err] = A_times_B_kronecker_C(ghxx,yhat1,yhat2,threads);
mexErrCheck('A_times_B_kronecker_C', err);
yhat3 = ghx*yhat3 + gyyy + guuu + 3*(gyyu + gyuu ...
+ gyy12 + ghxss*yhat1 + ghuss*u);
yhat2 = ghx*yhat2 + gyy + guu + 2*gyu + ghs2;
yhat1 = ghx*yhat1 + ghu*u;
y_(order_var,i) = constant + yhat1 + (1/2)*yhat2 + (1/6)*yhat3;
yhat1 = yhat1(ipred);
yhat2 = yhat2(ipred);
yhat3 = yhat3(ipred);
end
else
%construction follows Andreasen et al (2013), Technical
%Appendix, Formulas (65) and (66)
%split into first, second, and third order terms
yhat1 = y0(order_var(k2))-dr.ys(order_var(k2));
yhat2 = zeros(nspred,1);
yhat3 = zeros(nspred,1);
for i=2:iter+M_.maximum_lag
u = ex_(i-1,:)';
%construct terms of order 2 from second order part, based
%on linear component yhat1
[gyy, err] = A_times_B_kronecker_C(ghxx,yhat1,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[guu, err] = A_times_B_kronecker_C(ghuu,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[gyu, err] = A_times_B_kronecker_C(ghxu,yhat1,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
%construct terms of order 3 from second order part, based
%on order 2 component yhat2
[gyy12, err] = A_times_B_kronecker_C(ghxx,yhat1,yhat2,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[gy2u, err] = A_times_B_kronecker_C(ghxu,yhat2,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
%construct terms of order 3, all based on first order component yhat1
y2a = kron(yhat1,yhat1);
[gyyy, err] = A_times_B_kronecker_C(ghxxx,y2a,yhat1,threads);
mexErrCheck('A_times_B_kronecker_C', err);
u2a = kron(u,u);
[guuu, err] = A_times_B_kronecker_C(ghuuu,u2a,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
yu = kron(yhat1,u);
[gyyu, err] = A_times_B_kronecker_C(ghxxu,yhat1,yu,threads);
mexErrCheck('A_times_B_kronecker_C', err);
[gyuu, err] = A_times_B_kronecker_C(ghxuu,yu,u,threads);
mexErrCheck('A_times_B_kronecker_C', err);
%add all terms of order 3, linear component based on third
%order yhat3
yhat3 = ghx*yhat3 +gyy12 ... prefactor is 1/2*2=1, see (65) Appendix Andreasen et al.
+ gy2u ... prefactor is 1/2*2=1, see (65) Appendix Andreasen et al.
+ 1/6*(gyyy + guuu + 3*(gyyu + gyuu + ghxss*yhat1 + ghuss*u)); %note: s is treated as variable, thus xss and uss are third order
yhat2 = ghx*yhat2 + 1/2*(gyy + guu + 2*gyu + ghs2);
yhat1 = ghx*yhat1 + ghu*u;
y_(order_var,i) = dr.ys(order_var)+yhat1 + yhat2 + yhat3; %combine terms again
yhat1 = yhat1(ipred);
yhat2 = yhat2(ipred);
yhat3 = yhat3(ipred);
end
end
end
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment