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