From d5a3a8e16af56af23c8eb27a693b2e67579bfde0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 15 Jun 2023 16:59:37 +0200
Subject: [PATCH] Various improvements to perfect foresight homotopy
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
– new option “endval_steady” to pf_setup command to recompute terminal
steady state in the homotopy loop
– new options “homotopy_linearization_fallback” and
“homotopy_marginal_linearization_fallback” to pf_solver and pfwee_solver
commands, to get an approximate solution when homotopy fails to go to 100%
– new options “homotopy_initial_step_size”, “homotopy_min_step_size”,
“homotopy_step_size_increase_success_count” and “homotopy_max_completion_share”
to pf_solver and pfwee_solver commands to fine tune the homotopy behaviour
– removed option “homotopy_alt_starting_point” to pf_solver command, not really
useful
– new options “steady_solve_algo”, “steady_tolf”, “steady_tolx”,
“steady_maxit”, “steady_markowitz” to pf_solver and pfwee_solver commands, to
control the computation of the terminal steady state (and remove the
equivalent options which previously had different names in pfwee_solver command)
---
doc/manual/source/the-model-file.rst | 164 ++++---
matlab/default_option_values.m | 19 +-
matlab/evaluate_steady_state.m | 8 +-
.../perfect_foresight_setup.m | 8 +-
.../perfect_foresight_solver.m | 416 +++++++++++++-----
.../perfect_foresight_solver_core.m | 6 -
..._foresight_with_expectation_errors_setup.m | 15 -
...foresight_with_expectation_errors_solver.m | 37 +-
preprocessor | 2 +-
tests/Makefile.am | 4 +-
tests/deterministic_simulations/homotopy.mod | 7 +-
.../homotopy_alt_starting_point.mod | 50 ---
.../homotopy_endval_steady_linearization.mod | 43 ++
.../homotopy_histval.mod | 17 +-
.../homotopy_marginal_linearization.mod | 42 ++
tests/deterministic_simulations/pfwee.mod | 7 +-
.../pfwee_constant_sim_length.mod | 13 +-
.../pfwee_homotopy.mod | 40 ++
.../pfwee_learnt_in.mod | 7 +-
19 files changed, 628 insertions(+), 277 deletions(-)
delete mode 100644 tests/deterministic_simulations/homotopy_alt_starting_point.mod
create mode 100644 tests/deterministic_simulations/homotopy_endval_steady_linearization.mod
create mode 100644 tests/deterministic_simulations/homotopy_marginal_linearization.mod
create mode 100644 tests/deterministic_simulations/pfwee_homotopy.mod
diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index fc1b530553..999e689b21 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -3701,21 +3701,69 @@ speed-up on large models.
This option tells Dynare to not try a homotopy technique (as described
above) if the problem cannot be solved directly.
- .. option:: homotopy_alt_starting_point
-
- When the homotopy technique is tried (as described above), Dynare first
- tries to reduce the size of the shock in order to get a successful
- simulation on which to build upon for simulating a shock of the true
- size. However, if an ``endval`` block is present (*i.e.* if the terminal
- state differs from the initial state), there are two ways of reducing
- the size of the shock. By default, Dynare will perform this reduction by
- computing a simulation whose initial state is closer to the target
- terminal state; in other words, it will implicitly modify the contents
- of the ``initval`` and ``shocks`` blocks to make them closer to the the
- contents of the ``endval`` block. If this option is set, Dynare will do
- the opposite: it will implicitly modify the contents of the ``endval``
- and ``shocks`` blocks to make them closer to the contents of the
- ``initval`` block.
+ .. option:: homotopy_initial_step_size = DOUBLE
+
+ Specifies which share of the shock should be applied in the first
+ iteration of the homotopy procedure. This option is useful when it is
+ known that immediately trying 100% of the shock will fail, so as to save
+ computing time. Must be between ``0`` and ``1``. Default: ``1``.
+
+ .. option:: homotopy_min_step_size = DOUBLE
+
+ The homotopy procedure halves the size of the step whenever there is
+ a failure. This option specifies the minimum step size under which the
+ homotopy procedure is considered to have failed. Default: ``0.001``.
+
+ .. option:: homotopy_step_size_increase_success_count = INTEGER
+
+ Specifies after how many consecutive successful iterations the homotopy
+ procedure should double the size of the step. A zero value means that
+ the step size should never be increased. Default: ``3``.
+
+ .. option:: homotopy_linearization_fallback
+
+ Whenever the homotopy procedure is not able to find a solution for 100%
+ of the shock, but is able to find one for a smaller share, instructs
+ Dynare to compute an approximate solution by rescaling the solution
+ obtained for a fraction of the shock, as if the reaction of the model to
+ the shock was a linear function of the size of that shock. More
+ formally, if :math:`s` is the share of the shock applied (between
+ :math:`0` and :math:`1`), :math:`y(s)` is the value of a given
+ endogenous variable at a given period as a function of :math:`s` (in
+ particular, :math:`y(1)` corresponds to the exact solution of the
+ problem), and :math:`s^*` is the greatest share of the shock for which
+ the homotopy procedure has been able to find a solution, then the approximate
+ solution returned is :math:`\frac{y(s^*)-y(0)}{s^*}`.
+
+ .. option:: homotopy_marginal_linearization_fallback [= DOUBLE]
+
+ Whenever the homotopy procedure is not able to find a solution for 100%
+ of the shock, but is able to find one for a smaller share, instructs
+ Dynare to compute an approximate solution obtained by rescaling the
+ solution obtained for a fraction of the shock, obtained as if the
+ reaction of the model to the shock was, at the margin, a linear function
+ of the size of that shock. More formally, if :math:`s` is the share of
+ the shock applied (between :math:`0` and :math:`1`), :math:`y(s)` is the
+ value of a given endogenous variable at a given period as a function of
+ :math:`s` (in particular, :math:`y(1)` corresponds to the exact solution
+ of the problem), :math:`s^*` is the greatest share of the shock for
+ which the homotopy procedure has been able to find a solution, and
+ :math:`\epsilon` is a small step size, then the approximate solution
+ returned is
+ :math:`y(s^*)+(1-s^*)\frac{y(s^*)-y(s^*-\epsilon)}{\epsilon}`. The value
+ of :math:`\epsilon` is ``0.01`` by default, but can be modified by
+ passing some other value to the option.
+
+ .. option:: homotopy_max_completion_share = DOUBLE
+
+ Instructs Dynare, within the homotopy procedure, to not try to compute
+ the solution for a greater share than the one given as the option value.
+ This option only makes sense when used in conjunction with either the
+ ``homotopy_linearization_fallback`` or the
+ ``homotopy_marginal_linearization_fallback`` option. It is typically
+ used in situations where it is known that homotopy will fail to go
+ beyond a certain point, so as to save computing time, while at the same
+ time getting an approximate solution.
.. option:: markowitz = DOUBLE
@@ -3827,6 +3875,50 @@ speed-up on large models.
``endval`` or a subsequent ``steady``. Only available with option
``stack_solve_algo==0`` or ``stack_solve_algo==7``.
+ .. option:: endval_steady
+
+ In scenarios with a permanent shock, specifies that the terminal
+ condition is a steady state, even if the ``steady`` command has not been
+ called after the ``endval`` block. As a consequence, the
+ ``perfect_foresight_solver`` command will compute the terminal steady
+ state itself (given the value of the exogenous variables given in the
+ ``endval`` block). In practice, this option is useful when the permanent
+ shock is very large, in which case the homotopy procedure inside
+ ``perfect_foresight_solver`` will find both the terminal steady state
+ and the transitional dynamics within the same loop (which is less costly
+ than first computing the terminal steady state by homotopy, then
+ computing the transitional dynamics by homotopy).
+
+ .. option:: steady_solve_algo = INTEGER
+
+ See :ref:`solve_algo <solvalg>`. Used when computing the terminal steady
+ state when option ``endval_steady`` has been specified to the
+ ``perfect_foresight_setup`` command.
+
+ .. option:: steady_tolf = DOUBLE
+
+ See :ref:`tolf <steady_tolf>`. Used when computing the terminal steady
+ state when option ``endval_steady`` has been specified to the
+ ``perfect_foresight_setup`` command.
+
+ .. option:: steady_tolx = DOUBLE
+
+ See :ref:`tolx <steady_tolx>`. Used when computing the terminal steady
+ state when option ``endval_steady`` has been specified to the
+ ``perfect_foresight_setup`` command.
+
+ .. option:: steady_maxit = INTEGER
+
+ See :ref:`maxit <steady_maxit>`. Used when computing the terminal steady
+ state when option ``endval_steady`` has been specified to the
+ ``perfect_foresight_setup`` command.
+
+ .. option:: steady_markowitz = DOUBLE
+
+ See :ref:`markowitz <steady_markowitz>`. Used when computing the
+ terminal steady state when option ``endval_steady`` has been specified
+ to the ``perfect_foresight_setup`` command.
+
*Output*
The simulated endogenous variables are available in global matrix
@@ -4019,7 +4111,9 @@ and ``endval`` blocks which are given a special ``learnt_in`` option.
always a steady state. Hence, it will recompute the terminal steady state
as many times as the anticipation about the terminal condition changes. In
particular, the information about endogenous variables that may be given in
- the ``endval`` block is ignored.
+ the ``endval`` block is ignored. Said otherwise, the equivalent of option
+ ``endval_steady`` of the ``perfect_foresight_setup`` command is always
+ implicitly enabled.
*Options*
@@ -4057,33 +4151,12 @@ and ``endval`` blocks which are given a special ``learnt_in`` option.
the name of the variable on the first line and ``s`` on the second
line. Of course, values in cells corresponding to ``t<s`` are ignored.
- .. option:: solve_algo = INTEGER
-
- See :ref:`solve_algo <solvalg>`. Used when computing the terminal steady state.
-
- .. option:: tolf = DOUBLE
-
- See :ref:`tolf <steady_tolf>`. Used when computing the terminal steady state.
-
- .. option:: tolx = DOUBLE
-
- See :ref:`tolx <steady_tolx>`. Used when computing the terminal steady state.
-
- .. option:: maxit = INTEGER
-
- See :ref:`maxit <steady_maxit>`. Used when computing the terminal steady state.
-
- .. option:: markowitz = DOUBLE
-
- See :ref:`markowitz <steady_markowitz>`. Used when computing the terminal steady state.
-
*Output*
``oo_.exo_simul`` and ``oo_.endo_simul`` are initialized before the
simulation. Temporary shocks are stored in ``oo_.pfwee.shocks_info``,
terminal conditions for exogenous variables are stored in
- ``oo_.pfwee.terminal_info``, and terminal steady states are stored in
- ``oo_.pfwee.terminal_steady_state``.
+ ``oo_.pfwee.terminal_info``.
*Example*
@@ -4189,7 +4262,9 @@ and ``endval`` blocks which are given a special ``learnt_in`` option.
*Output*
The simulated paths of endogenous variables are available in
- ``oo_.endo_simul``.
+ ``oo_.endo_simul``. The terminal steady state values corresponding to the
+ last period of the information set are available in ``oo_.steady_state``
+ and ``oo_.exo_steady_state``.
.. matvar:: oo_.pfwee.shocks_info
@@ -4212,17 +4287,6 @@ and ``endval`` blocks which are given a special ``learnt_in`` option.
the terminal condition for exogenous indexed ``k``, as anticipated from
period ``s``, is stored in ``oo_.pfwee.terminal_info(k,s)``.
-.. matvar:: oo_.pfwee.terminal_steady_state
-
- |br| This variable stores the terminal steady states for endogenous
- variables used during perfect foresight simulations with expectation
- errors, after :comm:`perfect_foresight_with_expectation_errors_setup` has
- been run. It is a matrix, whose lines correspond to endogenous variables
- (in declaration order), and whose columns correspond to informational time.
- In other words, the terminal steady state for endogenous indexed ``k``, as
- anticipated from period ``s``, is stored in
- ``oo_.pfwee.terminal_steady_state(k,s)``.
-
.. _stoch-sol:
Stochastic solution and simulation
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 52bd5b6fea..c0a0bf67ec 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -49,6 +49,7 @@ options_.schur_vec_tol = 1e-11; % used to find nonstationary variables in Schur
options_.qz_criterium = [];
options_.qz_zero_threshold = 1e-6;
options_.lyapunov_complex_threshold = 1e-15;
+options_.solve_algo = 4;
options_.solve_tolf = eps^(1/3);
options_.solve_tolx = eps^(2/3);
options_.trust_region_initial_step_bound_factor = 1;
@@ -326,7 +327,22 @@ options_.markowitz = 0.5;
options_.minimal_solving_periods = 1;
options_.endogenous_terminal_period = false;
options_.no_homotopy = false;
-options_.homotopy_alt_starting_point = false;
+options_.simul.endval_steady = false;
+
+options_.simul.homotopy_max_completion_share = 1;
+options_.simul.homotopy_min_step_size = 1e-3;
+options_.simul.homotopy_step_size_increase_success_count = 3;
+options_.simul.homotopy_initial_step_size = 1;
+options_.simul.homotopy_linearization_fallback = false;
+options_.simul.homotopy_marginal_linearization_fallback = 0; % Size of the step used for the marginal linearization; 0 means disabled
+
+% Options used by perfect_foresight_* commands when they compute the steady
+% state corresponding to a terminal condition
+options_.simul.steady_solve_algo = options_.solve_algo;
+options_.simul.steady_maxit = options_.steady.maxit;
+options_.simul.steady_tolf = options_.solve_tolf;
+options_.simul.steady_tolx = options_.solve_tolx;
+options_.simul.steady_markowitz = options_.markowitz;
% Perfect foresight with expectation errors
options_.pfwee.constant_simulation_length = false;
@@ -334,7 +350,6 @@ options_.pfwee.constant_simulation_length = false;
% Solution
options_.order = 2;
options_.pruning = false;
-options_.solve_algo = 4;
options_.replic = 50;
options_.simul_replic = 1;
options_.drop = 100;
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index 4de31803f9..c23f562adc 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -394,7 +394,9 @@ elseif options.bytecode
ys = bytecode('static', ys_init, exo_ss, params);
end
catch ME
- disp(ME.message);
+ if options.verbosity >= 1
+ disp(ME.message);
+ end
ys = ys_init;
check = 1;
end
@@ -420,7 +422,9 @@ elseif options.bytecode
[~, ~, ys, T] = bytecode(ys, exo_ss, params, ys, 1, ys, T, 'evaluate', 'static', ...
'block_decomposed', ['block = ' int2str(b)]);
catch ME
- disp(ME.message);
+ if options.verbosity >= 1
+ disp(ME.message);
+ end
check = 1;
break
end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_setup.m b/matlab/perfect-foresight-models/perfect_foresight_setup.m
index db0e04fbf5..f187b9dc2e 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_setup.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_setup.m
@@ -12,7 +12,7 @@ function perfect_foresight_setup()
% SPECIAL REQUIREMENTS
% none
-% Copyright © 1996-2021 Dynare Team
+% Copyright © 1996-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -64,5 +64,9 @@ if ~isempty(M_.det_shocks) && options_.periods<max([M_.det_shocks.periods])
error('PERFECT_FORESIGHT_SETUP: Please check the declaration of the shocks or increase the value of the periods option.')
end
+if options_.simul.endval_steady && M_.maximum_lead == 0
+ error('PERFECT_FORESIGHT_SETUP: Option endval_steady cannot be used on a purely backward or static model.')
+end
+
oo_ = make_ex_(M_,options_,oo_);
-oo_ = make_y_(M_,options_,oo_);
\ No newline at end of file
+oo_ = make_y_(M_,options_,oo_);
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index bfbf63b0b9..d412afb301 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -51,162 +51,318 @@ if options_.debug
end
end
-initperiods = 1:M_.maximum_lag;
-lastperiods = (M_.maximum_lag+periods+1):(M_.maximum_lag+periods+M_.maximum_lead);
+% Various sanity checks
+
+if options_.no_homotopy && (options_.simul.homotopy_initial_step_size ~= 1 ...
+ || options_.simul.homotopy_max_completion_share ~= 1 ...
+ || options_.simul.homotopy_linearization_fallback ...
+ || options_.simul.homotopy_marginal_linearization_fallback ~= 0)
+ error('perfect_foresight_solver: the no_homotopy option is incompatible with homotopy_initial_step_size, homotopy_max_completion_share, homotopy_linearization_fallback and homotopy_marginal_linearization_fallback options')
+end
-[oo_.endo_simul, success, maxerror, solver_iter, per_block_status] = perfect_foresight_solver_core(M_,options_,oo_);
+if options_.simul.homotopy_initial_step_size > 1 || options_.simul.homotopy_initial_step_size < 0
+ error('perfect_foresight_solver: The value given to homotopy_initial_step_size option must be in the [0,1] interval')
+end
-% If simulation failed try homotopy.
-if ~success && ~options_.no_homotopy
+if options_.simul.homotopy_min_step_size > 1 || options_.simul.homotopy_min_step_size < 0
+ error('perfect_foresight_solver: The value given to homotopy_min_step_size option must be in the [0,1] interval')
+end
- if ~options_.noprint
- fprintf('\nSimulation of the perfect foresight model failed!\n')
- fprintf('Switching to a homotopy method...\n')
+if options_.simul.homotopy_max_completion_share > 1 || options_.simul.homotopy_max_completion_share < 0
+ error('perfect_foresight_solver: The value given to homotopy_max_completion_share option must be in the [0,1] interval')
+end
+
+if options_.simul.homotopy_marginal_linearization_fallback > 1 || options_.simul.homotopy_marginal_linearization_fallback < 0
+ error('perfect_foresight_solver: The value given to homotopy_marginal_linearization_fallback option must be in the [0,1] interval')
+end
+
+if options_.simul.homotopy_initial_step_size < options_.simul.homotopy_min_step_size
+ error('perfect_foresight_solver: The value given to homotopy_initial_step_size option must be greater or equal to that given to homotopy_min_step_size option')
+end
+
+if options_.simul.homotopy_linearization_fallback && options_.simul.homotopy_marginal_linearization_fallback > 0
+ error('perfect_foresight_solver: Options homotopy_linearization_fallback and homotopy_marginal_linearization_fallback cannot be used together')
+end
+
+if options_.simul.homotopy_max_completion_share < 1 && ~options_.simul.homotopy_linearization_fallback && options_.simul.homotopy_marginal_linearization_fallback == 0
+ error('perfect_foresight_solver: Option homotopy_max_completion_share has a value less than 1, so you must also specify either homotopy_linearization_fallback or homotopy_marginal_linearization_fallback')
+end
+
+initperiods = 1:M_.maximum_lag;
+simperiods = M_.maximum_lag+(1:periods);
+lastperiods = M_.maximum_lag+periods+(1:M_.maximum_lead);
+
+% Create base scenario for homotopy, which corresponds to the initial steady
+% state (i.e. a known solution to the perfect foresight problem, assuming that
+% oo_.steady_state/ys0_ effectively contains a steady state)
+if isempty(ys0_)
+ endobase = repmat(oo_.steady_state, 1,M_.maximum_lag+periods+M_.maximum_lead);
+ exobase = repmat(oo_.exo_steady_state',M_.maximum_lag+periods+M_.maximum_lead,1);
+else
+ endobase = repmat(ys0_, 1, M_.maximum_lag+periods+M_.maximum_lead);
+ exobase = repmat(ex0_', M_.maximum_lag+periods+M_.maximum_lead, 1);
+end
+
+% Determine whether the terminal condition is not a steady state (typically
+% because steady was not called after endval)
+if ~options_.simul.endval_steady && ~isempty(ys0_)
+ terminal_condition_is_a_steady_state = true;
+ for j = lastperiods
+ endval_resid = evaluate_static_model(oo_.endo_simul(:,j), oo_.exo_simul(j,:), M_.params, M_, options_);
+ if norm(endval_resid, 'Inf') > options_.simul.steady_tolf
+ terminal_condition_is_a_steady_state = false;
+ break
+ end
end
+end
- % Disable warnings if homotopy
- warning_old_state = warning;
- warning off all
- % Do not print anything
- oldverbositylevel = options_.verbosity;
- options_.verbosity = 0;
+% Copy the paths for the exogenous and endogenous variables, as given by perfect_foresight_setup
+exoorig = oo_.exo_simul;
+endoorig = oo_.endo_simul;
+
+current_share = 0; % Share of shock successfully completed so far
+step = min(options_.simul.homotopy_initial_step_size, options_.simul.homotopy_max_completion_share);
+success_counter = 0;
+iteration = 0;
+
+function local_success = create_scenario(share)
+ % For a given share, updates the exogenous path and also the initial and
+ % terminal conditions for the endogenous path (but do not modify the initial
+ % guess for endogenous)
+
+ % Compute convex combination for the path of exogenous
+ oo_.exo_simul = exoorig*share + exobase*(1-share);
+
+ % Compute convex combination for the initial condition
+ % In most cases, the initial condition is a steady state and this does nothing
+ % This is for cases when the initial condition is out of equilibrium
+ oo_.endo_simul(:, initperiods) = share*endoorig(:, initperiods)+(1-share)*endobase(:, initperiods);
+
+ % If there is a permanent shock, ensure that the rescaled terminal condition is
+ % a steady state (if the user asked for this recomputation, or if the original
+ % terminal condition is a steady state)
+ local_success = true;
+ if options_.simul.endval_steady || (~isempty(ys0_) && terminal_condition_is_a_steady_state)
+
+ % Set “local” options for steady state computation (after saving the global values)
+ saved_steady_solve_algo = options_.solve_algo;
+ options_.solve_algo = options_.simul.steady_solve_algo;
+ saved_steady_maxit = options_.steady.maxit;
+ options_.steady.maxit = options_.simul.steady_maxit;
+ saved_steady_tolf = options_.solve_tolf;
+ options_.solve_tolf = options_.simul.steady_tolf;
+ saved_steady_tolx = options_.solve_tolx;
+ options_.solve_tolx = options_.simul.steady_tolx;
+ saved_steady_markowitz = options_.markowitz;
+ options_.markowitz = options_.simul.steady_markowitz;
+
+ saved_ss = oo_.endo_simul(:, lastperiods);
+ % Effectively compute the terminal steady state
+ for j = lastperiods
+ % First use the terminal steady of the previous homotopy iteration as guess value (or the contents of the endval block if this is the first iteration)
+ [oo_.endo_simul(:, j), ~, info] = evaluate_steady_state(oo_.endo_simul(:, j), oo_.exo_simul(j, :), M_, options_, true);
+ if info(1)
+ % If this fails, then try again using the initial steady state as guess value
+ if isempty(ys0_)
+ guess_value = oo_.steady_state;
+ else
+ guess_value = ys0_;
+ end
+ [oo_.endo_simul(:, j), ~, info] = evaluate_steady_state(guess_value, oo_.exo_simul(j, :), M_, options_, true);
+ if info(1)
+ % If this fails again, give up and restore last periods in oo_.endo_simul
+ oo_.endo_simul(:, lastperiods) = saved_ss;
+ local_success = false;
+ break;
+ end
+ end
+ end
- % Set initial paths for the endogenous and exogenous variables.
- if ~options_.homotopy_alt_starting_point
- endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+periods+M_.maximum_lead);
- exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+periods+M_.maximum_lead,1);
+ % The following is needed for the STEADY_STATE() operator to work properly
+ oo_.steady_state = oo_.endo_simul(:, end);
+
+ options_.solve_algo = saved_steady_solve_algo;
+ options_.steady.maxit = saved_steady_maxit;
+ options_.solve_tolf = saved_steady_tolf;
+ options_.solve_tolx = saved_steady_tolx;
+ options_.markowitz = saved_steady_markowitz;
else
- if isempty(ys0_) || isempty(ex0_)
- error('The homotopy_alt_starting_point option cannot be used without an endval block');
- end
- endoinit = repmat(ys0_, 1, M_.maximum_lag+periods+M_.maximum_lead);
- exoinit = repmat(ex0_', M_.maximum_lag+periods+M_.maximum_lead, 1);
+ % The terminal condition is not a steady state, compute a convex combination
+ oo_.endo_simul(:, lastperiods) = share*endoorig(:, lastperiods)+(1-share)*endobase(:, lastperiods);
end
+end
- % Copy the current paths for the exogenous and endogenous variables.
- exosim = oo_.exo_simul;
- endosim = oo_.endo_simul;
+while step > options_.simul.homotopy_min_step_size
- current_weight = 0; % Current weight of target point in convex combination.
- step = .5; % Set default step size.
- success_counter = 0;
- iteration = 0;
+ iteration = iteration+1;
- if ~options_.noprint
- fprintf('Iter. \t | Lambda \t | status \t | Max. residual\n')
- fprintf('++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n')
+ saved_endo_simul = oo_.endo_simul;
+
+ new_share = current_share + step; % Try this share, and see if it succeeds
+
+ if new_share > options_.simul.homotopy_max_completion_share
+ new_share = options_.simul.homotopy_max_completion_share; % Don't go beyond target point
+ step = new_share - current_share;
end
- while (step > options_.dynatol.x)
- if ~isequal(step,1)
- options_.verbosity = 0;
+ steady_success = create_scenario(new_share);
+
+ if steady_success
+ % At the first iteration, use the initial guess given by
+ % perfect_foresight_setup or the user (but only if new_share=1, otherwise it
+ % does not make much sense). Afterwards, until a converging iteration has been obtained,
+ % use the rescaled terminal condition (or, if there is no lead, the base
+ % scenario / initial steady state).
+ if current_share == 0
+ if iteration == 1 && new_share == 1
+ oo_.endo_simul(:, simperiods) = endoorig(:, simperiods);
+ elseif M_.maximum_lead > 0
+ oo_.endo_simul(:, simperiods) = repmat(oo_.endo_simul(:, lastperiods(1)), 1, options_.periods);
+ else
+ oo_.endo_simul(:, simperiods) = endobase(:, simperiods);
+ end
end
- iteration = iteration+1;
- new_weight = current_weight + step; % Try this weight, and see if it succeeds
+ % Solve for the paths of the endogenous variables.
+ [oo_.endo_simul, success, maxerror, solver_iter, per_block_status] = perfect_foresight_solver_core(M_, options_, oo_);
+ else
+ success = false;
+ maxerror = NaN;
+ solver_iter = [];
+ per_block_status = [];
+ end
- if new_weight >= 1
- new_weight = 1; % Don't go beyond target point
- step = new_weight - current_weight;
+ if options_.no_homotopy || (iteration == 1 && success && new_share == 1)
+ % Skip homotopy
+ if success
+ current_share = new_share;
end
+ did_homotopy = false;
+ break
+ end
- % Compute convex combination for exo path and initial/terminal endo conditions
- % But take care of not overwriting the computed part of oo_.endo_simul
- oo_.exo_simul = exosim*new_weight + exoinit*(1-new_weight);
- oo_.endo_simul(:,[initperiods, lastperiods]) = new_weight*endosim(:,[initperiods, lastperiods])+(1-new_weight)*endoinit(:,[initperiods, lastperiods]);
-
- % Detect Nans or complex numbers in the solution.
- path_with_nans = any(any(isnan(oo_.endo_simul)));
- path_with_cplx = any(any(~isreal(oo_.endo_simul)));
-
- if isequal(iteration, 1)
- % First iteration, same initial guess as in the first call to perfect_foresight_solver_core routine.
- oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = endoinit(:,1:periods);
- elseif path_with_nans || path_with_cplx
- % If solver failed with NaNs or complex number, use previous solution as an initial guess.
- oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = saved_endo_simul(:,1+M_.maximum_lag:end-M_.maximum_lead);
- end
+ if iteration == 1
+ % First iteration failed, so we enter homotopy
+ did_homotopy = true;
- % Make a copy of the paths.
- saved_endo_simul = oo_.endo_simul;
+ if ~options_.noprint
+ fprintf('\nEntering the homotopy method iterations...\n')
+ fprintf('\nIter. \t | Share \t | Status \t | Max. residual\n')
+ fprintf('++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n')
+ end
- % Solve for the paths of the endogenous variables.
- [oo_.endo_simul, success, maxerror, solver_iter, per_block_status] = perfect_foresight_solver_core(M_,options_,oo_);
+ % Disable warnings if homotopy
+ warning_old_state = warning;
+ warning off all
+ % Do not print anything
+ oldverbositylevel = options_.verbosity;
+ options_.verbosity = 0;
+ end
- 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', maxerror)
- end
- break
- end
- success_counter = success_counter + 1;
- if success_counter >= 3
- success_counter = 0;
- step = step * 2;
- end
- if ~options_.noprint
- fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', maxerror)
- end
- else
- % If solver failed, then go back.
- oo_.endo_simul = saved_endo_simul;
+ if success
+ % Successful step
+ if ~options_.noprint
+ fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_share, 'succeeded', maxerror)
+ end
+ current_share = new_share;
+ if current_share >= options_.simul.homotopy_max_completion_share
+ break
+ end
+ success_counter = success_counter + 1;
+ if options_.simul.homotopy_step_size_increase_success_count > 0 ...
+ && success_counter >= options_.simul.homotopy_step_size_increase_success_count
success_counter = 0;
- step = step / 2;
- if ~options_.noprint
- 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
+ step = step * 2;
+ end
+ else
+ oo_.endo_simul = saved_endo_simul;
+ success_counter = 0;
+ step = step / 2;
+ if ~options_.noprint
+ if ~steady_success
+ fprintf('%i \t | %1.5f \t | %s\n', iteration, new_share, 'failed (in endval steady)')
+ elseif isreal(maxerror)
+ fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_share, 'failed', maxerror)
+ else
+ fprintf('%i \t | %1.5f \t | %s \t | %s\n', iteration, new_share, 'failed', 'Complex')
end
end
end
- if ~options_.noprint
- fprintf('++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n')
- end
- options_.verbosity = oldverbositylevel;
- warning(warning_old_state);
+end
+
+if did_homotopy && ~options_.noprint
+ fprintf('++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n')
end
%If simulated paths are complex, take real part and recompute the residuals to check whether this is actually a solution
if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_foresight_problem DLL
- ny = size(oo_.endo_simul, 1);
- if M_.maximum_lag > 0
- y0 = real(oo_.endo_simul(:, M_.maximum_lag));
+ real_simul = real(oo_.endo_simul);
+ real_maxerror = recompute_maxerror(real_simul, oo_, M_, options_);
+ if real_maxerror <= options_.dynatol.f
+ oo_.endo_simul = real_simul;
+ maxerror = real_maxerror;
else
- y0 = NaN(ny, 1);
- end
- if M_.maximum_lead > 0
- yT = real(oo_.endo_simul(:, M_.maximum_lag+periods+1));
- else
- yT = NaN(ny, 1);
- end
- 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_);
-
- maxerror = max(max(abs(residuals)))
- if maxerror < options_.dynatol.f
- success = true;
- oo_.endo_simul=real(oo_.endo_simul);
- else
- success = false;
+ current_share = 0;
disp('Simulation terminated with imaginary parts in the residuals or endogenous variables.')
end
end
-if success
+% Put solver status information in oo_, and do linearization if needed and requested
+if current_share == 1
if ~options_.noprint
fprintf('Perfect foresight solution found.\n\n')
end
+ oo_.deterministic_simulation.status = true;
+elseif options_.simul.homotopy_linearization_fallback && current_share > 0
+ oo_.endo_simul = endobase + (oo_.endo_simul - endobase)/current_share;
+ oo_.exo_simul = exoorig;
+ if options_.simul.endval_steady
+ % The following is needed for the STEADY_STATE() operator to work properly,
+ % and thus must come before computing the maximum error.
+ % This is not a true steady state, but it is the closest we can get to
+ oo_.steady_state = oo_.endo_simul(:, end);
+ end
+ maxerror = recompute_maxerror(oo_.endo_simul, oo_, M_, options_);
+
+ if ~options_.noprint
+ fprintf('Perfect foresight solution found for %.1f%% of the shock, then extrapolation was performed using linearization\n\n', current_share*100)
+ end
+ oo_.deterministic_simulation.status = true;
+ oo_.deterministic_simulation.homotopy_linearization = true;
+elseif options_.simul.homotopy_marginal_linearization_fallback > 0 && current_share > options_.simul.homotopy_marginal_linearization_fallback
+ saved_endo_simul = oo_.endo_simul;
+ new_share = current_share - options_.simul.homotopy_marginal_linearization_fallback;
+ new_success = create_scenario(new_share);
+ if new_success
+ [oo_.endo_simul, new_success] = perfect_foresight_solver_core(M_, options_, oo_);
+ end
+ if new_success
+ oo_.endo_simul = saved_endo_simul + (saved_endo_simul - oo_.endo_simul)*(1-current_share)/options_.simul.homotopy_marginal_linearization_fallback;
+ oo_.exo_simul = exoorig;
+ if options_.simul.endval_steady
+ % The following is needed for the STEADY_STATE() operator to work properly,
+ % and thus must come before computing the maximum error.
+ % This is not a true steady state, but it is the closest we can get to
+ oo_.steady_state = oo_.endo_simul(:, end);
+ end
+ maxerror = recompute_maxerror(oo_.endo_simul, oo_, M_, options_);
+
+ if ~options_.noprint
+ fprintf('Perfect foresight solution found for %.1f%% of the shock, then extrapolation was performed using marginal linearization\n\n', current_share*100)
+ end
+ oo_.deterministic_simulation.homotopy_marginal_linearization = true;
+ else
+ fprintf('perfect_foresight_solver: marginal linearization failed, unable to find solution for %.1f%% of the shock. Try to modify the value of homotopy_marginal_linearization_fallback option\n\n', new_share*100)
+ end
+ oo_.deterministic_simulation.status = new_success;
else
fprintf('Failed to solve perfect foresight model\n\n')
+ oo_.deterministic_simulation.status = false;
end
-% Put solver status information in oo_
-oo_.deterministic_simulation.status = success;
oo_.deterministic_simulation.error = maxerror;
+oo_.deterministic_simulation.homotopy_completion_share = current_share;
+oo_.deterministic_simulation.homotopy_iterations = iteration;
if ~isempty(solver_iter)
oo_.deterministic_simulation.iterations = solver_iter;
end
@@ -214,6 +370,12 @@ if ~isempty(per_block_status)
oo_.deterministic_simulation.block = per_block_status;
end
+% Must come after marginal linearization
+if did_homotopy
+ options_.verbosity = oldverbositylevel;
+ warning(warning_old_state);
+end
+
dyn2vec(M_, oo_, options_);
if isfield(oo_, 'initval_series') && ~isempty(oo_.initval_series)
@@ -236,3 +398,29 @@ assignin('base', 'Simulated_time_series', ts);
if success
oo_.gui.ran_perfect_foresight = true;
end
+
+end
+
+function maxerror = recompute_maxerror(endo_simul, oo_, M_, options_)
+ % Computes ∞-norm of residuals for a given path of endogenous,
+ % given the exogenous path and parameters in oo_ and M_
+ if options_.bytecode
+ residuals = bytecode('dynamic', 'evaluate', endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1);
+ else
+ ny = size(endo_simul, 1);
+ periods = size(endo_simul, 2) - M_.maximum_lag - M_.maximum_lead;
+ if M_.maximum_lag > 0
+ y0 = endo_simul(:, M_.maximum_lag);
+ else
+ y0 = NaN(ny, 1);
+ end
+ if M_.maximum_lead > 0
+ yT = endo_simul(:, M_.maximum_lag+periods+1);
+ else
+ yT = NaN(ny, 1);
+ end
+ yy = 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 = norm(vec(residuals), 'Inf');
+end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index e1a384b20b..a4e9ec5845 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -68,9 +68,6 @@ if options_.block
if options_.verbosity >= 1
disp(ME.message)
end
- if options_.no_homotopy
- error('Error in bytecode')
- end
y = oo_.endo_simul; % Set something for y, need for computing maxerror
success = false;
end
@@ -86,9 +83,6 @@ else
if options_.verbosity >= 1
disp(ME.message)
end
- if options_.no_homotopy
- error('Error in bytecode')
- end
y = oo_.endo_simul; % Set something for y, need for computing maxerror
success = false;
end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_setup.m b/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_setup.m
index 2a89e13efe..39ae5ede12 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_setup.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_setup.m
@@ -125,21 +125,6 @@ else
end
end
-%% Compute the terminal steady state for all informational periods
-oo_.pfwee.terminal_steady_state = NaN(M_.endo_nbr, periods);
-for p = 1:periods
- if p > 1 && all(oo_.pfwee.terminal_info(:, p) == oo_.pfwee.terminal_info(:, p-1))
- oo_.pfwee.terminal_steady_state(:, p) = oo_.pfwee.terminal_steady_state(:, p-1);
- else
- if p == 1
- init = oo_.steady_state;
- else
- init = oo_.pfwee.terminal_steady_state(:, p-1);
- end
- oo_.pfwee.terminal_steady_state(:, p) = evaluate_steady_state(init, oo_.pfwee.terminal_info(:, p), M_, options_, true);
- end
-end
-
% Build initial paths for endos and exos (only initial conditions are set, the rest is NaN)
if isempty(ys0_)
oo_.endo_simul = repmat(oo_.steady_state, 1, M_.maximum_lag+periods+M_.maximum_lead);
diff --git a/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m b/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m
index c2e7f0fed3..5c05102eba 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m
@@ -19,18 +19,23 @@ function perfect_foresight_with_expectation_errors_solver
global M_ oo_ options_ ys0_
-% Save original steady state, for restoring it at the end
-orig_steady_state = oo_.steady_state;
-orig_exo_steady_state = oo_.exo_steady_state;
-
% Same for periods (it will be modified before calling perfect_foresight_solver if constants_simulation_length option is false)
periods = options_.periods;
+% Save some options
+orig_homotopy_max_completion_share = options_.simul.homotopy_max_completion_share;
+orig_endval_steady = options_.simul.endval_steady;
+
% Retrieve initial paths built by pfwee_setup
% (the versions in oo_ will be truncated before calling perfect_foresight_solver)
endo_simul = oo_.endo_simul;
exo_simul = oo_.exo_simul;
+% Enforce the endval steady option in pf_solver
+if M_.maximum_lead > 0
+ options_.simul.endval_steady = true;
+end
+
% Start main loop around informational periods
info_period = 1;
increment = 0;
@@ -41,10 +46,11 @@ while info_period <= periods
% Compute terminal steady state as anticipated
oo_.exo_steady_state = oo_.pfwee.terminal_info(:, info_period);
- oo_.steady_state = oo_.pfwee.terminal_steady_state(:, info_period);
+ % oo_.steady_state will be updated by pf_solver (since endval_steady=true)
if options_.pfwee.constant_simulation_length && increment > 0
- endo_simul = [ endo_simul NaN(M_.endo_nbr, increment)];
+ % Use previous terminal steady state as guess value for: simulation periods that don’t yet have an initial guess (i.e. are NaNs at this point); and also for the terminal steady state
+ endo_simul = [ endo_simul repmat(oo_.steady_state, 1, increment)];
exo_simul = [ exo_simul; NaN(increment, M_.exo_nbr)];
end
@@ -54,11 +60,7 @@ while info_period <= periods
else
sim_length = periods - info_period + 1;
end
- if options_.pfwee.constant_simulation_length && increment > M_.maximum_lead
- % Use terminal steady state as guess value for simulation periods that don’t yet have an initial guess (i.e. are NaNs at this point)
- oo_.endo_simul(:, M_.maximum_lag+periods-(0:increment-M_.maximum_lead-1)) = repmat(oo_.steady_state, 1, increment-M_.maximum_lead);
- end
- oo_.endo_simul(:, end-M_.maximum_lead+1:end) = repmat(oo_.steady_state, 1, M_.maximum_lead);
+
oo_.exo_simul = exo_simul(info_period:end, :);
oo_.exo_simul(M_.maximum_lag+(1:periods-info_period+1), :) = oo_.pfwee.shocks_info(:, info_period:end, info_period)';
oo_.exo_simul(M_.maximum_lag+periods-info_period+2:end, :) = repmat(oo_.exo_steady_state', sim_length+M_.maximum_lead-(periods-info_period+1), 1);
@@ -71,6 +73,13 @@ while info_period <= periods
error('perfect_foresight_with_expectation_errors_solver: failed to compute solution for information available at period %d\n', info_period)
end
+ if info_period == 1
+ homotopy_completion_share = oo_.deterministic_simulation.homotopy_completion_share;
+ options_.simul.homotopy_max_completion_share = homotopy_completion_share;
+ elseif oo_.deterministic_simulation.homotopy_completion_share ~= homotopy_completion_share
+ error('perfect_foresight_solver_with_expectation_errors: could not find a solution for information available at period %d with the same homotopy completion share as period 1\n', info_period)
+ end
+
endo_simul(:, info_period:end) = oo_.endo_simul;
exo_simul(info_period:end, :) = oo_.exo_simul;
@@ -88,7 +97,7 @@ end
oo_.endo_simul = endo_simul;
oo_.exo_simul = exo_simul;
-% Restore some values
-oo_.steady_state = orig_steady_state;
-oo_.exo_steady_state = orig_exo_steady_state;
+% Restore some options
options_.periods = periods;
+options_.simul.homotopy_max_completion_share = orig_homotopy_max_completion_share;
+options_.simul.endval_steady = orig_endval_steady;
diff --git a/preprocessor b/preprocessor
index a1b8602760..fa5c4c5191 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit a1b8602760e77895ee065de4dadb0d532c4e926e
+Subproject commit fa5c4c5191600143f2203ae8cc1e87b709ff6921
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 9a055bfd0b..9d3704d09c 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -372,8 +372,9 @@ MODFILES = \
deterministic_simulations/rbc_det5.mod \
deterministic_simulations/rbc_det6.mod \
deterministic_simulations/homotopy.mod \
- deterministic_simulations/homotopy_alt_starting_point.mod \
deterministic_simulations/homotopy_histval.mod \
+ deterministic_simulations/homotopy_endval_steady_linearization.mod \
+ deterministic_simulations/homotopy_marginal_linearization.mod \
deterministic_simulations/rbc_det_exo_lag_2a.mod \
deterministic_simulations/rbc_det_exo_lag_2b.mod \
deterministic_simulations/rbc_det_exo_lag_2c.mod \
@@ -397,6 +398,7 @@ MODFILES = \
deterministic_simulations/pfwee_constant_sim_length.mod \
deterministic_simulations/pfwee_learnt_in.mod \
deterministic_simulations/pfwee_multiple_shocks.mod \
+ deterministic_simulations/pfwee_homotopy.mod \
lmmcp/rbc.mod \
lmmcp/sw_lmmcp.mod \
lmmcp/sw_newton.mod \
diff --git a/tests/deterministic_simulations/homotopy.mod b/tests/deterministic_simulations/homotopy.mod
index 64ccde0d84..a1043102e8 100644
--- a/tests/deterministic_simulations/homotopy.mod
+++ b/tests/deterministic_simulations/homotopy.mod
@@ -1,4 +1,5 @@
// Example that triggers homotopy in perfect foresight simulation.
+// Simulation with a permanent shock.
var Consumption, Capital, LoggedProductivity;
@@ -34,12 +35,12 @@ set_time(1Q1);
initval;
LoggedProductivityInnovation = 0;
- LoggedProductivity = 10;
- Capital = 1;
end;
+steady;
+
endval;
- LoggedProductivityInnovation = 0;
+ LoggedProductivityInnovation = 1;
end;
steady;
diff --git a/tests/deterministic_simulations/homotopy_alt_starting_point.mod b/tests/deterministic_simulations/homotopy_alt_starting_point.mod
deleted file mode 100644
index 81af0a728d..0000000000
--- a/tests/deterministic_simulations/homotopy_alt_starting_point.mod
+++ /dev/null
@@ -1,50 +0,0 @@
-// Test for the homotopy_alt_starting_point option of perfect_foresight_solver
-
-var Consumption, Capital, LoggedProductivity;
-
-varexo LoggedProductivityInnovation;
-
-parameters beta, alpha, delta, rho;
-
-beta = .985;
-alpha = 1/3;
-delta = alpha/10;
-rho = .9;
-
-model;
-
-[name='Euler equation']
-1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta);
-
-[name='Physical capital stock law of motion']
-Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption;
-
-[name='Logged productivity law of motion']
-LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
-
-end;
-
-steady_state_model;
- LoggedProductivity = LoggedProductivityInnovation/(1-rho);
- Capital = (exp(LoggedProductivity)*alpha/(1/beta-1+delta))^(1/(1-alpha));
- Consumption = exp(LoggedProductivity)*Capital^alpha-delta*Capital;
-end;
-
-initval;
- LoggedProductivityInnovation = 0;
-end;
-
-steady;
-
-endval;
- LoggedProductivityInnovation = 0.4;
-end;
-
-steady;
-
-perfect_foresight_setup(periods=200);
-perfect_foresight_solver(homotopy_alt_starting_point);
-
-if ~oo_.deterministic_simulation.status
- error('Perfect foresight simulation failed')
-end
diff --git a/tests/deterministic_simulations/homotopy_endval_steady_linearization.mod b/tests/deterministic_simulations/homotopy_endval_steady_linearization.mod
new file mode 100644
index 0000000000..97b8fbae01
--- /dev/null
+++ b/tests/deterministic_simulations/homotopy_endval_steady_linearization.mod
@@ -0,0 +1,43 @@
+// Example that triggers homotopy in perfect foresight simulation.
+// Tests the endval_steady, homotopy_linearization_fallback options, and a few more.
+
+var Consumption, Capital, LoggedProductivity;
+
+varexo LoggedProductivityInnovation;
+
+parameters beta, alpha, delta, rho;
+
+beta = .985;
+alpha = 1/3;
+delta = alpha/10;
+rho = .9;
+
+model;
+ 1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta);
+ Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption;
+ LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
+end;
+
+initval;
+ LoggedProductivityInnovation = 0;
+end;
+
+steady;
+
+endval;
+ LoggedProductivityInnovation = 1;
+ Consumption = 0.1;
+ Capital = 1;
+end;
+
+perfect_foresight_setup(periods=200, endval_steady);
+
+perfect_foresight_solver(homotopy_initial_step_size = 0.5,
+ homotopy_max_completion_share = 0.6,
+ homotopy_min_step_size = 0.0001,
+ homotopy_linearization_fallback,
+ steady_solve_algo = 13, steady_tolf=1e-7);
+
+if ~oo_.deterministic_simulation.status
+ error('Perfect foresight simulation failed')
+end
diff --git a/tests/deterministic_simulations/homotopy_histval.mod b/tests/deterministic_simulations/homotopy_histval.mod
index 8a1bba3013..551bab78ff 100644
--- a/tests/deterministic_simulations/homotopy_histval.mod
+++ b/tests/deterministic_simulations/homotopy_histval.mod
@@ -1,4 +1,5 @@
// Example that triggers homotopy in perfect foresight simulation.
+// Simulation starts out of steady state and returns to it.
var Consumption, Capital, LoggedProductivity;
@@ -24,11 +25,11 @@ LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
end;
-% steady_state_model;
-% LoggedProductivity = LoggedProductivityInnovation/(1-rho);
-% Capital = (exp(LoggedProductivity)*alpha/(1/beta-1+delta))^(1/(1-alpha));
-% Consumption = exp(LoggedProductivity)*Capital^alpha-delta*Capital;
-% end;
+steady_state_model;
+ LoggedProductivity = LoggedProductivityInnovation/(1-rho);
+ Capital = (exp(LoggedProductivity)*alpha/(1/beta-1+delta))^(1/(1-alpha));
+ Consumption = exp(LoggedProductivity)*Capital^alpha-delta*Capital;
+end;
set_time(1Q1);
@@ -45,12 +46,6 @@ histval;
LoggedProductivity(0)=10;
end;
-// endval;
-// LoggedProductivityInnovation = 0;
-// end;
-
-// steady;
-
perfect_foresight_setup(periods=200);
perfect_foresight_solver;
diff --git a/tests/deterministic_simulations/homotopy_marginal_linearization.mod b/tests/deterministic_simulations/homotopy_marginal_linearization.mod
new file mode 100644
index 0000000000..080e11b67c
--- /dev/null
+++ b/tests/deterministic_simulations/homotopy_marginal_linearization.mod
@@ -0,0 +1,42 @@
+// Example that triggers homotopy in perfect foresight simulation.
+// Tests the homotopy_marginal_linearization_fallback and homotopy_step_size_increase_success_count options
+
+var Consumption, Capital, LoggedProductivity;
+
+varexo LoggedProductivityInnovation;
+
+parameters beta, alpha, delta, rho;
+
+beta = .985;
+alpha = 1/3;
+delta = alpha/10;
+rho = .9;
+
+model;
+ 1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta);
+ Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption;
+ LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
+end;
+
+initval;
+ LoggedProductivityInnovation = 0;
+end;
+
+steady;
+
+endval;
+ LoggedProductivityInnovation = 1;
+ Consumption = 0.1;
+ Capital = 1;
+end;
+
+perfect_foresight_setup(periods=200, endval_steady);
+
+perfect_foresight_solver(homotopy_max_completion_share = 0.7,
+ homotopy_step_size_increase_success_count = 3,
+ homotopy_marginal_linearization_fallback,
+ steady_solve_algo = 13);
+
+if ~oo_.deterministic_simulation.status
+ error('Perfect foresight simulation failed')
+end
diff --git a/tests/deterministic_simulations/pfwee.mod b/tests/deterministic_simulations/pfwee.mod
index a42300cbbc..fb7f4cf015 100644
--- a/tests/deterministic_simulations/pfwee.mod
+++ b/tests/deterministic_simulations/pfwee.mod
@@ -27,13 +27,18 @@ steady;
check;
+// Save initial steady state (it will be modified by pfwee)
+orig_steady_state = oo_.steady_state;
+orig_exo_steady_state = oo_.exo_steady_state;
+
perfect_foresight_with_expectation_errors_setup(periods = 7, datafile = 'pfwee.csv');
-// First simulation with default options
perfect_foresight_with_expectation_errors_solver;
pfwee_simul = oo_.endo_simul;
// Now compute the solution by hand to verify the results
+oo_.steady_state = orig_steady_state;
+oo_.exo_steady_state = orig_exo_steady_state;
perfect_foresight_setup;
diff --git a/tests/deterministic_simulations/pfwee_constant_sim_length.mod b/tests/deterministic_simulations/pfwee_constant_sim_length.mod
index 82a5199452..efe0d2efa2 100644
--- a/tests/deterministic_simulations/pfwee_constant_sim_length.mod
+++ b/tests/deterministic_simulations/pfwee_constant_sim_length.mod
@@ -26,17 +26,21 @@ steady;
check;
-// First simulation with default options
+// Save initial steady state (it will be modified by pfwee)
+orig_steady_state = oo_.steady_state;
+orig_exo_steady_state = oo_.exo_steady_state;
+
perfect_foresight_with_expectation_errors_setup(periods = 7, datafile = 'pfwee.csv');
+
perfect_foresight_with_expectation_errors_solver(constant_simulation_length);
pfwee_simul = oo_.endo_simul;
// Now compute the solution by hand to verify the results
+oo_.steady_state = orig_steady_state;
+oo_.exo_steady_state = orig_exo_steady_state;
perfect_foresight_setup;
-initial_steady_state = oo_.steady_state;
-
// Information arriving in period 1 (temp shock now)
oo_.exo_simul(2,1) = 1.2;
perfect_foresight_solver;
@@ -74,8 +78,9 @@ oo_.exo_simul(7,1) = 1.1;
oo_.exo_simul(8,1) = 1.1;
oo_.exo_steady_state = 1.1;
oo_.exo_simul(9:14, 1) = repmat(oo_.exo_steady_state', 6, 1);
+oo_.endo_simul(:, 12:13) = repmat(oo_.steady_state, 1, 2);
oo_.steady_state = evaluate_steady_state(oo_.steady_state, oo_.exo_steady_state, M_, options_, true);
-oo_.endo_simul(:, 12:14) = repmat(oo_.steady_state, 1, 3);
+oo_.endo_simul(:, 14) = oo_.steady_state;
saved_endo = oo_.endo_simul(:, 1:5);
saved_exo = oo_.exo_simul(1:5, :);
oo_.endo_simul = oo_.endo_simul(:, 6:end);
diff --git a/tests/deterministic_simulations/pfwee_homotopy.mod b/tests/deterministic_simulations/pfwee_homotopy.mod
new file mode 100644
index 0000000000..c4175c0906
--- /dev/null
+++ b/tests/deterministic_simulations/pfwee_homotopy.mod
@@ -0,0 +1,40 @@
+/* Example that triggers homotopy in perfect foresight simulation with
+ expectation errors. */
+
+var Consumption, Capital, LoggedProductivity;
+
+varexo LoggedProductivityInnovation;
+
+parameters beta, alpha, delta, rho;
+
+beta = .985;
+alpha = 1/3;
+delta = alpha/10;
+rho = .9;
+
+model;
+ 1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta);
+ Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption;
+ LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
+end;
+
+initval;
+ LoggedProductivityInnovation = 0;
+end;
+
+steady;
+
+endval;
+ LoggedProductivityInnovation = 0.6;
+end;
+
+endval(learnt_in = 5);
+ LoggedProductivityInnovation = 1;
+end;
+
+perfect_foresight_with_expectation_errors_setup(periods=200);
+perfect_foresight_with_expectation_errors_solver(homotopy_max_completion_share = 0.8, homotopy_linearization_fallback, steady_solve_algo = 13);
+
+if ~oo_.deterministic_simulation.status
+ error('Perfect foresight simulation failed')
+end
diff --git a/tests/deterministic_simulations/pfwee_learnt_in.mod b/tests/deterministic_simulations/pfwee_learnt_in.mod
index 8c2a0dc148..bef2a04674 100644
--- a/tests/deterministic_simulations/pfwee_learnt_in.mod
+++ b/tests/deterministic_simulations/pfwee_learnt_in.mod
@@ -70,13 +70,18 @@ endval(learnt_in = 6);
x *= 0.75;
end;
+// Save initial steady state (it will be modified by pfwee)
+orig_steady_state = oo_.steady_state;
+orig_exo_steady_state = oo_.exo_steady_state;
+
perfect_foresight_with_expectation_errors_setup(periods = 7);
-// First simulation with default options
perfect_foresight_with_expectation_errors_solver;
pfwee_simul = oo_.endo_simul;
// Now compute the solution by hand to verify the results
+oo_.steady_state = orig_steady_state;
+oo_.exo_steady_state = orig_exo_steady_state;
perfect_foresight_setup;
--
GitLab