stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Perfect foresight: inner functions no longer return a modified oo_ 

They now only return what’s really their output (simulated paths, maximum
residual error…). This is a move towards a more functional programming style.
remove-submodule
Sébastien Villemot 2023-06-06 18:13:10 +02:00
parent 6ed90b3dbf
commit 7722e8e36b
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
16 changed files with 171 additions and 249 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
matlab/ep/extended_path_core.m
View File

@ -4,7 +4,7 @@ function [y, info_convergence, endogenousvariablespaths] = extended_path_core(pe
debug,order,M,pfm,algo,solve_algo,stack_solve_algo,...
olmmcp,options,oo,initialguess)
% Copyright © 2016-2020 Dynare Team
% Copyright © 2016-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -56,32 +56,25 @@ if order == 0
options.lmmcp = olmmcp;
options.solve_algo = solve_algo;
options.stack_solve_algo = stack_solve_algo;
tmp = perfect_foresight_solver_core(M, options, oo);
if ~tmp.deterministic_simulation.status
info_convergence = false;
else
info_convergence = true;
end
[endogenousvariablespaths, info_convergence] = perfect_foresight_solver_core(M, options, oo);
else
switch(algo)
case 0
[flag, tmp.endo_simul] = ...
[flag, endogenousvariablespaths] = ...
solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1
[flag, tmp.endo_simul] = ...
[flag, endogenousvariablespaths] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order);
end
info_convergence = ~flag;
end
if ~info_convergence && ~options.no_homotopy
[info_convergence, tmp.endo_simul] = extended_path_homotopy(endo_simul, exo_simul, M, options, oo, pfm, ep, order, algo, 2, debug);
[info_convergence, endogenousvariablespaths] = extended_path_homotopy(endo_simul, exo_simul, M, options, oo, pfm, ep, order, algo, 2, debug);
end
if info_convergence
y = tmp.endo_simul(:,2);
y = endogenousvariablespaths(:,2);
else
y = NaN(size(endo_nbr,1));
end
endogenousvariablespaths = tmp.endo_simul;

20
matlab/ep/extended_path_homotopy.m
View File

@ -35,20 +35,20 @@ if ismember(method, [1, 2])
oo.endo_simul(:,1) = oo.steady_state + weight*(endo_simul0(:,1) - oo.steady_state);
oo.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo.exo_steady_state));
if order==0
tmp = perfect_foresight_solver_core(M, options, oo);
[endo_simul_new, success] = perfect_foresight_solver_core(M, options, oo);
else
switch(algo)
case 0
[flag, tmp.endo_simul] = ...
[flag, endo_simul_new] = ...
solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1
[flag, tmp.endo_simul] = ...
[flag, endo_simul_new] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order);
end
end
if isequal(order, 0)
% Logical variable flag is false iff the solver fails.
flag = tmp.deterministic_simulation.status;
flag = success;
else
% Fix convention issue on the value of flag.
flag = ~flag;
@ -70,7 +70,7 @@ if ismember(method, [1, 2])
oldweight = weight;
weight = min(weight*increase_factor, 1);
increase_flag = true;
endo_simul = tmp.endo_simul;
endo_simul = endo_simul_new;
else
if increase_flag
weight = oldweight + (weight-oldweight)/100;
@ -102,20 +102,20 @@ if isequal(method, 3) || (isequal(method, 2) && noconvergence)
oo.endo_simul = endo_simul;
oo.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo.exo_steady_state));
if order==0
tmp = perfect_foresight_solver_core(M, options, oo);
[endo_simul_new, success] = perfect_foresight_solver_core(M, options, oo);
else
switch(algo)
case 0
[flag, tmp.endo_simul] = ...
[flag, endo_simul_new] = ...
solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1
[flag, tmp.endo_simul] = ...
[flag, endo_simul_new] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order);
end
end
if isequal(order, 0)
% Logical variable flag is false iff the solver fails.
flag = tmp.deterministic_simulation.status;
flag = success;
else
% Fix convention issue on the value of flag.
flag = ~flag;
@ -130,7 +130,7 @@ if isequal(method, 3) || (isequal(method, 2) && noconvergence)
continue
end
index = index+1;
endo_simul = tmp.endo_simul;
endo_simul = endo_simul_new;
else
break
end

10
matlab/evaluate_steady_state.m
View File

