2015-12-15 17:53:13 +01:00
function [endogenousvariables, info] = sim1 ( endogenousvariables, exogenousvariables, steadystate, M, options)
2019-06-24 17:52:09 +02:00
% Performs deterministic simulations with lead or lag of one period, using
% a basic Newton solver on sparse matrices.
% Uses perfect_foresight_problem DLL to construct the stacked problem.
2006-10-29 18:27:48 +01:00
%
% INPUTS
2019-06-24 17:52:09 +02:00
% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (initial condition + initial guess + terminal condition).
2016-07-05 19:41:24 +02:00
% - exogenousvariables [double] T*M array, paths for the exogenous variables.
2021-03-17 15:02:00 +01:00
% - steadystate [double] N*1 array, steady state for the endogenous variables.
2016-07-05 19:41:24 +02:00
% - M [struct] contains a description of the model.
% - options [struct] contains various options.
2006-10-29 18:27:48 +01:00
% OUTPUTS
2019-06-24 17:52:09 +02:00
% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (solution of the perfect foresight model).
2016-07-05 19:41:24 +02:00
% - info [struct] contains informations about the results.
2008-08-01 14:40:33 +02:00
2022-04-13 13:15:19 +02:00
% Copyright © 1996-2021 Dynare Team
2008-08-01 14:40:33 +02:00
%
% 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
2021-06-09 17:33:48 +02:00
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
2005-02-18 20:54:39 +01:00
2019-09-12 14:28:35 +02:00
verbose = options . verbosity && ~ options . noprint ;
2015-02-12 15:26:11 +01:00
2015-02-11 22:23:21 +01:00
ny = M . endo_nbr ;
2019-06-24 17:52:09 +02:00
periods = options . periods ;
vperiods = periods * ones ( 1 , options . simul . maxit ) ;
2012-06-04 17:23:14 +02:00
2019-06-24 17:52:09 +02:00
if M . maximum_lag > 0
y0 = endogenousvariables ( : , M . maximum_lag ) ;
else
y0 = NaN ( ny , 1 ) ;
end
2005-02-18 20:54:39 +01:00
2019-06-24 17:52:09 +02:00
if M . maximum_lead > 0
yT = endogenousvariables ( : , M . maximum_lag + periods + 1 ) ;
else
yT = NaN ( ny , 1 ) ;
end
2015-12-15 17:53:13 +01:00
2019-06-24 17:52:09 +02:00
y = reshape ( endogenousvariables ( : , M . maximum_lag + ( 1 : periods ) ) , ny * periods , 1 ) ;
2012-06-04 17:23:14 +02:00
2019-06-24 17:52:09 +02:00
stop = false ;
2012-06-04 17:23:14 +02:00
2015-02-12 15:26:11 +01:00
if verbose
skipline ( )
2015-02-18 23:58:37 +01:00
printline ( 56 )
disp ( ' MODEL SIMULATION:' )
skipline ( )
2015-02-12 15:26:11 +01:00
end
2005-02-18 20:54:39 +01:00
2019-06-24 17:52:09 +02:00
h1 = clock ;
2015-12-15 17:53:13 +01:00
2015-05-25 09:43:21 +02:00
for iter = 1 : options . simul . maxit
2019-06-24 17:52:09 +02:00
h2 = clock ;
2019-07-09 14:05:45 +02:00
[ res , A ] = perfect_foresight_problem ( y , y0 , yT , exogenousvariables , M . params , steadystate , periods , M , options ) ;
2021-03-17 15:02:00 +01:00
% A is the stacked Jacobian with period x equations alongs the rows and
% periods times variables (in declaration order) along the columns
if options . debug && iter == 1
row = find ( all ( A == 0 , 2 ) ) ;
column = find ( all ( A == 0 , 1 ) ) ;
if ~ isempty ( row ) || ~ isempty ( column )
fprintf ( ' The stacked Jacobian is singular. The problem derives from:\n' )
if ~ isempty ( row )
time_period = ceil ( row / ny ) ;
equation = row - ny * ( time_period - 1 ) ;
for eq_iter = 1 : length ( equation )
fprintf ( ' The derivative of equation %d at time %d is zero for all variables\n' , equation ( eq_iter ) , time_period ( eq_iter ) ) ;
end
end
if ~ isempty ( column )
time_period = ceil ( column / ny ) ;
variable = column - ny * ( time_period - 1 ) ;
for eq_iter = 1 : length ( variable )
fprintf ( ' The derivative with respect to variable %d at time %d is zero for all equations\n' , variable ( eq_iter ) , time_period ( eq_iter ) ) ;
end
end
end
end
2019-06-24 17:52:09 +02:00
if options . endogenous_terminal_period && iter > 1
for it = 1 : periods
if max ( abs ( res ( ( it - 1 ) * ny + ( 1 : ny ) ) ) ) < options . dynatol . f / 1e7
if it < periods
res = res ( 1 : ( it * ny ) ) ;
A = A ( 1 : ( it * ny ) , 1 : ( it * ny ) ) ;
yT = y ( it * ny + ( 1 : ny ) ) ;
endogenousvariables ( : , M . maximum_lag + ( ( it + 1 ) : periods ) ) = reshape ( y ( it * ny + ( 1 : ( ny * ( periods - it ) ) ) ) , ny , periods - it ) ;
y = y ( 1 : ( it * ny ) ) ;
periods = it ;
end
2013-12-27 18:35:53 +01:00
break
end
end
2019-06-24 17:52:09 +02:00
vperiods ( iter ) = periods ;
2008-10-21 17:33:51 +02:00
end
2019-06-24 17:52:09 +02:00
2012-06-04 17:23:14 +02:00
err = max ( abs ( res ) ) ;
2015-07-21 15:32:33 +02:00
if options . debug
2015-07-07 15:26:49 +02:00
fprintf ( ' \nLargest absolute residual at iteration %d: %10.3f\n' , iter , err ) ;
2019-06-24 17:52:09 +02:00
if any ( isnan ( res ) ) || any ( isinf ( res ) ) || any ( any ( isnan ( endogenousvariables ) ) ) || any ( any ( isinf ( endogenousvariables ) ) )
2015-07-07 15:26:49 +02:00
fprintf ( ' \nWARNING: NaN or Inf detected in the residuals or endogenous variables.\n' ) ;
2013-10-06 20:22:20 +02:00
end
skipline ( )
end
2015-02-18 23:58:37 +01:00
if verbose
2019-06-24 17:52:09 +02:00
fprintf ( ' Iter: %d,\t err. = %g,\t time = %g\n' , iter , err , etime ( clock , h2 ) ) ;
2015-02-18 23:58:37 +01:00
end
2015-02-11 22:23:21 +01:00
if err < options . dynatol . f
2019-06-24 17:52:09 +02:00
stop = true ;
2008-10-21 17:33:51 +02:00
break
end
2019-06-24 17:52:09 +02:00
if options . simul . robust_lin_solve
2020-01-06 09:06:18 +01:00
dy = - lin_solve_robust ( A , res , verbose , options ) ;
2013-12-27 18:35:53 +01:00
else
2019-06-24 17:52:09 +02:00
dy = - lin_solve ( A , res , verbose ) ;
2013-12-27 18:35:53 +01:00
end
2019-06-24 17:52:09 +02:00
if any ( isnan ( dy ) ) || any ( isinf ( dy ) )
2017-03-23 00:01:11 +01:00
if verbose
2019-09-12 14:28:35 +02:00
display_critical_variables ( reshape ( dy , [ ny periods ] ) ' , M , options . noprint ) ;
2017-03-23 00:01:11 +01:00
end
end
2019-06-24 17:52:09 +02:00
y = y + dy ;
2005-02-18 20:54:39 +01:00
end
2019-06-24 17:52:09 +02:00
endogenousvariables ( : , M . maximum_lag + ( 1 : periods ) ) = reshape ( y , ny , periods ) ;
2013-12-28 17:41:36 +01:00
2019-06-24 17:52:09 +02:00
if options . endogenous_terminal_period
periods = options . periods ;
err = evaluate_max_dynamic_residual ( str2func ( [ M . fname , ' .dynamic' ] ) , endogenousvariables , exogenousvariables , M . params , steadystate , periods , ny , M . maximum_endo_lag , M . lead_lag_incidence ) ;
end
2012-06-04 17:23:14 +02:00
2013-10-06 20:22:20 +02:00
if stop
2019-12-27 18:58:32 +01:00
% initial or terminal observations may contain
% harmless NaN or Inf. We test only values computed above
if any ( any ( isnan ( y ) ) ) || any ( any ( isinf ( y ) ) )
2015-12-15 17:53:13 +01:00
info . status = false ; % NaN or Inf occurred
info . error = err ;
info . iterations = iter ;
info . periods = vperiods ( 1 : iter ) ;
2015-02-12 15:26:11 +01:00
if verbose
2015-02-18 23:58:37 +01:00
skipline ( )
2019-06-24 17:52:09 +02:00
fprintf ( ' Total time of simulation: %g.\n' , etime ( clock , h1 ) )
disp ( ' Simulation terminated with NaN or Inf in the residuals or endogenous variables.' )
2019-09-12 14:28:35 +02:00
display_critical_variables ( reshape ( dy , [ ny periods ] ) ' , M , options . noprint ) ;
2015-02-18 23:58:37 +01:00
disp ( ' There is most likely something wrong with your model. Try model_diagnostics or another simulation method.' )
printline ( 105 )
2014-11-06 09:08:50 +01:00
end
2013-10-06 20:22:20 +02:00
else
2015-02-12 15:26:11 +01:00
if verbose
skipline ( ) ;
2019-06-24 17:52:09 +02:00
fprintf ( ' Total time of simulation: %g.\n' , etime ( clock , h1 ) )
2015-02-18 23:58:37 +01:00
printline ( 56 )
2015-02-12 15:26:11 +01:00
end
2015-12-15 17:53:13 +01:00
info . status = true ; % Convergency obtained.
info . error = err ;
info . iterations = iter ;
info . periods = vperiods ( 1 : iter ) ;
2013-10-06 20:22:20 +02:00
end
elseif ~ stop
2015-02-12 15:26:11 +01:00
if verbose
2013-10-06 20:22:20 +02:00
skipline ( ) ;
2019-06-24 17:52:09 +02:00
fprintf ( ' Total time of simulation: %g.\n' , etime ( clock , h1 ) )
2015-02-18 23:58:37 +01:00
disp ( ' Maximum number of iterations is reached (modify option maxit).' )
printline ( 62 )
2013-10-06 20:22:20 +02:00
end
2015-12-15 17:53:13 +01:00
info . status = false ; % more iterations are needed.
info . error = err ;
info . periods = vperiods ( 1 : iter ) ;
info . iterations = options . simul . maxit ;
2005-02-18 20:54:39 +01:00
end
2015-02-12 15:26:11 +01:00
if verbose
skipline ( ) ;
2015-02-16 09:08:02 +01:00
end
2015-09-06 11:41:49 +02:00
2019-06-24 17:52:09 +02:00
function x = lin_solve ( A, b, verbose)
if norm ( b ) < sqrt ( eps ) % then x = 0 is a solution
2017-05-16 15:10:20 +02:00
x = 0 ;
return
end
x = A \ b ;
2019-06-24 17:52:09 +02:00
x ( ~ isfinite ( x ) ) = 0 ;
relres = norm ( b - A * x ) / norm ( b ) ;
2017-05-16 15:10:20 +02:00
if relres > 1e-6 && verbose
2019-06-24 17:52:09 +02:00
fprintf ( ' WARNING : Failed to find a solution to the linear system.\n' ) ;
2017-05-16 15:10:20 +02:00
end
2020-01-06 09:06:18 +01:00
function [ x, flag, relres ] = lin_solve_robust ( A, b ,verbose, options)
2019-06-24 17:52:09 +02:00
if norm ( b ) < sqrt ( eps ) % then x = 0 is a solution
2017-05-16 15:10:20 +02:00
x = 0 ;
flag = 0 ;
relres = 0 ;
return
end
2015-09-06 11:41:49 +02:00
2017-05-16 15:10:20 +02:00
x = A \ b ;
2019-06-24 17:52:09 +02:00
x ( ~ isfinite ( x ) ) = 0 ;
[ x , flag , relres ] = bicgstab ( A , b , [ ] , [ ] , [ ] , [ ] , x ) ; % returns immediately if x is a solution
2017-05-16 15:10:20 +02:00
if flag == 0
return
end
2015-09-06 11:41:49 +02:00
2019-09-12 14:28:35 +02:00
if ~ options . noprint
disp ( relres ) ;
end
2015-09-06 11:41:49 +02:00
2017-05-16 15:10:20 +02:00
if verbose
2019-06-24 17:52:09 +02:00
fprintf ( ' Initial bicgstab failed, trying alternative start point.\n' ) ;
2017-05-16 15:10:20 +02:00
end
old_x = x ;
old_relres = relres ;
2019-06-24 17:52:09 +02:00
[ x , flag , relres ] = bicgstab ( A , b ) ;
2017-05-16 15:10:20 +02:00
if flag == 0
return
end
2015-09-06 11:41:49 +02:00
2017-05-16 15:10:20 +02:00
if verbose
2019-06-24 17:52:09 +02:00
fprintf ( ' Alternative start point also failed with bicgstab, trying gmres.\n' ) ;
2017-05-16 15:10:20 +02:00
end
if old_relres < relres
x = old_x ;
end
2019-06-24 17:52:09 +02:00
[ x , flag , relres ] = gmres ( A , b , [ ] , [ ] , [ ] , [ ] , [ ] , x ) ;
2017-05-16 15:10:20 +02:00
if flag == 0
return
end
2015-09-06 11:41:49 +02:00
2017-05-16 15:10:20 +02:00
if verbose
2019-06-24 17:52:09 +02:00
fprintf ( ' Initial gmres failed, trying alternative start point.\n' ) ;
2017-05-16 15:10:20 +02:00
end
old_x = x ;
old_relres = relres ;
2019-06-24 17:52:09 +02:00
[ x , flag , relres ] = gmres ( A , b ) ;
2017-05-16 15:10:20 +02:00
if flag == 0
return
end
2017-03-23 00:01:11 +01:00
2017-05-16 15:10:20 +02:00
if verbose
2019-06-24 17:52:09 +02:00
fprintf ( ' Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n' ) ;
2017-05-16 15:10:20 +02:00
end
if old_relres < relres
x = old_x ;
relres = old_relres ;
end
old_x = x ;
old_relres = relres ;
2019-06-24 17:52:09 +02:00
x = pinv ( full ( A ) ) * b ;
relres = norm ( b - A * x ) / norm ( b ) ;
2017-05-16 15:10:20 +02:00
if old_relres < relres
x = old_x ;
relres = old_relres ;
end
flag = relres > 1e-6 ;
if flag ~= 0 && verbose
2019-06-24 17:52:09 +02:00
fprintf ( ' WARNING : Failed to find a solution to the linear system\n' ) ;
2017-05-16 15:10:20 +02:00
end
2019-09-12 14:28:35 +02:00
function display_critical_variables ( dyy, M, noprint)
if noprint
return
end
2017-03-23 00:01:11 +01:00
2017-05-16 15:10:20 +02:00
if any ( isnan ( dyy ) )
indx = find ( any ( isnan ( dyy ) ) ) ;
2017-10-10 10:05:59 +02:00
endo_names = M . endo_names ( indx ) ;
2017-05-16 15:10:20 +02:00
disp ( ' Last iteration provided NaN for the following variables:' )
2017-10-10 10:05:59 +02:00
fprintf ( ' %s, ' , endo_names { : } ) ,
2017-05-16 15:10:20 +02:00
fprintf ( ' \n' ) ,
end
if any ( isinf ( dyy ) )
indx = find ( any ( isinf ( dyy ) ) ) ;
2017-10-10 10:05:59 +02:00
endo_names = M . endo_names ( indx ) ;
2017-05-16 15:10:20 +02:00
disp ( ' Last iteration diverged (Inf) for the following variables:' )
2017-10-10 10:05:59 +02:00
fprintf ( ' %s, ' , endo_names { : } ) ,
2017-05-16 15:10:20 +02:00
fprintf ( ' \n' ) ,
end