function likel = jpslikel(param,m,W,mats,Sigma_cP_OLS);
% Value of likelihood function for DTSM with unspanned macro risks
% Bauer, Rudebusch, Wu (2014, AER) Comment on ``Term Premia and
% Inflation Uncertainty: Empirical Evidence from an International Panel Dataset'' 
%
% This code was generously provided to us by Jonathan Wright.

lb=[-0.6;-0.6;-0.6;0;-inf*ones(15,1)];
ub=[-.0001;-.0001;-.0001;0.4;inf*ones(15,1)];
A = [[-1 1 0; 0 -1 1] zeros(2,16)];
b = [0;0];
if any(param>ub)||any(param<lb)||any(A*param>b)
    likel=100000;
else
    nf=size(Sigma_cP_OLS,1);
    [T,J] = size(m(2:end,:));
    N = size(W,1);
    ch=param(5:end);
    for j=1:nf-1; ch=[ch(1:nf*j);zeros(j,1);ch((nf*j)+1:end)]; end;
    ch=reshape(ch,nf,nf);
    Sigma_cP=ch*ch';
    K1Q_X=diag(param(1:3)); rinfQ=param(4);
    cP = m*W'; % (T+1)*N, cP stands for math caligraphic P.
    [BcP,AcP,AX,BX] = jszLoadings(W,K1Q_X,rinfQ,Sigma_cP(1:3,1:3),mats,1/4);
    [K0Q_cP,K1Q_cP,rho0_cP,rho1_cP] = jszRotation(W,K1Q_X,rinfQ,1/4,[],[],BX,AX);
    yields_m =  ones(T+1,1)*AcP + cP*BcP; % (T+1)*J, model-implied yields
    yield_errors = m(2:end,:) - yields_m(2:end,:); % 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
    sigma_e = sqrt(sum(square_orthogonal_yield_errors(:))/(T*(J-N)));
    llkQ = 0.5*sum(sum(square_orthogonal_yield_errors.')/sigma_e^2 + (J-N)*log(2*pi) + (J-N)*log(sigma_e^2));
    llkP=(0.5*N*T*log(2*pi))+(0.5*T*log(det(Sigma_cP)))+(0.5*T*sum(diag(pinv(Sigma_cP)*Sigma_cP_OLS)));
    likel=llkQ+llkP;
end