subroutine loglike11(n,x,ldat,dat, lnl)
!
!  Log likelihood function for the following VARMA(1,1),
!
!        Y[t] = Phi(x) Y[t-1] + e[t] - Theta(x) e[t-1] 
!
!  where Y[t] is a vector of length k, e[t] ~ N(0,Sig), 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   = [T,k,Y(:)], where Y is T x k
!
!  The vector x is 
! 
!     x     = A11,A12,A21,A22, B11,B12,B21,B22, C11,C21,C22
!
!  where Mij is element (i,j) of matrix M and
!
!     Phi   = inverse[chol(I+AA')'] * A
!     Theta = inverse[chol(I+BB')'] * B
!     Sig   = CC'
!
!  Note that these transformations ensure stationarity and invertibility.
!

!  Ellen McGrattan, 4-23-06
!  Revised, ERM, 5-12-06

  implicit none
  integer, parameter                 :: tmax=5000,kmax=5,lwork=5000
  integer, intent(in)                :: n,ldat
  integer, dimension(lwork)          :: iwork,ipiv
  integer                            :: i,j,l,t,k,info,kk
  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,at,tem1,tem2,kf,akc,       &
                                        phi,thet,ahat,bhat,atem
  real, dimension(1,1)               :: uiu
  real, dimension(lwork)             :: work
  real                               :: det

  !
  !  Set up data matrix and parameters 
  !

  t       = int(dat(1))
  k       = int(dat(2))
  kk      = 2*k

  allocate(ahat(k,k))
  allocate(phi(k,k))
  allocate(bhat(k,k))
  allocate(thet(k,k))
  allocate(cs(k,k))
  allocate(a(kk,kk))
  allocate(atem(kk,kk))
  allocate(b(kk,k))
  allocate(c(k,kk))
  allocate(sige(k,k))
  allocate(bsigb(kk,kk))
  allocate(s(kk,kk))
  allocate(as(kk,kk))
  allocate(xhat(kk,1))
  allocate(u(k,1))
  allocate(ome(k,k))
  allocate(iome(k,k))
  allocate(at(kk,kk))
  allocate(tem1(k,kk))
  allocate(tem2(1,k))
  allocate(kf(kk,k))
  allocate(akc(kk,kk))

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

  !
  ! Set up the state space system
  !
  !     X[t+1] = [Phi(x) I] X[t] + [     I     ] e[t+1]
  !              [0      0]        [ -Theta(x) ]
  ! 
  !       Y[t] = [I      0] X[t]
  !
  ! 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.
  l    = 1
  do i = 1,k
    a(i,i+k) = 1.
    b(i,i)   = 1.
    c(i,i)   = 1.
    do j=1,k
      ahat(i,j) = x(l)
      l         = l+1
    enddo
  enddo
  do i=1,k
    do j=1,k
      bhat(i,j) = x(l)
      l         = l+1
    enddo
  enddo
  l    = 2*k*k+1
  do j=1,k
    do i=j,k
      cs(i,j) = x(l)
      l       = l+1
    enddo
  enddo         
  sige    = matmul(cs,transpose(cs))

  call stationary(k,ahat,sige, phi)
  call stationary(k,bhat,sige, thet)
  a(1:k,1:k)    = phi
  b(k+1:kk,1:k) = -thet

  !
  ! 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 Sig B'
  !

  lnl     = 0.
  xhat    = 0.
  tem1    = matmul(sige,transpose(b))
  bsigb   = matmul(b,tem1)
  s       = bsigb
  atem    = a
  at      = -transpose(a)

  call sb04qd(kk,kk,atem,kk,at,kk,s,kk,as,kk,iwork,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in loglike11.f90: problem with sb04qd'
    write(*,*) 'info = ',info
    return
  endif

  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 loglike11.f90: problem with third 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 loglike11.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(ahat)
  deallocate(phi)
  deallocate(bhat)
  deallocate(thet)
  deallocate(cs)
  deallocate(a)
  deallocate(atem)
  deallocate(b)
  deallocate(c)
  deallocate(sige)
  deallocate(bsigb)
  deallocate(s)
  deallocate(as)
  deallocate(xhat)
  deallocate(u)
  deallocate(ome)
  deallocate(iome)
  deallocate(at)
  deallocate(tem1)
  deallocate(tem2)
  deallocate(kf)
  deallocate(akc)

end subroutine loglike11