function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,smoother_results,occbin_info)
% [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,smoother_results,occbin_info)
% Writes the smoother results into respective fields in oo_
%
% Inputs:
%   M_              [structure]     storing the model information
%   oo_             [structure]     storing the results
%   options_        [structure]     storing the options
%   bayestopt_      [structure]     storing information about priors
%   dataset_        [structure]     storing the dataset
%   dataset_info    [structure]     dataset descriptive statistics
%   smoother_results        [structure]     output from DsgeSmoother with fields:
%                   smoother_results.smoothed:
%                     .alphahat      [double]  (m*T) smoothed states alpha_{t|T} (decision-rule order)
%                     .etahat        [double]  (r*T) smoothed structural shocks eta_{t|T}
%                     .epsilonhat    [double]  (n*T) smoothed measurement errors epsilon_{t|T}
%                     .alphahat0     [double]  (m*1) smoothed initial state alpha_{0|T} (decision-rule order)
%                     .V_tgT         [double]  (m*m*T) Var(alpha_t|Y_{1:T}) state uncertainty (decision-rule order)
%                     .V_0gT         [double]  (m*m) Var(alpha_0|Y_{1:T}) state uncertainty at t=0 (decision-rule order)
%                     .V_eta         [double]  (r*r*T) Var(eta_t|Y_{1:T})
%                   smoother_results.filtered:
%                     .ahat          [double]  (m*T) updated variables a_{t|t} (decision-rule order)
%                     .P             [double]  (m*m*(T+1)) one-step ahead forecast error variance (decision-rule order)
%                     .Ptt           [double]  (m*m*T) updated covariance P_{t|t} (decision-rule order)
%                     .aK            [double]  (K*m*(T+K)) k-step ahead filtered variables (decision-rule order)
%                     .PK            [double]  (K*m*m*(T+K)) k-step ahead forecast error variance (decision-rule order)
%                     .decomp        [double]  (K*m*r*(T+K)) shock decomposition of k-step ahead (decision-rule order)
%                   smoother_results.model:
%                     .SteadyState   [double]  (m*1) steady state (declaration order)
%                     .trend_coeff   [double]  (n*1) trend slope for each observed variable
%                     .trend_addition [double] (n*T) deterministic trend component; stored in options_.varobs order
%                   smoother_results.d, smoother_results.info, smoother_results.regimes
%   occbin_info     [structure]     (optional) OccBin-specific output from occbin.DSGE_smoother.
%                   When provided, all oo_.occbin.smoother.* and oo_.occbin.linear_smoother.*
%                   fields are stored internally. Fields:
%                     .error_flag               [integer]    0=converged, 321=diverged, 322=failure
%                     .regime_history           [struct]     converged per-period regime sequence
%                     .realtime_regime_history  [struct]     first-pass 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
%                     .is_converged             [logical]   true if smoother converged
%                     .is_periodic              [logical]   true if regime sequence is periodic
%                     .regime_history_iter      [cell]      iteration-by-iteration regime log
%                     .simul                    [struct]    (optional) simulation paths from occbin.solver
%                       .piecewise              [double]    piecewise-linear simulation
%                       .linear                 [double]    linear simulation (unless piecewise_only)
%
% Outputs:
%   oo_             [structure] storing the results:
%                   oo_.Smoother.SteadyState: steady states (declaration order)
%                   oo_.Smoother.TrendCoeffs: trend coefficients, with zeros where no trend applies (declaration order)
%                   oo_.Smoother.Variance: one-step ahead forecast error variance (declaration order)
%                   oo_.Smoother.UpdatedVariance: updated covariance (declaration order)
%                   oo_.Smoother.Shock_variance: smoothed shock variance Var(eta_t|Y_{1:T})
%                   oo_.Smoother.Constant: structure storing the constant term of the smoother_results
%                   oo_.Smoother.Trend: structure storing the trend term of the smoother_results
%                   oo_.Smoother.Init_State: smoothed initial state by variable name
%                   oo_.Smoother.Init_State_Vector: smoothed initial state vector (declaration order)
%                   oo_.Smoother.State_uncertainty: smoothed state uncertainty (declaration order)
%                   oo_.Smoother.Init_State_uncertainty: smoothed state uncertainty at t=0 (declaration order)
%                   oo_.FilteredVariablesKStepAhead: k-step ahead filtered variables (declaration order)
%                   oo_.FilteredVariablesKStepAheadVariances: k-step ahead forecast error variance matrices (declaration order)
%                   oo_.FilteredVariablesShockDecomposition: shock decomposition of k-step ahead filtered variables (declaration order)
%                   oo_.SmoothedVariables: structure storing the smoothed variables
%                   oo_.FilteredVariables: structure storing the filtered variables
%                   oo_.UpdatedVariables: structure storing the updated variables
%                   oo_.SmoothedShocks: structure storing the smoothed shocks
%                   oo_.SmoothedShockVariances: structure storing the smoothed shock variances
%                   oo_.SmoothedMeasurementErrors: structure storing the smoothed measurement errors
%   yf              [double]    (nvarobs*T) matrix storing the smoothed observed variables (order of options_.varobs)
%
% Notes:
%   m:  number of endogenous variables (M_.endo_nbr)
%   T:  number of Time periods (options_.nobs)
%   r:  number of structural shocks (M_.exo_nbr)
%   n:  number of observables (length(options_.varobs))
%   K:  maximum forecast horizon (max(options_.nk))
%
%   First all smoothed variables are saved without trend and constant.
%       Then trend and constant are added for the observed variables.

