From 4caebd76c9872bdba7f57cf54e89f35c8e605d31 Mon Sep 17 00:00:00 2001
From: Marco Ratto <marco.ratto@ec.europa.eu>
Date: Wed, 21 Jul 2021 18:03:23 +0200
Subject: [PATCH] added 2 new occbin.simul options (defaults preserve current
 behavior): 1) algo_truncation: sets the number of iterations for a truncated
 guess-verify algorithm (i.e. if max_iter<=algo_truncation, no error is
 triggered on output, but the user is happy with the last regime in the
 algorithm) 2) reset_regime_in_new_period: is set to true, it resets the guess
 regimes to unconstrained when a new shock arrives, instead of starting with a
 guess regime consistent with the one identified in previous time periods.
 this sometimes allows more robust convergence/identification of the regimes

---
 matlab/+occbin/set_default_options.m   |  2 ++
 matlab/+occbin/solve_one_constraint.m  | 35 +++++++++++++++--------
 matlab/+occbin/solve_two_constraints.m | 39 ++++++++++++++++++--------
 3 files changed, 54 insertions(+), 22 deletions(-)

diff --git a/matlab/+occbin/set_default_options.m b/matlab/+occbin/set_default_options.m
index a3b437cd16..568fbc5418 100644
--- a/matlab/+occbin/set_default_options.m
+++ b/matlab/+occbin/set_default_options.m
@@ -160,6 +160,7 @@ if ismember(flag,{'shock_decomp','all'})
 end
 
 if ismember(flag,{'simul','all'})
+    options_occbin_.simul.algo_truncation = 1;
     options_occbin_.simul.debug = false;
     options_occbin_.simul.curb_retrench=false;
     options_occbin_.simul.endo_init=zeros(M_.endo_nbr,1);
@@ -174,6 +175,7 @@ if ismember(flag,{'simul','all'})
     options_occbin_.simul.check_ahead_periods=200;
     options_occbin_.simul.periodic_solution=true;
     options_occbin_.simul.piecewise_only = false;
+    options_occbin_.simul.reset_regime_in_new_period = false;
     options_occbin_.simul.restrict_state_space=false;
     options_occbin_.simul.SHOCKS=zeros(1,M_.exo_nbr);
     options_occbin_.simul.waitbar=true;
diff --git a/matlab/+occbin/solve_one_constraint.m b/matlab/+occbin/solve_one_constraint.m
index 34b2e25ba2..0801f8c1eb 100644
--- a/matlab/+occbin/solve_one_constraint.m
+++ b/matlab/+occbin/solve_one_constraint.m
@@ -163,6 +163,15 @@ for shock_period = 1:n_shocks_periods
         
         % get the hypothesized piece wise linear solution
         if shock_period==1 || shock_period>1 && any(data.shocks_sequence(shock_period,:))
+            if iter==1 && opts_simul_.reset_regime_in_new_period
+                binding_indicator=false(size(binding_indicator));
+                binding_indicator_history{iter}=binding_indicator;
+                % analyse violvec and isolate contiguous periods in the other regime.
+                [regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug);
+                regime_history(shock_period).regime = regime;
+                regime_history(shock_period).regimestart = regime_start;
+            end
+            
             [zdatalinear_, SS_out.T(:,:,shock_period), SS_out.R(:,:,shock_period), SS_out.C(:,shock_period), SS, update_flag]=occbin.mkdatap_anticipated_dyn(nperiods_0,DM,...
                 regime_start(end)-1,binding_indicator,...
                 data.exo_pos,data.shocks_sequence(shock_period,:),endo_init,update_flag);
@@ -243,22 +252,26 @@ for shock_period = 1:n_shocks_periods
         
     end
     
-    if regime_change_this_iteration ==1 && max_iter>1
-        disp_verbose(['occbin solver:: period ' int2str(shock_period) '::'],opts_simul_.debug)
-        if is_periodic
-            disp_verbose('Occbin solver loops between two regimes.',opts_simul_.debug)
-            if periodic_solution
-                disp_verbose(['Max error:' num2str(merr) '.'],opts_simul_.debug)
+    if regime_change_this_iteration ==1 
+        if max_iter>opts_simul_.algo_truncation
+            disp_verbose(['occbin solver:: period ' int2str(shock_period) '::'],opts_simul_.debug)
+            if is_periodic
+                disp_verbose('Occbin solver loops between two regimes.',opts_simul_.debug)
+                if periodic_solution
+                    disp_verbose(['Max error:' num2str(merr) '.'],opts_simul_.debug)
+                else
+                    if opts_simul_.waitbar; dyn_waitbar_close(hh); end
+                    error_flag = 310;
+                    return
+                end
             else
