function val = jsz_llk(yields_o, W, K1Q_X, kinfQ, K0P_cP, K1P_cP, Sigma_cP, mats, dt, sigma_e, restr, ind_restr,Lam1_rank,ev_max,Lam0,Lam1)
%     ## Compute the likelihood for a Gaussian term structure.
%     ## Source "A New Perspective on Gaussian Dynamic Term Structure Models" by Joslin, Singleton and Zhu
%     ##
%     ## INPUTS:
%     ## yields.o   : (T+1)*J,  matrix of observed yields (first row are t=0 observations, which likelihood conditions on)
%     ## mats       : 1*J,      maturities in years
%     ## dt         : scalar,   length of period in years
%     ##
%     ## W          : N*J,      vector of portfolio weights to fit without error.
%     ## K1Q,X      : N*N,      normalized latent-model matrix (does not have to be diagonal, see form below)
%     ## kinfQ      : scalar,  when the model is stationary, the long run mean of the annualized short rate under Q is -kinfQ/K1(m1,m1)
%     ## Sigma.cP   : N*N,      positive definite matrix that is the covariance of innovations to cP
%     ##
%     ## OPTIONAL INPUTS -- concentrated out if not supplied:
%     ## K0P.cP     : N*1       intercept in VAR for cP
%     ## K1P.cP     : N*N       mean reversion matrix in VAR for cP
%     ## sigma.e    : scalar    standard error of yield observation errors
%     ## restr      : scalar    what type of restrictions on market prices of risk or VAR params
%     ##                        0 -- no restrictions
%     ##                        1 -- zero restrictions indicated by vector ind.restr
%     ##                        2 -- reduced rank restriction (JSZ)
%     ##                        3 -- eigenvalue restr. lam_1^Q = lam_1^P (JPS)
%     ##                        4 -- restrictions on VAR as in Duffee-Forecasting
%     ##                        planned -- zero restrictions and eigenvalue restriction
%     ## ind.restr  : N^2*1     for case 1 -- vector with indicators for zero restrictions on Lam1
%     ## Lam1.rank  : scalar    for case 2 -- rank of Lamda.1
%     ## ev.max     : scalar    for cases 3 and 4 -- largest eigenvalue of Phi -- NA means Q-ev
%     ##
%     ## OUTPUT:
%     ## llk        : T*1       time series of -log likelihoods (includes 2-pi constants)
%     ## AcP        : 1*J       yt = AcP' + BcP'*Xt  (yt is J*1 vector)
%     ## BcP        : N*J       AcP, BcP satisfy internal consistency condition that AcP*W' = 0, BcP*W' = I_N
%     ## AX         : 1*J       yt = AX' + BX'*Xt
%     ## BX         : N*J       Xt is the 'jordan-normalized' latent state
%     ## ...
%     ##
%     ## The model takes the form:
%     ##   r(t) = rho0.cP + rho1.cP'*cPt
%     ##        = 1'*Xt  (Xt is the 'jordan-normalized' state
%     ##        = 1 period discount rate (annualized)
%     ##
%     ## Under Q:
%     ##   X(t+1) - X(t)   = K0Q.X  + K1Q.X*X(t)  + eps_X(t+1),   cov(eps_X(t+1)) = Sigma.X
%     ##   cP(t+1) - cP(t) = K0Q.cP + K1Q.cP*X(t) + eps_cP(t+1),  cov(eps_cP(t+1)) = Sigma.cP
%     ##   where Sigma.X is chosen to match Sigma.cP
%     ## and K0Q_X(m1) = kinfQ where m1 is the multiplicity of the highest eigenvalue (typically 1)
% 
%     ## Under P:
%     ##   cP(t+1) - cP(t) = K0P.cP + K1P.cP*X(t) + eps_cP(t+1),  cov(eps_cP(t+1)) = Sigma.cP
%     ##
%     ## Model yields are given by:
%     ##   yt^m = AcP' + BcP'*cPt  (J*1)
%     ## And observed yields are given by:
%     ##  yt^o = yt^m + epsilon.e(t)
%     ## where V*epsilon.e~N(0,sigma.e^2 I_(J-N))
%     ## and V is an (J-N)*J matrix which projects onto the span orthogonal to the
%     ## row span of W.  This means errors are orthogonal to cPt and cPt^o = cPt^m.
%     ##

% setup
T = size(yields_o,1)-1  ;
J = size(yields_o,2)  ;
N = size(W,1)   ;
cP = yields_o * W'  ;  % (T+1)*N

%     ## COMPUTE THE Q-LIKELIHOOD:
%     ## First find the loadings for the model:
%     ## yt = AcP' + BcP'*cPt, AcP is 1*J, BcP is N*J

