program call
!
! RUNLIKE   Use an example from CKM (with the baseline parameters) to
!           test the loglike.m program
!
!
!           Ellen McGrattan, 6-2-06
!           Revised, ERM, 10-2-09

  use random
  implicit none
  integer, parameter                  :: n=24,t=200,k=3,ldat=t*k+2,          &
                                         nsimmle=1000,nboot=500,nperkp=200,  &
                                         nsimkp=500,lwork=5000,              &
                                         lrandstat=k*(nperkp*nsimkp)
  integer                             :: i,j,l,iseed,s,nseed,info,ise,istat, &
                                         nwant,nstat
  integer, dimension(n)               :: ipiv
  integer, dimension(:), allocatable  :: seed
  real, dimension(k,k)                :: b1,m,a,b,phi,thet,csig,sig,ahat,    &
                                         bhat,cshat,b1hat,mhat,sighat,       &
                                         mkbyk,skbyk,phihat,thethat,         &
                                         mb1hat,sb1hat,mmhat,smhat,          &
                                         mb1mhat,sb1mhat,mcshat,scshat
  real, dimension(k,1)                :: e,elag,ylag,alp,mkby1,skby1
  real, dimension(k,t)                :: y,ehat
  real, dimension(ldat)               :: dat
  real, dimension(n)                  :: x0,x,se,x1,mx,sx,mse,sse,xhld
  real, dimension(lwork)              :: work
  real, dimension(nsimmle)            :: logl
  real, dimension(n,nsimmle)          :: est
  real, dimension(lrandstat)          :: randstat
  real, dimension(:), allocatable     :: stat0,stat,mstat,sstat
  real, dimension(:,:), allocatable   :: dwant
  real                                :: lnl,mlnl,slnl
  external loglike11

  open(unit=5,file='run11.inp')
  open(unit=7,file='call11.est')
  open(unit=8,file='call11.lnl')
  open(unit=9,file='call11.out')
  open(unit=10,file='call11.stat')

  read(5,*) b1
  read(5,*) m
  read(5,*) sig
  read(5,*) ise
  read(5,*) istat
  read(5,*) alp
  read(5,*) nwant


  !
  ! Allocate space to matrices used when computing statistics
  !
  if (istat==0) nwant = 1
  nstat  = (3*k+1)*k+(nwant+k+2)*nwant
  allocate(dwant(nwant,nwant))
  allocate(stat0(nstat))
  allocate(stat(nstat))
  allocate(mstat(nstat))
  allocate(sstat(nstat))
  read(5,*) dwant

  !
  ! Set up the initial parameter vector x
  !
  phi     = b1+m
  thet    = 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
  !
  ! Use general matrices A and B that will be mapped to Phi and Theta
  !
  call istationary(k,phi,sig, a)
  call istationary(k,thet,sig, b)
  l = 1
  do i=1,k
    do j=1,k
      x0(l)  = a(i,j)
      l      = l+1
    enddo
  enddo
  do i=1,k
    do j=1,k
      x0(l)  = b(i,j)
      l      = l+1
    enddo
  enddo  
  do j=1,k
    do i=j,k
      x0(l) = csig(i,j)
      l     = l+1
    enddo
  enddo

  !
  ! Initialize the random generator's seed
  !
  call random_seed(size=nseed)
  allocate(seed(nseed))
  call random_seed(get=seed)
  call random_seed(put=seed)

  !
  ! Draw random numbers for bootstrapping or for computing statistics
  ! that will be held fixed across draws of new data sets for MLE
  !
  if (istat==1) then
    do i=1,lrandstat
      randstat(i) = random_normal()
    enddo
    !
    ! Compute statistics for the benchmark parameters
    !
    call stat11(k,b1,m,sig,alp,nwant,dwant,lrandstat,randstat,nperkp,        &
                nsimkp,nstat,  stat0)
  endif

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

  dat(1:2) = (/ float(t), float(k) /)
  mlnl     = 0.
  slnl     = 0.
  mx       = 0.
  sx       = 0.
  mse      = 0.
  sse      = 0.
  mstat    = 0.
  sstat    = 0.
  mb1hat   = 0.
  sb1hat   = 0.
  mmhat    = 0.
  smhat    = 0.
  mb1mhat  = 0.
  sb1mhat  = 0.
  mcshat   = 0.
  scshat   = 0.
  do s=1,nsimmle
    !
    ! 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(phi,ylag)+e-matmul(thet,elag)
      y(:,j)    = ylag(:,1)
      elag      = e
    enddo  
    do i=1,k
      dat((i-1)*t+3:i*t+2) = y(i,:)
    enddo
        
    !
    ! Compute MLE estimates for B1,M,Sig
    !

    xhld  = x0
    call uncmnd(n,xhld,ldat,dat,loglike11,x,lnl,info,work,lwork)
    mlnl  = mlnl + lnl
    slnl  = slnl + lnl*lnl

    l = 1
    do i=1,k
      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  
    do j=1,k
      do i=j,k
        cshat(i,j) = x(l) 
        l          = l+1
      enddo
    enddo
    sighat     = matmul(cshat,transpose(cshat))
    call stationary(k,ahat,sighat, phihat)
    call stationary(k,bhat,sighat, thethat)
    b1hat      = phihat-thethat
    mhat       = thethat
    l          = 1
    do i=1,k
      do j=1,k
        est(l,s) = b1hat(i,j)
        l        = l+1
      enddo
    enddo
    do i=1,k
      do j=1,k
        est(l,s) = mhat(i,j)
        l        = l+1
      enddo
    enddo
    est(2*k*k+1:n,s) = x(2*k*k+1:n) 
    logl(s)          = lnl
    mx           = mx      + est(:,s)
    sx           = sx      + est(:,s)*est(:,s)
    mb1hat       = mb1hat  + b1hat
    sb1hat       = sb1hat  + b1hat*b1hat
    mmhat        = mmhat   + mhat
    smhat        = smhat   + mhat*mhat
    mb1mhat      = mb1mhat + b1hat+mhat
    sb1mhat      = sb1mhat + (b1hat+mhat)*(b1hat+mhat)
    mcshat       = mcshat  + cshat
    scshat       = scshat  + cshat*cshat

    !
    ! Compute statistics
    ! 
    if (istat==1) then
      call stat11(k,b1hat,mhat,sighat,alp,nwant,dwant,lrandstat,randstat,    &
                  nperkp,nsimkp,nstat,  stat)
      mstat = mstat + stat
      sstat = sstat + stat*stat
      write(10,'(60(1x,f13.4)') stat
    endif

    !
    ! If ise=1, compute standard errors for first MLE 
    !
    if ((ise==1).and.(s==1)) then
      x1  = x
      call stderr11(n,x,ldat,dat,nboot, se)
    endif
    if (ise>1) then
      call stderr11(n,x,ldat,dat,nboot, se)
      mse = mse + se
      sse = sse + se*se
    endif
  enddo
  mse      = mse/float(nsimmle)
  sse      = (sse-float(nsimmle)*mse*mse)/float(nsimmle-1)
  sse      = sqrt(sse)

  mx       = mx/float(nsimmle)
  sx       = (sx-float(nsimmle)*mx*mx)/float(nsimmle-1)
  sx       = sqrt(sx)

  mb1hat   = mb1hat/float(nsimmle)
  sb1hat   = (sb1hat-float(nsimmle)*mb1hat*mb1hat)/float(nsimmle-1)
  sb1hat   = sqrt(sb1hat)
  mmhat    = mmhat/float(nsimmle)
  smhat    = (smhat-float(nsimmle)*mmhat*mmhat)/float(nsimmle-1)
  smhat    = sqrt(smhat)
  mb1mhat  = mb1mhat/float(nsimmle)
  sb1mhat  = (sb1mhat-float(nsimmle)*mb1mhat*mb1mhat)/float(nsimmle-1)
  sb1mhat  = sqrt(sb1mhat)
  mcshat   = mcshat/float(nsimmle)
  scshat   = (scshat-float(nsimmle)*mcshat*mcshat)/float(nsimmle-1)
  scshat   = sqrt(scshat)

  mstat    = mstat/float(nsimmle)
  sstat    = (sstat-float(nsimmle)*mstat*mstat)/float(nsimmle-1)
  sstat    = sqrt(sstat)
  !
  ! Print all results
  ! 
  write(9,'(a41)') '========================================='
  write(9,'(a41)') ' Estimates and Statistics for VARMA(1,1) '
  write(9,'(a41)') '========================================='
  write(9,'(a41)') '  Model:                                 '
  write(9,'(a41)') '    Y[t] =(B1+M) Y[t-1] +e[t]-M e[t-1]   '
  write(9,'(a41)') '    Eee'' =Sig = A0*A0''                   '
  write(9,*)
  write(9,'(a41)') '  Possible observables (Y):              '
  write(9,'(a41)') '    (1) Delta log(y[t]/l[t])             '
  write(9,'(a41)') '    (2) log(l[t])                        '
  write(9,'(a41)') '    (3) (1-.99L) log(l[t])               '
  write(9,'(a41)') '    (4) log(x[t]/y[t])                   '
  write(9,'(a41)') '    (5) log(c[t]/y[t])                   '
  write(9,'(a41)') '    (6) Delta log(y[t])                  '
  write(9,*)
  write(9,'(a41)') '  Possible unobserved shocks:            '
  write(9,'(a41)') '    (1) log z[t]                         '
  write(9,'(a41)') '    (2) taul[t]                          '
  write(9,'(a41)') '    (3) taux[t]                          '
  write(9,'(a41)') '    (4) log(ghat[t])                     '
  write(9,*)
  write(9,'(a41)') '  Variables used for KP-statistics       '
  write(9,'(a41)') '    X[i,t] = alp[i] X[i,t-1] + Y[i,t]    '
  write(9,'(a41)') '    W[t]   = D X[t]                      '
  write(9,*)
  write(9,'(a41)') '  Options set for results below:         '
  write(9,'(a17,f7.3)')         '    alp  =       ',alp(1,1)
  do i=2,nwant
    write(9,'(a17,f7.3)')       '                 ',alp(i,1)
  enddo
  write(9,'(a17,6(f7.3,1x))')   '    D    =       ',dwant(1,:)
  do i=2,nwant
    write(9,'(a17,6(f7.3,1x))') '                 ',dwant(i,:)
  enddo
  write(9,*)
  write(9,'(a41)') '  Initial Argument of Log(L)             '
  write(9,*)
  do i=1,n
    write(9,'(a12,1x,f7.3)')     '            ',x0(i)
  enddo
  write(9,*)
  write(9,'(a41)') '  Initial matrices:                      '
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    B1    = ',b1(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',b1(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    M     = ',m(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',m(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '   B1+M   = ',b1(1,:)+m(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',b1(i,:)+m(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   ' chol(Sig)= ',csig(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',csig(i,:)
  enddo

  write(9,*)
  write(9,'(a41)') '  ---------------------------------------'
  write(9,'(a41)') '         Statistics for True DGP         '
  write(9,'(a41)') '  ---------------------------------------'
  write(9,*)
  write(9,'(a42)') '    Matrix A0 Measuring Impact of 1% Shock'
  l = 1
  do i=1,k
    do j=1,k
      mkbyk(i,j) = stat0(l)
      l          = l+1
    enddo
  enddo
  do i=1,k
    write(9,'(10x,10(1x,f7.3))') mkbyk(i,:)
  enddo
  write(9,*)
  write(9,*) '    Variance of Y '
  write(9,*) 
  do i=1,k
    write(9,'(10x,(1x,f7.3))') stat0(l)
    l          = l+1
  enddo
  do i=1,k
    do j=1,k
      mkbyk(i,j) = stat0(l)
      l          = l+1
    enddo
  enddo
  write(9,*)
  write(9,*) '    Decomposition of Variance(Y), Sig=A0*A0''      '
  write(9,*)
  do i=1,k
    write(9,'(10x,10(1x,f7.3))') mkbyk(i,:)
  enddo
  do i=1,k
    do j=1,k
      mkbyk(i,j) = stat0(l)
      l          = l+1
    enddo
  enddo
  write(9,*)
  write(9,*) '    Decomposition of Variance(Y), chol(Sig)      '
  write(9,*)
  do i=1,k
    write(9,'(10x,10(1x,f7.3))') mkbyk(i,:)
  enddo
  write(9,*)
  write(9,*) '    KP-statistics '
  write(9,*)
  do i=1,(nwant+k+2)*nwant
    write(9,'(10x,2(1x,f7.3))') stat0(l)
    l         = l+1
  enddo

  write(9,*)
  write(9,'(a41)') '  ---------------------------------------'
  write(9,'(a41)') '       Statistics for MLE Estimates      '
  write(9,'(a41)') '  ---------------------------------------'
  write(9,*)
  write(9,'(a19,i4,a16,i4,a8)') '  Laboratory data: ',nsimmle,               &
                                ' simulations of ',t,' periods'
  write(9,*)
  write(9,'(a41)') '  Argument of Log(L), Mean and Std.Dev.  '
  write(9,*)
  do i=1,n
    write(9,'(a12,2(1x,f7.3))')   '            ',mx(i),sx(i)
  enddo
  write(9,*)
  write(9,'(a41)') '  Means of Estimated Matrices:           '
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    B1    = ',mb1hat(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',mb1hat(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    M     = ',mmhat(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',mmhat(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '   B1+M   = ',mb1mhat(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',mb1mhat(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   ' chol(Sig)= ',mcshat(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',mcshat(i,:)
  enddo
  write(9,*)
  write(9,*)
  write(9,'(a41)') '  Std. Dev. of Estimated Matrices:       '
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    B1    = ',sb1hat(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',sb1hat(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    M     = ',smhat(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',smhat(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '   B1+M   = ',sb1mhat(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',sb1mhat(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   ' chol(Sig)= ',scshat(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',scshat(i,:)
  enddo
  write(9,*)
  write(9,'(a42)') '  Statistics averaged across MLE estimates'
  l = 1
  do i=1,k
    do j=1,k
      mkbyk(i,j) = mstat(l)
      skbyk(i,j) = sstat(l)
      l          = l+1
    enddo
  enddo
  write(9,*) '    Matrix A0 Measuring Impact of 1% Shock '
  write(9,*) ' '
  write(9,*) '      Mean:'
  do i=1,k
    write(9,'(10x,10(1x,f7.3))') mkbyk(i,:)
  enddo
  write(9,*)
  write(9,*) '      Std. Dev.:'
  do i=1,k
    write(9,'(10x,10(1x,f7.3))') skbyk(i,:)
  enddo
  write(9,*)
  write(9,*) '    Variance of Y '
  write(9,*)
  write(9,*) '      Mean and Std.Dev.:'
  do i=1,k
    write(9,'(10x,2(1x,f7.3))') mstat(l),sstat(l)
    l          = l+1
  enddo
  do i=1,k
    do j=1,k
      mkbyk(i,j) = mstat(l)
      skbyk(i,j) = sstat(l)
      l          = l+1
    enddo
  enddo
  write(9,*)
  write(9,*) '    Decomposition of Variance(Y), Sig=A0*A0''      '
  write(9,*)
  write(9,*) '      Mean:'
  do i=1,k
    write(9,'(10x,10(1x,f7.3))') mkbyk(i,:)
  enddo
  write(9,*)
  write(9,*) '      Std. Dev.:'
  do i=1,k
    write(9,'(10x,10(1x,f7.3))') skbyk(i,:)
  enddo
  do i=1,k
    do j=1,k
      mkbyk(i,j) = mstat(l)
      skbyk(i,j) = sstat(l)
      l          = l+1
    enddo
  enddo
  write(9,*)
  write(9,*) '    Decomposition of Variance(Y), chol(Sig)      '
  write(9,*)
  write(9,*) '      Mean:'
  do i=1,k
    write(9,'(10x,10(1x,f7.3))') mkbyk(i,:)
  enddo
  write(9,*)
  write(9,*) '      Std. Dev.:'
  do i=1,k
    write(9,'(10x,10(1x,f7.3))') skbyk(i,:)
  enddo
  write(9,*)
  write(9,*) '    KP-statistics '
  write(9,*)
  write(9,*) '      Mean and Std.Dev.:'
  do i=1,(nwant+k+2)*nwant
    write(9,'(10x,2(1x,f7.3))') mstat(l),sstat(l)
    l         = l+1
  enddo
                                                                                
  do s=1,nsimmle
    write(7,'(30(1x,f11.6)') est(:,s)
    write(8,*) logl(s)
  enddo

  deallocate(seed)
  deallocate(dwant)
  deallocate(stat)
  deallocate(mstat)
  deallocate(sstat)

end program call