From 17ebcf4363b701f641eefe33146d07009e2a3dfb Mon Sep 17 00:00:00 2001
From: Marco Ratto <marco.ratto@ec.europa.eu>
Date: Sat, 16 Dec 2023 12:06:38 +0100
Subject: [PATCH] re introduce use_relaxation, it was not fully implicitly
 embedded in the new kalman update engine. now use_relaxation is embedded in
 the new update engine and is extended to work with 2 constraints.

---
 matlab/+occbin/kalman_update_engine.m         | 43 ++++++++++++++++---
 matlab/+occbin/set_default_options.m          |  1 +
 matlab/+occbin/setup.m                        |  1 +
 .../dynrbc_common.inc                         |  3 ++
 4 files changed, 43 insertions(+), 5 deletions(-)

diff --git a/matlab/+occbin/kalman_update_engine.m b/matlab/+occbin/kalman_update_engine.m
index 2f117977d5..7dac3dfb51 100644
--- a/matlab/+occbin/kalman_update_engine.m
+++ b/matlab/+occbin/kalman_update_engine.m
@@ -12,6 +12,8 @@ if nargin>26
     is_multivariate = false;
 end
 
+use_relaxation = false;
+
 if is_multivariate
     [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,struct(),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
 else
@@ -73,6 +75,31 @@ else
     end
 end
 
+diffstart=0;
+if info==0
+if M_.occbin.constraint_nbr==1
+    oldstart = regimes_(1).regimestart(end);
+    newstart = regx(1).regimestart(end);
+    diffstart = newstart-oldstart;
+else
+    newstart1 = regx(1).regimestart1(end);
+    newstart2 = regx(1).regimestart2(end);
+    oldstart1 = regimes_(1).regimestart1(end);
+    oldstart2 = regimes_(1).regimestart2(end);
+    diffstart = max(newstart1-oldstart1,newstart2-oldstart2);
+end
+end
+
+if options_.occbin.filter.use_relaxation && diffstart>2
+    if info0==0
+        % make sure we match criteria to enter further solution attempts
+        info1=1;
+    end
+    options_.occbin.likelihood.brute_force_regime_guess   = true;
+    options_.occbin.likelihood.loss_function_regime_guess = true;
+    use_relaxation = true;
+end
+
 if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (info==0 &&  ~isequal(regx(1),base_regime))
 
     guess_regime = [base_regime base_regime];
@@ -90,8 +117,9 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
             end
             if info2==0
                 use_index= 1;
-                if not(info==0 && isequal(regx2{1},regx))
-                    % found a solution, different from previous one 
+                if not(info==0 && isequal(regx2{1},regx)) && not(use_relaxation && likx2{1}>=likx)
+                    % found a solution, different from previous or
+                    % use_relaxation and likelihood is better
                     break
                 end
             end
@@ -105,6 +133,8 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
             end
             if info2==0
                 use_index = 2;
+                % if use_relaxation and we are here, previous guess did not
+                % improve solution, so we test for this one
             end
 
             if use_index
@@ -166,16 +196,18 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
                             end
                             if info2==0
                                 use_index= gindex;
-                                if not(info==0 && isequal(regx2{gindex},regx))
+                                if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx)
                                     % found a solution, different from previous one
+                                    % use_relaxation and likelihood improves
                                     break
                                 end
                             end
                         end % loop over other regime slack, binding in expectation or binding in current period
 
                         if info2==0
-                            if not(info==0 && isequal(regx2{gindex},regx))
+                            if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx)
                                 % found a solution, different from previous one
+                                % use_relaxation and likelihood improves
                                 break
                             end
                         end
@@ -183,8 +215,9 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
                     end % loop over current regime binding in expectation vs binding in current period
 
                     if info2==0
-                        if not(info==0 && isequal(regx2{gindex},regx))
+                        if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx)
                             % found a solution, different from previous one
+                            % use_relaxation and likelihood improves
                             break
                         end
                     end
diff --git a/matlab/+occbin/set_default_options.m b/matlab/+occbin/set_default_options.m
index 7f6c299353..d3dda39f3c 100644
--- a/matlab/+occbin/set_default_options.m
+++ b/matlab/+occbin/set_default_options.m
@@ -44,6 +44,7 @@ end
 if ismember(flag,{'filter','all'})
     options_occbin_.filter.state_covariance = false;
     options_occbin_.filter.guess_regime = false;
+    options_occbin_.filter.use_relaxation = false;
 end
 
 if ismember(flag,{'forecast','all'})
diff --git a/matlab/+occbin/setup.m b/matlab/+occbin/setup.m
index 4c239102d4..a050ee787b 100644
--- a/matlab/+occbin/setup.m
+++ b/matlab/+occbin/setup.m
@@ -47,6 +47,7 @@ options_ = occbin.set_option(options_,options_occbin_,'likelihood.maxit');
 options_ = occbin.set_option(options_,options_occbin_,'likelihood.check_ahead_periods');
 options_ = occbin.set_option(options_,options_occbin_,'likelihood.periodic_solution');
 options_ = occbin.set_option(options_,options_occbin_,'likelihood.max_number_of_iterations');
+options_ = occbin.set_option(options_,options_occbin_,'filter.use_relaxation');
 options_ = occbin.set_option(options_,options_occbin_,'likelihood.inversion_filter');
 
 if isfield(M_,'surprise_shocks') && ~isempty(M_.surprise_shocks)
diff --git a/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc b/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc
index 9f504605ef..51e22ec23e 100644
--- a/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc
+++ b/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc
@@ -89,6 +89,9 @@ occbin_graph(noconstant) c erra lambdak k i a k;
     occbin_solver(simul_periods=200,simul_maxit=200,simul_curb_retrench,simul_check_ahead_periods=200);
     occbin_setup(smoother_periods=200,smoother_maxit=200,smoother_curb_retrench,smoother_check_ahead_periods=200);
     calib_smoother(datafile=datasim);
+    oo0=oo_;
+    occbin_setup(filter_use_relaxation);
+    calib_smoother(datafile=datasim);
     oo_= occbin.unpack_simulations(M_,oo_,options_);
 
     titlelist = char('c','lambdak','k','i','a','erra');
-- 
GitLab