@ -370,11 +370,11 @@ elseif ~options.bytecode && options.block
end
else
nze = length(M.block_structure_stat.block(b).g1_sparse_rowval);
[ys, T, ~, info2] = solve_one_boundary(fh_static, ys, exo_ss, ...
params, [], T, mfs_idx, nze, 1, false, b, 0, options.simul.maxit, ...
options.solve_tolf, ...
0, options.solve_algo, true, false, false, M, options, []);
if info2
[ys, T, success] = solve_one_boundary(fh_static, ys, exo_ss, ...
params, [], T, mfs_idx, nze, 1, false, b, 0, options.simul.maxit, ...
options.solve_tolf, ...
0, options.solve_algo, true, false, false, M, options);
if ~success
check = 1;
break
end

37
matlab/perfect-foresight-models/perfect_foresight_solver.m
View File

@ -54,10 +54,10 @@ end
initperiods = 1:M_.maximum_lag;
lastperiods = (M_.maximum_lag+periods+1):(M_.maximum_lag+periods+M_.maximum_lead);
oo_ = perfect_foresight_solver_core(M_,options_,oo_);
[oo_.endo_simul, success, maxerror, solver_iter, per_block_status] = perfect_foresight_solver_core(M_,options_,oo_);
% If simulation failed try homotopy.
if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
if ~success && ~options_.no_homotopy
if ~options_.noprint
fprintf('\nSimulation of the perfect foresight model failed!\n')
@ -131,13 +131,13 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
saved_endo_simul = oo_.endo_simul;
% Solve for the paths of the endogenous variables.
[oo_,me] = perfect_foresight_solver_core(M_,options_,oo_);
[oo_.endo_simul, success, maxerror, solver_iter, per_block_status] = perfect_foresight_solver_core(M_,options_,oo_);
if oo_.deterministic_simulation.status
if success
current_weight = new_weight;
if current_weight >= 1
if ~options_.noprint
fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', me)
fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', maxerror)
end
break
end
@ -147,7 +147,7 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
step = step * 2;
end
if ~options_.noprint
fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', me)
fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', maxerror)
end
else
% If solver failed, then go back.
@ -155,8 +155,8 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
success_counter = 0;
step = step / 2;
if ~options_.noprint
if isreal(me)
fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'failed', me)
if isreal(maxerror)
fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'failed', maxerror)
else
fprintf('%i \t | %1.5f \t | %s \t | %s\n', iteration, new_weight, 'failed', 'Complex')
end
@ -186,16 +186,17 @@ if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_fores
yy = real(oo_.endo_simul(:,M_.maximum_lag+(1:periods)));
residuals = perfect_foresight_problem(yy(:), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_, options_);
if max(abs(residuals))< options_.dynatol.f
oo_.deterministic_simulation.status = true;
maxerror = max(max(abs(residuals)))
if maxerror < options_.dynatol.f
success = true;
oo_.endo_simul=real(oo_.endo_simul);
else
oo_.deterministic_simulation.status = false;
success = false;
disp('Simulation terminated with imaginary parts in the residuals or endogenous variables.')
end
end
if oo_.deterministic_simulation.status
if success
if ~options_.noprint
fprintf('Perfect foresight solution found.\n\n')
end
@ -203,6 +204,16 @@ else
fprintf('Failed to solve perfect foresight model\n\n')
end
% Put solver status information in oo_
oo_.deterministic_simulation.status = success;
oo_.deterministic_simulation.error = maxerror;
if ~isempty(solver_iter)
oo_.deterministic_simulation.iterations = solver_iter;
end
if ~isempty(per_block_status)
oo_.deterministic_simulation.block = per_block_status;
end
dyn2vec(M_, oo_, options_);
if isfield(oo_, 'initval_series') && ~isempty(oo_.initval_series)
@ -222,6 +233,6 @@ end
assignin('base', 'Simulated_time_series', ts);
if oo_.deterministic_simulation.status
if success
oo_.gui.ran_perfect_foresight = true;
end

87
matlab/perfect-foresight-models/perfect_foresight_solver_core.m
View File

