subroutine loglikess(n,x,ldat,dat, lnl)
!
!  Log likelihood function for the following State Space
!
!      S[t+1] = A S[t] + B eps[t],  E eps[t]eps[t]' = I
!        Y[t] = C S[t] 
!
!        S[t] = [log k[t], s[t], s[t-1]]'
!        s[t] = P s[t-1] + Q eps[t]
!
!  and
!      A = [ gammak  , gamma, 0_{1 x k}]   B = [0_{1 x k}]
!          [0_{k x 1},   P  , 0_{k x k}]       [  Q      ]
!          [0_{k x 1},   I  , 0_{k x k}]       [0_{k x k}]
!
!  where Y[t], s[t], and eps[t] are vectors of length k, and 
!  the sample size for observables in 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,ind,Y(:)], where Y is T x k and
!               ind is a vector of indices indicating
!               which parameters are to be estimated 
!
!  The order of parameters that could be estimated are:
!
!              Order             Parameter
!              -----             ---------
!                1                gammak  
!              2 : k+1            gamma  
!            k+2 : (k+1)k+1       P(1,:),P(2,:),...,P(k,:)
!       (k+1)k+2 : (2k+1)k+1      Q(1,:),Q(2,:),...,Q(k,:)
!      (2k+1)k+2 : 2(2k+1)k+1     C(1,:),C(2,:),....C(k,:)
!

!  Ellen McGrattan, 6-3-06
!  Revised, ERM, 6-4-06

  implicit none
  integer, parameter                 :: tmax=5000,kmax=5,lwork=5000
  integer, intent(in)                :: n,ldat
  integer, dimension(lwork)          :: iwork,ipiv
  integer, dimension(:), allocatable :: ind
  integer                            :: i,j,l,t,k,info,ns,lind
  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,bb,s,as,xhat,u,ome,iome,at,    &  
                                        tem1,tem2,kf,akc,atem,imat,p,psta
  real, dimension(:), allocatable    :: par
  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))
  ns      = 2*k+1
  lind    = 2*k*ns+1

  allocate(a(ns,ns))
  allocate(atem(ns,ns))
  allocate(b(ns,k))
  allocate(c(k,ns))
  allocate(bb(ns,ns))
  allocate(s(ns,ns))
  allocate(as(ns,ns))
  allocate(xhat(ns,1))
  allocate(u(k,1))
  allocate(ome(k,k))
  allocate(iome(k,k))
  allocate(at(ns,ns))
  allocate(tem1(k,ns))
  allocate(tem2(1,k))
  allocate(kf(ns,k))
  allocate(akc(ns,ns))
  allocate(ind(lind))
  allocate(par(lind))
  allocate(imat(k,k))
  allocate(p(k,k))
  allocate(psta(k,k))

  ind     = int(dat(3:lind+2))

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

  par        = 0.
  par(1)     = .96
  par(2)     = 1.
  l          = (k+1)*k+2
  do j=1,k
    par(l)   = 1.
    l        = l+k+1
  enddo

  !
  ! Update default parameters in case of those to be estimated
  !
  l          = 1
  do j=1,lind
    if (ind(j).eq.1) then
      par(j) = x(l)
      l      = l+1
    endif
  enddo

  !
  ! Set up the state space system
  !
  !     S[t+1] = A S[t] + B e[t+1]
  !       Y[t] = C X[t]
  !
  a              = 0.
  b              = 0.
  c              = 0.
  a(1,1)         = par(1)/sqrt(1.+par(1)*par(1))
  a(1,2:k+1)     = par(2:k+1)
  l              = k+2
  imat           = 0.
  do i=1,k
    imat(i,i)    = 1.
    do j=1,k
      p(i,j)     = par(l)
      l          = l+1
    enddo
  enddo
  call stationary(k,p,imat, psta)
  a(2:k+1,2:k+1) = psta
  do i=2,k+1
    a(i+k,i)     = 1.
  enddo 
  l              = (k+1)*k+2
  do i=2,k+1
    do j=1,k
      b(i,j)     = par(l)
      l          = l+1
    enddo
  enddo
  l              = (2*k+1)*k+2
  do i=1,k
    do j=1,ns
      c(i,j)     = par(l)
      l          = l+1
    enddo
  enddo

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

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

  call sb04qd(ns,ns,atem,ns,at,ns,s,ns,as,ns,iwork,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in loglikess.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 loglikess.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 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))+bb
    xhat        = matmul(a,xhat)+matmul(kf,u)
  enddo
  lnl   = -lnl +0.91893853320467*t*k
  deallocate(a)
  deallocate(atem)
  deallocate(b)
  deallocate(c)
  deallocate(bb)
  deallocate(s)
  deallocate(as)
  deallocate(xhat)
  deallocate(u)
  deallocate(ome)
  deallocate(iome)
  deallocate(at)
  deallocate(tem1)
  deallocate(tem2)
  deallocate(kf)
  deallocate(akc)
  deallocate(ind)
  deallocate(par)
  deallocate(imat)
  deallocate(p)
  deallocate(psta)

end subroutine loglikess