diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index fc1b5305536f417bca417aaa5d79453c77188e62..999e689b21c26c36bfceb630719d3c8bd7c722fe 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 52bd5b6feabaa16b42e4b1c8c85b26c6abb7810d..c0a0bf67ec95cbed267007f14d1b861a266508bc 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 4de31803f9b9515ac81ce48fdd6211eb4041ed57..c23f562adc7aab1c8f4fbc369dca0281e86df8b2 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 db0e04fbf5483e383a97d51896621c15a75bd25a..f187b9dc2e38ddd6e7a7449b581c592d7d20c527 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 bfbf63b0b9d167f9ebc6b663c13f8c51e2071dee..d412afb30122659873708df3efbefbb9166657db 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 e1a384b20bfcbda14829154415d9d3d0045ae835..a4e9ec584534ae323872795c271f38dd65896e5c 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 2a89e13efe5846e7c63bfcdc40e0aa8929cd2d6a..39ae5ede12e351ca4b96f925948a167e55311117 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 c2e7f0fed3226dd54c8cf391d565a334ca11e4d1..5c05102ebae46623e2f32b679b22656cb714357f 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 a1b8602760e77895ee065de4dadb0d532c4e926e..fa5c4c5191600143f2203ae8cc1e87b709ff6921 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit a1b8602760e77895ee065de4dadb0d532c4e926e
+Subproject commit fa5c4c5191600143f2203ae8cc1e87b709ff6921
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 9a055bfd0ba9d677df30b21ec94f8a61062ec910..9d3704d09c55465f570b94d0226a1a1e37c84294 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 64ccde0d8480a8bb9f5bca6d6c261588bcdd8981..a1043102e83e815452538496cf165af529771fe4 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 81af0a728d0e979237e70a80b60caa8dd4b71139..0000000000000000000000000000000000000000
--- 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 0000000000000000000000000000000000000000..97b8fbae011c6baa30a41663c12c2c8e0383287c
--- /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 8a1bba3013930b3f9d2ad88c040f37f1f74226d7..551bab78ffa04d70214cea615317a2b50bf5b5a1 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 0000000000000000000000000000000000000000..080e11b67ceb7da602cff95c1df6ca5036c05a2f
--- /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 a42300cbbc53ab6dc985829df297730953e338be..fb7f4cf0151b444b95a544d4302767ebfc38fa27 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 82a5199452bd9ad9b247aca182089ca47f9a4954..efe0d2efa2d217ffcba5b94c1da99b708992fc8e 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 0000000000000000000000000000000000000000..c4175c0906284d1c91ba0966b08ee1d234919527
--- /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 8c2a0dc148b94073c73e78227809d7699b9a0706..bef2a04674b3ca8acb1b082614d63f48fd162e60 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;