function [smoother_results, occbin_info] = DSGE_smoother(xparam1,M_,dr,steady_state,exo_steady_state,exo_det_steady_state,options_,estim_params_,dataset_,dataset_info,smoother_info,measurement_info)
% [smoother_results, occbin_info] = DSGE_smoother(xparam1,M_,dr,steady_state,exo_steady_state,exo_det_steady_state,options_,estim_params_,dataset_,dataset_info)
% Runs a DSGE smoother with occasionally binding constraints
%
% INPUTS
% - xparam1              [double]     (p*1) vector of (estimated) parameters.
% - M_                   [structure]  MATLAB's structure describing the model (M_).
% - dr                   [structure]  Decision rules structure.
% - steady_state         [double]     (m*1) vector of steady state values (declaration order).
% - exo_steady_state     [double]     (q*1) vector of exogenous steady state values.
% - exo_det_steady_state [double]     Vector of deterministic exogenous steady state values.
% - options_             [structure]  MATLAB's structure describing the current options (options_).
% - estim_params_        [structure]  characterizing parameters to be estimated
% - dataset_             [structure]  the dataset after required transformation
% - dataset_info         [structure]  Various information about the dataset (descriptive statistics and missing observations)
%
% OUTPUTS
% - smoother_results  [structure]  Struct with same layout as DsgeSmoother output:
%                   smoother_results.smoothed:
%                     .alphahat      (m*T)         smoothed states (a_{t|T}, decision-rule order, deviations from SS)
%                     .etahat        (r*T)         smoothed structural shocks
%                     .epsilonhat    (n*T)         smoothed measurement errors
%                     .alphahat0     (m*1)         smoothed initial state (a_{0|T})
%                     .V_tgT         (m*m*T)       smoothed state uncertainty (Var[alpha_t|Y_T])
%                     .V_0gT         (m*m)         initial state uncertainty (Var[alpha_0|Y_T])
%                     .V_eta         []            smoothed shock covariance (if requested)
%                   smoother_results.filtered:
%                     .ahat          (m*T)         updated (filtered) variables (a_{t|t})
%                     .P             (m*m*(T+1))   one-step ahead forecast error variance
%                     .Ptt           (m*m*T)       updated variance (if requested)
%                     .aK            (K*m*(T+K))   k-step ahead filtered variables
%                     .PK            (K*m*m*(T+K)) k-step ahead forecast error variance
%                     .decomp        (K*m*r*(T+K)) shock decomposition of k-step filtered
%                   smoother_results.model:
%                     .T             (m*m)         baseline transition matrix
%                     .R             (m*r)         baseline shock impact matrix
%                     .SteadyState   (m*1)         steady state (declaration order)
%                     .trend_coeff   (n*1)         trend slopes
%                     .trend_addition (n*T)        pure trend component (varobs order)
%                     .mf            (n*1)         observed variable indices (decision-rule order)
% - occbin_info       [structure]  OccBin-specific outputs. Pass to store_smoother_results
%                     as 8th argument to store all oo_.occbin.smoother.* and
%                     oo_.occbin.linear_smoother.* fields.
%                     Fields:
%                     .error_flag               [integer]    0=converged, 321=diverged, 322=failure
%                     .regime_history           [struct]     converged per-period regime sequence
%                     .realtime_regime_history  [struct]     first-pass (initial) regime history
%                     .TT                       {T}(m*m)    time-varying transition matrices
%                     .RR                       {T}(m*r)    time-varying shock impact matrices
%                     .CC                       {T}(m*1)    time-varying constant vectors
%                     .linear                   [struct]    linear smoother results (DsgeSmoother layout), or []
%                     .warning_flag             [integer]   convergence warning code (0/320/321/322/323)
%                     .is_converged             [boolean]   smoother iteration convergence flag
%                     .is_periodic              [boolean]   regime oscillation detected
%                     .regime_history_iter      [cell]      iteration-by-iteration regime log
%
%   Producers/consumers for smoothed states and shocks stored via store_smoother_results:
%     - smoother_results.smoothed.alphahat/etahat are raw (m*T)/(r*T) arrays in
%       decision-rule order (deviations from steady state). These get stored:
%       (a) by store_smoother_results into oo_.SmoothedVariables/oo_.SmoothedShocks
%           (converted to named structs in declaration order with trend+constant added), and
%       (b) into oo_.occbin.smoother.alphahat/etahat (raw DR-order arrays, deviations from SS).
%     - oo_.occbin.smoother.alphahat/etahat are consumed by:
%         occbin.shock_decomposition  — for piecewise-linear shock decomposition
%         occbin.irf                  — for OccBin impulse response functions
%     - oo_.occbin.smoother.T0/R0/C0 (time-varying matrices) are consumed by:
%         occbin.shock_decomposition  — for reconstructing state paths

