1function [dr, info] = stochastic_solvers(dr, task, M_, options_, oo_) 2 3% Computes the reduced form solution of a rational expectations model (first, second or third 4% order approximation of the stochastic model around the deterministic steady state). 5% 6% INPUTS 7% - dr [struct] Decision rules for stochastic simulations. 8% - task [integer] scalar, if task = 0 then decision rules are computed and if task = 1 then only eigenvales are computed. 9% - M_ [struct] Definition of the model. 10% - options_ [struct] Options. 11% - oo_ [struct] Results 12% 13% OUTPUTS 14% - dr [struct] Decision rules for stochastic simulations. 15% - info [integer] scalar, error code: 16% 17% info=1 -> the model doesn't define current variables uniquely 18% info=2 -> problem in mjdgges.dll info(2) contains error code. 19% info=3 -> BK order condition not satisfied info(2) contains "distance" 20% absence of stable trajectory. 21% info=4 -> BK order condition not satisfied info(2) contains "distance" 22% indeterminacy. 23% info=5 -> BK rank condition not satisfied. 24% info=6 -> The jacobian matrix evaluated at the steady state is complex. 25% info=9 -> k_order_pert was unable to compute the solution 26 27% Copyright (C) 1996-2019 Dynare Team 28% 29% This file is part of Dynare. 30% 31% Dynare is free software: you can redistribute it and/or modify 32% it under the terms of the GNU General Public License as published by 33% the Free Software Foundation, either version 3 of the License, or 34% (at your option) any later version. 35% 36% Dynare is distributed in the hope that it will be useful, 37% but WITHOUT ANY WARRANTY; without even the implied warranty of 38% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 39% GNU General Public License for more details. 40% 41% You should have received a copy of the GNU General Public License 42% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 43 44info = 0; 45 46if options_.linear 47 options_.order = 1; 48end 49 50local_order = options_.order; 51if local_order~=1 && M_.hessian_eq_zero 52 local_order = 1; 53 warning('stochastic_solvers: using order = 1 because Hessian is equal to zero'); 54end 55if options_.order>2 && ~options_.k_order_solver 56 error('You need to set k_order_solver for order>2') 57end 58 59if options_.aim_solver && (local_order > 1) 60 error('Option "aim_solver" is incompatible with order >= 2') 61end 62 63if M_.maximum_endo_lag == 0 64 if local_order >= 2 65 fprintf('\nSTOCHASTIC_SOLVER: Dynare does not solve purely forward models at higher order.\n') 66 fprintf('STOCHASTIC_SOLVER: To circumvent this restriction, you can add a backward-looking dummy equation of the form:\n') 67 fprintf('STOCHASTIC_SOLVER: junk=0.9*junk(-1);\n') 68 error(['2nd and 3rd order approximation not implemented for purely ' ... 69 'forward models']) 70 end 71 if M_.exo_det_nbr~=0 72 fprintf('\nSTOCHASTIC_SOLVER: Dynare does not solve purely forward models with var_exo_det.\n') 73 fprintf('STOCHASTIC_SOLVER: To circumvent this restriction, you can add a backward-looking dummy equation of the form:\n') 74 fprintf('STOCHASTIC_SOLVER: junk=0.9*junk(-1);\n') 75 error(['var_exo_det not implemented for purely forward models']) 76 end 77end 78 79if M_.maximum_endo_lead==0 && M_.exo_det_nbr~=0 80 fprintf('\nSTOCHASTIC_SOLVER: Dynare does not solve purely backward models with var_exo_det.\n') 81 fprintf('STOCHASTIC_SOLVER: To circumvent this restriction, you can add a foward-looking dummy equation of the form:\n') 82 fprintf('STOCHASTIC_SOLVER: junk=0.9*junk(+1);\n') 83 error(['var_exo_det not implemented for purely backwards models']) 84end 85 86if options_.k_order_solver 87 if options_.risky_steadystate 88 [dr,info] = dyn_risky_steadystate_solver(oo_.steady_state,M_,dr, ... 89 options_,oo_); 90 else 91 orig_order = options_.order; 92 options_.order = local_order; 93 dr = set_state_space(dr,M_,options_); 94 [dr,info] = k_order_pert(dr,M_,options_); 95 options_.order = orig_order; 96 end 97 return 98end 99 100klen = M_.maximum_lag + M_.maximum_lead + 1; 101exo_simul = [repmat(oo_.exo_steady_state',klen,1) repmat(oo_.exo_det_steady_state',klen,1)]; 102iyv = M_.lead_lag_incidence'; 103iyv = iyv(:); 104iyr0 = find(iyv) ; 105 106if M_.exo_nbr == 0 107 oo_.exo_steady_state = [] ; 108end 109 110it_ = M_.maximum_lag + 1; 111z = repmat(dr.ys,1,klen); 112if local_order == 1 113 if (options_.bytecode) 114 [chck, ~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... 115 M_.params, dr.ys, 1); 116 jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd]; 117 else 118 [~,jacobia_] = feval([M_.fname '.dynamic'],z(iyr0),exo_simul, ... 119 M_.params, dr.ys, it_); 120 end 121elseif local_order == 2 122 if (options_.bytecode) 123 [chck, ~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... 124 M_.params, dr.ys, 1); 125 jacobia_ = [loc_dr.g1 loc_dr.g1_x]; 126 else 127 [~,jacobia_,hessian1] = feval([M_.fname '.dynamic'],z(iyr0),... 128 exo_simul, ... 129 M_.params, dr.ys, it_); 130 end 131 if options_.use_dll 132 % In USE_DLL mode, the hessian is in the 3-column sparse representation 133 hessian1 = sparse(hessian1(:,1), hessian1(:,2), hessian1(:,3), ... 134 size(jacobia_, 1), size(jacobia_, 2)*size(jacobia_, 2)); 135 end 136 [infrow,infcol]=find(isinf(hessian1)); 137 if options_.debug 138 if ~isempty(infrow) 139 fprintf('\nSTOCHASTIC_SOLVER: The Hessian of the dynamic model contains Inf.\n') 140 fprintf('STOCHASTIC_SOLVER: Try running model_diagnostics to find the source of the problem.\n') 141 save([M_.fname '_debug.mat'],'hessian1') 142 end 143 end 144 if ~isempty(infrow) 145 info(1)=11; 146 return 147 end 148 [nanrow,nancol]=find(isnan(hessian1)); 149 if options_.debug 150 if ~isempty(nanrow) 151 fprintf('\nSTOCHASTIC_SOLVER: The Hessian of the dynamic model contains NaN.\n') 152 fprintf('STOCHASTIC_SOLVER: Try running model_diagnostics to find the source of the problem.\n') 153 save([M_.fname '_debug.mat'],'hessian1') 154 end 155 end 156 if ~isempty(nanrow) 157 info(1)=12; 158 return 159 end 160end 161 162[infrow,infcol]=find(isinf(jacobia_)); 163 164if options_.debug 165 if ~isempty(infrow) 166 fprintf('\nSTOCHASTIC_SOLVER: The Jacobian of the dynamic model contains Inf. The problem is associated with:\n\n') 167 display_problematic_vars_Jacobian(infrow,infcol,M_,dr.ys,'dynamic','STOCHASTIC_SOLVER: ') 168 save([M_.fname '_debug.mat'],'jacobia_') 169 end 170end 171 172if ~isempty(infrow) 173 info(1)=10; 174 return 175end 176 177if ~isreal(jacobia_) 178 if max(max(abs(imag(jacobia_)))) < 1e-15 179 jacobia_ = real(jacobia_); 180 else 181 if options_.debug 182 [imagrow,imagcol]=find(abs(imag(jacobia_))>1e-15); 183 fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n') 184 display_problematic_vars_Jacobian(imagrow,imagcol,M_,dr.ys,'dynamic','STOCHASTIC_SOLVER: ') 185 end 186 info(1) = 6; 187 info(2) = sum(sum(imag(jacobia_).^2)); 188 return 189 end 190end 191 192[nanrow,nancol]=find(isnan(jacobia_)); 193if options_.debug 194 if ~isempty(nanrow) 195 fprintf('\nSTOCHASTIC_SOLVER: The Jacobian of the dynamic model contains NaN. The problem is associated with:\n\n') 196 display_problematic_vars_Jacobian(nanrow,nancol,M_,dr.ys,'dynamic','STOCHASTIC_SOLVER: ') 197 save([M_.fname '_debug.mat'],'jacobia_') 198 end 199end 200 201if ~isempty(nanrow) 202 info(1) = 8; 203 NaN_params=find(isnan(M_.params)); 204 info(2:length(NaN_params)+1) = NaN_params; 205 return 206end 207 208kstate = dr.kstate; 209nstatic = M_.nstatic; 210nfwrd = M_.nfwrd; 211nspred = M_.nspred; 212nboth = M_.nboth; 213nsfwrd = M_.nsfwrd; 214order_var = dr.order_var; 215nd = size(kstate,1); 216nz = nnz(M_.lead_lag_incidence); 217 218sdyn = M_.endo_nbr - nstatic; 219 220[~,cols_b,cols_j] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, ... 221 order_var)); 222b = zeros(M_.endo_nbr,M_.endo_nbr); 223b(:,cols_b) = jacobia_(:,cols_j); 224 225if M_.maximum_endo_lead == 0 226 % backward models: simplified code exist only at order == 1 227 if local_order == 1 228 [k1,~,k2] = find(kstate(:,4)); 229 dr.ghx(:,k1) = -b\jacobia_(:,k2); 230 if M_.exo_nbr 231 dr.ghu = -b\jacobia_(:,nz+1:end); 232 end 233 dr.eigval = eig(kalman_transition_matrix(dr,nstatic+(1:nspred),1:nspred,M_.exo_nbr)); 234 dr.full_rank = 1; 235 dr.edim = nnz(abs(dr.eigval) > options_.qz_criterium); 236 dr.sdim = nd-dr.edim; 237 if dr.edim 238 temp = sort(abs(dr.eigval)); 239 temp = temp(dr.sdim+1:nd)-1-options_.qz_criterium; 240 info(1) = 3; 241 info(2) = temp'*temp; 242 end 243 else 244 fprintf('\nSTOCHASTIC_SOLVER: Dynare does not solve purely backward models at higher order.\n') 245 fprintf('STOCHASTIC_SOLVER: To circumvent this restriction, you can add a forward-looking dummy equation of the form:\n') 246 fprintf('STOCHASTIC_SOLVER: junk=0.9*junk(+1);\n') 247 error(['2nd and 3rd order approximation not implemented for purely ' ... 248 'backward models']) 249 end 250elseif options_.risky_steadystate 251 orig_order = options_.order; 252 options_.order = local_order; 253 [dr,info] = dyn_risky_steadystate_solver(oo_.steady_state,M_,dr, ... 254 options_,oo_); 255 options_.order = orig_order; 256else 257 % If required, use AIM solver if not check only 258 if options_.aim_solver && (task == 0) 259 [dr,info] = AIM_first_order_solver(jacobia_,M_,dr,options_.qz_criterium); 260 261 else % use original Dynare solver 262 [dr,info] = dyn_first_order_solver(jacobia_,M_,dr,options_,task); 263 if info(1) || task 264 return 265 end 266 end 267 268 if local_order > 1 269 % Second order 270 dr = dyn_second_order_solver(jacobia_,hessian1,dr,M_,... 271 options_.threads.kronecker.sparse_hessian_times_B_kronecker_C); 272 273 % reordering second order derivatives, used for deterministic 274 % variables below 275 k1 = nonzeros(M_.lead_lag_incidence(:,order_var)'); 276 kk = [k1; length(k1)+(1:M_.exo_nbr+M_.exo_det_nbr)']; 277 nk = size(kk,1); 278 kk1 = reshape([1:nk^2],nk,nk); 279 kk1 = kk1(kk,kk); 280 hessian1 = hessian1(:,kk1(:)); 281 end 282end 283 284 285%exogenous deterministic variables 286if M_.exo_det_nbr > 0 287 gx = dr.gx; 288 f1 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+2:end,order_var)))); 289 f0 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var)))); 290 fudet = sparse(jacobia_(:,nz+M_.exo_nbr+1:end)); 291 M1 = inv(f0+[zeros(M_.endo_nbr,nstatic) f1*gx zeros(M_.endo_nbr,nsfwrd-nboth)]); 292 M2 = M1*f1; 293 dr.ghud = cell(M_.exo_det_length,1); 294 dr.ghud{1} = -M1*fudet; 295 for i = 2:M_.exo_det_length 296 dr.ghud{i} = -M2*dr.ghud{i-1}(end-nsfwrd+1:end,:); 297 end 298 299 if local_order > 1 300 lead_lag_incidence = M_.lead_lag_incidence; 301 k0 = find(lead_lag_incidence(M_.maximum_endo_lag+1,order_var)'); 302 k1 = find(lead_lag_incidence(M_.maximum_endo_lag+2,order_var)'); 303 hu = dr.ghu(nstatic+[1:nspred],:); 304 hud = dr.ghud{1}(nstatic+1:nstatic+nspred,:); 305 zx = [eye(nspred);dr.ghx(k0,:);gx*dr.Gy;zeros(M_.exo_nbr+M_.exo_det_nbr, ... 306 nspred)]; 307 zu = [zeros(nspred,M_.exo_nbr); dr.ghu(k0,:); gx*hu; zeros(M_.exo_nbr+M_.exo_det_nbr, ... 308 M_.exo_nbr)]; 309 zud=[zeros(nspred,M_.exo_det_nbr);dr.ghud{1};gx(:,1:nspred)*hud;zeros(M_.exo_nbr,M_.exo_det_nbr);eye(M_.exo_det_nbr)]; 310 R1 = hessian1*kron(zx,zud); 311 dr.ghxud = cell(M_.exo_det_length,1); 312 kf = [M_.endo_nbr-nfwrd-nboth+1:M_.endo_nbr]; 313 kp = nstatic+[1:nspred]; 314 dr.ghxud{1} = -M1*(R1+f1*dr.ghxx(kf,:)*kron(dr.ghx(kp,:),dr.ghud{1}(kp,:))); 315 Eud = eye(M_.exo_det_nbr); 316 for i = 2:M_.exo_det_length 317 hudi = dr.ghud{i}(kp,:); 318 zudi=[zeros(nspred,M_.exo_det_nbr);dr.ghud{i};gx(:,1:nspred)*hudi;zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)]; 319 R2 = hessian1*kron(zx,zudi); 320 dr.ghxud{i} = -M2*(dr.ghxud{i-1}(kf,:)*kron(dr.Gy,Eud)+dr.ghxx(kf,:)*kron(dr.ghx(kp,:),dr.ghud{i}(kp,:)))-M1*R2; 321 end 322 R1 = hessian1*kron(zu,zud); 323 dr.ghudud = cell(M_.exo_det_length,1); 324 dr.ghuud{1} = -M1*(R1+f1*dr.ghxx(kf,:)*kron(dr.ghu(kp,:),dr.ghud{1}(kp,:))); 325 Eud = eye(M_.exo_det_nbr); 326 for i = 2:M_.exo_det_length 327 hudi = dr.ghud{i}(kp,:); 328 zudi=[zeros(nspred,M_.exo_det_nbr);dr.ghud{i};gx(:,1:nspred)*hudi;zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)]; 329 R2 = hessian1*kron(zu,zudi); 330 dr.ghuud{i} = -M2*dr.ghxud{i-1}(kf,:)*kron(hu,Eud)-M1*R2; 331 end 332 R1 = hessian1*kron(zud,zud); 333 dr.ghudud = cell(M_.exo_det_length,M_.exo_det_length); 334 dr.ghudud{1,1} = -M1*R1-M2*dr.ghxx(kf,:)*kron(hud,hud); 335 for i = 2:M_.exo_det_length 336 hudi = dr.ghud{i}(nstatic+1:nstatic+nspred,:); 337 zudi=[zeros(nspred,M_.exo_det_nbr);dr.ghud{i};gx(:,1:nspred)*hudi+dr.ghud{i-1}(kf,:);zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)]; 338 R2 = hessian1*kron(zudi,zudi); 339 dr.ghudud{i,i} = -M2*(dr.ghudud{i-1,i-1}(kf,:)+... 340 2*dr.ghxud{i-1}(kf,:)*kron(hudi,Eud) ... 341 +dr.ghxx(kf,:)*kron(hudi,hudi))-M1*R2; 342 R2 = hessian1*kron(zud,zudi); 343 dr.ghudud{1,i} = -M2*(dr.ghxud{i-1}(kf,:)*kron(hud,Eud)+... 344 dr.ghxx(kf,:)*kron(hud,hudi))... 345 -M1*R2; 346 for j=2:i-1 347 hudj = dr.ghud{j}(kp,:); 348 zudj=[zeros(nspred,M_.exo_det_nbr);dr.ghud{j};gx(:,1:nspred)*hudj;zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)]; 349 R2 = hessian1*kron(zudj,zudi); 350 dr.ghudud{j,i} = -M2*(dr.ghudud{j-1,i-1}(kf,:)+dr.ghxud{j-1}(kf,:)* ... 351 kron(hudi,Eud)+dr.ghxud{i-1}(kf,:)* ... 352 kron(hudj,Eud)+dr.ghxx(kf,:)*kron(hudj,hudi))-M1*R2; 353 end 354 end 355 end 356end 357 358if options_.loglinear 359 % this needs to be extended for order=2,3 360 [il,il1,ik,k1] = indices_lagged_leaded_exogenous_variables(dr.order_var,M_); 361 [illag,illag1,iklag,klag1] = indices_lagged_leaded_exogenous_variables(dr.order_var(M_.nstatic+(1:M_.nspred)),M_); 362 if ~isempty(ik) 363 if M_.nspred > 0 364 dr.ghx(ik,iklag) = repmat(1./dr.ys(k1),1,length(klag1)).*dr.ghx(ik,iklag).* ... 365 repmat(dr.ys(klag1)',length(ik),1); 366 dr.ghx(ik,illag) = repmat(1./dr.ys(k1),1,length(illag)).*dr.ghx(ik,illag); 367 end 368 if M_.exo_nbr > 0 369 dr.ghu(ik,:) = repmat(1./dr.ys(k1),1,M_.exo_nbr).*dr.ghu(ik,:); 370 end 371 end 372 if ~isempty(il) && M_.nspred > 0 373 dr.ghx(il,iklag) = dr.ghx(il,iklag).*repmat(dr.ys(klag1)', ... 374 length(il),1); 375 end 376 if local_order > 1 377 error('Loglinear options currently only works at order 1') 378 end 379end 380end 381