function y = llk_sdtsm(theta, Y, prior, rlb_in, mats, dt)
% log-likelihood of our Shadow rate DTSM

T = size(Y,2);
N = 3;
lenJ = length(mats);

Phi = reshape(theta(1:9),3,3);
% dlamQ = -exp(theta(10:12))./(1 + exp(theta(10:12)));
% lamQ = 1 + cumsum(dlamQ) ;
lamQ = theta(10:12);
% sigma = [theta(13) 0 0 ; theta(14) theta(16) 0; theta(15) theta(17) theta(18)]; 
Sigma_tilde = [theta(13) theta(14) theta(15); theta(14) theta(16) theta(17); theta(15) theta(17) theta(18)]; 
% omega = sigma*sigma'; 
kinfQ = theta(19);
% sigtau2 = exp(theta(25));
sigtau2 = theta(25);
a = theta(21:22);
p = theta(23:24);
gamma = [1-sum(a); a];
Xbar = [-sum(p); p];
% sige2 = exp(theta(20));
sige2 = theta(20);
% Sigma_tilde = omega - (gamma * gamma')*sigtau2 ;
omega = Sigma_tilde + (gamma * gamma')*sigtau2 ;

if any(abs(eig(Phi))>0.9999)
    y = 1e10;
    return
end
if any(lamQ<0) || any(isnan(lamQ)) || any(lamQ>0.99999)
    y = 1e10;
    return
end
if invalidCovMat(Sigma_tilde)
    y = 1e10 ;
    return
end

% 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)*(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,1) = 6;
errvec = nan(lenJ,T);
loads = SRTSM_loadings(kinfQ, diag(lamQ), omega, mats, dt) ;
aJ = loads.an;
bJ = loads.bn;
sigmasJ = loads.SigmaQ;

for t = 1:1:T
    rlb = rlb_in(t);
    musJ = aJ + bJ*X1(2:N+1,t);  % (K-1)x1
    z1_temp = (musJ-rlb)./sigmasJ;
    z2_temp = rlb + (musJ-rlb).*normcdf(z1_temp) + sigmasJ.*normpdf(z1_temp);
    h = normcdf(z1_temp) ;  
    h = [h h h].*bJ ;         % Kx3 
    H = [zeros(lenJ,1) h];    % Kx4
    err = Y(:,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;
end
llf = sum(loglikvec(3:end)); 
% lr_prior = log(diwish(Sigma_tilde, prior.kappa, prior.Psi))+ ...
%     log(mvnpdf([a;p], prior.ap_mean, prior.ap_cov))+...
%     log(dinvgamma(sigtau2, prior.alpha_eta/2, prior.delta_eta/2))+...
%     logprior_Phi(Phi,prior.Phi, prior.Omega)+...
%     log(dinvgamma(sige2, prior.alpha_e/2, prior.delta_e/2))+...
%     log(normpdf(kinfQ, prior.kmean, prior.kvar));

lr_prior = log(logprior_Phi(Sigma_tilde, prior.Psi, prior.PsiCov))+ ...
    log(mvnpdf([a;p], prior.ap_mean, prior.ap_cov))+...
    log(normpdf(log(sigtau2), prior.mean_sigtau2, prior.delta_eta))+...
    logprior_Phi(Phi, prior.Phi, prior.Omega)+...
    log(normpdf(log(sige2), prior.mean_sige2, prior.delta_e))+...
    log(normpdf(kinfQ, prior.kmean, prior.kvar));

y = -(llf+lr_prior) ;

end


