1function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariables, steadystate, M, options) 2% Performs deterministic simulations with lead or lag of one period, using 3% a basic Newton solver on sparse matrices. 4% Uses perfect_foresight_problem DLL to construct the stacked problem. 5% 6% INPUTS 7% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (initial condition + initial guess + terminal condition). 8% - exogenousvariables [double] T*M array, paths for the exogenous variables. 9% - steadystate [double] N*1 array, steady state for the endogenous variables. 10% - M [struct] contains a description of the model. 11% - options [struct] contains various options. 12% OUTPUTS 13% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (solution of the perfect foresight model). 14% - info [struct] contains informations about the results. 15 16% Copyright (C) 1996-2019 Dynare Team 17% 18% This file is part of Dynare. 19% 20% Dynare is free software: you can redistribute it and/or modify 21% it under the terms of the GNU General Public License as published by 22% the Free Software Foundation, either version 3 of the License, or 23% (at your option) any later version. 24% 25% Dynare is distributed in the hope that it will be useful, 26% but WITHOUT ANY WARRANTY; without even the implied warranty of 27% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 28% GNU General Public License for more details. 29% 30% You should have received a copy of the GNU General Public License 31% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 32 33verbose = options.verbosity && ~options.noprint; 34 35ny = M.endo_nbr; 36periods = options.periods; 37vperiods = periods*ones(1,options.simul.maxit); 38 39if M.maximum_lag > 0 40 y0 = endogenousvariables(:, M.maximum_lag); 41else 42 y0 = NaN(ny, 1); 43end 44 45if M.maximum_lead > 0 46 yT = endogenousvariables(:, M.maximum_lag+periods+1); 47else 48 yT = NaN(ny, 1); 49end 50 51y = reshape(endogenousvariables(:, M.maximum_lag+(1:periods)), ny*periods, 1); 52 53stop = false; 54 55if verbose 56 skipline() 57 printline(56) 58 disp('MODEL SIMULATION:') 59 skipline() 60end 61 62h1 = clock; 63 64for iter = 1:options.simul.maxit 65 h2 = clock; 66 67 [res, A] = perfect_foresight_problem(y, y0, yT, exogenousvariables, M.params, steadystate, periods, M, options); 68 69 if options.endogenous_terminal_period && iter > 1 70 for it = 1:periods 71 if max(abs(res((it-1)*ny+(1:ny)))) < options.dynatol.f/1e7 72 if it < periods 73 res = res(1:(it*ny)); 74 A = A(1:(it*ny), 1:(it*ny)); 75 yT = y(it*ny+(1:ny)); 76 endogenousvariables(:, M.maximum_lag+((it+1):periods)) = reshape(y(it*ny+(1:(ny*(periods-it)))), ny, periods-it); 77 y = y(1:(it*ny)); 78 periods = it; 79 end 80 break 81 end 82 end 83 vperiods(iter) = periods; 84 end 85 86 err = max(abs(res)); 87 if options.debug 88 fprintf('\nLargest absolute residual at iteration %d: %10.3f\n',iter,err); 89 if any(isnan(res)) || any(isinf(res)) || any(any(isnan(endogenousvariables))) || any(any(isinf(endogenousvariables))) 90 fprintf('\nWARNING: NaN or Inf detected in the residuals or endogenous variables.\n'); 91 end 92 skipline() 93 end 94 if verbose 95 fprintf('Iter: %d,\t err. = %g,\t time = %g\n', iter, err, etime(clock,h2)); 96 end 97 if err < options.dynatol.f 98 stop = true; 99 break 100 end 101 if options.simul.robust_lin_solve 102 dy = -lin_solve_robust(A, res, verbose, options); 103 else 104 dy = -lin_solve(A, res, verbose); 105 end 106 if any(isnan(dy)) || any(isinf(dy)) 107 if verbose 108 display_critical_variables(reshape(dy,[ny periods])', M, options.noprint); 109 end 110 end 111 y = y + dy; 112end 113 114endogenousvariables(:, M.maximum_lag+(1:periods)) = reshape(y, ny, periods); 115 116if options.endogenous_terminal_period 117 periods = options.periods; 118 err = evaluate_max_dynamic_residual(str2func([M.fname,'.dynamic']), endogenousvariables, exogenousvariables, M.params, steadystate, periods, ny, M.maximum_endo_lag, M.lead_lag_incidence); 119end 120 121if stop 122 if any(any(isnan(endogenousvariables))) || any(any(isinf(endogenousvariables))) 123 info.status = false;% NaN or Inf occurred 124 info.error = err; 125 info.iterations = iter; 126 info.periods = vperiods(1:iter); 127 if verbose 128 skipline() 129 fprintf('Total time of simulation: %g.\n', etime(clock,h1)) 130 disp('Simulation terminated with NaN or Inf in the residuals or endogenous variables.') 131 display_critical_variables(reshape(dy,[ny periods])', M, options.noprint); 132 disp('There is most likely something wrong with your model. Try model_diagnostics or another simulation method.') 133 printline(105) 134 end 135 else 136 if verbose 137 skipline(); 138 fprintf('Total time of simulation: %g.\n', etime(clock,h1)) 139 printline(56) 140 end 141 info.status = true;% Convergency obtained. 142 info.error = err; 143 info.iterations = iter; 144 info.periods = vperiods(1:iter); 145 end 146elseif ~stop 147 if verbose 148 skipline(); 149 fprintf('Total time of simulation: %g.\n', etime(clock,h1)) 150 disp('Maximum number of iterations is reached (modify option maxit).') 151 printline(62) 152 end 153 info.status = false;% more iterations are needed. 154 info.error = err; 155 info.periods = vperiods(1:iter); 156 info.iterations = options.simul.maxit; 157end 158 159if verbose 160 skipline(); 161end 162 163function x = lin_solve(A, b, verbose) 164if norm(b) < sqrt(eps) % then x = 0 is a solution 165 x = 0; 166 return 167end 168 169x = A\b; 170x(~isfinite(x)) = 0; 171relres = norm(b - A*x) / norm(b); 172if relres > 1e-6 && verbose 173 fprintf('WARNING : Failed to find a solution to the linear system.\n'); 174end 175 176function [ x, flag, relres ] = lin_solve_robust(A, b ,verbose, options) 177if norm(b) < sqrt(eps) % then x = 0 is a solution 178 x = 0; 179 flag = 0; 180 relres = 0; 181 return 182end 183 184x = A\b; 185x(~isfinite(x)) = 0; 186[ x, flag, relres ] = bicgstab(A, b, [], [], [], [], x); % returns immediately if x is a solution 187if flag == 0 188 return 189end 190 191if ~options.noprint 192 disp( relres ); 193end 194 195if verbose 196 fprintf('Initial bicgstab failed, trying alternative start point.\n'); 197end 198old_x = x; 199old_relres = relres; 200[ x, flag, relres ] = bicgstab(A, b); 201if flag == 0 202 return 203end 204 205if verbose 206 fprintf('Alternative start point also failed with bicgstab, trying gmres.\n'); 207end 208if old_relres < relres 209 x = old_x; 210end 211[ x, flag, relres ] = gmres(A, b, [], [], [], [], [], x); 212if flag == 0 213 return 214end 215 216if verbose 217 fprintf('Initial gmres failed, trying alternative start point.\n'); 218end 219old_x = x; 220old_relres = relres; 221[ x, flag, relres ] = gmres(A, b); 222if flag == 0 223 return 224end 225 226if verbose 227 fprintf('Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n'); 228end 229if old_relres < relres 230 x = old_x; 231 relres = old_relres; 232end 233old_x = x; 234old_relres = relres; 235x = pinv(full(A)) * b; 236relres = norm(b - A*x) / norm(b); 237if old_relres < relres 238 x = old_x; 239 relres = old_relres; 240end 241flag = relres > 1e-6; 242if flag ~= 0 && verbose 243 fprintf('WARNING : Failed to find a solution to the linear system\n'); 244end 245 246function display_critical_variables(dyy, M, noprint) 247 248if noprint 249 return 250end 251 252if any(isnan(dyy)) 253 indx = find(any(isnan(dyy))); 254 endo_names= M.endo_names(indx); 255 disp('Last iteration provided NaN for the following variables:') 256 fprintf('%s, ', endo_names{:}), 257 fprintf('\n'), 258end 259if any(isinf(dyy)) 260 indx = find(any(isinf(dyy))); 261 endo_names = M.endo_names(indx); 262 disp('Last iteration diverged (Inf) for the following variables:') 263 fprintf('%s, ', endo_names{:}), 264 fprintf('\n'), 265end 266