% Copyright © 2021-2026 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

% ALGORITHM OVERVIEW
% The OccBin smoother iterates between Kalman smoothing and piecewise-linear
% simulation to find a self-consistent regime sequence under occasionally
% binding constraints:
%   1. Linear smoother (optional): baseline DsgeSmoother without constraints
%   2. Initial (realtime) smoother: DsgeSmoother with initial regime guess;
%      results snapshotted in 'realtime' struct for potential fallback
%   3. Regime discovery: simulate with smoothed shocks; compare implied regimes
%   4. conditional smoother (if needed): re-smooth conditioning on initial state
%      from step 2; results stored in 'conditional' struct
%   5. Iteration loop: alternate smoothing and simulation until convergence
%   6. Fallback: use realtime results if no convergence achieved
%
% State-space matrix naming:
%   T, R        - baseline (unconstrained, time-invariant) from DsgeSmoother
%   TT, RR, CC  - time-varying (piecewise-linear) from occbin.solver

Ptt = [];
V_eta = [];
linear_smoother = [];

if not(isempty(xparam1))
    M_ = set_all_parameters(xparam1,estim_params_,M_);
end
if options_.occbin.smoother.linear_smoother
    %% linear smoother
    options_.occbin.smoother.status=false;
    linear_smoother = DsgeSmoother(xparam1,dataset_,dataset_info,M_,dr,steady_state,exo_steady_state,exo_det_steady_state,options_,estim_params_,smoother_info);
    dr = linear_smoother.model.dr;
    disp_verbose('OccBin: linear smoother done.',options_.verbosity)
    options_.occbin.smoother.status=true;
end
% if init_mode
%% keep these commented lines for the moment, should some nasty problem in the future call back for this option
%     shocks1     = etahat(:,2:end)';
%     clear regime_history ;
%
%     for ii = 1:gend-1
%
%         fprintf('Period number %d out of %d \n',ii,gend-1)
%         opts_simul.SHOCKS = shocks1(ii,:);
%         opts_simul.endo_init = alphahat(dr.inv_order_var,ii);
%         out =  occbin.runsim_fn(opts_simul, M_, options_);
%         regime_history(ii) = out.regime_history;
%     end
%
%     [TT,RR] = dynare_resolve(M_,options_);
%     CC = zeros([size(RR,1),gend ]);
% else

opts_simul = options_.occbin.simul;
list_of_options_to_copy = {'curb_retrench','maxit','waitbar','periods', ...
    'check_ahead_periods','max_check_ahead_periods','periodic_solution', ...
    'full_output','piecewise_only'};
for option_iter=1:length(list_of_options_to_copy)
    opts_simul.(list_of_options_to_copy{option_iter}) = options_.occbin.smoother.(list_of_options_to_copy{option_iter});
end
occbin_options = struct();

occbin_options.first_period_occbin_update = options_.occbin.smoother.first_period_occbin_update;
occbin_options.opts_simul = opts_simul; % this builds the opts_simul options field needed by occbin.solver
occbin_options.opts_regime.binding_indicator = options_.occbin.smoother.init_binding_indicator;
occbin_options.opts_regime.regime_history=options_.occbin.smoother.init_regime_history;

options_.noprint = true;

is_realtime_smoother_successful = true;

%% Phase 2: Initial (realtime) smoother run with regime guess
smoother_opts = struct('occbin_options',occbin_options);
smoother_results = DsgeSmoother(xparam1,dataset_,dataset_info,M_,dr,steady_state,exo_steady_state,exo_det_steady_state,options_,estim_params_,smoother_info,smoother_opts);
% Unpack only variables used in expressions below; filtered.{ahat,aK,PK},
% filtered.{P,Ptt,decomp}, smoothed.{V_tgT,V_0gT,V_eta}, and
% model.{trend_coeff,trend_addition} stay in smoother_results.
alphahat   = smoother_results.smoothed.alphahat;
etahat     = smoother_results.smoothed.etahat;
epsilonhat = smoother_results.smoothed.epsilonhat;
alphahat0  = smoother_results.smoothed.alphahat0;
T = smoother_results.model.T;  R = smoother_results.model.R;
SteadyState = smoother_results.model.SteadyState;
dr = smoother_results.model.dr;
mf = smoother_results.model.mf;
error_indicator = smoother_results.info;
regime_history = smoother_results.regimes;

