program call
!
! RUNLIKE   Use an example from CKM (with the baseline parameters) to
!           test the loglike.m program
!
!
!           Ellen McGrattan

  use random
  implicit none
  integer, parameter          :: p=1,q=1,t=300,k=2,n=11,nimp=13,           &
                                 nsmax=5000,lwork=1000,ldat=t*k+4,bmax=1000
  integer                     :: i,j,nsim,iseed,s,seedsize,info,ise,l,     &
                                 nboot,iu
  integer, dimension(n)       :: ipiv
  integer, allocatable        :: seed(:)
  real, dimension(k,k)        :: b1,m,a,csig,sig,ahat,mhat,cshat
  real, dimension(k,1)        :: e,elag,ylag
  real, dimension(k,t)        :: y,ehat
  real, dimension(ldat)       :: dat
  real, dimension(n)          :: x0,x,xf,xb,xff,xfb,xbf,xbb,dl,se1,se2,    &
                                 xboot,bse,sumx,sumx2,xnew  
  real, dimension(n,t)        :: dlt
  real, dimension(n,1)        :: dltj
  real, dimension(n,bmax)     :: xball
  real, dimension(n,n)        :: imat,d2l
  real, dimension(lwork)      :: work
  real, dimension(nsmax)      :: logl
  real, dimension(n,nsmax)    :: est
  real, dimension(nimp,nsmax) :: imp
  real                        :: alpha,lnl,lnlt,lnlf,lnlb,lnltf,lnltb,     &
                                 lnlff,lnlfb,lnlbf,lnlbb,deli,delj,u,lnlnew
  external loglike

  open(unit=5,file='runlike.inp')
  open(unit=7,file='seeds.inp')

  read(5,*) iseed
  read(5,*) nsim
  read(5,*) alpha
  read(5,*) b1
  read(5,*) m
  read(5,*) sig
  read(5,*) ise
  read(5,*) nboot
  if (nsim>nsmax) write(*,*) 'Increase nsmax in runlike.f90'

  a       = b1+m
  csig    = sig
  call dpotrf('L',k,csig,k,info)
  do j=2,k
    do i=1,j-1
      csig(i,j) = 0.
    enddo
  enddo

  l = 1
  do j=1,k
    do i=1,k
      x0(l)  = a(i,j)
      l      = l+1
    enddo
    do i=1,k
      x0(l)  = m(i,j)
      l      = l+1
    enddo
  enddo  
  l = 2*k*k+1
  do j=1,k
    do i=j,k
      x0(l) = csig(i,j)
      l     = l+1
    enddo
  enddo

  !
  ! If iseed=1, initialize the random generator's seed with 
  ! vector in seeds.inp
  !
  if (iseed==1) then
    call random_seed(size=seedsize)
    allocate(seed(seedsize))
    read(7,*) seed
    call random_seed(put=seed)
    deallocate(seed)
  endif

  !
  ! Simulate data, estimate parameters, and compute statistics nsim times
  !

  write(*,*) '---------------------------------------'
  write(*,*) ' MLE ESTIMATION OF Y=A Y(-1)+e-B e(-1) '
  write(*,*) '---------------------------------------'
  write(*,'(a25)') '   Parameters of the True Model '
  write(*,*) '  -------------------------------------------------'
  write(*,*) '         A              B             chol(Eee'')'
  write(*,*) '  -------------------------------------------------'
  write(*,'(3x,2F7.3,3x,2F7.3,3x,F7.3)') x0(1),x0(5),x0(3),x0(7),x0(9)
  write(*,'(3x,2F7.3,3x,2F7.3,3x,2F7.3)') x0(2),x0(6),x0(4),x0(8),x0(10),x0(11)
  write(*,*) 
  write(*,*) 

  dat(1:4) = (/ p,q,t,k /)
  do s=1,nsim

    !
    ! Generate a Txk matrix of observables y
    !

    ylag   = 0.
    elag   = 0.
    do j=1,t
      do i=1,k 
        e(i,1)  = random_normal()
      enddo
      e         = matmul(csig,e)
      ylag      = matmul(a,ylag)+e-matmul(m,elag)
      y(:,j)    = ylag(:,1)
      elag      = e
    enddo  
    do i=1,k
      dat((i-1)*t+5:i*t+4) = y(i,:)
    enddo
    !
    ! Compute MLE estimates for B1,M,Sig
    !

    call uncmnd(n,x0,ldat,dat,loglike,x,lnl,info,work,lwork)
    est(1:n,s) = x
    logl(s)    = lnl

    !
    ! If ise=1, compute standard errors using two estimates 
    ! of the information matrix
    !
    if (ise==1) then
      d2l  = 0.
      imat = 0.
      do j=1,n
        delj  = max(abs(x(j))*.0001,.00000001)
        xf    = x
        xb    = x
        xf(j) = x(j)+delj
        xb(j) = x(j)-delj
        call loglike(n,xf,ldat,dat, lnlf)
        call loglike(n,xb,ldat,dat, lnlb)
        dl(j) = .5*(lnlf-lnlb)/delj
        do i=1,t
          call loglike_term(n,xf,ldat,dat,i, lnltf)
          call loglike_term(n,xb,ldat,dat,i, lnltb)
          dlt(j,i) = .5*(lnltf-lnltb)/delj
        enddo
        do i=1,j
          deli   = max(abs(x(i))*.0001,.00000001)
          xff    = x
          xfb    = x
          xbf    = x
          xbb    = x
          xff(i) = xff(i) + deli
          xff(j) = xff(j) + delj
          xfb(i) = xfb(i) + deli
          xfb(j) = xfb(j) - delj
          xbf(i) = xbf(i) - deli
          xbf(j) = xbf(j) + delj
          xbb(i) = xbb(i) - deli
          xbb(j) = xbb(j) - delj
          call loglike(n,xff,ldat,dat, lnlff)
          call loglike(n,xfb,ldat,dat, lnlfb)
          call loglike(n,xbf,ldat,dat, lnlbf)
          call loglike(n,xbb,ldat,dat, lnlbb)
          d2l(j,i) = .25*(lnlff-lnlfb-lnlbf+lnlbb)/(delj*deli)
          if (i/=j) d2l(i,j) = d2l(j,i)
        enddo
      enddo
      do j=1,t
        dltj(1:n,1) = dlt(1:n,j)
        imat        = imat + matmul(dltj,transpose(dltj))
      enddo
      call dgetrf(n,n,imat,n,ipiv,info)
      call dgetri(n,imat,n,ipiv,work,lwork,info)
      call dgetrf(n,n,d2l,n,ipiv,info)
      call dgetri(n,d2l,n,ipiv,work,lwork,info)
      do i=1,n
        se1(i) = sqrt(imat(i,i))
        se2(i) = sqrt(d2l(i,i))
      enddo 
      bse       = 0.
      sumx      = 0.
      sumx2     = 0.

      l = 1
      do j=1,k
        do i=1,k
          ahat(i,j) = x(l)  
          l         = l+1
        enddo
        do i=1,k
          mhat(i,j) = x(l) 
          l         = l+1
        enddo
      enddo  
      !l = 2*k*k+1
      !do j=1,k
      !  do i=j,k
      !    cshat(i,j) = x(l) 
      !    l         = l+1
      !  enddo
      !enddo
      ylag   = 0.
      elag   = 0.
      do j=1,t
        e(:,1)    = y(:,j)
        e         = e-matmul(ahat,ylag)+matmul(mhat,elag)
        ehat(:,j) = e(:,1)
        ylag(:,1) = y(:,j)
        elag      = e
      enddo
      do l=1,nboot
        ylag   = 0.
        elag   = 0.
        do j=1,t
          call random_number(u)
          iu       = ceiling(u*float(t)+.0000000001)
          e(:,1)   = ehat(:,iu)
          ylag     = matmul(ahat,ylag)+e-matmul(mhat,elag)
          y(:,j)   = ylag(:,1)
          elag     = e
        enddo  
        do i=1,k
          dat((i-1)*t+5:i*t+4) = y(i,:)
        enddo
        call uncmnd(n,x,ldat,dat,loglike,xboot,lnl,info,work,lwork)
        xball(:,l) = xboot
        sumx       = sumx+xboot
        sumx2      = sumx2+xboot*xboot
      enddo
      do l=1,nboot
        write(20,'(11f11.4)') (xball(i,l),i=1,11)
      enddo
      write(*,'(11f11.4)') sumx
      write(*,*)
      write(*,'(11f11.4)') sumx2
      bse   = sqrt((sumx2-sumx*sumx/float(nboot))/float(nboot-1))
      do j=1,n
        write(*,*) se1(j),se2(j),bse(j)
      enddo
    endif
  enddo
  !
  ! print results
  ! 
  x  = est(:,1)
  write(*,'(a25,i4)') '   Results for Simulation ',s 
  write(*,*) '  -------------------------------------------------'
  write(*,*) '         A              B             chol(Eee'')'
  write(*,*) '  -------------------------------------------------'
  write(*,'(3x,2F8.4,1x,2F8.4,1x,F7.3)') x(1),x(5),x(3),x(7),x(9)
  write(*,'(3x,2F8.4,1x,2F8.4,1x,2F7.3)') x(2),x(6),x(4),x(8),x(10),x(11)
  write(*,*) 
  write(*,'(A9,F10.3)') '   Ln(L):',lnl
  write(*,*) 
  !call varmastat(n,x, stat)
  !write(*,*) '   Statistics:'

end program call