From 85f7af913341d45f59d893c9dd09cb461478152c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Thu, 10 Apr 2014 16:38:39 +0200 Subject: [PATCH] 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 | 6 ++ matlab/global_initialization.m | 3 +- matlab/perfect_foresight_solver.m | 74 ++++++++++++++++++-- preprocessor/DynareBison.yy | 4 +- preprocessor/DynareFlex.ll | 1 + tests/Makefile.am | 1 + tests/deterministic_simulations/homotopy.mod | 49 +++++++++++++ tests/deterministic_simulations/rbc_det1.mod | 2 +- 8 files changed, 133 insertions(+), 7 deletions(-) create mode 100644 tests/deterministic_simulations/homotopy.mod diff --git a/doc/dynare.texi b/doc/dynare.texi index 453794632..9152838ea 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 1414c5de7..cedf42ceb 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 ce9d94306..5c6321462 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 c986eac70..c2cb723e7 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 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 cc1fffd7d..6ab6faab0 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 planner_discount {return token::PLANNER_DISCOUNT;} calibration {return token::CALIBRATION;} irf_plot_threshold {return token::IRF_PLOT_THRESHOLD;} +no_homotopy {return token::NO_HOMOTOPY;} equation {return token::EQUATION;} exclusion {return token::EXCLUSION;} diff --git a/tests/Makefile.am b/tests/Makefile.am index 9fac03949..b53d1f5a6 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 000000000..cd17a513e --- /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 e0a6877b6..a5fe6f5b3 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