1function y_=simult_(M_,options_,y0,dr,ex_,iorder)
2% Simulates the model using a perturbation approach, given the path for the exogenous variables and the
3% decision rules.
4%
5% INPUTS
6%    M_       [struct]   model
7%    options_ [struct]   options
8%    y0       [double]   n*1 vector, initial value (n is the number of declared endogenous variables plus the number
9%                        of auxilliary variables for lags and leads); must be in declaration order, i.e. as in M_.endo_names
10%    dr       [struct]   matlab's structure where the reduced form solution of the model is stored.
11%    ex_      [double]   T*q matrix of innovations.
12%    iorder   [integer]  order of the taylor approximation.
13%
14% OUTPUTS
15%    y_       [double]   n*(T+1) time series for the endogenous variables.
16%
17% SPECIAL REQUIREMENTS
18%    none
19
20% Copyright (C) 2001-2019 Dynare Team
21%
22% This file is part of Dynare.
23%
24% Dynare is free software: you can redistribute it and/or modify
25% it under the terms of the GNU General Public License as published by
26% the Free Software Foundation, either version 3 of the License, or
27% (at your option) any later version.
28%
29% Dynare is distributed in the hope that it will be useful,
30% but WITHOUT ANY WARRANTY; without even the implied warranty of
31% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
32% GNU General Public License for more details.
33%
34% You should have received a copy of the GNU General Public License
35% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
36
37iter = size(ex_,1);
38endo_nbr = M_.endo_nbr;
39exo_nbr = M_.exo_nbr;
40
41y_ = zeros(size(y0,1),iter+M_.maximum_lag);
42y_(:,1) = y0;
43
44if options_.loglinear && ~options_.logged_steady_state
45    dr.ys=log(dr.ys);
46end
47
48if ~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
49    if iorder==1
50        y_(:,1) = y_(:,1)-dr.ys;
51    end
52end
53
54if options_.k_order_solver && ~options_.pruning % Call dynare++ routines.
55    ex_ = [zeros(M_.maximum_lag,M_.exo_nbr); ex_];
56    [err, y_] = dynare_simul_(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,exo_nbr, ...
57                              y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed, ...
58                              dr.ys(dr.order_var),dr);
59    mexErrCheck('dynare_simul_', err);
60    y_(dr.order_var,:) = y_;
61else
62    if options_.block
63        if M_.maximum_lag > 0
64            k2 = dr.state_var;
65        else
66            k2 = [];
67        end
68        order_var = 1:endo_nbr;
69        dr.order_var = order_var;
70    else
71        k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
72        k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr;
73        order_var = dr.order_var;
74    end
75
76    switch iorder
77      case 1
78        if isempty(dr.ghu)% For (linearized) deterministic models.
79            for i = 2:iter+M_.maximum_lag
80                yhat = y_(order_var(k2),i-1);
81                y_(order_var,i) = dr.ghx*yhat;
82            end
83        elseif isempty(dr.ghx)% For (linearized) purely forward variables (no state variables).
84            y_(dr.order_var,:) = dr.ghu*transpose(ex_);
85        else
86            epsilon = dr.ghu*transpose(ex_);
87            for i = 2:iter+M_.maximum_lag
88                yhat = y_(order_var(k2),i-1);
89                y_(order_var,i) = dr.ghx*yhat + epsilon(:,i-1);
90            end
91        end
92        y_ = bsxfun(@plus,y_,dr.ys);
93      case 2
94        constant = dr.ys(order_var)+.5*dr.ghs2;
95        if options_.pruning
96            y__ = y0;
97            for i = 2:iter+M_.maximum_lag
98                yhat1 = y__(order_var(k2))-dr.ys(order_var(k2));
99                yhat2 = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
100                epsilon = ex_(i-1,:)';
101                [abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1);
102                mexErrCheck('A_times_B_kronecker_C', err);
103                [abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
104                mexErrCheck('A_times_B_kronecker_C', err);
105                [abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon);
106                mexErrCheck('A_times_B_kronecker_C', err);
107                y_(order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
108                    + abcOut1 + abcOut2 + abcOut3;
109                y__(order_var) = dr.ys(order_var) + dr.ghx*yhat1 + dr.ghu*epsilon;
110            end
111        else
112            for i = 2:iter+M_.maximum_lag
113                yhat = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
114                epsilon = ex_(i-1,:)';
115                [abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat);
116                mexErrCheck('A_times_B_kronecker_C', err);
117                [abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
118                mexErrCheck('A_times_B_kronecker_C', err);
119                [abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon);
120                mexErrCheck('A_times_B_kronecker_C', err);
121                y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
122                    + abcOut1 + abcOut2 + abcOut3;
123            end
124        end
125      case 3
126        % only with pruning
127        % the third moments of the shocks are assumed null. We don't have
128        % an interface for specifying them
129        ghx = dr.ghx;
130        ghu = dr.ghu;
131        ghxx = dr.ghxx;
132        ghxu = dr.ghxu;
133        ghuu = dr.ghuu;
134        ghs2 = dr.ghs2;
135        ghxxx = dr.ghxxx;
136        ghxxu = dr.ghxxu;
137        ghxuu = dr.ghxuu;
138        ghuuu = dr.ghuuu;
139        ghxss = dr.ghxss;
140        ghuss = dr.ghuss;
141        nspred = M_.nspred;
142        ipred = M_.nstatic+(1:nspred);
143        %construction follows Andreasen et al (2013), Technical
144        %Appendix, Formulas (65) and (66)
145        %split into first, second, and third order terms
146        yhat1 = y0(order_var(k2))-dr.ys(order_var(k2));
147        yhat2 = zeros(nspred,1);
148        yhat3 = zeros(nspred,1);
149        for i=2:iter+M_.maximum_lag
150            u = ex_(i-1,:)';
151            %construct terms of order 2 from second order part, based
152            %on linear component yhat1
153            [gyy, err] = A_times_B_kronecker_C(ghxx,yhat1);
154            mexErrCheck('A_times_B_kronecker_C', err);
155            [guu, err] = A_times_B_kronecker_C(ghuu,u);
156            mexErrCheck('A_times_B_kronecker_C', err);
157            [gyu, err] = A_times_B_kronecker_C(ghxu,yhat1,u);
158            mexErrCheck('A_times_B_kronecker_C', err);
159            %construct terms of order 3 from second order part, based
160            %on order 2 component yhat2
161            [gyy12, err] = A_times_B_kronecker_C(ghxx,yhat1,yhat2);
162            mexErrCheck('A_times_B_kronecker_C', err);
163            [gy2u, err] = A_times_B_kronecker_C(ghxu,yhat2,u);
164            mexErrCheck('A_times_B_kronecker_C', err);
165            %construct terms of order 3, all based on first order component yhat1
166            y2a = kron(yhat1,yhat1);
167            [gyyy, err] = A_times_B_kronecker_C(ghxxx,y2a,yhat1);
168            mexErrCheck('A_times_B_kronecker_C', err);
169            u2a = kron(u,u);
170            [guuu, err] = A_times_B_kronecker_C(ghuuu,u2a,u);
171            mexErrCheck('A_times_B_kronecker_C', err);
172            yu = kron(yhat1,u);
173            [gyyu, err] = A_times_B_kronecker_C(ghxxu,yhat1,yu);
174            mexErrCheck('A_times_B_kronecker_C', err);
175            [gyuu, err] = A_times_B_kronecker_C(ghxuu,yu,u);
176            mexErrCheck('A_times_B_kronecker_C', err);
177            %add all terms of order 3, linear component based on third
178            %order yhat3
179            yhat3 = ghx*yhat3 +gyy12 ... % prefactor is 1/2*2=1, see (65) Appendix Andreasen et al.
180                    + gy2u ... % prefactor is 1/2*2=1, see (65) Appendix Andreasen et al.
181                    + 1/6*(gyyy + guuu + 3*(gyyu + gyuu +  ghxss*yhat1 + ghuss*u)); %note: s is treated as variable, thus xss and uss are third order
182            yhat2 = ghx*yhat2 + 1/2*(gyy + guu + 2*gyu + ghs2);
183            yhat1 = ghx*yhat1 + ghu*u;
184            y_(order_var,i) = dr.ys(order_var)+yhat1 + yhat2 + yhat3; %combine terms again
185            yhat1 = yhat1(ipred);
186            yhat2 = yhat2(ipred);
187            yhat3 = yhat3(ipred);
188        end
189      otherwise
190          error(['pruning not available for order = ' int2str(options_.order)])
191    end
192end
193