1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
|
function y_=simult_(M_,options_,y0,dr,ex_,iorder)
% Simulates the model using a perturbation approach, given the path for the exogenous variables and the
% decision rules.
%
% INPUTS
% M_ [struct] model
% options_ [struct] options
% 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); must be in declaration order, i.e. as in M_.endo_names
% 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.
%
% OUTPUTS
% y_ [double] n*(T+1) time series for the endogenous variables.
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2001-2019 Dynare Team
%
% 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/>.
iter = size(ex_,1);
endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
y_ = zeros(size(y0,1),iter+M_.maximum_lag);
y_(:,1) = y0;
if options_.loglinear && ~options_.logged_steady_state
dr.ys=log(dr.ys);
end
if ~options_.k_order_solver || (options_.k_order_solver && options_.pruning) %if k_order_pert is not used or if we do not use Dynare++ with k_order_pert
if iorder==1
y_(:,1) = y_(:,1)-dr.ys;
end
end
if options_.k_order_solver && ~options_.pruning % Call dynare++ routines.
ex_ = [zeros(M_.maximum_lag,M_.exo_nbr); ex_];
[err, y_] = dynare_simul_(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,exo_nbr, ...
y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed, ...
dr.ys(dr.order_var),dr);
mexErrCheck('dynare_simul_', err);
y_(dr.order_var,:) = y_;
else
if options_.block
if M_.maximum_lag > 0
k2 = dr.state_var;
else
k2 = [];
end
order_var = 1:endo_nbr;
dr.order_var = order_var;
else
k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr;
order_var = dr.order_var;
end
switch iorder
case 1
if isempty(dr.ghu)% For (linearized) deterministic models.
for i = 2:iter+M_.maximum_lag
yhat = y_(order_var(k2),i-1);
y_(order_var,i) = dr.ghx*yhat;
end
elseif isempty(dr.ghx)% For (linearized) purely forward variables (no state variables).
y_(dr.order_var,:) = dr.ghu*transpose(ex_);
else
epsilon = dr.ghu*transpose(ex_);
for i = 2:iter+M_.maximum_lag
yhat = y_(order_var(k2),i-1);
y_(order_var,i) = dr.ghx*yhat + epsilon(:,i-1);
end
end
y_ = bsxfun(@plus,y_,dr.ys);
case 2
constant = dr.ys(order_var)+.5*dr.ghs2;
if options_.pruning
y__ = y0;
for i = 2:iter+M_.maximum_lag
yhat1 = y__(order_var(k2))-dr.ys(order_var(k2));
yhat2 = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
epsilon = ex_(i-1,:)';
[abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1);
mexErrCheck('A_times_B_kronecker_C', err);
[abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
mexErrCheck('A_times_B_kronecker_C', err);
[abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon);
mexErrCheck('A_times_B_kronecker_C', err);
y_(order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
+ abcOut1 + abcOut2 + abcOut3;
y__(order_var) = dr.ys(order_var) + dr.ghx*yhat1 + dr.ghu*epsilon;
end
else
for i = 2:iter+M_.maximum_lag
yhat = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
epsilon = ex_(i-1,:)';
[abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat);
mexErrCheck('A_times_B_kronecker_C', err);
[abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
mexErrCheck('A_times_B_kronecker_C', err);
[abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon);
mexErrCheck('A_times_B_kronecker_C', err);
y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
+ abcOut1 + abcOut2 + abcOut3;
end
end
case 3
% only with pruning
% the third moments of the shocks are assumed null. We don't have
% an interface for specifying them
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;
nspred = M_.nspred;
ipred = M_.nstatic+(1:nspred);
%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);
mexErrCheck('A_times_B_kronecker_C', err);
[guu, err] = A_times_B_kronecker_C(ghuu,u);
mexErrCheck('A_times_B_kronecker_C', err);
[gyu, err] = A_times_B_kronecker_C(ghxu,yhat1,u);
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);
mexErrCheck('A_times_B_kronecker_C', err);
[gy2u, err] = A_times_B_kronecker_C(ghxu,yhat2,u);
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);
mexErrCheck('A_times_B_kronecker_C', err);
u2a = kron(u,u);
[guuu, err] = A_times_B_kronecker_C(ghuuu,u2a,u);
mexErrCheck('A_times_B_kronecker_C', err);
yu = kron(yhat1,u);
[gyyu, err] = A_times_B_kronecker_C(ghxxu,yhat1,yu);
mexErrCheck('A_times_B_kronecker_C', err);
[gyuu, err] = A_times_B_kronecker_C(ghxuu,yu,u);
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
otherwise
error(['pruning not available for order = ' int2str(options_.order)])
end
end
|