if isempty(kinfQ) 
%      ## concentrate out kinfQ
%      ## AcP = alpha0_cP*kinf + alpha1_cP    
    rho0_cP = 0 ;
    % AcP0, AX0 will be the loadings with rho0_cP = 0, which won't be correct
    loads = jsz_loadings_rho0cP(W, K1Q_X, rho0_cP, Sigma_cP, mats, dt)  ;
    BcP = loads.BcP ; 
    BX = loads.BX  ;
    alpha0_X = loads.alpha0_X  ; 
    alpha1_X = loads.alpha1_X  ;
    alpha0_cP = loads.alpha0_cP  ; 
    alpha1_cP = loads.alpha1_cP  ;
    K1Q_X = loads.K1Q_X ;
    m1 = loads.m1;
    % back out kinfQ that fits average yields
    V = transpose(null(W)) ;
    kinfQ = transpose(mean(yields_o(2:(T+1),:),1)' - alpha1_cP' - BcP'* (mean(cP(2:(T+1),:),1))') * (V'*V*alpha0_cP') /(alpha0_cP* (V)' *V*alpha0_cP') ;
    
    % get correct loadings
    AX = alpha0_X*kinfQ + alpha1_X ;
    AcP = alpha0_cP*kinfQ + alpha1_cP ;
    
    % get these to return to caller (not used in jsz.llk)
    K0Q_X = zeros(N,1) ;
    K0Q_X(m1) = kinfQ;
    params = jsz_rotation(W, K1Q_X, K0Q_X, dt, BX, AX) ;
    K0Q_cP = params.K0Q_cP ;
    K1Q_cP = params.K1Q_cP  ;
    rho0_cP = params.rho0_cP  ; 
    rho1_cP = params.rho1_cP  ;
    
else
    loads = jsz_loadings(WR, K1Q_X, kinfQ, Omega_cP, mats, dt) ;
    BcP = loads.BcP   ; 
    AcP = loads.AcP   ; 
    K0Q_cP=loads.K0Q_cP  ; 
    K1Q_cP=loads.K1Q_cP  ;
    rho0_cP = loads.rho0_cP   ; 
    rho1_cP = loads.rho1_cP   ;
    K0Q_X = loads.K0Q_X ; 
    K1Q_X=LOADS.K1Q_X  ;
    AX = loads.AX  ; 
    BX = loads.BX  ;
end

yields_m = ones(T+1,1) * AcP + cP * BcP ; % (T+1)*J, model-implied yields
yield_errors = yields_o(2:(T+1),:) - yields_m(2:(T+1),:) ;   % T*J
square_orthogonal_yield_errors = yield_errors.^2 ;     % T*J  but N-dimensional projection onto W is always 0, so effectively (J-N) dimensional

% Compute optimal sigma.e if it is not supplied
if isempty(sigma_e)
    sigma_e = sqrt(sum(square_orthogonal_yield_errors(:))/(T*(J-N))) ;
end

term1 = 0.5* transpose(sum(square_orthogonal_yield_errors,2)/sigma_e^2) ;
term2 = (J-N)*0.5*log(2*pi) ;
term3 = 0.5*(J-N)*log(sigma_e^2) ;

llkQ = term1 + term2 + term3 ;  %1*T


