function y = llk_P(tau, X, Xbar, Phi, gamma, Omega, sigtau2, nobs)
% P-likelihood
% tau_{t+1} = tau_t + eta_{t+1},  Var(eta_{t+1}) = sigtau2
% Xtil_{t+1} = Phi * Xtil_t + util_{t+1},  Var(util_{t+1}) = Sigma_tilde

Xtilde = X - repmat(Xbar',nobs,1) - tau.*repmat(gamma',nobs,1) ; % Tx3
Sigma_tilde = Omega - sigtau2*(gamma*gamma');
e = Xtilde(2:end,:)' - Phi*Xtilde(1:end-1,:)' ; % 3x(T-1) VAR errors (Xtilde only)
e = e(:,2:end); % 3x(T-2)
llk_var = -.5*(nobs-1)*log(det(Sigma_tilde)) - .5*sum(e.*(Sigma_tilde\e),'all') ;
e2 = diff(tau); % T-1
e2 = e2(2:end); % T-2
llk_rw = -.5*(nobs-1)*log(sigtau2) - .5*sum(e2.^2/sigtau2) ;
y= llk_var + llk_rw ;

% Z = [tau,X];  % Tx4
% mu_Z = [0; (eye(3)-Phi)*Xbar];
% Phi_Z = [1, zeros(1,3); (eye(3)-Phi)*gamma, Phi];
% e = Z(2:end,:)' - mu_Z - Phi_Z*Z(1:end-1,:)'; % 3x(T-1)
% Omega_Z = [sigtau2 gamma'*sigtau2; gamma*sigtau2 Omega];
% if invalidCovMat(Omega_Z)
%     y = 1e-10;
%     return
% end
% y = -.5*(nobs-1)*log(det(Omega_Z)) - 0.5*sum(e.*(Omega_Z\e),'all');