% Cache decision rules from the first occbin DsgeSmoother call.
% M_.params does not change after set_all_parameters, so all subsequent
% DsgeSmoother calls would recompute identical decision rules via
% dynare_resolve.  Passing them through smoother_opts.precomputed avoids
% this redundant work.
precomputed = struct('T',T,'R',R,'SteadyState',SteadyState,'dr',dr,'params',M_.params,'Sigma_e',M_.Sigma_e,'H',M_.H);
if options_.smoother_redux
    [~,~,~,~,~,~,~,~,~,T0_cached,R0_cached] = ...
        occbin.dynare_resolve(M_,options_,dr,steady_state,exo_steady_state,exo_det_steady_state,[],'restrict');
    precomputed.T0 = T0_cached;
    precomputed.R0 = R0_cached;
end

if error_indicator(1) || isempty(alphahat0)
    is_realtime_smoother_successful = false;
    if ~options_.occbin.smoother.linear_smoother
        options_.occbin.smoother.status=false;
        smoother_fallback = DsgeSmoother(xparam1,dataset_,dataset_info,M_,dr,steady_state,exo_steady_state,exo_det_steady_state,options_,estim_params_,smoother_info);
        etahat = smoother_fallback.smoothed.etahat;
        alphahat0 = smoother_fallback.smoothed.alphahat0;
        options_.occbin.smoother.status=true;
    else
        etahat = linear_smoother.smoothed.etahat;
        alphahat0 = linear_smoother.smoothed.alphahat0;
        smoother_results.smoothed.V_0gT = linear_smoother.smoothed.V_0gT;
    end
    base_regime = struct();
    if M_.occbin.constraint_nbr==1
        base_regime.regime = 0;
        base_regime.regimestart = 1;
    else
        base_regime.regime1 = 0;
        base_regime.regimestart1 = 1;
        base_regime.regime2 = 0;
        base_regime.regimestart2 = 1;
    end
    regime_history = [];
    for jper=1:size(alphahat,2)+1
        if jper == 1
            regime_history = base_regime;
        else
            regime_history(jper) = base_regime;
        end
    end
end

realtime_regime_history = regime_history;
opts_regime.regime_history = regime_history;

% ahat, aK, and PK are not overwritten by later iterations (their output
% positions use ~), so they naturally persist as output variables.

opts_regime.binding_indicator=[];
regime_history0 = regime_history;

disp_verbose('OccBin smoother iteration 1.',options_.verbosity)
opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(1,M_.exo_nbr)];
opts_simul.exo_pos = 1:M_.exo_nbr;
opts_simul.endo_init = alphahat0(dr.inv_order_var,1);
opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid multiple solution issues!
opts_simul.periods = size(opts_simul.SHOCKS,1);
options_.occbin.simul=opts_simul;
occbin_smoother_debug=options_.occbin.smoother.debug;
[~, out, ss] = occbin.solver(M_,options_,dr,steady_state,exo_steady_state,exo_det_steady_state);
error_flag=0;
is_last_simulation_converged = true;
if out.error_flag
    is_last_simulation_converged = false;
    disp_verbose('OccBin smoother:: simulation within smoother did not converge.',options_.verbosity)
    error_flag=321;
    % store regimes consistent with the last smoother run
    out.regime_history = regime_history ;
elseif not(isequal(out.regime_history(1:options_.occbin.likelihood.first_period_binding_regime_allowed-1),regime_history(1:options_.occbin.likelihood.first_period_binding_regime_allowed-1)))
    disp_verbose('Occbin smoother:: simulation violates first_period_binding_regime_allowed.',~options_.noprint)
    error_flag=322;
    % store regimes consistent with the last smoother run
    out.regime_history = regime_history ;
end
regime_history = out.regime_history;
if options_.smoother_redux
    occbin_options.opts_simul.restrict_state_space =1;
    [T,R] = dynare_resolve(M_,options_,dr,steady_state,exo_steady_state,exo_det_steady_state);
end
TT = ss.T(dr.order_var,dr.order_var,:);
RR = ss.R(dr.order_var,:,:);
CC = ss.C(dr.order_var,:);

opts_regime.regime_history = regime_history;
opts_regime.binding_indicator = [];
[TT, RR, CC, regime_history] = occbin.check_regimes(TT, RR, CC, opts_regime, M_, options_ , dr, steady_state, exo_steady_state, exo_det_steady_state);
is_changed = ~isequal(regime_history0,regime_history);
if isempty(regime_history0)
    regime_history0 = regime_history;
end
iter=1;

is_periodic = 0;
is_changed_start = 0;
maxiter = options_.occbin.smoother.max_number_of_iterations;
occbin_smoother_fast = options_.occbin.smoother.fast;

