function [ytT, n_T] = smoother(x_,T_,S_)

%% Data 
Zobs = S_.data;

[nobs, T_sample] = size(Zobs);

%% Options
Opt_ = S_.Options;

penalty = 1e10; % Penalty for non-uniqueness of a stable solution

%% Break dates

Tz   = T_(1); % date break in z
Tsig = T_(2); % date break in volatility

%% Parameters

S_.param = param_func(x_);

%% Variable and dimensions

v_    = var_dim ; 

n_    = v_.dim.n_    ; 
ns    = v_.dim.ns    ; 

%% Initial guess for steady state

    S_.initvals = [0.834755321 % ln c -- 1
              -0.095030773 % ln invest -- 2
              -0.973976034 % ln lambda -- 3
               0.005224883 % ln r_F -- 4
               4.631985754 % ln k -- 5
              -0.020702261 % ln l -- 6
               1.024943397 % ln w -- 7
              -4.461247468 % ln r_k -- 8
               1.220529165 % ln y --9
               0           % nx (can be negative) -- 10
               0           % ca (can be negative) -- 11
               0.003924883 % r_star -- 12
               0           % b_F (can be negative) -- 13
               0.014022286 % ln r -- 14
               0.0         % b  (can be negative) -- 15
              -0.2         % ln g -- 16
               0.18        % tr (can be negative) -- 17
               3.1845];    % p_b -- 18

%% Structure 1 (z_ss high, gy_ss)

S_.options.gy = 1 ;  % sets gy_ss

st_1 = structural_matrices(S_);
%
[   S0_1, ...
    S1_1, ...
    S2_1, ...
    ~, ...
    ~, ...
    ~, ...
    ~, ...
    unique_1 ] = Smats(st_1.C_, st_1.GAM0, st_1.GAM1, st_1.PSI, st_1.PPI);

if unique_1 ~= 1
   out = penalty;
   return
end

if isempty(unique_1)
   out = penalty;
   return
end


% BP Solution matrices
C_1 = S0_1(1:n_,:);
Q_1 = S1_1(1:n_,1:n_);
G_1 = S2_1(1:n_,:);

% For numerical precision of discrete-time Lyapunov equation
tole = 1e-14;
C_1(abs(C_1) < tole) = 0;
Q_1(abs(Q_1) < tole) = 0;
G_1(abs(G_1) < tole) = 0;

% y_bar_1 = (eye(n_)-Q_1) \ C_1; 

%% Structure 2 (z_ss low, gy_ss low , g_ss at structure 1)

S_.param.z_ss = S_.param.z_ss + S_.param.Delta_z ; 
S_.param.tau_w_ss = S_.param.tau_w_ss + S_.param.Delta_tauw;
S_.param.tau_K_ss = S_.param.tau_K_ss + S_.param.Delta_tauK;


S_.param.g_ss = st_1.y_steady(16);

S_.initvals = [0.834755321 % ln c -- 1
              -0.095030773 % ln invest -- 2
              -0.973976034 % ln lambda -- 3
               0.005224883 % ln r_F -- 4
               4.631985754 % ln k -- 5
              -0.020702261 % ln l -- 6
               1.024943397 % ln w -- 7
              -4.461247468 % ln r_k -- 8
               1.220529165 % ln y --9
               0           % nx (can be negative) -- 10
               0           % ca (can be negative) -- 11
               0.003924883 % r_star -- 12
               0           % b_F (can be negative) -- 13
               0.014022286 % ln r -- 14
               0.0         % b  (can be negative) -- 15
              -1.6         % ln gy -- 16
               0.18      % tr (can be negative) -- 17
              3.1845];    % p_b -- 18

S_.options.gy = 0 ;  % sets g_ss

st_2 = structural_matrices(S_);

[   S0_2, ...
    S1_2, ...
    S2_2, ...
    ~, ...
    ~, ...
    ~, ...
    ~, ...
    unique_2 ] = Smats(st_2.C_, st_2.GAM0, st_2.GAM1, st_2.PSI, st_2.PPI);

if unique_2 ~= 1
   out = penalty;
   return
end

if isempty(unique_2)
   out = penalty;
   return
end

% BP Solution matrices
C_2 = S0_2(1:n_,:);
Q_2 = S1_2(1:n_,1:n_);
G_2 = S2_2(1:n_,:);

% For numerical precision of discrete-time Lyapunov

C_2(abs(C_2) < tole) = 0;
Q_2(abs(Q_2) < tole) = 0;
G_2(abs(G_2) < tole) = 0;

%% Kulish-Pagan matrices over sample

Omegat = zeros(ns,ns,T_sample);

Omegat(:,:,1:Tsig-1) = repmat(S_.param.mu_sig*st_1.OMEGA,[1,1,Tsig-1]);

Omegat(:,:,Tsig:T_sample) = repmat(st_1.OMEGA,[1,1,T_sample-Tsig+1]);

Ct = zeros(n_,1, T_sample+1);
Pt = zeros(n_,n_,T_sample+1);
Dt = zeros(n_,ns,T_sample+1);

