subroutine loglikefiml(n,x,ldat,dat, lnl)
!
!  Log likelihood function for the following model
!
!      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], 1]'
!        s[t] = P s[t-1] + Q eps[t]
!
!  and
!      A = [ gammak  ,   gamma  , 0_{1 x k+1}]  B = [0_{1 x k}]
!          [0_{k x 1},     P    , 0_{k x k+1}]      [   Q     ]
!          [0_{k x 1},     I    , 0_{k x k+1}]      [0_{k x k}]
!          [   0     , 0_{1 x k}, 0_{1 x k},1]      [0_{1 x k}]
!
!  where k<=4 and gammak, gamma, P, Q, and C are functions of the 
!  parameters in x. The vectors Y[t], s[t], and eps[t] are 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,ishk,iobs,ind,icon,ub,lb,Y(:)], where
!               ishk = k x 1 vector listing shocks
!               iobs = k x 1 vector listing observables
!               ind  = 1 for parameters to be estimated
!                      0 otherwise
!               icon = 1 for parameters to be transformed
!                      0 otherwise
!               ub   = upper bound for parameter 
!               lb   = lower bound for parameter 
!               Y    = T x k matrix of observations
!
!  The order of parameters that could be estimated is:
!
!              Order             Parameter
!              -----             ---------
!                1                int. rate
!                2                delta
!                3                gn
!                4                psi
!                5                sigma
!                6                theta
!               7-10              sbar
!              11-26              P (by rows)
!              27-42              Q (by rows)
!