sto_alphahat=alphahat;
sto_etahat={etahat};
sto_CC = CC;
sto_RR = RR;
sto_TT = TT;
sto_eee=NaN(size(TT,1),size(TT,3));
for k=1:size(TT,3)
    sto_eee(:,k) = eig(TT(:,:,k));
end

is_conditional_smoother_converged = false;
is_realtime_smoother_converged = true;

if (is_changed || error_flag)
    is_realtime_smoother_converged = false;
end
if is_realtime_smoother_successful
    % Snapshot of the initial (realtime) smoother results for potential fallback.
    % Working locals:
    realtime.alphahat=alphahat;
    realtime.etahat=etahat;
    realtime.epsilonhat=epsilonhat;
    realtime.alphahat0=alphahat0;
    realtime.SteadyState=SteadyState;
    realtime.T=T;
    realtime.R=R;
    realtime.dr=dr;
    % Pass-through fields from smoother_results:
    realtime.P=smoother_results.filtered.P;
    realtime.decomp=smoother_results.filtered.decomp;
    realtime.trend_coeff=smoother_results.model.trend_coeff;
    realtime.trend_addition=smoother_results.model.trend_addition;
    realtime.state_uncertainty=smoother_results.smoothed.V_tgT;
    realtime.state_uncertainty0=smoother_results.smoothed.V_0gT;
    % Other:
    realtime.opts_simul=opts_simul;
    realtime.regime_history=regime_history;
end