% Copyright © 2014-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/>.

% Unpack struct into local variables used throughout this function
atT                = smoother_results.smoothed.alphahat;
innov              = smoother_results.smoothed.etahat;
measurement_error  = smoother_results.smoothed.epsilonhat;
a0T                = smoother_results.smoothed.alphahat0;
state_uncertainty  = smoother_results.smoothed.V_tgT;
state_uncertainty0 = smoother_results.smoothed.V_0gT;
V_eta              = smoother_results.smoothed.V_eta;
updated_variables  = smoother_results.filtered.ahat;
P                  = smoother_results.filtered.P;
Ptt                = smoother_results.filtered.Ptt;
aK                 = smoother_results.filtered.aK;
PK                 = smoother_results.filtered.PK;
decomp             = smoother_results.filtered.decomp;
ys                 = smoother_results.model.SteadyState;
trend_coeff        = smoother_results.model.trend_coeff;
Trend              = smoother_results.model.trend_addition;

gend=dataset_.nobs;
if isempty(Trend)
    Trend=zeros(length(options_.varobs),gend);
end


%make sure there are no stale results
field_names={'Smoother','SmoothedVariables','UpdatedVariables','FilteredVariables','FilteredVariablesKStepAhead','FilteredVariablesShockDecomposition','FilteredVariablesKStepAheadVariances','SmoothedShocks','SmoothedMeasurementErrors','SmoothedShockVariances'};
for field_iter=1:length(field_names)
    if isfield(oo_,field_names(field_iter))
        oo_=rmfield(oo_,field_names(field_iter));
    end
end

if options_.occbin.smoother.status
    oo_.Smoother.occbin = true;
else
    oo_.Smoother.occbin = false;
end


if options_.loglinear
    oo_.Smoother.loglinear = true;
    if ~isempty(M_.aux_vars) % deal with lead/lag of exogenous variables
        exo_lead_lag_index = M_.orig_endo_nbr + find(([M_.aux_vars.type] == 2) | ([M_.aux_vars.type]  == 3));
    else
        exo_lead_lag_index =[];
    end

else
    oo_.Smoother.loglinear = false;
end
%% write variable steady state
oo_.Smoother.SteadyState = ys;

%% write trend coefficients and trend
oo_.Smoother.TrendCoeffs = zeros(size(ys));

all_varobs_id = options_.varobs_id;
if ~isempty(M_.filter_tunes)
    all_varobs_id = [all_varobs_id unique([M_.filter_tunes.endo_id])];
end
oo_.Smoother.TrendCoeffs(all_varobs_id) = trend_coeff; %are in order of options_.varobs followed by filter tunes

if ~isempty(Trend)
    for var_iter=1:M_.endo_nbr
        if isempty(strmatch(M_.endo_names{var_iter}, options_.varobs, 'exact'))
            oo_.Smoother.Trend.(M_.endo_names{var_iter}) = zeros(gend,1);
        end
    end
    for var_iter=1:length(options_.varobs)
        oo_.Smoother.Trend.(options_.varobs{var_iter}) = Trend(var_iter,:)';
    end