@ -1,4 +1,4 @@
function [oo_, maxerror] = perfect_foresight_solver_core(M_, options_, oo_)
function [y, success, maxerror, iter, per_block_status] = perfect_foresight_solver_core(M_, options_, oo_)
% Core function calling solvers for perfect foresight model
%
@ -8,8 +8,11 @@ function [oo_, maxerror] = perfect_foresight_solver_core(M_, options_, oo_)
% - oo_ [struct] contains results
%
% OUTPUTS
% - oo_ [struct] contains results
% - y [double array] path for the endogenous variables (solution)
% - success [logical] Whether a solution was found
% - maxerror [double] contains the maximum absolute error
% - iter [integer] Number of iterations of the underlying nonlinear solver (empty for non-iterative methods)
% - per_block_status [struct] In the case of block decomposition, provides per-block solver status information (empty if no block decomposition)
% Copyright © 2015-2023 Dynare Team
%
@ -49,62 +52,60 @@ if options_.linear && (isequal(options_.stack_solve_algo, 0) || isequal(options_
options_.linear_approximation = true;
end
maxerror = [];
iter = [];
per_block_status = [];
if options_.block
if M_.block_structure.time_recursive
error('Internal error: can''t perform stacked perfect foresight simulation with time-recursive block decomposition')
end
if options_.bytecode
try
oo_.endo_simul = bytecode('dynamic', 'block_decomposed', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
oo_.deterministic_simulation.status = true;
y = bytecode('dynamic', 'block_decomposed', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
success = true;
catch ME
disp(ME.message)
if options_.no_homotopy
error('Error in bytecode')
end
oo_.deterministic_simulation.status = false;
success = false;
end
else
oo_ = solve_block_decomposed_problem(options_, M_, oo_);
[y, success, maxerror, per_block_status] = solve_block_decomposed_problem(options_, M_, oo_);
end
else
if options_.bytecode
try
oo_.endo_simul = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods);
oo_.deterministic_simulation.status = true;
y = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods);
success = true;
catch ME
disp(ME.message)
if options_.no_homotopy
error('Error in bytecode')
end
oo_.deterministic_simulation.status = false;
success = false;
end
else
if M_.maximum_endo_lead == 0 && M_.maximum_endo_lag>0 && ~options_.lmmcp.status % Purely backward model
[oo_.endo_simul, oo_.deterministic_simulation] = ...
sim1_purely_backward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
[y, success] = sim1_purely_backward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
elseif M_.maximum_endo_lag == 0 && M_.maximum_endo_lead>0 && ~options_.lmmcp.status % Purely forward model
[oo_.endo_simul, oo_.deterministic_simulation] = ...
sim1_purely_forward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
[y, success] = sim1_purely_forward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
elseif M_.maximum_endo_lag == 0 && M_.maximum_endo_lead == 0 && ~options_.lmmcp.status % Purely static model
[oo_.endo_simul, oo_.deterministic_simulation] = ...
sim1_purely_static(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
[y, success] = sim1_purely_static(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
else % General case
switch options_.stack_solve_algo
case 0
if options_.linear_approximation
[oo_.endo_simul, oo_.deterministic_simulation] = ...
sim1_linear(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
[y, success, maxerror] = sim1_linear(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
else
[oo_.endo_simul, oo_.deterministic_simulation] = ...
sim1(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
[y, success, maxerror, iter] = sim1(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
end
case {1 6}
if options_.linear_approximation
error('Invalid value of stack_solve_algo option!')
end
[oo_.endo_simul, oo_.deterministic_simulation] = ...
sim1_lbj(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
[y, success, maxerror, iter] = sim1_lbj(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
case 7
if options_.linear_approximation
if isequal(options_.solve_algo, 10)
@ -116,11 +117,9 @@ else
warning('It would be more efficient to set option solve_algo equal to 0!')
end
end
[oo_.endo_simul, oo_.deterministic_simulation] = ...
solve_stacked_linear_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
[y, success] = solve_stacked_linear_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
else
[oo_.endo_simul, oo_.deterministic_simulation, residuals] = ...
solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
[y, success, maxerror] = solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
end
otherwise
error('Invalid value of stack_solve_algo option!')
@ -129,29 +128,21 @@ else
end
end
if nargout>1
if options_.lmmcp.status
maxerror = max(max(abs(residuals)));
elseif options_.block && ~options_.bytecode
maxerror = oo_.deterministic_simulation.error;
% Some solvers do not compute the maximum error, so do it here if needed
if nargout > 2 && isempty(maxerror)
ny = size(oo_.endo_simul, 1);
if M_.maximum_lag > 0
y0 = y(:, M_.maximum_lag);
else
if options_.bytecode
residuals = bytecode('dynamic','evaluate', oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1);
else
if M_.maximum_lag > 0
y0 = oo_.endo_simul(:, M_.maximum_lag);
else
y0 = NaN(ny, 1);
end
if M_.maximum_lead > 0
yT = oo_.endo_simul(:, M_.maximum_lag+periods+1);
else
yT = NaN(ny, 1);
end
yy = oo_.endo_simul(:,M_.maximum_lag+(1:periods));
residuals = perfect_foresight_problem(yy(:), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_, options_);
end
maxerror = max(max(abs(residuals)));
y0 = NaN(ny, 1);
end
if M_.maximum_lead > 0
yT = y(:, M_.maximum_lag+periods+1);
else
yT = NaN(ny, 1);
end
yy = y(:,M_.maximum_lag+(1:periods));
residuals = perfect_foresight_problem(yy(:), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_, options_);
maxerror = max(max(abs(residuals)));
end

23
matlab/perfect-foresight-models/sim1.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, exogenousvariables, steadystate, M, options)
% 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.
@ -11,9 +11,11 @@ function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariab
% - options [struct] contains various options.
% OUTPUTS
% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (solution of the perfect foresight model).
% - info [struct] contains informations about the results.
% - success [logical] Whether a solution was found
% - err [double] ∞-norm of the residual
% - iter [integer] Number of iterations
% Copyright © 1996-2021 Dynare Team
% Copyright © 1996-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -145,10 +147,7 @@ if stop
% 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)))
info.status = false;% NaN or Inf occurred
info.error = err;
info.iterations = iter;
info.periods = vperiods(1:iter);
success = false; % NaN or Inf occurred
if verbose
skipline()
fprintf('Total time of simulation: %g.\n', etime(clock,h1))
@ -163,10 +162,7 @@ if stop
fprintf('Total time of simulation: %g.\n', etime(clock,h1))
printline(56)
end
info.status = true;% Convergency obtained.
info.error = err;
info.iterations = iter;
info.periods = vperiods(1:iter);
success = true; % Convergency obtained.
end
elseif ~stop
if verbose
@ -175,10 +171,7 @@ elseif ~stop
disp('Maximum number of iterations is reached (modify option maxit).')
printline(62)
end
info.status = false;% more iterations are needed.
info.error = err;
info.periods = vperiods(1:iter);
info.iterations = options.simul.maxit;
success = false; % more iterations are needed.
end
if verbose

11
matlab/perfect-foresight-models/sim1_lbj.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, info] = sim1_lbj(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success, err, iter] = sim1_lbj(endogenousvariables, exogenousvariables, steadystate, M, options)
% Performs deterministic simulations with lead or lag on one period using the historical LBJ algorithm
%
@ -94,9 +94,7 @@ for iter = 1:options.simul.maxit
skipline()
disp(sprintf('Total time of simulation: %s', num2str(etime(clock,h1))))
end
info.status = 1;% Convergency obtained.
info.error = err;
info.iterations = iter;
success = true; % Convergency obtained.
break
end
end
@ -106,10 +104,7 @@ if ~stop
disp(sprintf('Total time of simulation: %s.', num2str(etime(clock,h1))))
disp('Maximum number of iterations is reached (modify option maxit).')
end
info.status = 0;% more iterations are needed.
info.error = err;
info.errors = c/abs(err);
info.iterations = options.simul.maxit;
success = false; % More iterations are needed.
end
if verbose

20
matlab/perfect-foresight-models/sim1_linear.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, info] = sim1_linear(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options)
function [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options)
% Solves a linear approximation of a perfect foresight model using sparse matrix.
%
@ -12,7 +12,8 @@ function [endogenousvariables, info] = sim1_linear(endogenousvariables, exogenou
%
% OUTPUTS
% - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
% - info [struct] contains informations about the results.
% - success [logical] Whether a solution was found
% - ERR [double] ∞-norm of the residual
%
% NOTATIONS
% - N is the number of endogenous variables.
@ -38,7 +39,7 @@ function [endogenousvariables, info] = sim1_linear(endogenousvariables, exogenou
% to center the variables around the deterministic steady state to solve the
% perfect foresight model.
% Copyright © 2015-2020 Dynare Team
% Copyright © 2015-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -209,9 +210,8 @@ try
catch
% Normally, because the model is linear, the solution of the perfect foresight model should
% be obtained in one Newton step. This is not the case if the model is singular.
info.status = false;
info.error = NaN;
info.iterations = 1;
success = false;
ERR = [];
if verbose
skipline()
disp('Singularity problem! The jacobian matrix of the stacked model cannot be inverted.')
@ -238,9 +238,7 @@ if verbose
end
if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y)) || ~isreal(res) || ~isreal(Y)
info.status = false;% NaN or Inf occurred
info.error = ERR;
info.iterations = 1;
success = false; % NaN or Inf occurred
endogenousvariables = reshape(Y, ny, periods+maximum_lag+M.maximum_lead);
if verbose
skipline()
@ -252,9 +250,7 @@ if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y)) || ~isre
disp('There is most likely something wrong with your model. Try model_diagnostics or another simulation method.')
end
else
info.status = true;% Convergency obtained.
info.error = ERR;
info.iterations = 1;
success = true; % Convergency obtained.
endogenousvariables = bsxfun(@plus, reshape(Y, ny, periods+maximum_lag+M.maximum_lead), steadystate_y);
end

8
matlab/perfect-foresight-models/sim1_purely_backward.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, info] = sim1_purely_backward(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success] = sim1_purely_backward(endogenousvariables, exogenousvariables, steadystate, M, options)
% Performs deterministic simulation of a purely backward model
@ -34,7 +34,7 @@ function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse
sparse_rowval, sparse_colval, sparse_colptr, T);
end
info.status = true;
success = true;
for it = M.maximum_lag + (1:options.periods)
y = endogenousvariables(:,it-1); % Values at previous period, also used as guess value for current period
@ -52,7 +52,7 @@ for it = M.maximum_lag + (1:options.periods)
options.dynatol.x, options, ...
feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
if check
info.status = false;
success = false;
if options.debug
dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
end
@ -71,7 +71,7 @@ for it = M.maximum_lag + (1:options.periods)
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, dynamic_resid, dynamic_g1, y, x, M.params, steadystate, M.dynamic_g1_sparse_rowval, M.dynamic_g1_sparse_colval, M.dynamic_g1_sparse_colptr);
if check
info.status = false;
success = false;
dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i in period %i', errorcode, it)
break
end