if ~is_realtime_smoother_converged && not(options_.lik_init==2 && options_.Harvey_scale_factor==0)
    % try conditional smoother, where state uncertainty is smaller [unless too many shocks in excess]
    disp_verbose(sprintf('OccBin: try conditional smoother iteration.'),options_.verbosity)
    occbin_options.opts_regime.regime_history=regime_history;
    %%%% resort to conditional smoother
    alphahat1 = alphahat0(dr.inv_order_var)+SteadyState;
    alphahat1 = alphahat1(dr.state_var);
    M_local = M_;
    % set direct assignment of initial states
    M_local.endo_initial_state.status = true;
    M_local.endo_initial_state.values = zeros(M_.endo_nbr,1);
    opts_local1 = options_;
    opts_local1.lik_init = 2;
    opts_local1.Harvey_scale_factor = 0;
    M_local.endo_initial_state.values(dr.state_var) = alphahat1;
    occbin_options.first_period_occbin_update = 1;
    %% Phase 4: conditional smoother - re-smooth conditioning on initial state
    smoother_opts = struct('precomputed',precomputed,'occbin_options',occbin_options);
    smoother_cond = DsgeSmoother(xparam1,dataset_,dataset_info,M_local,dr,steady_state,exo_steady_state,exo_det_steady_state,opts_local1,estim_params_,smoother_info,smoother_opts);
    conditional.alphahat = smoother_cond.smoothed.alphahat;
    conditional.etahat = smoother_cond.smoothed.etahat;
    conditional.epsilonhat = smoother_cond.smoothed.epsilonhat;
    conditional.alphahat0 = smoother_cond.smoothed.alphahat0;
    conditional.T = smoother_cond.model.T;
    conditional.R = smoother_cond.model.R;
    conditional.P = smoother_cond.filtered.P;
    conditional.error_indicator = smoother_cond.info;
    conditional.regime_history = smoother_cond.regimes;
    if conditional.error_indicator(1)
        disp_verbose('OccBin smoother:: there was an error in running conditional smoother.',options_.verbosity)
        error_flag=322;
    else
        %%%% compute state uncertainty consistent with conditional smoother
        %%%% regime sequence
        regime_history = conditional.regime_history;
        opts_simul.SHOCKS = [conditional.etahat(:,1:end)'; zeros(1,M_.exo_nbr)];
        opts_simul.endo_init = alphahat0(dr.inv_order_var,1);
        options_.occbin.simul=opts_simul;
        out0 = out;
        ss0=ss;
        [~, out, ss] = occbin.solver(M_,options_,dr,steady_state,exo_steady_state,exo_det_steady_state);
        error_flag=0;
        if out.error_flag
            disp_verbose('OccBin smoother:: simulation within conditional smoother did not converge.',options_.verbosity)
            error_flag=321;
            % use regimes consistent with the last smoother run
            out = out0;
            out.regime_history = regime_history;
            ss=ss0;
        elseif not(isequal(out.regime_history(1:options_.occbin.likelihood.first_period_binding_regime_allowed-1),regime_history(1:options_.occbin.likelihood.first_period_binding_regime_allowed-1)))
            disp_verbose('Occbin smoother:: simulation violates first_period_binding_regime_allowed.',~options_.noprint)
            error_flag=322;
            % use regimes consistent with the last smoother run
            out.regime_history = regime_history ;
        end
        conditional.TT = ss.T(dr.order_var,dr.order_var,:);
        conditional.RR = ss.R(dr.order_var,:,:);
        conditional.CC = ss.C(dr.order_var,:);

        opts_regime.regime_history = out.regime_history;
        [conditional.TT, conditional.RR, conditional.CC, regime_history] = occbin.check_regimes(conditional.TT, conditional.RR, conditional.CC, opts_regime, M_, options_ , dr, steady_state, exo_steady_state, exo_det_steady_state);
        is_changed = ~isequal(conditional.regime_history,regime_history);
        if not(is_changed || error_flag)
            is_conditional_smoother_converged = true;
            occbin_options.opts_regime.regime_history=conditional.regime_history;
            occbin_options.first_period_occbin_update = inf;
            % compute uncertainty consistent with sequence of regimes of conditional smoother
            CC = conditional.CC;
            TT = conditional.TT;
            RR = conditional.RR;
            smoother_opts = struct('precomputed',precomputed,'occbin_options',occbin_options,'TT',TT,'RR',RR,'CC',CC);
            smoother_unc = DsgeSmoother(xparam1,dataset_,dataset_info,M_,dr,steady_state,exo_steady_state,exo_det_steady_state,options_,estim_params_,smoother_info,smoother_opts);
            % Update pass-through fields from uncertainty re-run
            smoother_results.filtered.P     = smoother_unc.filtered.P;
            smoother_results.filtered.Ptt   = smoother_unc.filtered.Ptt;
            smoother_results.filtered.decomp = smoother_unc.filtered.decomp;
            smoother_results.smoothed.V_tgT = smoother_unc.smoothed.V_tgT;
            smoother_results.smoothed.V_0gT = smoother_unc.smoothed.V_0gT;
            smoother_results.smoothed.V_eta = smoother_unc.smoothed.V_eta;
            % Working locals from conditional smoother
            alphahat = conditional.alphahat;
            etahat = conditional.etahat;
            epsilonhat = conditional.epsilonhat;
            T = conditional.T;
            R = conditional.R;
            alphahat0 = conditional.alphahat0;
            regime_history = conditional.regime_history;
        else
            disp_verbose(sprintf('OccBin: conditional smoother did not converge.'),options_.verbosity)
        end
    end
end


occbin_options.first_period_occbin_update = inf;

% try consistent estimation of shocks and regime sequence
while is_changed && maxiter>iter && ~is_periodic && is_last_simulation_converged
    if iter==1
        if is_realtime_smoother_successful
            regime_history = realtime.regime_history;
        end
        error_flag=0;
    end
    iter=iter+1;
    disp_verbose(sprintf('OccBin smoother iteration %u.', iter),options_.verbosity)
    occbin_options.opts_regime.regime_history=regime_history;
    smoother_opts = struct('precomputed',precomputed,'occbin_options',occbin_options,'TT',TT,'RR',RR,'CC',CC);
    smoother_iter = DsgeSmoother(xparam1,dataset_,dataset_info,M_,dr,steady_state,exo_steady_state,exo_det_steady_state,options_,estim_params_,smoother_info,smoother_opts);
    % Working locals
    alphahat   = smoother_iter.smoothed.alphahat;
    etahat     = smoother_iter.smoothed.etahat;
    epsilonhat = smoother_iter.smoothed.epsilonhat;
    alphahat0  = smoother_iter.smoothed.alphahat0;
    T = smoother_iter.model.T;  R = smoother_iter.model.R;
    SteadyState = smoother_iter.model.SteadyState;
    dr = smoother_iter.model.dr;
    % Update pass-through fields
    smoother_results.smoothed.V_tgT = smoother_iter.smoothed.V_tgT;
    smoother_results.smoothed.V_0gT = smoother_iter.smoothed.V_0gT;
    smoother_results.smoothed.V_eta = smoother_iter.smoothed.V_eta;
    smoother_results.filtered.P     = smoother_iter.filtered.P;
    smoother_results.filtered.Ptt   = smoother_iter.filtered.Ptt;
    smoother_results.filtered.decomp = smoother_iter.filtered.decomp;
    smoother_results.model.trend_coeff    = smoother_iter.model.trend_coeff;
    smoother_results.model.trend_addition = smoother_iter.model.trend_addition;
    sto_etahat(iter)={etahat};
    regime_history0(iter,:) = regime_history;
    if occbin_smoother_debug
        save('Occbin_smoother_debug_regime_history','regime_history0');
    end

    sto_CC = CC;
    sto_RR = RR;
    sto_TT = TT;

    opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(1,M_.exo_nbr)];
    opts_simul.endo_init = alphahat0(dr.inv_order_var,1);
    options_.occbin.simul=opts_simul;
    ss0=ss;
    [~, out, ss] = occbin.solver(M_,options_,dr,steady_state,exo_steady_state,exo_det_steady_state);
    if out.error_flag
        disp_verbose('OccBin smoother:: simulation within smoother did not converge.',options_.verbosity)
        error_flag=321;
        % use regimes consistent with the last smoother run
        out.regime_history = regime_history;
        break
    elseif not(isequal(out.regime_history(1:options_.occbin.likelihood.first_period_binding_regime_allowed-1),regime_history(1:options_.occbin.likelihood.first_period_binding_regime_allowed-1)))
        disp_verbose('Occbin smoother:: simulation violates first_period_binding_regime_allowed.',~options_.noprint)
        error_flag=322;
        % use regimes consistent with the last smoother run
        out.regime_history = regime_history ;
        break
    end
    regime_history = out.regime_history;
    TT = ss.T(dr.order_var,dr.order_var,:);
    RR = ss.R(dr.order_var,:,:);
    CC = ss.C(dr.order_var,:);

    opts_regime.regime_history = regime_history;
    [TT, RR, CC, regime_history] = occbin.check_regimes(TT, RR, CC, opts_regime, M_, options_ , dr, steady_state, exo_steady_state, exo_det_steady_state);
    is_changed = ~isequal(regime_history0(iter,:),regime_history);
    isdiff_regime=NaN(size(regime_history0,2),M_.occbin.constraint_nbr);
    isdiff_start=NaN(size(isdiff_regime));
    isdiff_=NaN(size(isdiff_regime));
    if M_.occbin.constraint_nbr==2
        for k=1:size(regime_history0,2)
            isdiff_regime(k,1) = ~isequal(regime_history0(end,k).regime1,regime_history(k).regime1);
            isdiff_start(k,1) = ~isequal(regime_history0(end,k).regimestart1,regime_history(k).regimestart1);
            isdiff_(k,1) =  isdiff_regime(k,1) ||  isdiff_start(k,1);
            isdiff_regime(k,2) = ~isequal(regime_history0(end,k).regime2,regime_history(k).regime2);
            isdiff_start(k,2) = ~isequal(regime_history0(end,k).regimestart2,regime_history(k).regimestart2);
            isdiff_(k,2) =  isdiff_regime(k,2) ||  isdiff_start(k,2);
        end
        is_changed_regime = any(isdiff_regime(:,1)) || any(isdiff_regime(:,2));
        is_changed_start = any(isdiff_start(:,1)) || any(isdiff_start(:,2));
    else
        for k=1:size(regime_history0,2)
            isdiff_regime(k,1) = ~isequal(regime_history0(end,k).regime,regime_history(k).regime);
            isdiff_start(k,1) = ~isequal(regime_history0(end,k).regimestart,regime_history(k).regimestart);
            isdiff_(k,1) =  isdiff_regime(k,1) ||  isdiff_start(k,1);
        end
        is_changed_regime = any(isdiff_regime(:,1));
        is_changed_start = any(isdiff_start(:,1));
    end
    if occbin_smoother_fast
        is_changed = is_changed_regime;
    end
    if iter>1
        for kiter=1:iter-1
            is_tmp = isequal(regime_history0(kiter,:),regime_history);
            if is_tmp
                break
            end
        end
        is_periodic = (is_changed && is_tmp);
    end

    if is_changed
        eee=NaN(size(TT,1),size(TT,3));
        for k=1:size(TT,3)
            eee(:,k) = eig(TT(:,:,k));
        end
    end

    if occbin_smoother_debug || (is_periodic && ~options_.noprint)
        regime_ = cell(0);
        regime_new = regime_;
        start_ = regime_;
        start_new = regime_;
        if M_.occbin.constraint_nbr==2
            indx_init_1 = find(isdiff_(:,1));
            if ~isempty(indx_init_1)
                qq={regime_history0(end,indx_init_1).regime1}';
                for j=1:length(qq), regime_(j,1) = {int2str(qq{j})}; end
                qq={regime_history0(end,indx_init_1).regimestart1}';
                for j=1:length(qq), start_(j,1) = {int2str(qq{j})}; end
                qq={regime_history(indx_init_1).regime1}';
                for j=1:length(qq), regime_new(j,1) = {int2str(qq{j})}; end
                qq={regime_history(indx_init_1).regimestart1}';
                for j=1:length(qq), start_new(j,1) = {int2str(qq{j})}; end
                disp('Time points where regime 1 differs')
                if ~isoctave
                    disp(table(indx_init_1, regime_, start_, regime_new, start_new))
                else % The table() function is not implemented in Octave, print something more or less equivalent (though much less readable)
                    disp(vertcat({'indx_init_1', 'regime_', 'start_', 'regime_new', 'start_new'}, ...
                                 horzcat(num2cell(indx_init_1), regime_, start_, regime_new, start_new)))
                end
            end

            indx_init_2 = find(isdiff_(:,2));
            if ~isempty(indx_init_2)
                regime_ = cell(0);
                regime_new = regime_;
                start_ = regime_;
                start_new = regime_;
                qq={regime_history0(end,indx_init_2).regime2}';
                for j=1:length(qq), regime_(j,1) = {int2str(qq{j})}; end
                qq={regime_history0(end,indx_init_2).regimestart2}';
                for j=1:length(qq), start_(j,1) = {int2str(qq{j})}; end
                qq={regime_history(indx_init_2).regime2}';
                for j=1:length(qq), regime_new(j,1) = {int2str(qq{j})}; end
                qq={regime_history(indx_init_2).regimestart2}';
                for j=1:length(qq), start_new(j,1) = {int2str(qq{j})}; end
                disp('Time points where regime 2 differs ')
                if ~isoctave
                    disp(table(indx_init_2, regime_, start_, regime_new, start_new))
                else % The table() function is not implemented in Octave, print something more or less equivalent (though much less readable)
                    disp(vertcat({'indx_init_2', 'regime_', 'start_', 'regime_new', 'start_new'}, ...
                                 horzcat(num2cell(indx_init_2), regime_, start_, regime_new, start_new)))
                end
            end
        else
            indx_init_1 = find(isdiff_(:,1));
            if ~isempty(indx_init_1)
                qq={regime_history0(end,indx_init_1).regime}';
                for j=1:length(qq), regime_(j,1) = {int2str(qq{j})}; end
                qq={regime_history0(end,indx_init_1).regimestart}';
                for j=1:length(qq), start_(j,1) = {int2str(qq{j})}; end
                qq={regime_history(indx_init_1).regime}';
                for j=1:length(qq), regime_new(j,1) = {int2str(qq{j})}; end
                qq={regime_history(indx_init_1).regimestart}';
                for j=1:length(qq), start_new(j,1) = {int2str(qq{j})}; end
                disp('Time points where regime differs')
                if ~isoctave
                    disp(table(indx_init_1, regime_, start_, regime_new, start_new))
                else % The table() function is not implemented in Octave, print something more or less equivalent (though much less readable)
                    disp(vertcat({'indx_init_1', 'regime_', 'start_', 'regime_new', 'start_new'}, ...
                                 horzcat(num2cell(indx_init_1), regime_, start_, regime_new, start_new)))
                end
            end

        end
    end
    if is_changed
        sto_alphahat=alphahat;
        sto_eee = eee;
    end
