From e714dc1a9f36fcd34f8a3da9e1de70bda814063c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=28Charybdis=29?= Date: Thu, 24 Mar 2016 22:42:44 +0100 Subject: [PATCH] Fixed bug in homotopy. Convex combination for initial/terminal endogenous variables paths was wrong. --- .../perfect_foresight_solver.m | 62 +++++++++++-------- 1 file changed, 35 insertions(+), 27 deletions(-) diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m index e171c749d..464bfff24 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m @@ -1,18 +1,18 @@ function perfect_foresight_solver() % Computes deterministic simulations -% +% % INPUTS % None -% +% % OUTPUTS % none -% +% % ALGORITHM -% +% % SPECIAL REQUIREMENTS % none -% Copyright (C) 1996-2015 Dynare Team +% Copyright (C) 1996-2016 Dynare Team % % This file is part of Dynare. % @@ -45,12 +45,12 @@ if options_.debug [residual(:,ii)] = model_static(oo_.steady_state, oo_.exo_simul(ii,:),M_.params); end problematic_periods=find(any(isinf(residual)) | any(isnan(residual)))-M_.maximum_endo_lag; - if ~isempty(problematic_periods) + if ~isempty(problematic_periods) period_string=num2str(problematic_periods(1)); for ii=2:length(problematic_periods) period_string=[period_string, ', ', num2str(problematic_periods(ii))]; end - fprintf('\n\nWARNING: Value for the exogenous variable(s) in period(s) %s inconsistent with the static model.\n',period_string); + fprintf('\n\nWARNING: Value for the exogenous variable(s) in period(s) %s inconsistent with the static model.\n',period_string); fprintf('WARNING: Check for division by 0.\n') end end @@ -62,15 +62,17 @@ oo_ = perfect_foresight_solver_core(M_,options_,oo_); % If simulation failed try homotopy. if ~oo_.deterministic_simulation.status && ~options_.no_homotopy + skipline() disp('Simulation of the perfect foresight model failed!') disp('Switching to a homotopy method...') + skipline() + if ~M_.maximum_lag disp('Homotopy not implemented for purely forward models!') disp('Failed to solve the model!') disp('Return with empty oo_.endo_simul.') oo_.endo_simul = []; - skipline() return end if ~M_.maximum_lead @@ -78,25 +80,26 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy disp('Failed to solve the model!') disp('Return with empty oo_.endo_simul.') oo_.endo_simul = []; - skipline() return end - skipline() + % Disable warnings if homotopy warning_old_state = warning; warning off all % Do not print anything oldverbositylevel = options_.verbosity; options_.verbosity = 0; - - 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_lag+options_.periods+M_.maximum_lead); - current_weight = 0; % Current weight of target point in convex combination - step = .5; + % Set initial paths for the endogenous and exogenous variables. + endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+options_.periods+M_.maximum_lead); + exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1); + + % Copy the current paths for the exogenous and endogenous variables. + exosim = oo_.exo_simul; + endosim = oo_.endo_simul; + + current_weight = 0; % Current weight of target point in convex combination. + step = .5; % Set default step size. success_counter = 0; iteration = 0; @@ -120,21 +123,24 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy % 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*oo_.endo_simul(:,[initperiods, lastperiods])+(1-new_weight)*endoinit(:,[initperiods, lastperiods]); - + 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) + + 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:options_.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_,me] = perfect_foresight_solver_core(M_,options_,oo_); if oo_.deterministic_simulation.status == 1 @@ -150,6 +156,7 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy end fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', me) else + % If solver failed, then go back. oo_.endo_simul = saved_endo_simul; success_counter = 0; step = step / 2; @@ -168,11 +175,12 @@ end if oo_.deterministic_simulation.status == 1 disp('Perfect foresight solution found.') - skipline() else warning('Failed to solve perfect foresight model') end +skipline() + dyn2vec; if ~isdates(options_.initial_period) && isnan(options_.initial_period) @@ -182,4 +190,4 @@ else end ts = dseries(transpose(oo_.endo_simul),initial_period,cellstr(M_.endo_names)); -assignin('base', 'Simulated_time_series', ts); +assignin('base', 'Simulated_time_series', ts); \ No newline at end of file