Various improvements to perfect foresight homotopy

– 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

@@ -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
+    .. option:: homotopy_initial_step_size = DOUBLE
 
-       When the homotopy technique is tried (as described above), Dynare first
+       Specifies which share of the shock should be applied in the first
-       tries to reduce the size of the shock in order to get a successful
+       iteration of the homotopy procedure. This option is useful when it is
-       simulation on which to build upon for simulating a shock of the true
+       known that immediately trying 100% of the shock will fail, so as to save
-       size. However, if an ``endval`` block is present (*i.e.* if the terminal
+       computing time. Must be between ``0`` and ``1``. Default: ``1``.
-       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
+    .. option:: homotopy_min_step_size = DOUBLE
-       computing a simulation whose initial state is closer to the target
+
-       terminal state; in other words, it will implicitly modify the contents
+       The homotopy procedure halves the size of the step whenever there is
-       of the ``initval`` and ``shocks`` blocks to make them closer to the the
+       a failure. This option specifies the minimum step size under which the
-       contents of the ``endval`` block. If this option is set, Dynare will do
+       homotopy procedure is considered to have failed. Default: ``0.001``.
-       the opposite: it will implicitly modify the contents of the ``endval``
+
-       and ``shocks`` blocks to make them closer to the contents of the
+    .. option:: homotopy_step_size_increase_success_count = INTEGER
-       ``initval`` block.
+
+       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_info``.
-    ``oo_.pfwee.terminal_steady_state``.
 
     *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

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;

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

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_);
+oo_ = make_y_(M_,options_,oo_);

matlab/perfect-foresight-models/perfect_foresight_solver.m

@@ -51,162 +51,318 @@ if options_.debug
     end
 end
 
+% 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
+
+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 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_.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;
-lastperiods = (M_.maximum_lag+periods+1):(M_.maximum_lag+periods+M_.maximum_lead);
+simperiods = M_.maximum_lag+(1:periods);
+lastperiods = M_.maximum_lag+periods+(1:M_.maximum_lead);
 
-[oo_.endo_simul, success, maxerror, solver_iter, per_block_status] = perfect_foresight_solver_core(M_,options_,oo_);
+% 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
 
-% If simulation failed try homotopy.
+% Determine whether the terminal condition is not a steady state (typically
-if ~success && ~options_.no_homotopy
+% because steady was not called after endval)
-
+if ~options_.simul.endval_steady && ~isempty(ys0_)
-    if ~options_.noprint
+    terminal_condition_is_a_steady_state = true;
-        fprintf('\nSimulation of the perfect foresight model failed!\n')
+    for j = lastperiods
-        fprintf('Switching to a homotopy method...\n')
+        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
+% Copy the paths for the exogenous and endogenous variables, as given by perfect_foresight_setup
-    warning_old_state = warning;
+exoorig = oo_.exo_simul;
-    warning off all
+endoorig = oo_.endo_simul;
-    % Do not print anything
-    oldverbositylevel = options_.verbosity;
-    options_.verbosity = 0;
 
-    % Set initial paths for the endogenous and exogenous variables.
+current_share = 0;  % Share of shock successfully completed so far
-    if ~options_.homotopy_alt_starting_point
+step = min(options_.simul.homotopy_initial_step_size, options_.simul.homotopy_max_completion_share);
-        endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+periods+M_.maximum_lead);
+success_counter = 0;
-        exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+periods+M_.maximum_lead,1);
+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
+
+        % 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_)
+        % The terminal condition is not a steady state, compute a convex combination
-            error('The homotopy_alt_starting_point option cannot be used without an endval block');
+        oo_.endo_simul(:, lastperiods) = share*endoorig(:, lastperiods)+(1-share)*endobase(:, lastperiods);
-        end
+    end
-        endoinit = repmat(ys0_, 1, M_.maximum_lag+periods+M_.maximum_lead);
+end
-        exoinit = repmat(ex0_', M_.maximum_lag+periods+M_.maximum_lead, 1);
+
+while step > options_.simul.homotopy_min_step_size
+
+    iteration = iteration+1;
+
+    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
 
-    % Copy the current paths for the exogenous and endogenous variables.
+    steady_success = create_scenario(new_share);
-    exosim = oo_.exo_simul;
-    endosim = oo_.endo_simul;
 
-    current_weight = 0;    % Current weight of target point in convex combination.
+    if steady_success
-    step = .5;             % Set default step size.
+        % At the first iteration, use the initial guess given by
-    success_counter = 0;
+        % perfect_foresight_setup or the user (but only if new_share=1, otherwise it
-    iteration = 0;
+        % 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
-    if ~options_.noprint
+        % scenario / initial steady state).
-        fprintf('Iter. \t | Lambda \t | status \t | Max. residual\n')
+        if current_share == 0
-        fprintf('++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n')
+            if iteration == 1 && new_share == 1
-    end
+                oo_.endo_simul(:, simperiods) = endoorig(:, simperiods);
-    while (step > options_.dynatol.x)
+            elseif M_.maximum_lead > 0
-
+                oo_.endo_simul(:, simperiods) = repmat(oo_.endo_simul(:, lastperiods(1)), 1, options_.periods);
-        if ~isequal(step,1)
+            else
-            options_.verbosity = 0;
+                oo_.endo_simul(:, simperiods) = endobase(:, simperiods);
+            end
         end
 
-        iteration = iteration+1;
-        new_weight = current_weight + step; % Try this weight, and see if it succeeds
-
-        if new_weight >= 1
-            new_weight = 1; % Don't go beyond target point
-            step = new_weight - current_weight;
-        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
-
-        % Make a copy of the paths.
-        saved_endo_simul = oo_.endo_simul;
-
         % 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_);
