Fixed bug in homotopy.

Convex combination for initial/terminal endogenous variables paths was wrong.

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
+    % Set initial paths for the endogenous and exogenous variables.
-    step = .5;
+    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*endosim(:,[initperiods, lastperiods])+(1-new_weight)*endoinit(:,[initperiods, lastperiods]);
-        oo_.endo_simul(:,[initperiods, lastperiods]) = ...
+
-            new_weight*oo_.endo_simul(:,[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);