end
%% Compute constant for observables
if options_.prefilter == 1 %as mean is taken after log transformation, no distinction is needed here
    constant_part=repmat(dataset_info.descriptive.mean',1,gend);
elseif options_.prefilter == 0 && options_.loglinear %logged steady state must be used
    constant_part=repmat(log(ys(bayestopt_.mfys)),1,gend);
elseif options_.prefilter == 0 && ~options_.loglinear %unlogged steady state must be used
    constant_part=repmat(ys(bayestopt_.mfys),1,gend);
end

%% get observed variables including trend and constant
trend_constant_observables=constant_part+Trend;
yf = atT(bayestopt_.mf,:)+trend_constant_observables;

if options_.nk > 0
    %filtered variable E_t(y_t+k) requires to shift trend by k periods
    filter_steps_required=union(1,options_.filter_step_ahead); % 1 is required for standard filtered variables
    for filter_iter=1:length(filter_steps_required)
        filter_step=filter_steps_required(filter_iter);
        trend_constant_observables_filtered.(['filter_ahead_' num2str(filter_step)])=constant_part+[Trend+repmat(filter_step*trend_coeff,1,gend)];
    end
end
%% write Smoother variance
if options_.filter_covariance
    oo_.Smoother.Variance = P;
end
if ~isempty(Ptt) && options_.updated_covariance
    oo_.Smoother.UpdatedVariance = Ptt;
end
if ~isempty(V_eta) && options_.smoothed_shock_variance
    oo_.Smoother.Shock_variance = V_eta;
end

%get indices of smoothed variables
i_endo_in_bayestopt_smoother_varlist = bayestopt_.smoother_saved_var_list;
i_endo_in_dr_matrices=bayestopt_.smoother_var_list(i_endo_in_bayestopt_smoother_varlist);
if ~isempty(options_.nk) && options_.nk ~= 0
    %% Compute constant
    i_endo_declaration_order = oo_.dr.order_var(i_endo_in_dr_matrices); %get indices of smoothed variables in name vector
    if  options_.loglinear %logged steady state must be used
        constant_all_variables=repmat(log(ys(i_endo_declaration_order))',[length(options_.filter_step_ahead),1,gend+max(options_.filter_step_ahead)]);
        if ~isempty(exo_lead_lag_index) && any(ismember(i_endo_declaration_order,exo_lead_lag_index)) % deal with lead/lag of exogenous variables
            constant_all_variables(:,ismember(i_endo_declaration_order,exo_lead_lag_index),:)=repmat(ys(i_endo_declaration_order(ismember(i_endo_declaration_order,exo_lead_lag_index)))',[length(options_.filter_step_ahead),1,gend+max(options_.filter_step_ahead)]);
        end
    elseif ~options_.loglinear %unlogged steady state must be used
        constant_all_variables=repmat((ys(i_endo_declaration_order))',[length(options_.filter_step_ahead),1,gend+max(options_.filter_step_ahead)]);
    end
    % add constant
    oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead,i_endo_in_dr_matrices,:)+constant_all_variables;
    if ~isempty(PK) && options_.filter_covariance %get K-step ahead variances
        oo_.FilteredVariablesKStepAheadVariances = ...
            PK(options_.filter_step_ahead,i_endo_in_dr_matrices,i_endo_in_dr_matrices,:);
    end
    if ~isempty(decomp) %get decomposition
        oo_.FilteredVariablesShockDecomposition = ...
            decomp(options_.filter_step_ahead,i_endo_in_dr_matrices,:,:);
    end
end

for i_endo_in_bayestopt_smoother_varlist=bayestopt_.smoother_saved_var_list'
    i_endo_in_dr=bayestopt_.smoother_var_list(i_endo_in_bayestopt_smoother_varlist);
    i_endo_declaration_order = oo_.dr.order_var(i_endo_in_dr); %get indices of smoothed variables in name vector
    %% Compute constant
    if  options_.loglinear == 1 %logged steady state must be used
        if ~isempty(exo_lead_lag_index) && ismember(i_endo_declaration_order,exo_lead_lag_index)% deal with lead/lag of exogenous variables
            constant_current_variable=repmat(ys(i_endo_declaration_order),gend,1);
        else
            constant_current_variable=repmat(log(ys(i_endo_declaration_order)),gend,1);
        end
    elseif options_.loglinear == 0 %unlogged steady state must be used
        constant_current_variable=repmat((ys(i_endo_declaration_order)),gend,1);
    end
    oo_.Smoother.Constant.(M_.endo_names{i_endo_declaration_order})=constant_current_variable;
    if ~isempty(a0T)
        oo_.Smoother.Init_State.(M_.endo_names{i_endo_declaration_order})=a0T(i_endo_in_dr)+constant_current_variable(1);
    end
    oo_.SmoothedVariables.(M_.endo_names{i_endo_declaration_order})=atT(i_endo_in_dr,:)'+constant_current_variable;
    if ~isempty(options_.nk) && options_.nk > 0 % && ~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file))
        oo_.FilteredVariables.(M_.endo_names{i_endo_declaration_order})=squeeze(aK(1,i_endo_in_dr,2:end-(options_.nk-1)))+constant_current_variable;
    end
    oo_.UpdatedVariables.(M_.endo_names{i_endo_declaration_order})=updated_variables(i_endo_in_dr,:)'+constant_current_variable;
end
if ~isempty(a0T)
    oo_.Smoother.Init_State_Vector(:,1)=a0T(oo_.dr.inv_order_var)+ys;
end

%% Add trend and constant for observed variables
for pos_iter=1:length(bayestopt_.mf)
    oo_.Smoother.Constant.(M_.endo_names{bayestopt_.mfys(pos_iter)})=constant_part(pos_iter,:)';
    if ismember(bayestopt_.mf(pos_iter),bayestopt_.smoother_var_list(bayestopt_.smoother_saved_var_list))
        oo_.SmoothedVariables.(M_.endo_names{bayestopt_.mfys(pos_iter)})=yf(pos_iter,:)';
        if ~isempty(options_.nk) && options_.nk > 0
            %filtered variable E_t(y_t+1) requires to shift trend by 1 period
            oo_.FilteredVariables.(M_.endo_names{bayestopt_.mfys(pos_iter)})=...
                squeeze(aK(1,bayestopt_.mf(pos_iter),2:end-(options_.nk-1)))...
                +trend_constant_observables_filtered.filter_ahead_1(pos_iter,:)';
            for filter_iter=1:length(options_.filter_step_ahead)
                filter_step=options_.filter_step_ahead(filter_iter);
                oo_.FilteredVariablesKStepAhead(filter_iter,find(i_endo_in_dr_matrices==bayestopt_.mf(pos_iter)),1+filter_step:end-(max(options_.filter_step_ahead)-filter_step)) = ...
                    squeeze(aK(filter_step,bayestopt_.mf(pos_iter),1+filter_step:end-(max(options_.filter_step_ahead)-filter_step)))...
                    +trend_constant_observables_filtered.(['filter_ahead_' num2str(filter_step)])(pos_iter,:)';
            end
        end
        %updated variables are E_t(y_t) so no trend shift is required
        oo_.UpdatedVariables.(M_.endo_names{bayestopt_.mfys(pos_iter)})=...
            updated_variables(bayestopt_.mf(pos_iter),:)'+trend_constant_observables(pos_iter,:)';
    end
end

%% resort fields that are in decision rule order to declaration order
if ~isempty(options_.nk) && options_.nk ~= 0
    positions_in_declaration_order=oo_.dr.order_var(bayestopt_.smoother_var_list(bayestopt_.smoother_saved_var_list));
    if ~(options_.selected_variables_only && ~(options_.forecast > 0)) %happens only when selected_variables_only is not used
        oo_.FilteredVariablesKStepAhead(:,positions_in_declaration_order,:)=oo_.FilteredVariablesKStepAhead;
        if ~isempty(PK) && options_.filter_covariance %get K-step ahead variances
            oo_.FilteredVariablesKStepAheadVariances(:,positions_in_declaration_order,positions_in_declaration_order,:)=oo_.FilteredVariablesKStepAheadVariances;
        end
        if ~isempty(decomp)
            oo_.FilteredVariablesShockDecomposition(:,positions_in_declaration_order,:,:)=oo_.FilteredVariablesShockDecomposition;
        end
    else
        positions_in_declaration_order=oo_.dr.order_var(bayestopt_.smoother_var_list(bayestopt_.smoother_saved_var_list));
        [~,sorted_index_declaration_order]=sort(positions_in_declaration_order);
        oo_.FilteredVariablesKStepAhead(:,sorted_index_declaration_order,:)=oo_.FilteredVariablesKStepAhead;
        if ~isempty(PK) && options_.filter_covariance %get K-step ahead variances
            oo_.FilteredVariablesKStepAheadVariances(:,sorted_index_declaration_order,sorted_index_declaration_order,:)=oo_.FilteredVariablesKStepAheadVariances;
        end
        if ~isempty(decomp)
            oo_.FilteredVariablesShockDecomposition(:,sorted_index_declaration_order,:,:)=oo_.FilteredVariablesShockDecomposition;
        end
    end
end

if options_.filter_covariance
    oo_.Smoother.Variance(oo_.dr.order_var,oo_.dr.order_var,:)=oo_.Smoother.Variance;
end
if ~isempty(Ptt) && options_.updated_covariance
    oo_.Smoother.UpdatedVariance(oo_.dr.order_var,oo_.dr.order_var,:) = oo_.Smoother.UpdatedVariance;
end
if options_.smoothed_state_uncertainty
    oo_.Smoother.State_uncertainty(oo_.dr.order_var,oo_.dr.order_var,:)=state_uncertainty;
    oo_.Smoother.Init_State_uncertainty(oo_.dr.order_var,oo_.dr.order_var)=state_uncertainty0;
end

%% get smoothed shocks
for exo_iter=1:M_.exo_nbr
    oo_.SmoothedShocks.(M_.exo_names{exo_iter})=innov(exo_iter,:)';
end

%% store per-shock smoothed variance V_eta(j,j,t) = Var(eta_j(t)|Y_{1:T})
if ~isempty(V_eta) && options_.smoothed_shock_variance
    for exo_iter=1:M_.exo_nbr
        oo_.SmoothedShockVariances.(M_.exo_names{exo_iter}) = squeeze(V_eta(exo_iter,exo_iter,:));
    end
end

%%  Smoothed measurement errors
if ~isequal(M_.H,0)
    %     measurement_error_indices=find(diag(M_.H)~=0);
    for meas_error_iter=1:length(options_.varobs)
        oo_.SmoothedMeasurementErrors.(options_.varobs{meas_error_iter})= measurement_error(meas_error_iter,:)';
    end
end

%% Store OccBin-specific results (when occbin_info is provided)
if nargin > 7 && ~isempty(occbin_info)
    % Always store error flag and regime history
    oo_.occbin.smoother.error_flag = occbin_info.error_flag;
    oo_.occbin.smoother.regime_history = occbin_info.regime_history;
    % Store remaining OccBin fields
    oo_.occbin.smoother.warning_flag = occbin_info.warning_flag;
    oo_.occbin.smoother.realtime_regime_history = occbin_info.realtime_regime_history;
    oo_.occbin.smoother.regime_history_iter = occbin_info.regime_history_iter;
    % Time-varying state space matrices (consumed by occbin.shock_decomposition)
    oo_.occbin.smoother.T0 = occbin_info.TT;
    oo_.occbin.smoother.R0 = occbin_info.RR;
    oo_.occbin.smoother.C0 = occbin_info.CC;
    % Raw smoothed states/shocks in DR order (consumed by occbin.shock_decomposition, occbin.irf)
    oo_.occbin.smoother.alphahat = atT;
    oo_.occbin.smoother.etahat = innov;
    oo_.occbin.smoother.epsilonhat = measurement_error;
    oo_.occbin.smoother.ahat = updated_variables;
    oo_.occbin.smoother.SteadyState = ys;
    oo_.occbin.smoother.trend_coeff = trend_coeff;
    oo_.occbin.smoother.aK = aK;
    oo_.occbin.smoother.is_converged = occbin_info.is_converged;
    oo_.occbin.smoother.is_periodic = occbin_info.is_periodic;
    % Simulation paths from occbin.solver
    if isfield(occbin_info,'simul')
        oo_.occbin.smoother.simul = occbin_info.simul;
    end
    % Linear smoother results
    if ~isempty(occbin_info.linear)
        linear_smoother_tmp = store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,occbin_info.linear);
        for jf = {'SmoothedVariables', 'UpdatedVariables', 'SmoothedShocks'}
            oo_.occbin.linear_smoother.(jf{1}) = linear_smoother_tmp.(jf{1});
        end
        oo_.occbin.linear_smoother.alphahat = occbin_info.linear.smoothed.alphahat;
        oo_.occbin.linear_smoother.etahat = occbin_info.linear.smoothed.etahat;
        oo_.occbin.linear_smoother.epsilonhat = occbin_info.linear.smoothed.epsilonhat;
        oo_.occbin.linear_smoother.ahat = occbin_info.linear.filtered.ahat;
        oo_.occbin.linear_smoother.SteadyState = occbin_info.linear.model.SteadyState;
        oo_.occbin.linear_smoother.trend_coeff = occbin_info.linear.model.trend_coeff;
        oo_.occbin.linear_smoother.aK = occbin_info.linear.filtered.aK;
        oo_.occbin.linear_smoother.T0 = occbin_info.linear.model.T;
        oo_.occbin.linear_smoother.R0 = occbin_info.linear.model.R;
        oo_.occbin.linear_smoother.decomp = occbin_info.linear.filtered.decomp;
        oo_.occbin.linear_smoother.alphahat0 = occbin_info.linear.smoothed.alphahat0;
        oo_.occbin.linear_smoother.state_uncertainty0 = occbin_info.linear.smoothed.V_0gT;
    end
end
