function y = esefit(Y, Phi, omega, sige2, gamma, Xbar, sigtau2, aJ, bJ, sigmasJ, rlb, lenJ)
% fit ese model given parameters 
% Output: fitted factor X_(t|t)

T = size(Y,1);
N = 3;
% lenJ = length(mats);
% pars_scale = [ones(9,1);  0.1;0.1;0.1;  ones(6,1);  10; 0.1;   1;1;0.1;0.1;  0.1];
% parameters = parameters./pars_scale;
% parameterization
% Phi = parameters(1:9); Phi = reshape(Phi,3,3); 
% if any(abs(eig(Phi))>1)
%     y = nan;
%     return
% end
% dlamQ = -exp(parameters(10:12))./(1 + exp(parameters(10:12)));
% lamQ = 1 + cumsum(dlamQ) ;
% if any(lamQ<0) || any(isnan(lamQ))
%     y = 1e8;
%     return
% end
% sigma = [parameters(13) 0 0 ; parameters(14) parameters(16) 0; parameters(15) parameters(17) parameters(18)]; 
% omega = sigma*sigma'; 
% kinfQ = parameters(19); 
% sige2 = exp(parameters(20));
% if sige2 > 10 || sige2 < 0.001
%     y = 1e9 ;
%     return
% end
% a = parameters(21:22);
% p = parameters(23:24);
% gamma = [1-sum(a); a];
% Xbar = [-sum(p); p];
% sigtau2 = exp(parameters(25));
% if sigtau2 > 10 || sigtau2<0.01
%     y = 1e9;
%     return
% end
% Sigma_tilde = omega - (gamma * gamma')*sigtau2 ;
% if invalidCovMat(Sigma_tilde)
%     y = nan ;
%     return
% end
fwdr = Y' ;
% parameterization
R = sige2 * eye(lenJ) ;
mu_Z = [0; (eye(N)-Phi)*Xbar] ;
Phi_Z = [1 zeros(1,N); (eye(N)-Phi)*gamma Phi] ;
Q = [sigtau2 sigtau2*gamma'; sigtau2*gamma omega] ;
if invalidCovMat(Q)
    y = nan ;
    return
end
% Initialization
X1 = NaN(N+1,T+1); % E_t-1[X_t]; t = 1,...,T+1
X2 = NaN(N+1,T+1); % E_t-1[X_t-1]; t = 1,...,T+1; only meanningful after t>=2
V1 = NaN((N+1)*(T+1),N+1); % COV_t-1[X_t]; t = 1,...,T+1
V2 = NaN((N+1)*(T+1),N+1); % COV_t-1[X_t-1]; t = 1,...,T+1; only meanningful after t>=2
% loglikvec = NaN(1,T); % lf(t); t = 1,...,T
% X1(:,1) = zeros(4,1); 
% V1(1:N+1,:) = 100*eye(N+1) ;
X1(:,1) = [5; 5; 1.5; 0.5];
V1(1:N+1,:) = 1*eye(N+1);
% Compute aJ, bJ and sigmasJ
% loads = SRTSM_loadings2(kinfQ, diag(lamQ), omega, mats, dt) ;
% aJ = loads.an;
% bJ = loads.bn;
% sigmasJ = loads.SigmaQ;
% errvec = nan(lenJ,T);
% z2 = nan(lenJ,T);
% Filtering the shadow rate
for t = 1:1:T
    musJ = aJ + bJ*X1(2:N+1,t);  % Jx1
    z1_temp = (musJ-rlb)./sigmasJ;
    z2_temp = rlb + (musJ-rlb).*normcdf(z1_temp) + sigmasJ.*normpdf(z1_temp);
    h = normcdf(z1_temp);   % Jx1
    h = [h h h].*bJ;        % Jx3  
    H = [zeros(lenJ,1), h];    % Jx4
    err = fwdr(:,t) - z2_temp;    % jx1
    S = H*V1((t-1)*4+1:t*4,:)*H' + R;   % jxj
    invS = S\eye(size(S,1)) ;
    K = V1((t-1)*4+1:t*4,:)*H'*invS ;   % 4xj
    % Calculate the likelihood
%     loglikvec(t) = -1/2 * (lenJ*log(2*pi) + log(det(S)) + err'*invS*err) ;
    % Update
    X2(:,t+1) = X1(:,t) + K*err;
    V2(4*t+1:4*(t+1),:) = (eye(N+1)-K*H)*V1((t-1)*4+1:t*4,:);
    % Predict
    X1(:,t+1) = mu_Z + Phi_Z * X2(:,t+1);
    V1(4*t+1:4*(t+1),:) = Phi_Z * V2(4*t+1:4*(t+1),:) * Phi_Z' + Q;
%     errvec(:,t) = err;
%     z2(:,t) = z2_temp;
end


y = X2(:,end);


end