end
regime_history0(max(iter+1,1),:) = regime_history;
regime_history = regime_history0(end,:);
regime_history_iter = regime_history0;
if occbin_smoother_debug
    save('Occbin_smoother_debug_regime_history','regime_history0')
end

is_smoother_converged=false;
warning_flag = 0;

if ( (maxiter==iter || ~is_last_simulation_converged) && is_changed) || is_periodic || error_flag
    disp_verbose('occbin.DSGE_smoother: smoother did not converge.',options_.verbosity)
    disp_verbose('occbin.DSGE_smoother: The algorithm did not reach a fixed point for the smoothed regimes.',options_.verbosity)
    if is_periodic
        is_smoother_converged=true;
        error_flag=0;
        warning_flag = 321;
        disp_verbose('occbin.DSGE_smoother: For the periods indicated above, regimes loops between the "regime_" and the "regime_new_" pattern displayed above.',options_.verbosity)
        disp_verbose('occbin.DSGE_smoother: We provide smoothed shocks consistent with "regime_".',options_.verbosity)
    elseif ~is_realtime_smoother_successful
        disp_verbose('occbin.DSGE_smoother: The realtime smoother was not successful either.',options_.verbosity)
        disp_verbose('occbin.DSGE_smoother: The respective fields will be left empty.',options_.verbosity)
        regime_history = [];
        error_flag=322;
    end