8
matlab/perfect-foresight-models/sim1_purely_forward.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, info] = sim1_purely_forward(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success] = sim1_purely_forward(endogenousvariables, exogenousvariables, steadystate, M, options)
% Performs deterministic simulation of a purely forward model
% Copyright © 2012-2023 Dynare Team
@ -33,7 +33,7 @@ function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse
sparse_rowval, sparse_colval, sparse_colptr, T);
end
info.status = true;
success = true;
for it = options.periods:-1:1
yf = endogenousvariables(:,it+1); % Values at next period, also used as guess value for current period
@ -51,7 +51,7 @@ for it = options.periods:-1:1
options.dynatol.x, options, ...
feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
if check
info.status = false;
success = false;
if options.debug
dprintf('sim1_purely_forward: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
end
@ -70,7 +70,7 @@ for it = options.periods:-1:1
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, dynamic_resid, dynamic_g1, yf, x, M.params, steadystate, M.dynamic_g1_sparse_rowval, M.dynamic_g1_sparse_colval, M.dynamic_g1_sparse_colptr);
if check
info.status = false;
success = false;
dprintf('sim1_purely_forward: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
break
end

8
matlab/perfect-foresight-models/sim1_purely_static.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, info] = sim1_purely_static(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success] = sim1_purely_static(endogenousvariables, exogenousvariables, steadystate, M, options)
% Performs deterministic simulation of a purely static model
@ -34,7 +34,7 @@ function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse
sparse_rowval, sparse_colval, sparse_colptr, T);
end
info.status = true;
success = true;
y = endogenousvariables(:,1);
@ -53,7 +53,7 @@ for it = 1:options.periods
options.dynatol.x, options, ...
feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
if check
info.status = false;
success = false;
if options.debug
dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
end
@ -72,7 +72,7 @@ for it = 1:options.periods
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, dynamic_resid, dynamic_g1, x, M.params, steadystate, M.dynamic_g1_sparse_rowval, M.dynamic_g1_sparse_colval, M.dynamic_g1_sparse_colptr);
if check
info.status = false;
success = false;
if options.debug
dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
end

