simult_.m 7.98 KB
Newer Older
sebastien's avatar
sebastien committed
1
function y_=simult_(y0,dr,ex_,iorder)
2
% Simulates the model using a perturbation approach, given the path for the exogenous variables and the
sebastien's avatar
sebastien committed
3
% decision rules.
michel's avatar
michel committed
4
%
assia's avatar
assia committed
5
% INPUTS
6
7
8
9
10
%    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.
assia's avatar
assia committed
11
12
%
% OUTPUTS
13
%    y_       [double]   n*(T+1) time series for the endogenous variables.
assia's avatar
assia committed
14
15
16
%
% SPECIAL REQUIREMENTS
%    none
17

18
% Copyright (C) 2001-2011 Dynare Team
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
%
% 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/>.
michel's avatar
michel committed
34

35
36
global M_ options_

37
iter = size(ex_,1);
38
39
endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
40

41
y_ = zeros(size(y0,1),iter+M_.maximum_lag);
42
y_(:,1) = y0;
43

44
45
46
47
48
49
% 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

50
51
if ~options_.k_order_solver
    if iorder==1
52
        y_(:,1) = y_(:,1)-dr.ys;
53
54
55
    end
end

56
57
if options_.k_order_solver && ~options_.pruning % Call dynare++ routines.
    ex_ = [zeros(1,exo_nbr); ex_];
58
59
    switch options_.order
      case 1
60
        [err, y_] = dynare_simul_(1,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,exo_nbr, ...
61
                                  y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed,dr.ys(dr.order_var),...
62
                                  zeros(endo_nbr,1),dr.g_1);
63
      case 2
64
        [err, y_] = dynare_simul_(2,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,exo_nbr, ...
65
66
                                  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);
67
      case 3
68
        [err, y_] = dynare_simul_(3,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,exo_nbr, ...
69
70
                                  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);
71
72
73
      otherwise
        error(['order = ' int2str(order) ' isn''t supported'])
    end
74
    mexErrCheck('dynare_simul_', err);
75
76
    y_(dr.order_var,:) = y_;
else
77
    if options_.block
78
79
80
81
82
        if M_.maximum_lag > 0
            k2 = dr.state_var;
        else
            k2 = [];
        end;
83
        order_var = 1:endo_nbr;
84
        dr.order_var = order_var;
85
86
    else
        k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
87
        k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr;
88
89
90
        order_var = dr.order_var;
    end;
    
91
92
    switch iorder
      case 1
93
        if isempty(dr.ghu)% For (linearized) deterministic models.
94
            for i = 2:iter+M_.maximum_lag
95
96
                yhat = y_(order_var(k2),i-1);
                y_(order_var,i) = dr.ghx*yhat;
97
            end
98
99
        elseif isempty(dr.ghx)% For (linearized) purely forward variables (no state variables).
            y_(dr.order_var,:) = dr.ghu*transpose(ex_);
100
        else
101
            epsilon = dr.ghu*transpose(ex_);
102
            for i = 2:iter+M_.maximum_lag
103
104
                yhat = y_(order_var(k2),i-1);
                y_(order_var,i) = dr.ghx*yhat + epsilon(:,i-1);
105
106
            end
        end
107
        y_ = bsxfun(@plus,y_,dr.ys);
108
      case 2
109
        constant = dr.ys(order_var)+.5*dr.ghs2;
110
111
        if options_.pruning
            y__ = y0;
112
            for i = 2:iter+M_.maximum_lag
113
114
                yhat1 = y__(order_var(k2))-dr.ys(order_var(k2));
                yhat2 = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
115
                epsilon = ex_(i-1,:)';
116
                [abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1,options_.threads.kronecker.A_times_B_kronecker_C);
117
                mexErrCheck('A_times_B_kronecker_C', err);
118
                [abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
119
                mexErrCheck('A_times_B_kronecker_C', err);
120
                [abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
121
                mexErrCheck('A_times_B_kronecker_C', err);
122
                y_(order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
123
                    + abcOut1 + abcOut2 + abcOut3;
124
                y__(order_var) = dr.ys(order_var) + dr.ghx*yhat1 + dr.ghu*epsilon;
125
126
            end
        else
127
            for i = 2:iter+M_.maximum_lag
128
                yhat = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
129
                epsilon = ex_(i-1,:)';
130
                [abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat,options_.threads.kronecker.A_times_B_kronecker_C);
131
                mexErrCheck('A_times_B_kronecker_C', err);
132
                [abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
133
                mexErrCheck('A_times_B_kronecker_C', err);
134
                [abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
135
                mexErrCheck('A_times_B_kronecker_C', err);
136
                y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
137
                    + abcOut1 + abcOut2 + abcOut3;
138
            end
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
         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);
189
190
        end
    end
191
end