else
    disp_verbose('occbin.DSGE_smoother: smoother converged.',options_.verbosity)
    is_smoother_converged=true;
    error_flag=0;
    if ~is_conditional_smoother_converged && ~is_realtime_smoother_converged
        warning_flag = 320;
        disp_verbose('occbin.DSGE_smoother: WARNING: algorithm converged to a different regime wrt realtime regime',options_.verbosity)
        disp_verbose('occbin.DSGE_smoother: WARNING: this usually indicates that for some period there is no unique regime dominance.',options_.verbosity)
    elseif occbin_smoother_fast && is_changed_start
        warning_flag = 323;
        disp_verbose('occbin.DSGE_smoother: WARNING: fast algo is used, regime duration was not forced to converge',options_.verbosity)
        disp_verbose('occbin.DSGE_smoother: WARNING: this usually indicates that for some period there is no unique regime dominance.',options_.verbosity)
    end
end

if is_realtime_smoother_successful && ~is_conditional_smoother_converged && ~is_realtime_smoother_converged && ~is_smoother_converged
    % fallback solution
    disp_verbose('occbin.DSGE_smoother: We provide results consistent with realtime regime',options_.verbosity)
    disp_verbose('occbin.DSGE_smoother: WARNING: this usually indicates that for some period there is no unique regime dominance.',options_.verbosity)
    warning_flag = 322;
    error_flag = 0;
    options_.occbin.simul=realtime.opts_simul;
    options_.occbin.simul.maxit=1;
    [~, out] = occbin.solver(M_,options_,dr,steady_state,exo_steady_state,exo_det_steady_state);
    regime_history = realtime_regime_history;
    is_changed=false;
    % Restore working locals
    alphahat=realtime.alphahat;
    etahat=realtime.etahat;
    epsilonhat=realtime.epsilonhat;
    alphahat0=realtime.alphahat0;
    SteadyState=realtime.SteadyState;
    T=realtime.T;
    R=realtime.R;
    dr=realtime.dr;
    % Restore pass-through fields
    smoother_results.filtered.P              = realtime.P;
    smoother_results.filtered.decomp         = realtime.decomp;
    smoother_results.model.trend_coeff       = realtime.trend_coeff;
    smoother_results.model.trend_addition    = realtime.trend_addition;
    smoother_results.smoothed.V_tgT          = realtime.state_uncertainty;
    smoother_results.smoothed.V_0gT          = realtime.state_uncertainty0;
