Perfect foresight solver now uses a homotopy technique by default.

This commit introduces a "no_homotopy" option to restore the old behavior.

Ref #220

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}.

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;

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

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

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

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 \

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

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;