Ct(:,:,1:Tz-1) = repmat(C_1,[1,1,Tz-1]);
Pt(:,:,1:Tz-1) = repmat(Q_1,[1,1,Tz-1]);
Dt(:,:,1:Tz-1) = repmat(G_1,[1,1,Tz-1]);

Ct(:,:,Tz:T_sample) = repmat(C_2,[1,1,T_sample-Tz+1]);
Pt(:,:,Tz:T_sample) = repmat(Q_2,[1,1,T_sample-Tz+1]);
Dt(:,:,Tz:T_sample) = repmat(G_2,[1,1,T_sample-Tz+1]);

%% Kalman filter

% Build matrix H for measurement equation


H = zeros(nobs,n_);

j_ = 1;

H(j_,v_.dy_obs)      = 1; j_ = j_ + 1;
H(j_,v_.dc_obs)      = 1; j_ = j_ + 1;
H(j_,v_.nx_obs)      = 1; j_ = j_ + 1;
H(j_,v_.dw_obs)      = 1; j_ = j_ + 1;
H(j_,v_.r_obs)       = 1; j_ = j_ + 1;
H(j_,v_.r_star_obs)       = 1; j_ = j_ + 1;
H(j_,v_.gy_obs)      = 1; j_ = j_ + 1;
H(j_,v_.by_obs)      = 1; j_ = j_ + 1;
H(j_,v_.rev_c_y_obs) = 1; j_ = j_ + 1;
H(j_,v_.rev_w_y_obs) = 1; j_ = j_ + 1;
H(j_,v_.rev_K_y_obs) = 1; 

V = diag(0.00*diag(cov(Zobs'))); % Var-Cov matrix of measurement error

DOmega = Dt(:,:,1)*Omegat(:,:,1)*Dt(:,:,1)';

    Sigma = (eye((n_)^2)-(kron(Pt(:,:,1),Pt(:,:,1))))\DOmega(:); % eye(n_) because I remove the constant
    Sigma = reshape(Sigma,n_,n_);

% Construct tensor to store the l x l matrices P(t+1/t)

    St = zeros(n_,n_,T_sample);
    St(:,:,1) = Sigma;  % Initial condition
    
 
 for t=1:(T_sample-1)
     St(:,:,t+1) = Pt(:,:,t+1)*( St(:,:,t)- St(:,:,t)* H' /( H*St(:,:,t)*H' + V ) * H*St(:,:,t))* Pt(:,:,t+1)' + Dt(:,:,t+1)*Omegat(:,:,t+1)*Dt(:,:,t+1)';
 end

 Kt = zeros(n_,nobs,T_sample);

 for t=1:T_sample
     Kt(:,:,t) = (Pt(:,:,t+1) * St(:,:,t) * H') / (H * St(:,:,t) * H' + V);
 end

% Construct series for the forecasts y(t+1/t) 

 y_hat = zeros(n_,T_sample+1);
 y_hat(:,1) = ((eye(n_) - Pt(:,:,1)))\Ct(:,:,1); 

 
 for t=1:T_sample
     y_hat(:,t+1) = Ct(:,:,t+1) + Pt(:,:,t+1)* y_hat(:,t) + Kt(:,:,t)*( Zobs(:,t) - H * y_hat(:,t));
 end

%% Loglikelihood 
U = zeros(nobs,T_sample);
out = 0;

for t = 1:T_sample
    U(:,t) = Zobs(:,t) - H * y_hat(:,t);    
    out = out + (nobs/2)*log(2*pi)+(1/2)*log(det(H*St(:,:,t)*H'+V))+(1/2)*(U(:,t)'/(H*St(:,:,t)*H'+V)*U(:,t));
end


if  abs(imag(out)) > 0
    out = penalty;
end 


%% Kalman Smoother

% Construct y(t|t)
 ytt = zeros(n_,T_sample);
    
for t=1:T_sample
    ytt(:,t) = y_hat(:,t) + St(:,:,t)* H'*(inv(H*St(:,:,t)*H'+V))*U(:,t);
end
    
% Construct y(t|T)

ytT = zeros(n_,T_sample);
ytT(:,T_sample) = ytt(:,T_sample);
r_T = zeros(n_,T_sample+1);
n_T = zeros(ns,T_sample);

r_T(:,T_sample) = H' / (H*St(:,:,T_sample)*H'+ V ) * (Zobs(:,T_sample) - H*y_hat(:,T_sample));
ytT(:,T_sample) = y_hat(:,T_sample) + St(:,:,T_sample) * r_T(:,T_sample);
n_T(:,T_sample) = Omegat(:,:,T_sample) * Dt(:,:,T_sample)'* r_T(:,T_sample);

for t = T_sample-1:-1:1
    
    r_T(:,t) = H' / (H*St(:,:,t)*H'+ V ) * (Zobs(:,t) - H*y_hat(:,t)) + ...
            (eye(n_) - St(:,:,t)*H' / (H*St(:,:,t)*H'+ V)* H)' * Pt(:,:,t+1)' * r_T(:,t+1);
        
    ytT(:,t) = y_hat(:,t) + St(:,:,t) * r_T(:,t);
        
    n_T(:,t) = Omegat(:,:,t) * Dt(:,:,t)' * r_T(:,t);
    
end