!  Ellen McGrattan, 6-10-06
!  Revised, ERM, 10-2-09

  implicit none
  integer, parameter                 :: tmax=10000,kmax=5,lwork=5000,lpar=42
  integer, intent(in)                :: n,ldat
  integer, dimension(lwork)          :: iwork,ipiv
  integer, dimension(lpar)           :: ind,icon
  integer, dimension(:), allocatable :: ishk,iobs,istate
  integer                            :: i,j,l,t,k,info,ns,ns1,inbadregion
  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,s0,u,ome,iome,at, &  
                                        vs,tem1,tem2,tem3,kf,akc,atem
  real, dimension(:), allocatable    :: er,ei
  real, dimension(1,1)               :: uiu,vl,vr
  real, dimension(lwork)             :: work
  real, dimension(lpar)              :: par,lb,ub
  real, dimension(6)                 :: gamma,ss
  real, dimension(6)                 :: psil,psiy,psix,psic
  real, dimension(4,4)               :: p,q,imat
  real, dimension(4,1)               :: sbar,p0
  real, dimension(6,10)              :: allc
  real, dimension(10)                :: cdlogprd,cqdlogl,clogl,clogxy,     &
                                        clogcy,cdlogy,alls0
  real                               :: det,alpha,eiga,tem,largenum


  largenum       = 100000.
  !
  ! Read in matrix dimensions from dat vector
  !

  t              = int(dat(1))
  k              = int(dat(2))
  ns             = 2*k+2
  ns1            = ns-1
  if ((t>tmax).or.(k>kmax)) then
    write(*,*) 'Edit loglikefiml and increase dimensions of y'
    return
  endif

  !
  ! Allocate space for matrices 
  !
  allocate(a(ns,ns))
  allocate(atem(ns1,ns1))
  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(s0(ns,1))
  allocate(u(k,1))
  allocate(ome(k,k))
  allocate(iome(k,k))
  allocate(at(ns1,ns1))
  allocate(vs(ns1,ns1))
  allocate(tem1(k,ns))
  allocate(tem2(1,k))
  allocate(tem3(ns1,ns1))
  allocate(kf(ns,k))
  allocate(akc(ns,ns))
  allocate(ishk(k))
  allocate(iobs(k))
  allocate(istate(ns))
  allocate(er(ns1))
  allocate(ei(ns1))

  !
  ! Read off estimation options and data from dat vector
  !
  ishk           = int(dat(3:k+2))
  iobs           = int(dat(k+3:2*k+2))
  ind            = int(dat(2*k+3:2*k+lpar+2))
  icon           = int(dat(2*k+lpar+3:2*k+2*lpar+2))
  lb             = dat(2*k+2*lpar+3:2*k+3*lpar+2)
  ub             = dat(2*k+3*lpar+3:2*k+4*lpar+2)
  l              = 2*k+4*lpar+2
  do j=1,k
    y(:,j)       = dat(l+(j-1)*t+1:l+j*t)
  enddo
  istate         = (/ 1, 1+ishk, 5+ishk, 10 /)

  !
  ! Default parameters
  !
  par            = 0.
  par(1)         = 0.01010135957165
  par(2)         = 0.01534982278818
  par(3)         = 0.00249067931432
  par(4)         = 1.8
  par(5)         = 1.000001
  par(6)         = 0.33
  par(7)         = 0.005
  par(8)         = .25
  par(10)        = log(.00001)
  par(16)        = .95
  par(21)        = .95
  par(27)        = 1.0
  par(32)        = 1.0
  par(37)        = 1.0
  alpha          = .99

  !
  ! Update default parameters in case of those to be estimated
  !

  l  = 1
  do i=1,lpar
    if (ind(i)==1) then
      if (icon(i)==0) then
        par(i)   = x(l)
      else
        par(i)   = .5*(ub(i)+lb(i))+.5*(ub(i)-lb(i))*x(l)/              &
                   sqrt(1.+x(l)*x(l))
      endif
      l          = l+1
    endif
  enddo
  inbadregion = 0
  do i=1,lpar
    if ((icon(i)==0).and.((par(i)>ub(i)).or.(par(i)<lb(i)))) then
      write(*,*) 'in bad region',i,x(i)
      inbadregion = 1
    endif
  enddo
  if (inbadregion==1) then
    write(*,*) 'going out'
    lnl = largenum
    return
  endif

  !
  ! Fill in the matrices for the stochastic processes
  !

  sbar(:,1)      = par(7:10)
  imat           = 0.
  p              = 0.
  q              = 0.
  l              = 11
  do i=1,4
    imat(i,i)    = 1.
    do j=1,4
      p(i,j)     = par(l)
      q(i,j)     = par(l+16)
      l          = l+1
    enddo
  enddo

  !
  ! Compute an equilibrium for the model 
  !

  call model(lpar,par, gamma,ss,psil,psiy,psix,psic,info)
  if (info.ne.0) then
    write(*,*) 'thrown out with equilibrium'
    write(*,*) x
    lnl          = largenum
    return
  endif

  !
  ! 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.
  a(1,1)         = gamma(1)
  do i=2,k+1
    l            = ishk(i-1)
    a(1,i)       = gamma(l+1)
    a(i+k,i)     = 1.
  enddo
  a(1,ns)        = gamma(6)
  p0             = matmul(imat-p,sbar)
  a(2:k+1,ns)    = p0(ishk,1)
  a(2:k+1,2:k+1) = p(ishk,ishk)
  a(ns,ns)       = 1.
  b(2:k+1,:)     = q(ishk,ishk)
  cdlogprd(1)    = (psiy(1)-psil(1))*(1.-1./gamma(1))
  cdlogprd(2:5)  = psiy(2:5)-psil(2:5)
  cdlogprd(6:9)  =-psiy(2:5)+psil(2:5)+(psiy(1)-psil(1))*gamma(2:5)/    &
                   gamma(1)
  cdlogprd(10)   = (psiy(1)-psil(1))*gamma(6)/gamma(1)
  cdlogprd(2)    = cdlogprd(2)+1.
  clogl          = 0.
  clogl(1:5)     = psil(1:5)
  clogl(10)      = psil(6)
  cqdlogl(1)     = psil(1)*(1.-alpha/gamma(1))
  cqdlogl(2:5)   = psil(2:5)
  cqdlogl(6:9)   = -alpha*(psil(2:5)-psil(1)*gamma(2:5)/gamma(1))
  cqdlogl(10)    = (1.-alpha)*psil(1)+alpha*psil(1)*gamma(6)/gamma(1)
  clogxy         = 0.
  clogxy(1:5)    = psix(1:5)-psiy(1:5)
  clogxy(10)     = psix(6)-psiy(6)
  clogcy         = 0.
  clogcy(1:5)    = psic(1:5)-psiy(1:5)
  clogcy(10)     = psic(6)-psiy(6)
  cdlogy(1)      = psiy(1)*(1.-1./gamma(1))
  cdlogy(2:5)    = psiy(2:5)
  cdlogy(6:9)    =-psiy(2:5)+psiy(1)*gamma(2:5)/gamma(1)
  cdlogy(10)     = psiy(1)*gamma(6)/gamma(1)
  cdlogy(2)      = cdlogy(2)+1.
  allc(1,:)      = cdlogprd 
  allc(2,:)      = clogl
  allc(3,:)      = cqdlogl
  allc(4,:)      = clogxy
  allc(5,:)      = clogcy
  allc(6,:)      = cdlogy
  alls0(1:5)     = ss(1:5)
  alls0(6:9)     = ss(2:5)
  alls0(10)      = ss(6)
  c              = allc(iobs,istate)
  s0(:,1)        = alls0(istate)

  !
  ! Check that maximum eigenvalue of A is less than 1
  !
  at      = a(1:ns1,1:ns1)
  call dgeev('N','N',ns1,at,ns1,er,ei,vl,1,vr,1,work,lwork,info)
  eiga  = sqrt(er(1)*er(1)+ei(1)*ei(1))
  do i=2,ns1
    tem = sqrt(er(i)*er(i)+ei(i)*ei(i))
    if (tem>eiga) eiga = tem
  enddo
  if (eiga>=1) then
    write(*,*) 'eiga>=1'
    write(*,*) 'x=',x
    lnl    = largenum
    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 B'
  !

  lnl            = 0.
  xhat           = s0
  bb             = matmul(b,transpose(b))
  vs             = bb(1:ns1,1:ns1)
  atem           = a(1:ns1,1:ns1)
  at             = -transpose(atem)
  call sb04qd(ns1,ns1,atem,ns1,at,ns1,vs,ns1,tem3,ns1,iwork,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in loglikefiml.f90: problem with sb04qd'
    write(*,*) 'info = ',info
    lnl = largenum
    return
  endif

  s              = 0.
  s(1:ns1,1:ns1) = vs

  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 dgetrf'
      write(*,*) 'info = ',info
      lnl = largenum
      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
      lnl = largenum
      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(s0)
  deallocate(u)
  deallocate(ome)
  deallocate(iome)
  deallocate(at)
  deallocate(vs)
  deallocate(tem1)
  deallocate(tem2)
  deallocate(tem3)
  deallocate(kf)
  deallocate(akc)
  deallocate(ishk)
  deallocate(iobs)
  deallocate(istate)
  deallocate(er)
  deallocate(ei)

end subroutine loglikefiml