diff --git a/doc/dynare.texi b/doc/dynare.texi
index 453794632dd692beb49886f575e13f1104890ea7..9152838eacd3a258a914f8661b6d0e0e83d3368b 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -3405,6 +3405,12 @@ options).
@end table
+@item 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.
+
@item markowitz = @var{DOUBLE}
Value of the Markowitz criterion, used to select the pivot. Only used
when @code{stack_solve_algo = 5}. Default: @code{0.5}.
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 1414c5de7e07f92d50ad03d5fbda98a304f476ca..cedf42ceb08152c0aad84adc5cb6023695d61062 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -272,6 +272,7 @@ options_.stack_solve_algo = 0;
options_.markowitz = 0.5;
options_.minimal_solving_periods = 1;
options_.endogenous_terminal_period = 0;
+options_.no_homotopy = 0;
% Solution
options_.order = 2;
@@ -444,7 +445,7 @@ M_.endo_histval = [];
M_.Correlation_matrix = [];
M_.Correlation_matrix_ME = [];
-% homotopy
+% homotopy for steady state
options_.homotopy_mode = 0;
options_.homotopy_steps = 1;
options_.homotopy_force_continue = 0;
diff --git a/matlab/perfect_foresight_solver.m b/matlab/perfect_foresight_solver.m
index ce9d94306d4d8bf31d282766a3673cf898ccc5dd..5c6321462913e45d4d25982d0e17ed562c053984 100644
--- a/matlab/perfect_foresight_solver.m
+++ b/matlab/perfect_foresight_solver.m
@@ -84,6 +84,75 @@ if options_.debug
end
end
+% Effectively compute simulation, possibly with homotopy
+if options_.no_homotopy
+ simulation_core;
+else
+ exosim = oo_.exo_simul;
+ exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1);
+ endosim = oo_.endo_simul;
+ endoinit = repmat(oo_.steady_state, 1,M_.maximum_endo_lag+options_.periods+M_.maximum_endo_lead);
+
+ current_weight = 0; % Current weight of target point in convex combination
+ step = 1;
+ success_counter = 0;
+
+ while (step > options_.dynatol.x)
+
+ 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);
+ endocombi = endosim*new_weight + endoinit*(1-new_weight);
+ oo_.endo_simul(:,1:M_.maximum_endo_lag) = endocombi(:,1:M_.maximum_endo_lag);
+ oo_.endo_simul(:,(end-M_.maximum_endo_lead):end) = endocombi(:,(end-M_.maximum_endo_lead):end);
+
+ simulation_core;
+
+ if oo_.deterministic_simulation.status == 1
+ current_weight = new_weight;
+ if current_weight >= 1
+ break
+ end
+ success_counter = success_counter + 1;
+ if success_counter >= 3
+ success_counter = 0;
+ step = step * 2;
+ disp('Homotopy step succeeded, doubling step size')
+ else
+ disp('Homotopy step succeeded')
+ end
+ else
+ success_counter = 0;
+ step = step / 2;
+ disp('Homotopy step failed, halving step size')
+ end
+ end
+end
+
+if oo_.deterministic_simulation.status == 1
+ disp('Perfect foresight solution found.')
+else
+ disp('Failed to solve perfect foresight model')
+end
+
+dyn2vec;
+
+ts = dseries(transpose(oo_.endo_simul),options_.initial_period,cellstr(M_.endo_names));
+assignin('base', 'Simulated_time_series', ts);
+
+end
+
+
+function simulation_core()
+
+global M_ oo_ options_
+
if(options_.block)
if(options_.bytecode)
[info, oo_.endo_simul] = bytecode('dynamic');
@@ -120,7 +189,4 @@ else
end
end
-dyn2vec;
-
-ts = dseries(transpose(oo_.endo_simul),options_.initial_period,cellstr(M_.endo_names));
-assignin('base', 'Simulated_time_series', ts);
+end
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index c986eac70a77cfd8b1e2dc47cb44f8ff4001b08a..c2cb723e718060e13898aea04f5a0472d5e7fa3e 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -108,7 +108,7 @@ class ParsingDriver;
%token MODE_CHECK MODE_CHECK_NEIGHBOURHOOD_SIZE MODE_CHECK_SYMMETRIC_PLOTS MODE_CHECK_NUMBER_OF_POINTS MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN
%token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL MCMC_JUMPING_COVARIANCE MOMENT_CALIBRATION
%token <string_val> NAME
-%token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS
+%token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY
%token NOGRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS
%token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS MAX_DIM_COVA_GROUP ADVANCED OUTFILE OUTVARS OVERWRITE
%token PARALLEL_LOCAL_FILES PARAMETERS PARAMETER_SET PARTIAL_INFORMATION PERFECT_FORESIGHT PERIODS PERIOD PLANNER_OBJECTIVE PLOT_CONDITIONAL_FORECAST PLOT_PRIORS PREFILTER PRESAMPLE
@@ -1006,6 +1006,7 @@ perfect_foresight_solver_options : o_stack_solve_algo
| o_minimal_solving_periods
| o_simul_maxit
| o_endogenous_terminal_period
+ | o_no_homotopy
;
simul : SIMUL ';'
@@ -2896,6 +2897,7 @@ o_mcmc_jumping_covariance : MCMC_JUMPING_COVARIANCE EQUAL HESSIAN
o_irf_plot_threshold : IRF_PLOT_THRESHOLD EQUAL non_negative_number { driver.option_num("impulse_responses.plot_threshold", $3); };
o_consider_all_endogenous : CONSIDER_ALL_ENDOGENOUS { driver.option_str("endo_vars_for_moment_computations_in_estimation", "all_endogenous_variables"); };
o_consider_only_observed : CONSIDER_ONLY_OBSERVED { driver.option_str("endo_vars_for_moment_computations_in_estimation", "only_observed_variables"); };
+o_no_homotopy : NO_HOMOTOPY { driver.option_num("no_homotopy", "1"); };
o_infile : INFILE EQUAL filename { driver.option_str("infile", $3); };
o_invars : INVARS EQUAL '(' symbol_list ')' { driver.option_symbol_list("invars"); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index cc1fffd7da14b2c12da07d0eda9442559b870cac..6ab6faab0ce5ed866986d466fc4a9141ceb415ae 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -608,6 +608,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
<DYNARE_STATEMENT>planner_discount {return token::PLANNER_DISCOUNT;}
<DYNARE_STATEMENT>calibration {return token::CALIBRATION;}
<DYNARE_STATEMENT>irf_plot_threshold {return token::IRF_PLOT_THRESHOLD;}
+<DYNARE_STATEMENT>no_homotopy {return token::NO_HOMOTOPY;}
<DYNARE_BLOCK>equation {return token::EQUATION;}
<DYNARE_BLOCK>exclusion {return token::EXCLUSION;}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 9fac03949d74cd9c2c23849e647173d8fddcaae3..b53d1f5a616f6e9871979f845927787fa3bc4f3f 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -156,6 +156,7 @@ MODFILES = \
deterministic_simulations/rbc_det4.mod \
deterministic_simulations/rbc_det5.mod \
deterministic_simulations/rbc_det6.mod \
+ deterministic_simulations/homotopy.mod \
walsh.mod \
measurement_errors/fs2000_corr_me_ml_mcmc/fs2000_corr_ME.mod \
trend_var/fs2000_nonstationary.mod \
diff --git a/tests/deterministic_simulations/homotopy.mod b/tests/deterministic_simulations/homotopy.mod
new file mode 100644
index 0000000000000000000000000000000000000000..cd17a513e1c07ce5a288bb5991a7addc05fb8116
--- /dev/null
+++ b/tests/deterministic_simulations/homotopy.mod
@@ -0,0 +1,49 @@
+// Example that triggers homotopy in perfect foresight simulation.
+
+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'] // This is an equation tag!
+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;
+
+set_time(1Q1);
+
+initval;
+ LoggedProductivityInnovation = 0;
+ LoggedProductivity = 10;
+ Capital = 1;
+end;
+
+endval;
+ LoggedProductivityInnovation = 0;
+end;
+
+steady;
+
+simul(periods=200);
+
+plot(Simulated_time_series.Capital(1Q1:25Q4));
diff --git a/tests/deterministic_simulations/rbc_det1.mod b/tests/deterministic_simulations/rbc_det1.mod
index e0a6877b6db659b789d73d5d58d3c94ce7c2b371..a5fe6f5b3e857705bdb54261f38e0063c032979b 100644
--- a/tests/deterministic_simulations/rbc_det1.mod
+++ b/tests/deterministic_simulations/rbc_det1.mod
@@ -69,7 +69,7 @@ histval;
Capital(0) = CapitalSS/2;
end;
-simul(periods=300);
+simul(periods=20);
rplot Consumption;
rplot Capital;
\ No newline at end of file