function y = startval_ese(theta, Y, rlb, mats, dt)
% log-likelihood of srtsm combined with falling stars
T = size(Y,1);
% J = mats/dt;
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];
theta = theta./pars_scale;
% parameterization
Phi = theta(1:9); Phi = reshape(Phi,3,3); 
if any(abs(eig(Phi))>1)
    y = 1e10;
    return
end
dlamQ = -exp(theta(10:12))./(1 + exp(theta(10:12)));
lamQ = 1 + cumsum(dlamQ) ;
if any(lamQ<0) || any(isnan(lamQ)) || any(lamQ>0.999)
    y = 1e8;
    return
end
sigma = [theta(13) 0 0 ; theta(14) theta(16) 0; theta(15) theta(17) theta(18)]; 
omega = sigma*sigma'; 
kinfQ = theta(19); 
sige2 = exp(theta(20));
if sige2 > 10 || sige2 < 0.001
    y = 1e9 ;
    return
end
a = theta(21:22);
p = theta(23:24);
gamma = [1-sum(a); a];
Xbar = [-sum(p); p];
sigtau2 = exp(theta(25));
if sigtau2 > 10 || sigtau2<0.01
    y = 1e9;
    return
end
Sigma_tilde = omega - (gamma * gamma')*sigtau2 ;
if invalidCovMat(Sigma_tilde)
    y = 1e10 ;
    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 = 1e6 ;
    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,N+1,T+1); % COV_t-1[X_t]; t = 1,...,T+1
V2 = NaN(N+1,N+1,T+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); 
X1(1,1) = 9;
V1(:,:,1) = 100*eye(N+1) ; 

% Compute aJ, bJ and sigmasJ
loads = SRTSM_loadings(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)*H' + R;       % jxj
    invS = S\eye(size(S,1));
    K = V1(:,:,t)*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(:,:,t+1) = (eye(N+1)-K*H)*V1(:,:,t);
    % Predict
    X1(:,t+1) = mu_Z + Phi_Z * X2(:,t+1);
    V1(:,:,t+1) = Phi_Z * V2(:,:,t+1) * Phi_Z' + Q;
    errvec(:,t) = err;
    z2(:,t) = z2_temp;
end
llf = sum(loglikvec(3:end)); 
y = -llf;