+                disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
                 if opts_simul_.waitbar; dyn_waitbar_close(hh); end
-                error_flag = 310;
+                error_flag = 311;
                 return
             end
         else
-            disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
-            if opts_simul_.waitbar; dyn_waitbar_close(hh); end
-            error_flag = 311;
-            return
+            binding_indicator = binding_indicator_history{end};
         end
     end
     if any(error_code_period)
diff --git a/matlab/+occbin/solve_two_constraints.m b/matlab/+occbin/solve_two_constraints.m
index 957551d8d4..234aaa78b5 100644
--- a/matlab/+occbin/solve_two_constraints.m
+++ b/matlab/+occbin/solve_two_constraints.m
@@ -177,6 +177,17 @@ for shock_period = 1:n_shocks_periods
         regime_history(shock_period).regime2 = regime_2;
         regime_history(shock_period).regimestart2 = regime_start_2;
         if shock_period==1 || shock_period>1 && any(data.shocks_sequence(shock_period,:)) % first period or shock happening
+            if iter==1 && opts_simul_.reset_regime_in_new_period
+                binding_indicator=false(size(binding_indicator));
+                binding_indicator_history{iter}=binding_indicator;
+                % analyse violvec and isolate contiguous periods in the other regime.
+                [regime_1, regime_start_1, error_code_period(1)]=occbin.map_regime(binding_indicator(:,1),opts_simul_.debug);
+                regime_history(shock_period).regime1 = regime_1;
+                regime_history(shock_period).regimestart1 = regime_start_1;
+                [regime_2, regime_start_2, error_code_period(2)]=occbin.map_regime(binding_indicator(:,2),opts_simul_.debug);
+                regime_history(shock_period).regime2 = regime_2;
+                regime_history(shock_period).regimestart2 = regime_start_2;
+            end
             Tmax=max([regime_start_1(end) regime_start_2(end)])-1;
             [zdatalinear_, SS_out.T(:,:,shock_period), SS_out.R(:,:,shock_period), SS_out.C(:,shock_period), SS, update_flag]=occbin.mkdatap_anticipated_2constraints_dyn(nperiods_0,...
                 DM,Tmax,...
@@ -271,22 +282,28 @@ for shock_period = 1:n_shocks_periods
             binding_indicator_history{iter}=binding_indicator;
         end
     end
-    if regime_change_this_iteration && max_iter>1
-        disp_verbose(['occbin solver: period ' int2str(shock_period) ':'],opts_simul_.debug)
-        if is_periodic
-            disp_verbose('Occbin solver loops between two regimes.',opts_simul_.debug)
-            if periodic_solution
-                disp_verbose(['Max error:' num2str(min_err) '.'],opts_simul_.debug)
+    if regime_change_this_iteration
+        if max_iter>opts_simul_.algo_truncation
+            disp_verbose(['occbin solver: period ' int2str(shock_period) ':'],opts_simul_.debug)
+            if is_periodic
+                disp_verbose('Occbin solver loops between two regimes.',opts_simul_.debug)
+                if periodic_solution
+                    disp_verbose(['Max error:' num2str(min_err) '.'],opts_simul_.debug)
+                else
+                    error_flag = 310;
+                    if opts_simul_.waitbar; dyn_waitbar_close(hh); end
+                    return;
+                end
             else
-                error_flag = 310;
+                disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
+                error_flag = 311;
                 if opts_simul_.waitbar; dyn_waitbar_close(hh); end
                 return;
             end
         else
-            disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
-            error_flag = 311;
-            if opts_simul_.waitbar; dyn_waitbar_close(hh); end
-            return;
+            % if max_iter <= truncation, we force indicator to equal the
+            % last guess
+            binding_indicator = binding_indicator_history{end};
         end
     end
     if any(error_code_period)
-- 
GitLab