function [lnl,U] = loglike(Gamma,data);
%LOGLIKE  Log likelihood function for the following VARMA(p,q),
%
%            Y[t] = A1 Y[t-1] + A2 Y[t-2]+ ... Ap Y[t-p] 
%                               + e[t] - B1 e[t-1] - ... Bq e[t-q],
%
%         where Y[t] is a vector of length k and e[t] ~ N(0,Sige).
%         lnl = loglike(Gamma,data) has an input vector Gamma which 
%         contains the parameters to be estimated, i.e.,
%
%            Gamma = [vec([A1;
%                          A2;
%                          :
%                          Ap;
%                          B1;
%                          :
%                          Bq]);
%                     S(find(tril(ones(k)))]
%    
%         where the vec term contains the elements of the A's and 
%         B's and S=chol(Sige)'.  Note that only the nonzero elements
%         of S should be included.  The second input of loglike is 
%         data.  The first four elements of data are p, q, T,
%         and k.  Elements 5 through T*k+4 are the columns of Y,
%         vec(Y), where Y is T x k. Element T*k+5 of data is 1 if
%         the steady state is to be used to initialize the variance-
%         covariance of the hidden state vector of the state-space
%         system associated with the VARMA(p,q).  By default, a
%         sequence of variance-covariance matrices are computed 
%         starting with the variance-covariance derived by assuming 
%         there are no observations.
%           

%         Ellen McGrattan, 4-23-06

%
% Set up data matrix and parameters 
%

p       = data(1);
q       = data(2);
T       = data(3);
k       = data(4);
Y       = zeros(T,k);
Y(:)    = data(5:T*k+4);
vcss    = 0;
if length(data)>T*k+4;
  vcss  = data(T*k+5);
end;

AB      = zeros((p+q)*k,k);
AB(:)   = Gamma(1:(p+q)*k*k);
Sige    = tril(ones(k));
i       = find(Sige);
Sige(i) = Gamma(k*k*(p+q)+1:length(Gamma));
Sige    = Sige*Sige';

%
% Set up the state space system
%
%     X[t+1] = [A1   I  0  0 ...   0] X[t] + [  I  ] e[t+1]
%              [A2   0  I  0 ...   0]        [ -B1 ]
%              [:    :  :  :       :]        [ :   ]
%              [Ar-1 0  0  0       I]        [     ]
%              [Ar   0  0  0       0]        [-Br-1]
%
% 
%       Y[t] = [I    0 ...         0] X[t]
%
% where r = max(p,q+1), or more succinctly, 
%
%     X[t+1] = A X[t] + B e[t+1]
%       Y[t] = C X[t]
%

r       = max(p,q+1);
A       = [[AB(1:p*k,:);zeros((r-p)*k,k)],[eye((r-1)*k);zeros(k,(r-1)*k)]];
B       = [eye(k);-AB(p*k+1:(p+q)*k,:);zeros((r-1-q)*k,k)];
C       = [eye(k),zeros(k,(r-1)*k)];
lam     = max(abs(eig(A)));
if lam>=1;
 lnl    = 1e+10*(lam-1);
 return
end;
%
% Use the Kalman Filter to derive innovations u[t]:
%
%        
%     Xhat[t+1] = A Xhat[t] + K[t] u[t]
%          Y[t] = C Xhat[t] + u[t]
%     
% where 
%
%     K[t]   = A S[t] C' *inv(C S[t] C')
%     S[t+1] = (A-K[t]C') S[t] (A-K[t]C')' + B Sige B'
%

lnl     = 0;
innov   = zeros(T,k);
xhat    = zeros(r*k,1);
if vcss;
  [K,S]   = kfilter(A,C,B*Sige*B',1e-15*eye(k));
  Ome     = C*S*C';
  for t=1:T;
    u     = Y(t,:)'-C*xhat;
    xhat  = A*xhat+K*u;
    U(t,:)= u';
  end;
else
  S       = dlyap(A,B*Sige*B');
  for t=1:T;
    Ome   = C*S*C';
    u     = Y(t,:)'-C*xhat;
    lnl   = lnl -.5*k*log(2*pi)-.5*log(det(Ome))-.5*u'*inv(Ome)*u;
    K     = A*S*C'/Ome;    
    S     = (A-K*C)*S*(A-K*C)'+B*Sige*B';
    xhat  = A*xhat+K*u;
    U(t,:)= u';
  end;
end;

lnl     = -lnl;
