subroutine loglike(n,x,ldat,dat, lnl)
!
!  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, e[t] ~ N(0,Sige), and the 
!  sample size is for Y[t] is T, e.g., t=1,...T.
!
!  Inputs: 
!
!     n    = length of parameter vector (x)
!     x    = parameter vector
!     ldat = length of dat
!     dat  = [p,q,T,k,Y(:)]
!
!  The vector x is 
! 
!     x    = [vec([A1;
!                  A2;
!                  : 
!                  Ap;
!                  B1;
!                  :
!                  Bq]);
!             S(find(tril(ones(k)))]
!    
!  where vec(M) stacks columns of M into a vector and the argument
!  find(tril(ones(k))) takes elements, column by column, that are 
!  in the lower triangular part of the matrix.  The matrix S is
!  found by taking a Cholesky decomposition of Sige and is assumed
!  to be lower triangular (ie, S*S'=Sige where S is lower triangular.
!  The second input of loglike is dat.  The first four elements of 
!  dat 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. 
!           

!  Ellen McGrattan, 4-23-06

  implicit none
  integer, parameter                 :: tmax=5000,kmax=5,lwork=1000
  integer, intent(in)                :: n,ldat
  integer, dimension(lwork)          :: iwork,ipiv
  integer                            :: i,j,l,p,q,t,k,r,info,kr
  real, dimension(n), intent(in)     :: x
  real, dimension(ldat), intent(in)  :: dat
  real, intent(out)                  :: lnl
  real, dimension(tmax,kmax)         :: y
  real, dimension(:,:), allocatable  :: a,b,c,cs,sige,bsigb,s,as,xhat,u,    &
                                        ome,iome,atmp,tem1,tem2,kf,akc
  real, dimension(1,1)               :: vl,vr,uiu
  real, dimension(:), allocatable    :: er,ei
  real, dimension(lwork)             :: work
  real                               :: eiga,det,tem

  !
  !  Set up data matrix and parameters 
  !

  p       = dat(1)
  q       = dat(2)
  t       = dat(3)
  k       = dat(4)
  r       = max(p,q+1)
  kr      = k*r

  allocate(a(kr,kr))
  allocate(b(kr,k))
  allocate(c(k,kr))
  allocate(cs(k,k))
  allocate(sige(k,k))
  allocate(bsigb(kr,kr))
  allocate(s(kr,kr))
  allocate(as(kr,kr))
  allocate(xhat(kr,1))
  allocate(u(k,1))
  allocate(ome(k,k))
  allocate(iome(k,k))
  allocate(atmp(kr,kr))
  allocate(er(kr))
  allocate(ei(kr))
  allocate(tem1(k,kr))
  allocate(tem2(1,k))
  allocate(kf(kr,k))
  allocate(akc(kr,kr))

  if ((t>tmax).or.(k>kmax)) then
    write(*,*) 'Edit loglike and increase dimensions of y'
    return
  endif
  do j=1,k
    y(:,j) = dat((j-1)*t+5:j*t+4)
  enddo

  !
  ! 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]
  !
  a   = 0.
  b   = 0.
  c   = 0.
  cs  = 0.
  do j=1,k
    l         = (j-1)*(p+q)*k
    do i=1,k*p
      a(i,j)  = x(l+i)
    enddo
    b(j,j)    = 1.
    c(j,j)    = 1.
  enddo
  do j=1,(r-1)*k
    a(j,j+k)  = 1.
  enddo
  do j=1,k
    l         = j*p*k+(j-1)*q*k-k
    do i=k+1,(q+1)*k
      b(i,j)  = -x(l+i)
    enddo
  enddo
  l           = 1
  do j=1,k
    do i=j,k
      cs(i,j) = x((p+q)*k*k+l)
      l       = l+1
    enddo
  enddo         
  sige    = matmul(cs,transpose(cs))

  !
  ! Check that maximum eigenvalue of A is less than 1
  !

  atmp  = a
  call dgeev('N','N',kr,atmp,kr,er,ei,vl,1,vr,1,work,lwork,info)
  eiga  = sqrt(er(1)*er(1)+ei(1)*ei(1))  
  do i=2,k*r
    tem = sqrt(er(i)*er(i)+ei(i)*ei(i))
    if (tem>eiga) eiga = tem
  enddo
  if (eiga>=1) then
    lnl    = 1000.+1e+10*(eiga-1.)
    return
  endif
  !
  ! 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.
  xhat    = 0.
  tem1    = matmul(sige,transpose(b))
  bsigb   = matmul(b,tem1)
  s       = bsigb
  atmp    = -transpose(a)

  call sb04qd(kr,kr,a,kr,atmp,kr,s,kr,as,kr,iwork,work,lwork,info)

  do i=1,t
    tem2(1,1:k) = y(i,1:k)
    tem1        = matmul(c,s)
    ome         = matmul(tem1,transpose(c))
    u           = tem2-matmul(c,xhat)
    iome        = ome
    !
    ! Invert variance-covariance matrix ome
    !
    call dgetrf(k,k,iome,k,ipiv,info)
    if (info.ne.0) then
      write(*,*) 'ERROR in loglike.f90: problem with dgetrf'
      write(*,*) 'info = ',info
      return
    endif
    det         = 1.
    do j=1,k
      if (ipiv(j).ne.j) then
        det     = -det*iome(j,j)
      else
        det     = det*iome(j,j)
      endif
    enddo
    call dgetri(k,iome,k,ipiv,work,lwork,info)
    if (info.ne.0) then
      write(*,*) 'ERROR in loglike.f90: problem with dgetri'
      write(*,*) 'info = ',info
      return
    endif
    tem2        = matmul(transpose(u),iome) 
    uiu         = matmul(tem2,u)
    lnl         = lnl-.5*log(det)-.5*uiu(1,1)
    as          = matmul(a,s)
    kf          = matmul(as,transpose(c))
    kf          = matmul(kf,iome)
    akc         = a-matmul(kf,c)
    s           = matmul(akc,s)
    s           = matmul(s,transpose(akc))+bsigb
    xhat        = matmul(a,xhat)+matmul(kf,u)
  enddo
  lnl   = -lnl +0.91893853320467*t*k
  deallocate(a)
  deallocate(b)
  deallocate(c)
  deallocate(cs)
  deallocate(sige)
  deallocate(bsigb)
  deallocate(s)
  deallocate(as)
  deallocate(xhat)
  deallocate(u)
  deallocate(ome)
  deallocate(iome)
  deallocate(atmp)
  deallocate(er)
  deallocate(ei)
  deallocate(tem1)
  deallocate(tem2)
  deallocate(kf)
  deallocate(akc)

end subroutine loglike