+        [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 options_.no_homotopy || (iteration == 1 && success && new_share == 1)
+        % Skip homotopy
         if success
-            current_weight = new_weight;
+            current_share = new_share;
-            if current_weight >= 1
+        end
-                if ~options_.noprint
+        did_homotopy = false;
-                    fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', maxerror)
+        break
-                end
+    end
-                break
+
-            end
+    if iteration == 1
-            success_counter = success_counter + 1;
+        % First iteration failed, so we enter homotopy
-            if success_counter >= 3
+        did_homotopy = true;
-                success_counter = 0;
+
-                step = step * 2;
+        if ~options_.noprint
-            end
+            fprintf('\nEntering the homotopy method iterations...\n')
-            if ~options_.noprint
+            fprintf('\nIter. \t | Share \t | Status \t | Max. residual\n')
-                fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', maxerror)
+            fprintf('++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n')
-            end
+        end
-        else
+
-            % If solver failed, then go back.
+        % Disable warnings if homotopy
-            oo_.endo_simul = saved_endo_simul;
+        warning_old_state = warning;
+        warning off all
+        % Do not print anything
+        oldverbositylevel = options_.verbosity;
+        options_.verbosity = 0;
+    end
+
+    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;
+            step = step * 2;
-            if ~options_.noprint
+        end
-                if isreal(maxerror)
+    else
-                    fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'failed', maxerror)
+        oo_.endo_simul = saved_endo_simul;
-                else
+        success_counter = 0;
-                    fprintf('%i \t | %1.5f \t | %s \t | %s\n', iteration, new_weight, 'failed', 'Complex')
+        step = step / 2;
-                end
+        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
+end
-        fprintf('++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n')
+
-    end
+if did_homotopy && ~options_.noprint
-    options_.verbosity = oldverbositylevel;
+    fprintf('++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n')
-    warning(warning_old_state);
 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);
+    real_simul = real(oo_.endo_simul);
-    if M_.maximum_lag > 0
+    real_maxerror = recompute_maxerror(real_simul, oo_, M_, options_);
-        y0 = real(oo_.endo_simul(:, M_.maximum_lag));
+    if real_maxerror <= options_.dynatol.f
+        oo_.endo_simul = real_simul;
+        maxerror = real_maxerror;
     else
-        y0 = NaN(ny, 1);
+        current_share = 0;
-    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;
         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

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

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);

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
+% Restore some options
-oo_.steady_state = orig_steady_state;
-oo_.exo_steady_state = orig_exo_steady_state;
 options_.periods = periods;
+options_.simul.homotopy_max_completion_share = orig_homotopy_max_completion_share;
+options_.simul.endval_steady = orig_endval_steady;

preprocessor

@@ -1 +1 @@
-Subproject commit a1b8602760e77895ee065de4dadb0d532c4e926e
+Subproject commit fa5c4c5191600143f2203ae8cc1e87b709ff6921

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 \

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;

tests/deterministic_simulations/homotopy_alt_starting_point.mod

@@ -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

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

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;
+steady_state_model;
-%   LoggedProductivity = LoggedProductivityInnovation/(1-rho);
+  LoggedProductivity = LoggedProductivityInnovation/(1-rho);
-%   Capital = (exp(LoggedProductivity)*alpha/(1/beta-1+delta))^(1/(1-alpha));
+  Capital = (exp(LoggedProductivity)*alpha/(1/beta-1+delta))^(1/(1-alpha));
-%   Consumption = exp(LoggedProductivity)*Capital^alpha-delta*Capital;
+  Consumption = exp(LoggedProductivity)*Capital^alpha-delta*Capital;
-% end;
+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;
 

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

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;
 

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);

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

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;