diff --git a/matlab/+occbin/IVF_core.m b/matlab/+occbin/IVF_core.m
index 30666cd3220ee8bcd99793c54a9f9442514e4f45..a4c26fd0ad0a082f8e3e4651f84769fa5832649f 100644
--- a/matlab/+occbin/IVF_core.m
+++ b/matlab/+occbin/IVF_core.m
@@ -1,4 +1,4 @@
-function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
+function [filtered_errs, resids, Emat, stateval, error_code, regime_history] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
 % [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
 % Calls the solver to get the shocks explaining the data for the inversion filter
 %
@@ -99,9 +99,10 @@ for this_period=1:sample_length
 
     opts_simul.SHOCKS = err_vals_out;
 
-    [ resids(this_period,inan), ~, stateval(this_period,:), Emat(:,inan,this_period), M_] = occbin.match_function(...
+    [ resids(this_period,inan), ~, stateval(this_period,:), Emat(:,inan,this_period), M_, out] = occbin.match_function(...
         err_vals_out,obs_list,current_obs,opts_simul, M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_);
     init_val = stateval(this_period,:); %update
+    regime_history(this_period) = out.regime_history(1);
     if max(abs(err_vals_out))>1e8
         error_code(1) = 306;
         error_code(4) = max(abs(err_vals_out))/1000;
diff --git a/matlab/+occbin/IVF_posterior.m b/matlab/+occbin/IVF_posterior.m
index 69c6c110ddb34815d31c8df8701ae1e24c5fe29a..77329eb24ebc5503b91e147af8dc62bfcedc1fdf 100644
--- a/matlab/+occbin/IVF_posterior.m
+++ b/matlab/+occbin/IVF_posterior.m
@@ -1,4 +1,4 @@
-function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov] = IVF_posterior(xparam1,...
+function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov, regime_history] = IVF_posterior(xparam1,...
     dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
 % [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov] = IVF_posterior(xparam1,...
 %     dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
@@ -107,7 +107,7 @@ end
 sample_length = size(obs,1);
 filtered_errs_init = zeros(sample_length,sum(err_index));
 
-[filtered_errs, resids, Emat, stateval, info] = occbin.IVF_core(M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,err_index,filtered_errs_init,obs_list,obs);
+[filtered_errs, resids, Emat, stateval, info, regime_history] = occbin.IVF_core(M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,err_index,filtered_errs_init,obs_list,obs);
 if info(1)
     fval = Inf;
     exit_flag = 0;
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index 52a51a4cf5862903152be09b763afef77b5c07d7..a1a01039280d585d8009bea5f34d4244804ef580 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -189,7 +189,7 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.
     if options_.order==1 && ~options_.particle.status
         if options_.smoother
             if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
-                [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+                [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov, oo_.occbin.smoother.regime_history] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
                 if ismember(info(1),[303,304,306])
                     fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n')
                 else
diff --git a/matlab/kalman/save_display_classical_smoother_results.m b/matlab/kalman/save_display_classical_smoother_results.m
index 48dcad0b5b31a0e7eae91bbe043acf70331dac18..f7f2f981e710722c51b9f8ff466d92773f743cac 100644
--- a/matlab/kalman/save_display_classical_smoother_results.m
+++ b/matlab/kalman/save_display_classical_smoother_results.m
@@ -36,7 +36,7 @@ gend= dataset_.nobs;
 n_varobs = length(options_.varobs);
 
 if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
-    [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+    [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov, oo_.occbin.smoother.regime_history] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
     if ismember(info(1),[303,304,306])
         fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n')
     else