% COMPUTE THE P-LIKELIHOOD:
if isempty(K0P_cP) || isempty(K1P_cP)
    % VAR parameters not supplied
    if restr == 0
        % unrestricted MPR -- run OLS for unconstrained VAR
        % Run OLS to obtain maximum likelihood estimates of K0P, K1P
        mdl = varm(3,1) ;
        ar_ols = estimate(mdl,cP) ;
        K0P_cP = ar_ols.Constant  ;
        K1P_cP = cell2mat(ar_ols.AR) - eye(N)  ;
    elseif  restr == 1
        % zero restrictions on MPR
        if isempty(Lam0) || isempty(Lam1)
            % Lam0/Lam1 not provided -- concentrated out of LLK
            if size(ind_restr,1) ~= N + N^2
                disp('jsz.llk: number of indicators not equal to elements of [Lam0, Lam1]');
                return
            end
            if sum(ind_restr) > 0
                mat = getLambdaRVAR(ind_restr, cP, K0Q_cP, K1Q_cP, inv(Sigma_cP)) ;
                Lam0 = mat.Lam0 ;
                Lam1 = mat.Lam1 ;
            else
                Lam0 = zeros(N, 1) ;
                Lam1 = zeros(N, N) ;
            end
        else
            if isempty(ind_restr)
                warning('jsz.llk: restr=1, Lam0/Lam1 and ind.restr provided -> ignoring ind.restr')
            end
        end
        K0P_cP = K0Q_cP + Lam0 ;
        K1P_cP = K1Q_cP + Lam1 ;
    elseif restr ==2
        % reduced rank regression ---- JSZ equation (29)
        if isempty(Lam0) || isempty(Lam1)
            % Lam0,Lam1 concentrated out, obtained via reduced rank regression
            % dep. var. is P_{t+1} - E_t^Q P_t
            Ydat = cP(2:T,:) - ( (K0Q_cP * ones(1,T-1))' + cP(1:T-1,:)*(K1Q_cP + eye(N))' ) ;
            % ind. var. is P_t
            Xdat = cP(1:(T-1),:) ;
            Ybar = mean(Ydat,1)  ;
            Xbar = mean(Xdat,1)  ;
            Ydem = Ydat - Ybar ;
            Xdem = Xdat - Xbar ;
            % first perform unrestricted regression
            beta_ols = Ydem' * Xdem / ( Xdem'*Xdem) ;
            if isempty(Lam1.rank) 
                disp('need to provide Lam1.rank') ;
                return ;
            end
            if Lam1_rank<0 || Lam1_rank>N 
                disp('invalid value of Lam1.rank') ;
                return ;
            end
            if lam1_rank == N
                Lam1 = beta_ols;
            else
                % use SVD
                P = (chol(Xdem'*Xdem))' ;
                L = (chol(Sigma_cP))' ;
                [U,S,V] = svd( P*(beta_ols)'/ L) ;
                D  = diag(S) ;
                % ## set singular values equal to zero to get matrix with rank r
                D((Lam1_rank+1):N) = 0  ;
                Lam1 = ( P\  U * diag(D) * V' * L)'  ;
            end
            Lam0 = (Ybar - Xbar*Lam1')' ;   %%  this is the MLE
        else
            % Lam0, Lam1 provided -- check rank restriction
            if rank(Lam1) ~= Lam1_rank
                disp("Lam1 does not satisfy rank restriction")
                return
            end
        end
        
        K0P_cP = K0Q_cP + Lam0  ;
        K1P_cP = K1Q_cP + Lam1  ;
    
    elseif restr == 3
        %  largest P-ev restricted to largest Q eigenvalue or value of ev
        %  dep. var. is P_{t+1} - lamQ.max*I_N*P_t
        if isempty(ev_max)
            ev_max = K1Q_X(1,1)+1 ;
        end
        Ydat = cP(2:T,:) - cP(1:(T-1),:) *( ev_max*eye(N) )' ;
        %  ind. var. is P_t
        Xdat = cP(1:(T-1),:) ;
        Ybar = mean(Ydat,1)  ;
        Xbar = mean(Xdat,1)  ;
        Ydem = Ydat - Ybar ;
        Xdem = Xdat - Xbar ;
        % first perform unrestricted regression
        beta_ols = Ydem' * Xdem / ( Xdem'*Xdem) ;
        % use SVD
        P = (chol(Xdem'*Xdem))' ;
        L = (chol(Omega_hat))' ;
        [U,S,V] = svd( P*(beta_ols)'/ L) ;
        D  = diag(S) ;
        % set singular value equal to zero to get matrix with rank N-1
        D(N) = 0  ;
        beta_restr = ( P\  U * diag(D) * V' * L)'  ;
        Phi_restr = beta_restr + ev_max * eye(N) ;
        K1P_cP = Phi_restr - eye(N) ;
        % intercept
        K0P_cP = (Ybar - Xbar*beta_restr')' ;
    elseif restr == 4
         % restricted VAR like in Duffee
         % PC1_t+1 = ev.max*PC1_t + Phi_12 PC2_t + Phi_13 PC3_t + error
         % Phi_21 = 0, Phi_31 = 0
         if isempty(ev_max)
             ev_max = K1Q_X(1,1)+1 ;
         end
         % Run OLS to obtain maximum likelihood estimates of K0P, K1P, given restrictions
         % VAR for PC2..PCN
         mdl = varm(2,1) ;
         var23 = estimate(mdl,cP(:,2:N)) ;
         % forecast changes in PC1
         ydat = cP(2:T,1)-ev_max*cP(1:(T-1),1);  % changes/PC1(t)-ev*PC1(t-1) in first PC
         xdat = cP(1:(T-1),2:N) ;
         lm1 = fitlm(xdat,ydat,'Intercept',false) ;
         Phi_hat = [[ev_max,zeros(N-1,1)] , [(lm1.Coefficients.Estimate)';cell2mat(var23.AR)]] ;
         mu_hat = [ 0 ; var23.Constant ]  ;
         K1P_cP = Phi_hat - eye(N) ;
         K0P_cP = mu_hat ;
    else
        disp("not yet implemented");
        return
    end
else
    if restr ~= 0
        disp("can't provide VAR params AND have restrictions on MPR!")
        return
    end
end

innovations = transpose(cP(2:(T+1),:)) - (K0P_cP * ones(1,T) + (K1P_cP+eye(N)) * transpose(cP(1:T,:)))  ; % N*T

llkP = .5*N*log(2*pi) + .5*log(det(Sigma_cP)) + .5*sum((innovations .* (Sigma_cP \innovations)),1) ;   % 1*T

val.llk=(llkQ + llkP)'  ; 
val.AcP = AcP  ; 
val.BcP = BcP  ;
val.AX = AX    ; 
val.BX = BX    ;
val.K0P_cP = K0P_cP ;
val.K1P_cP = K1P_cP ;
val.sigma_e = sigma_e ;
val.K0Q_cP = K0Q_cP ;
val.K1Q_cP = K1Q_cP ;
val.rho0_cP = rho0_cP ; 
val.rho1_cP = rho1_cP ;
val.K0Q_X = K0Q_X  ;
val.K1Q_X = K1Q_X  ;
val.cP = cP  ; 
val.llkP = llkP  ; 
val.llkQ = llkQ  ;

end