38
matlab/perfect-foresight-models/solve_block_decomposed_problem.m
View File

@ -1,4 +1,4 @@
function oo_ = solve_block_decomposed_problem(options_, M_, oo_)
function [y, success, maxerror, per_block_status] = solve_block_decomposed_problem(options_, M_, oo_)
% Computes deterministic simulation with block option without bytecode
% Copyright © 2020-2023 Dynare Team
@ -41,23 +41,20 @@ end
y=oo_.endo_simul;
T=NaN(M_.block_structure.dyn_tmp_nbr, options_.periods+M_.maximum_lag+M_.maximum_lead);
oo_.deterministic_simulation.status = false;
maxerror = 0;
nblocks = length(M_.block_structure.block);
per_block_status = struct('success', cell(1, nblocks), 'error', cell(1, nblocks), 'iterations', cell(1, nblocks));
funcname = [ M_.fname '.dynamic'];
for blk = 1:length(M_.block_structure.block)
for blk = 1:nblocks
recursive_size = M_.block_structure.block(blk).endo_nbr - M_.block_structure.block(blk).mfs;
y_index = M_.block_structure.block(blk).variable((recursive_size+1):end);
fh_dynamic = str2func(sprintf('%s.sparse.block.dynamic_%d', M_.fname, blk));
if M_.block_structure.block(blk).Simulation_Type == 1 || ... % evaluateForward
M_.block_structure.block(blk).Simulation_Type == 2 % evaluateBackward
oo_.deterministic_simulation.status = true;
oo_.deterministic_simulation.error = 0;
oo_.deterministic_simulation.iterations = 0;
oo_.deterministic_simulation.block(blk).status = true;
oo_.deterministic_simulation.block(blk).error = 0;
oo_.deterministic_simulation.block(blk).iterations = 0;
if M_.block_structure.block(blk).Simulation_Type == 1
range = M_.maximum_lag+1:M_.maximum_lag+options_.periods;
else
@ -79,28 +76,33 @@ for blk = 1:length(M_.block_structure.block)
M_.block_structure.block(blk).g1_sparse_colptr, T(:, it_));
y(:, it_) = y3n(M_.endo_nbr+(1:M_.endo_nbr));
end
success = true;
maxblkerror = 0;
iter = [];
elseif M_.block_structure.block(blk).Simulation_Type == 3 || ... % solveForwardSimple
M_.block_structure.block(blk).Simulation_Type == 4 || ... % solveBackwardSimple
M_.block_structure.block(blk).Simulation_Type == 6 || ... % solveForwardComplete
M_.block_structure.block(blk).Simulation_Type == 7 % solveBackwardComplete
is_forward = M_.block_structure.block(blk).Simulation_Type == 3 || M_.block_structure.block(blk).Simulation_Type == 6;
[y, T, oo_] = solve_one_boundary(fh_dynamic, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, is_forward, true, false, M_, options_, oo_);
[y, T, success, maxblkerror, iter] = solve_one_boundary(fh_dynamic, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, is_forward, true, false, M_, options_);
elseif M_.block_structure.block(blk).Simulation_Type == 5 || ... % solveTwoBoundariesSimple
M_.block_structure.block(blk).Simulation_Type == 8 % solveTwoBoundariesComplete
[y, T, oo_] = solve_two_boundaries(fh_dynamic, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, options_, M_, oo_);
[y, T, success, maxblkerror, iter] = solve_two_boundaries(fh_dynamic, y, oo_.exo_simul, M_.params, oo_.steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, options_, M_);
end
tmp = y(M_.block_structure.block(blk).variable, :);
if any(isnan(tmp) | isinf(tmp))
disp(['Inf or Nan value during the resolution of block ' num2str(blk)]);
oo_.deterministic_simulation.status = false;
oo_.deterministic_simulation.error = 100;
oo_.deterministic_simulation.block(blk).status = false;
oo_.deterministic_simulation.block(blk).error = 100;
success = false;
end
if ~oo_.deterministic_simulation.status
per_block_status(blk).success = success;
per_block_status(blk).error = maxblkerror;
per_block_status(blk).iter = iter;
maxerror = max(maxblkerror, maxerror);
if ~success
return
end
end
oo_.endo_simul = y;

