function [l , g] = arimalik(phitheta,z,p,d)
% ARIMALIK(Z,PHI,THETA,D) (concentrated) -2 log likelihood of ARIMA model
% INPUTS:
% z = data vector; n by 1; missing values allowed
% phitheta = vector of AR coeffs. concatenated with MA coeffs. (p+q by 1)
% p = order of AR
% d = order of differencing
% RETURNS: l = n log (x-muhat)^T R^{-1} (x-muhat) + log |R|
%  where n = length(z)-d
%  x = d-fold differenced process
%  R = correlation matrix corresponding to phi, theta
%  muhat = generalized least squares estimate of mean mu
%  Except for constant terms, l is -2 log likelihood
% g = 2 by 1 vector of [ r1 ; r2] where r1 is 1 - minimum modulus of
%   roots of AR equation and r2 is for MA equation.
phi = phitheta(1:p);
[m n] = size(phitheta);
theta = phitheta((p+1):m);
x = z;
if d > 0
   for k=1:d
      x = diff(x);
   end
end
[n m] = size(x);
acf = armaacf(phi,theta,n-1);
R = toeplitz(acf);
T = chol(R);
dt = prod(diag(T))^2;
y = backsolve(T,[x ones(n,1)],1,1);
muhat = (y(:,1)'*y(:,2))/(y(:,2)'*y(:,2));
y1 = y(:,1)-muhat*y(:,2);
s2 = y1'*y1/n;
l = n*log(s2) + log(dt);
g = [ statcheck(phi) ; statcheck(theta) ];