end

% When debug mode and iterations didn't converge, revert to second-to-last
if is_changed && occbin_smoother_debug && iter>1 && oo_.occbin.smoother.error_flag == 0
    CC = sto_CC;
    RR = sto_RR;
    TT = sto_TT;
    regime_history = regime_history0(end-1,:);
end

% Write working locals back into the output struct
smoother_results.smoothed.alphahat  = alphahat;
smoother_results.smoothed.etahat    = etahat;
smoother_results.smoothed.epsilonhat = epsilonhat;
smoother_results.smoothed.alphahat0 = alphahat0;
smoother_results.model.T  = T;
smoother_results.model.R  = R;
smoother_results.model.SteadyState = SteadyState;
smoother_results.model.mf = mf;
% filtered.{ahat,aK,PK,P,Ptt,decomp}, smoothed.{V_tgT,V_0gT,V_eta},
% model.{trend_coeff,trend_addition} are already current in smoother_results.

% OccBin-specific comparison plot (piecewise vs linear smoothed shocks)
if options_.occbin.smoother.plot
    linear_etahat = [];
    if ~isempty(linear_smoother)
        linear_etahat = linear_smoother.smoothed.etahat;
    end
    occbin.smoother_plots(M_, options_, etahat, linear_etahat);
end

% Pack occbin_info (all OccBin-specific outputs)
occbin_info.error_flag = error_flag;
occbin_info.regime_history = regime_history;
occbin_info.realtime_regime_history = realtime_regime_history;
occbin_info.TT = TT;
occbin_info.RR = RR;
occbin_info.CC = CC;
occbin_info.linear = linear_smoother;
occbin_info.warning_flag = warning_flag;
occbin_info.is_converged = is_smoother_converged;
occbin_info.is_periodic = logical(is_periodic);
occbin_info.regime_history_iter = regime_history_iter;
% Simulation results from occbin.solver (piecewise/linear paths)
if isfield(out,'piecewise')
    occbin_info.simul.piecewise = out.piecewise(1:end-1,:);
end
if ~options_.occbin.simul.piecewise_only && isfield(out,'linear')
    occbin_info.simul.linear = out.linear(1:end-1,:);
end

