From 7722e8e36b1318584dd5fb44846d5617e7278dda Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Tue, 6 Jun 2023 18:13:10 +0200 Subject: [PATCH] Perfect foresight: inner functions no longer return a modified oo_ MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit They now only return what’s really their output (simulated paths, maximum residual error…). This is a move towards a more functional programming style. --- matlab/ep/extended_path_core.m | 19 ++-- matlab/ep/extended_path_homotopy.m | 20 ++--- matlab/evaluate_steady_state.m | 10 +-- .../perfect_foresight_solver.m | 37 +++++--- .../perfect_foresight_solver_core.m | 87 +++++++++---------- matlab/perfect-foresight-models/sim1.m | 23 ++--- matlab/perfect-foresight-models/sim1_lbj.m | 11 +-- matlab/perfect-foresight-models/sim1_linear.m | 20 ++--- .../sim1_purely_backward.m | 8 +- .../sim1_purely_forward.m | 8 +- .../sim1_purely_static.m | 8 +- .../solve_block_decomposed_problem.m | 38 ++++---- .../solve_stacked_linear_problem.m | 15 ++-- .../solve_stacked_problem.m | 19 ++-- matlab/solve_one_boundary.m | 61 +++---------- matlab/solve_two_boundaries.m | 36 ++------ 16 files changed, 171 insertions(+), 249 deletions(-) diff --git a/matlab/ep/extended_path_core.m b/matlab/ep/extended_path_core.m index 945b008cb..6cc5f8db9 100644 --- a/matlab/ep/extended_path_core.m +++ b/matlab/ep/extended_path_core.m @@ -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; diff --git a/matlab/ep/extended_path_homotopy.m b/matlab/ep/extended_path_homotopy.m index dca98a178..77fbf0e9f 100644 --- a/matlab/ep/extended_path_homotopy.m +++ b/matlab/ep/extended_path_homotopy.m @@ -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 diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m index 74a8284b0..dd3a368fe 100644 --- a/matlab/evaluate_steady_state.m +++ b/matlab/evaluate_steady_state.m @@ -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 diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m index b15135af2..bfbf63b0b 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m @@ -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 diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m index ecb3abba3..35db87b4e 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m @@ -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 diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m index 8f0f65a90..f9dafab9a 100644 --- a/matlab/perfect-foresight-models/sim1.m +++ b/matlab/perfect-foresight-models/sim1.m @@ -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 diff --git a/matlab/perfect-foresight-models/sim1_lbj.m b/matlab/perfect-foresight-models/sim1_lbj.m index 7ffc6eb0e..ce3e582f4 100644 --- a/matlab/perfect-foresight-models/sim1_lbj.m +++ b/matlab/perfect-foresight-models/sim1_lbj.m @@ -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 diff --git a/matlab/perfect-foresight-models/sim1_linear.m b/matlab/perfect-foresight-models/sim1_linear.m index 3bffa6d42..f6a67759b 100644 --- a/matlab/perfect-foresight-models/sim1_linear.m +++ b/matlab/perfect-foresight-models/sim1_linear.m @@ -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 diff --git a/matlab/perfect-foresight-models/sim1_purely_backward.m b/matlab/perfect-foresight-models/sim1_purely_backward.m index 35a4f9ff4..904b6e647 100644 --- a/matlab/perfect-foresight-models/sim1_purely_backward.m +++ b/matlab/perfect-foresight-models/sim1_purely_backward.m @@ -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 diff --git a/matlab/perfect-foresight-models/sim1_purely_forward.m b/matlab/perfect-foresight-models/sim1_purely_forward.m index 888d68d10..58778ec57 100644 --- a/matlab/perfect-foresight-models/sim1_purely_forward.m +++ b/matlab/perfect-foresight-models/sim1_purely_forward.m @@ -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 diff --git a/matlab/perfect-foresight-models/sim1_purely_static.m b/matlab/perfect-foresight-models/sim1_purely_static.m index 0c1774670..76913d740 100644 --- a/matlab/perfect-foresight-models/sim1_purely_static.m +++ b/matlab/perfect-foresight-models/sim1_purely_static.m @@ -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 diff --git a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m index 85989f652..ed6637f11 100644 --- a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m +++ b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m @@ -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; diff --git a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m index 80c876d13..c937fdfa8 100644 --- a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m +++ b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m @@ -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 diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m index b93ec7614..f9c2c4fcb 100644 --- a/matlab/perfect-foresight-models/solve_stacked_problem.m +++ b/matlab/perfect-foresight-models/solve_stacked_problem.m @@ -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 diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m index 932d50ff8..2c265e109 100644 --- a/matlab/solve_one_boundary.m +++ b/matlab/solve_one_boundary.m @@ -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) diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m index 7ad7a4f40..248b87d87 100644 --- a/matlab/solve_two_boundaries.m +++ b/matlab/solve_two_boundaries.m @@ -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);