New “homotopy_alt_starting_point” option to “perfect_foresight_solver” command

doc/manual/source/the-model-file.rst

@@ -3581,6 +3581,14 @@ speed-up on large models.
     Note that ``perfect_foresight_setup`` must be called before this
     command, in order to setup the environment for the simulation.
 
+    If the perfect foresight solver cannot directly find the solution of the
+    problem, it subsequently tries a homotopy technique (unless the
+    ``no_homotopy`` option is given). Concretely, this technique consists in
+    dividing the problem into smaller steps by diminishing the size of shocks
+    and increasing them progressively until the problem converges. Note that
+    the homotopy technique is not implemented for purely forward or backward
+    models.
+
     *Options*
 
     .. option:: maxit = INTEGER
@@ -3690,13 +3698,24 @@ speed-up on large models.
 
     .. option:: no_homotopy
 
-       By default, the perfect foresight solver uses a homotopy
-       technique if it cannot solve the problem. Concretely, it
-       divides the problem into smaller steps by diminishing the size
-       of shocks and increasing them progressively until the problem
-       converges. This option tells Dynare to disable that
-       behavior. Note that the homotopy is not implemented for purely
-       forward or backward 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:: markowitz = DOUBLE
 

matlab/default_option_values.m

@@ -12,7 +12,7 @@ function options_ = default_option_values(M_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2018-2022 Dynare Team
+% Copyright © 2018-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -328,6 +328,7 @@ options_.markowitz = 0.5;
 options_.minimal_solving_periods = 1;
 options_.endogenous_terminal_period = false;
 options_.no_homotopy = false;
+options_.homotopy_alt_starting_point = false;
 
 % Perfect foresight with expectation errors
 options_.pfwee.terminal_steady_state_as_guess_value = false;
@@ -781,4 +782,4 @@ options_.pac.estimation.ols.share_of_optimizing_agents.lb = 0.0;
 options_.pac.estimation.ols.share_of_optimizing_agents.ub = 1.0;
 
 options_.conditional_likelihood.status = false;
-options_.conditional_likelihood.order = 1;
+options_.conditional_likelihood.order = 1;

matlab/perfect-foresight-models/perfect_foresight_solver.m

@@ -29,7 +29,7 @@ function perfect_foresight_solver()
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-global M_ options_ oo_
+global M_ options_ oo_ ys0_ ex0_
 
 check_input_arguments(options_, M_, oo_);
 
@@ -102,8 +102,16 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
     options_.verbosity = 0;
 
     % Set initial paths for the endogenous and exogenous variables.
-    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);
+    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);
+    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);
+    end
 
     % Copy the current paths for the exogenous and endogenous variables.
     exosim = oo_.exo_simul;

preprocessor

@@ -1 +1 @@
-Subproject commit 389a2647d3b1067d3af54596ad3bd076fab8a4b4
+Subproject commit 114d8eadfb3dee24ab65f291253076c4d0ea83c3

tests/Makefile.am

@@ -375,6 +375,7 @@ 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/rbc_det_exo_lag_2a.mod \
 	deterministic_simulations/rbc_det_exo_lag_2b.mod \

tests/deterministic_simulations/homotopy_alt_starting_point.mod

@@ -0,0 +1,50 @@
+// 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