15
matlab/perfect-foresight-models/solve_stacked_linear_problem.m
View File

@ -1,6 +1,6 @@
function [endogenousvariables, info] = solve_stacked_linear_problem(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options)
function [endogenousvariables, success] = solve_stacked_linear_problem(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options)
% Copyright © 2015-2022 Dynare Team
% Copyright © 2015-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -60,11 +60,8 @@ end
endogenousvariables = [y0 bsxfun(@plus,reshape(y,M.endo_nbr,options.periods), steadystate_y) yT];
if check
info.status = false;
if options.debug
dprintf('solve_stacked_linear_problem: Nonlinear solver routine failed with errorcode=%i.', errorcode)
end
else
info.status = true;
success = ~check;
if ~success && options.debug
dprintf('solve_stacked_linear_problem: Nonlinear solver routine failed with errorcode=%i.', errorcode)
end

19
matlab/perfect-foresight-models/solve_stacked_problem.m
View File

@ -1,4 +1,4 @@
function [endogenousvariables, info, residuals] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options)
function [endogenousvariables, success, maxerror] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options)
% Solves the perfect foresight model using dynare_solve
%
@ -11,10 +11,10 @@ function [endogenousvariables, info, residuals] = solve_stacked_problem(endogeno
%
% OUTPUTS
% - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
% - info [struct] contains informations about the results.
% - residuals [double] N*T array, residuals of the equations (with 0 for initial condition)
% - success [logical] Whether a solution was found
% - maxerror [double] 1-norm of the residual
% Copyright © 2015-2022 Dynare Team
% Copyright © 2015-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -77,12 +77,9 @@ residuals(:, M.maximum_lag+(1:options.periods)) = reshape(res, M.endo_nbr, optio
if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity problem
residuals(eq_to_ignore,bsxfun(@le, endogenousvariables(eq_to_ignore,:), lb(eq_to_ignore)+eps) | bsxfun(@ge,endogenousvariables(eq_to_ignore,:),ub(eq_to_ignore)-eps))=0;
end
maxerror = max(max(abs(residuals)));
success = ~check;
if check
info.status = false;
if options.debug
dprintf('solve_stacked_problem: Nonlinear solver routine failed with errorcode=%i.', errorcode)
end
else
info.status = true;
if ~success && options.debug
dprintf('solve_stacked_problem: Nonlinear solver routine failed with errorcode=%i.', errorcode)
end

61
matlab/solve_one_boundary.m
View File

@ -1,5 +1,5 @@
function [y, T, oo_, info] = solve_one_boundary(fh, y, x, params, steady_state, T, ...
y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo, is_forward, is_dynamic, verbose, M, options, oo_)
function [y, T, success, max_res, iter] = solve_one_boundary(fh, y, x, params, steady_state, T, ...
y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo, is_forward, is_dynamic, verbose, M, options)
% Computes the deterministic simulation or the steady state for a block of equations containing
% only lags or only leads (but not both).
%
@ -25,17 +25,15 @@ function [y, T, oo_, info] = solve_one_boundary(fh, y, x, params, steady_state,
% stack_solve_algo [integer] linear solver method used in the Newton algorithm
% is_forward [logical] Whether the block has to be solved forward
% If false, the block is solved backward
% is_dynamic [logical] If true, this is a deterministic simulation
% and the oo_.deterministic_simulation field is updated.
% If false, this is a steady state computation
% (oo_.detereministic_simulation remains unchanged).
% is_dynamic [logical] Whether this is a deterministic simulation
% verbose [logical] Whether iterations are to be printed
%
% OUTPUTS
% y [matrix] All endogenous variables of the model
% T [matrix] Temporary terms
% info [integer] >=0 no error
% <0 error
% y [matrix] All endogenous variables of the model
% T [matrix] Temporary terms
% success [logical] Whether a solution was found
% max_res [double] ∞-norm of the residual
% iter [integer] Number of iterations
%
% ALGORITHM
% Newton with LU or GMRES or BicGstab for dynamic block
@ -135,15 +133,7 @@ for it_=start:incr:finish
if verbose
disp('The singularity of the jacobian matrix could not be corrected')
end
if is_dynamic
oo_.deterministic_simulation.status = false;
oo_.deterministic_simulation.error = max_res;
oo_.deterministic_simulation.iterations = iter;
oo_.deterministic_simulation.block(Block_Num).status = false;% Convergency failed.
oo_.deterministic_simulation.block(Block_Num).error = max_res;
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
end
info = -Block_Num*10;
success = false;
return
end
end
@ -167,15 +157,7 @@ for it_=start:incr:finish
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, it_, iter);
end
end
if is_dynamic
oo_.deterministic_simulation.status = false;
oo_.deterministic_simulation.error = max_res;
oo_.deterministic_simulation.iterations = iter;
oo_.deterministic_simulation.block(Block_Num).status = false;% Convergency failed.
oo_.deterministic_simulation.block(Block_Num).error = max_res;
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
end
info = -Block_Num*10;
success = false;
return
end
else
@ -296,30 +278,13 @@ for it_=start:incr:finish
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, it_,iter);
end
end
if is_dynamic
oo_.deterministic_simulation.status = false;
oo_.deterministic_simulation.error = max_res;
oo_.deterministic_simulation.iterations = iter;
oo_.deterministic_simulation.block(Block_Num).status = false;% Convergency failed.
oo_.deterministic_simulation.block(Block_Num).error = max_res;
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
end
info = -Block_Num*10;
success = false;
return
end
end
if is_dynamic
info = 1;
oo_.deterministic_simulation.status = true;
oo_.deterministic_simulation.error = max_res;
oo_.deterministic_simulation.iterations = iter;
oo_.deterministic_simulation.block(Block_Num).status = true;
oo_.deterministic_simulation.block(Block_Num).error = max_res;
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
else
info = 0;
end
success = true;
function y3n = dynendo(y, it_, M)
if it_ > 1 && it_ < size(y, 2)

36
matlab/solve_two_boundaries.m
View File

@ -1,4 +1,4 @@
function [y, T, oo]= solve_two_boundaries(fh, y, x, params, steady_state, T, y_index, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo,options,M, oo)
function [y, T, success, max_res, iter] = solve_two_boundaries(fh, y, x, params, steady_state, T, y_index, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, cutoff, stack_solve_algo,options,M)
% Computes the deterministic simulation of a block of equation containing
% both lead and lag variables using relaxation methods
%
@ -28,12 +28,13 @@ function [y, T, oo]= solve_two_boundaries(fh, y, x, params, steady_state, T, y_i
% - 3 BicGStab
% - 4 Optimal path length
% M [structure] Model description
% oo [structure] Results
%
% OUTPUTS
% y [matrix] All endogenous variables of the model
% T [matrix] Temporary terms
% oo [structure] Results
% success [logical] Whether a solution was found
% max_res [double] ∞-norm of the residual
% iter [integer] Number of iterations
%
% ALGORITHM
% Newton with LU or GMRES or BicGstab
@ -131,12 +132,7 @@ while ~(cvg || iter>maxit_)
continue
else
disp('The singularity of the jacobian matrix could not be corrected');
oo.deterministic_simulation.status = false;
oo.deterministic_simulation.error = max_res;
oo.deterministic_simulation.iterations = iter;
oo.deterministic_simulation.block(Block_Num).status = false;% Convergency failed.
oo.deterministic_simulation.block(Block_Num).error = max_res;
oo.deterministic_simulation.block(Block_Num).iterations = iter;
success = false;
return
end
end
@ -156,12 +152,7 @@ while ~(cvg || iter>maxit_)
fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, iter);
end
end
oo.deterministic_simulation.status = false;
oo.deterministic_simulation.error = max_res;
oo.deterministic_simulation.iterations = iter;
oo.deterministic_simulation.block(Block_Num).status = false;% Convergency failed.
oo.deterministic_simulation.block(Block_Num).error = max_res;
oo.deterministic_simulation.block(Block_Num).iterations = iter;
success = false;
return
end
else
@ -333,21 +324,12 @@ if (iter>maxit_)
printline(41)
%disp(['No convergence after ' num2str(iter,'%4d') ' iterations in Block ' num2str(Block_Num,'%d')])
end
oo.deterministic_simulation.status = false;
oo.deterministic_simulation.error = max_res;
oo.deterministic_simulation.iterations = iter;
oo.deterministic_simulation.block(Block_Num).status = false;% Convergency failed.
oo.deterministic_simulation.block(Block_Num).error = max_res;
oo.deterministic_simulation.block(Block_Num).iterations = iter;
success = false;
return
end
oo.deterministic_simulation.status = true;
oo.deterministic_simulation.error = max_res;
oo.deterministic_simulation.iterations = iter;
oo.deterministic_simulation.block(Block_Num).status = true;% Convergency obtained.
oo.deterministic_simulation.block(Block_Num).error = max_res;
oo.deterministic_simulation.block(Block_Num).iterations = iter;
success = true;
function y3n = dynendo(y, it_, M)
y3n = reshape(y(:, it_+(-1:1)), 3*M.endo_nbr, 1);