2008-12-19 17:49:52 +01:00
|
|
|
function y_=simult_(y0,dr,ex_,iorder)
|
2010-05-25 14:00:08 +02:00
|
|
|
% Simulates the model using a perturbation approach, given the path for the exogenous variables and the
|
2008-12-19 17:49:52 +01:00
|
|
|
% decision rules.
|
2005-02-18 20:54:39 +01:00
|
|
|
%
|
2007-12-21 18:13:52 +01:00
|
|
|
% INPUTS
|
2010-05-25 14:00:08 +02:00
|
|
|
% y0 [double] n*1 vector, initial value (n is the number of declared endogenous variables plus the number
|
|
|
|
% of auxilliary variables for lags and leads)
|
|
|
|
% dr [struct] matlab's structure where the reduced form solution of the model is stored.
|
|
|
|
% ex_ [double] T*q matrix of innovations.
|
|
|
|
% iorder [integer] order of the taylor approximation.
|
2007-12-21 18:13:52 +01:00
|
|
|
%
|
|
|
|
% OUTPUTS
|
2010-05-25 14:00:08 +02:00
|
|
|
% y_ [double] n*(T+1) time series for the endogenous variables.
|
2007-12-21 18:13:52 +01:00
|
|
|
%
|
|
|
|
% SPECIAL REQUIREMENTS
|
|
|
|
% none
|
2008-08-01 14:40:33 +02:00
|
|
|
|
2011-02-04 17:27:33 +01:00
|
|
|
% Copyright (C) 2001-2011 Dynare Team
|
2008-08-01 14:40:33 +02:00
|
|
|
%
|
|
|
|
% This file is part of Dynare.
|
|
|
|
%
|
|
|
|
% Dynare is free software: you can redistribute it and/or modify
|
|
|
|
% it under the terms of the GNU General Public License as published by
|
|
|
|
% the Free Software Foundation, either version 3 of the License, or
|
|
|
|
% (at your option) any later version.
|
|
|
|
%
|
|
|
|
% Dynare is distributed in the hope that it will be useful,
|
|
|
|
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
|
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
|
% GNU General Public License for more details.
|
|
|
|
%
|
|
|
|
% You should have received a copy of the GNU General Public License
|
|
|
|
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
2005-02-18 20:54:39 +01:00
|
|
|
|
2010-05-25 14:00:08 +02:00
|
|
|
global M_ options_
|
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
iter = size(ex_,1);
|
2012-07-31 21:50:28 +02:00
|
|
|
endo_nbr = M_.endo_nbr;
|
|
|
|
exo_nbr = M_.exo_nbr;
|
2010-05-25 14:00:08 +02:00
|
|
|
|
2010-11-19 16:30:55 +01:00
|
|
|
y_ = zeros(size(y0,1),iter+M_.maximum_lag);
|
2010-05-25 14:00:08 +02:00
|
|
|
y_(:,1) = y0;
|
2009-12-16 18:17:34 +01:00
|
|
|
|
2011-10-27 16:45:17 +02:00
|
|
|
% stoch_simul sets k_order_solver=1 if order=3, but does so only locally, so we
|
|
|
|
% have to do it here also
|
|
|
|
if options_.order == 3
|
|
|
|
options_.k_order_solver = 1;
|
|
|
|
end
|
|
|
|
|
2010-11-18 14:55:07 +01:00
|
|
|
if ~options_.k_order_solver
|
|
|
|
if iorder==1
|
2010-11-20 11:34:18 +01:00
|
|
|
y_(:,1) = y_(:,1)-dr.ys;
|
2010-11-18 14:55:07 +01:00
|
|
|
end
|
|
|
|
end
|
|
|
|
|
2012-07-31 21:50:28 +02:00
|
|
|
if options_.k_order_solver && ~options_.pruning % Call dynare++ routines.
|
|
|
|
ex_ = [zeros(1,exo_nbr); ex_];
|
2009-12-16 18:17:34 +01:00
|
|
|
switch options_.order
|
|
|
|
case 1
|
2012-07-31 21:50:28 +02:00
|
|
|
[err, y_] = dynare_simul_(1,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,exo_nbr, ...
|
2011-02-04 17:17:48 +01:00
|
|
|
y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed,dr.ys(dr.order_var),...
|
2012-07-31 21:50:28 +02:00
|
|
|
zeros(endo_nbr,1),dr.g_1);
|
2009-12-16 18:17:34 +01:00
|
|
|
case 2
|
2012-07-31 21:50:28 +02:00
|
|
|
[err, y_] = dynare_simul_(2,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,exo_nbr, ...
|
2011-02-04 17:17:48 +01:00
|
|
|
y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed,dr.ys(dr.order_var),dr.g_0, ...
|
|
|
|
dr.g_1,dr.g_2);
|
2009-12-16 18:17:34 +01:00
|
|
|
case 3
|
2012-07-31 21:50:28 +02:00
|
|
|
[err, y_] = dynare_simul_(3,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,exo_nbr, ...
|
2011-02-04 17:17:48 +01:00
|
|
|
y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed,dr.ys(dr.order_var),dr.g_0, ...
|
|
|
|
dr.g_1,dr.g_2,dr.g_3);
|
2009-12-16 18:17:34 +01:00
|
|
|
otherwise
|
|
|
|
error(['order = ' int2str(order) ' isn''t supported'])
|
|
|
|
end
|
2010-09-22 15:43:36 +02:00
|
|
|
mexErrCheck('dynare_simul_', err);
|
2009-12-16 18:17:34 +01:00
|
|
|
y_(dr.order_var,:) = y_;
|
|
|
|
else
|
2010-11-20 16:11:42 +01:00
|
|
|
if options_.block
|
2011-06-18 17:43:45 +02:00
|
|
|
if M_.maximum_lag > 0
|
|
|
|
k2 = dr.state_var;
|
|
|
|
else
|
|
|
|
k2 = [];
|
|
|
|
end;
|
2012-07-31 21:50:28 +02:00
|
|
|
order_var = 1:endo_nbr;
|
2011-06-18 17:43:45 +02:00
|
|
|
dr.order_var = order_var;
|
2010-11-20 16:11:42 +01:00
|
|
|
else
|
|
|
|
k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
|
2012-07-31 21:50:28 +02:00
|
|
|
k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr;
|
2010-11-20 16:11:42 +01:00
|
|
|
order_var = dr.order_var;
|
|
|
|
end;
|
|
|
|
|
2010-05-25 14:00:08 +02:00
|
|
|
switch iorder
|
|
|
|
case 1
|
2010-11-20 11:34:18 +01:00
|
|
|
if isempty(dr.ghu)% For (linearized) deterministic models.
|
2010-11-19 16:30:55 +01:00
|
|
|
for i = 2:iter+M_.maximum_lag
|
2010-11-20 16:11:42 +01:00
|
|
|
yhat = y_(order_var(k2),i-1);
|
|
|
|
y_(order_var,i) = dr.ghx*yhat;
|
2009-12-16 18:17:34 +01:00
|
|
|
end
|
2010-11-20 11:34:18 +01:00
|
|
|
elseif isempty(dr.ghx)% For (linearized) purely forward variables (no state variables).
|
|
|
|
y_(dr.order_var,:) = dr.ghu*transpose(ex_);
|
2009-12-16 18:17:34 +01:00
|
|
|
else
|
2010-11-20 11:34:18 +01:00
|
|
|
epsilon = dr.ghu*transpose(ex_);
|
2010-11-19 16:30:55 +01:00
|
|
|
for i = 2:iter+M_.maximum_lag
|
2010-11-20 16:11:42 +01:00
|
|
|
yhat = y_(order_var(k2),i-1);
|
|
|
|
y_(order_var,i) = dr.ghx*yhat + epsilon(:,i-1);
|
2009-12-16 18:17:34 +01:00
|
|
|
end
|
|
|
|
end
|
2010-11-18 14:55:07 +01:00
|
|
|
y_ = bsxfun(@plus,y_,dr.ys);
|
2010-05-25 14:00:08 +02:00
|
|
|
case 2
|
2010-11-20 16:11:42 +01:00
|
|
|
constant = dr.ys(order_var)+.5*dr.ghs2;
|
2010-05-25 14:00:08 +02:00
|
|
|
if options_.pruning
|
|
|
|
y__ = y0;
|
2010-11-19 16:30:55 +01:00
|
|
|
for i = 2:iter+M_.maximum_lag
|
2010-11-20 16:11:42 +01:00
|
|
|
yhat1 = y__(order_var(k2))-dr.ys(order_var(k2));
|
|
|
|
yhat2 = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
|
2010-05-25 14:00:08 +02:00
|
|
|
epsilon = ex_(i-1,:)';
|
2011-09-02 15:26:46 +02:00
|
|
|
[abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1,options_.threads.kronecker.A_times_B_kronecker_C);
|
2010-09-21 10:40:18 +02:00
|
|
|
mexErrCheck('A_times_B_kronecker_C', err);
|
2011-09-02 15:26:46 +02:00
|
|
|
[abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
2010-09-21 10:40:18 +02:00
|
|
|
mexErrCheck('A_times_B_kronecker_C', err);
|
2011-09-02 15:26:46 +02:00
|
|
|
[abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
2010-09-21 10:40:18 +02:00
|
|
|
mexErrCheck('A_times_B_kronecker_C', err);
|
2010-11-20 16:11:42 +01:00
|
|
|
y_(order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
|
2010-09-21 10:40:18 +02:00
|
|
|
+ abcOut1 + abcOut2 + abcOut3;
|
2010-11-20 16:11:42 +01:00
|
|
|
y__(order_var) = dr.ys(order_var) + dr.ghx*yhat1 + dr.ghu*epsilon;
|
2010-05-25 14:00:08 +02:00
|
|
|
end
|
|
|
|
else
|
2010-11-19 16:30:55 +01:00
|
|
|
for i = 2:iter+M_.maximum_lag
|
2010-11-20 16:11:42 +01:00
|
|
|
yhat = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
|
2010-05-25 14:00:08 +02:00
|
|
|
epsilon = ex_(i-1,:)';
|
2011-09-02 15:26:46 +02:00
|
|
|
[abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat,options_.threads.kronecker.A_times_B_kronecker_C);
|
2010-09-21 10:40:18 +02:00
|
|
|
mexErrCheck('A_times_B_kronecker_C', err);
|
2011-09-02 15:26:46 +02:00
|
|
|
[abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
2010-09-21 10:40:18 +02:00
|
|
|
mexErrCheck('A_times_B_kronecker_C', err);
|
2011-09-02 15:26:46 +02:00
|
|
|
[abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
2010-09-21 10:40:18 +02:00
|
|
|
mexErrCheck('A_times_B_kronecker_C', err);
|
2010-05-25 14:00:08 +02:00
|
|
|
y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
|
2010-09-21 10:40:18 +02:00
|
|
|
+ abcOut1 + abcOut2 + abcOut3;
|
2010-05-25 14:00:08 +02:00
|
|
|
end
|
2012-07-31 21:50:28 +02:00
|
|
|
end
|
|
|
|
case 3
|
|
|
|
% only with pruning
|
|
|
|
ghx = dr.ghx;
|
|
|
|
ghu = dr.ghu;
|
|
|
|
ghxx = dr.ghxx;
|
|
|
|
ghxu = dr.ghxu;
|
|
|
|
ghuu = dr.ghuu;
|
|
|
|
ghs2 = dr.ghs2;
|
|
|
|
ghxxx = dr.ghxxx;
|
|
|
|
ghxxu = dr.ghxxu;
|
|
|
|
ghxuu = dr.ghxuu;
|
|
|
|
ghuuu = dr.ghuuu;
|
|
|
|
ghxss = dr.ghxss;
|
|
|
|
ghuss = dr.ghuss;
|
|
|
|
threads = options_.threads.kronecker.A_times_B_kronecker_C;
|
|
|
|
npred = dr.npred;
|
|
|
|
ipred = dr.nstatic+(1:npred);
|
|
|
|
yhat1 = y0(order_var(k2))-dr.ys(order_var(k2));
|
|
|
|
yhat2 = zeros(npred,1);
|
|
|
|
yhat3 = zeros(npred,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) = yhat1 + (1/2)*yhat2 + (1/6)*yhat3;
|
|
|
|
yhat1 = yhat1(ipred);
|
|
|
|
yhat2 = yhat2(ipred);
|
|
|
|
yhat3 = yhat3(ipred);
|
2009-12-16 18:17:34 +01:00
|
|
|
end
|
|
|
|
end
|
2010-05-25 14:00:08 +02:00
|
|
|
end
|