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

  use random
  implicit none
  integer, parameter                  :: t=200,k=3,ns=2*k+1,n=2*k*k+3*k-7,   &
                                         nsimmle=1000,nboot=500,nperkp=200,  &
                                         nsimkp=500,lwork=5000,              &
                                         lrandstat=k*(nperkp*nsimkp),        &
                                         lind=2*k*ns+1,ldat=t*k+2+lind
  integer                             :: i,j,l,iseed,s,nseed,info,ise,istat, &
                                         nwant,nstat
  integer, dimension(lind)            :: ind
  integer, dimension(:), allocatable  :: seed
  real, dimension(ns,ns)              :: a,ahat,mahat,sahat
  real, dimension(ns,k)               :: b
  real, dimension(k,ns)               :: c,chat,mchat,schat
  real, dimension(ns,1)               :: state
  real, dimension(k,k)                :: mkbyk,skbyk
  real, dimension(k,1)                :: eps,alp,mkby1,skby1
  real, dimension(k-1)                :: rho,gamma
  real, dimension(k,t)                :: y
  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,gamk
  external loglikess

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

  read(5,*) gamk
  read(5,*) gamma
  do i=1,k-1
    read(5,*) rho(i)
  enddo
  read(5,*) c
  read(5,*) ise
  read(5,*) istat
  read(5,*) alp
  read(5,*) nwant
  a          = 0.
  a(1,1)     = gamk
  a(1,2)     = 1.
  do i=1,k-1
    a(1,i+2) = gamma(i)
  enddo
  do i=3,k+1
    a(i,i)   = rho(i-2)
  enddo
  do i=2,k+1
    a(i+k,i) = 1.
  enddo 
  do i=1,k
    b(i+1,i) = 1.
  enddo

  !
  ! Allocate space to matrices used when computing statistics
  !
  if (istat==0) nwant = 1
  nstat  = (2*k+1)*k+(nwant+2)*nwant+nwant*k
  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
  !
  x0(1)         = gamk/sqrt(1-gamk*gamk)
  x0(2:k)       = gamma
  l             = k+1
  do i=1,k-1    
    x0(l)       = rho(i)/sqrt(1-rho(i)*rho(i))
    l           = l+1
  enddo
  l = 2*k
  do j=1,ns
    x0(l)       = c(1,j)
    l           = l+1
  enddo
  do i=2,k
    do j=1,k+1
      x0(l)     = c(i,j)
      l         = l+1
    enddo
  enddo

  ind           = 0
  ind(1)        = 1
  ind(3:k+1)    = 1
  l             = 2*k+3
  do i=2,k
    ind(l)      = 1
    l           = l+k+1
  enddo
  l             = (2*k+1)*k+2
  do i=1,ns
    ind(l)      = 1
    l           = l+1
  enddo
  do i=1,k-1
    do j=1,k+1
      ind(l)      = 1
      l           = l+1
    enddo
    do j=k+2,ns
      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
    call statss(ns,k,a,b,c,alp,nwant,dwant,lrandstat,randstat,nperkp,        &
                nsimkp,nstat,  stat0)
  endif

  !
  ! Simulate data, estimate parameters, and compute statistics 
  !
  l           = 2*k*(2*k+1)+3
  dat(1:l)    = (/ float(t), float(k), float(ind) /)
  mlnl        = 0.
  slnl        = 0.
  mx          = 0.
  sx          = 0.
  mse         = 0.
  sse         = 0.
  mstat       = 0.
  sstat       = 0.
  mahat       = 0.
  sahat       = 0.
  mchat       = 0.
  schat       = 0.
  do s=1,nsimmle
    !
    ! Generate a Txk matrix of observables y
    !
    state      = 0.
    do j=1,t
      do i=1,k 
        eps(i,1) = random_normal()
      enddo
      state    = matmul(a,state)+matmul(b,eps)
      y(:,j)   = matmul(c,state(:,1))
    enddo  
    l          = 2*k*(2*k+1)+3
    do i=1,k
      dat(l+(i-1)*t+1:l+i*t) = y(i,:)
    enddo
                                                                                
    !
    ! Compute MLE estimates for elements of A,C
    !

    xhld      = x0
    call uncmnd(n,xhld,ldat,dat,loglikess,x,lnl,info,work,lwork)
    logl(s)       = lnl
    mlnl          = mlnl + lnl
    slnl          = slnl + lnl*lnl
    est(:,s)      = x
    est(1,s)      = x(1)/sqrt(1.+x(1)*x(1))
    do i=3,k+1
      est(i,s)    = x(i)/sqrt(1.+x(i)*x(i))
    enddo
    mx            = mx + est(:,s)
    sx            = sx + est(:,s)*est(:,s)

    !
    ! Fill estimates into matrices for computing statistics
    !
    ahat          = 0.
    chat          = 0.
    ahat(1,1)     = x(1)/sqrt(1.+x(1)*x(1))
    ahat(1,2)     = 1.
    do i=2,k
      ahat(1,i+1) = x(i)
    enddo
    l             = k+1
    do i=3,k+1
      ahat(i,i)   = x(l)/sqrt(1.+x(l)*x(l))
      l           = l+1
    enddo
    do i=2,k+1
      ahat(i+k,i) = 1.
    enddo 
    l             = 2*k
    do j=1,ns
      chat(1,j)   = x(l)
      l           = l+1
    enddo
    do i=2,k
      do j=1,k+1
        chat(i,j) = x(l)
        l         = l+1
      enddo
    enddo
 
    mahat         = mahat + ahat
    sahat         = sahat + ahat*ahat
    mchat         = mchat + chat
    schat         = schat + chat*chat

    !
    ! Compute statistics
    ! 
    if (istat==1) then
      call statss(ns,k,ahat,b,chat,alp,nwant,dwant,lrandstat,randstat,       &
                  nperkp,nsimkp,nstat,  stat)
      mstat     = mstat + stat
      sstat     = sstat + stat*stat
      write(10,'(50(1x,f11.6)') stat
    endif

    !
    ! If ise=1, compute standard errors for first MLE 
    !
    !if ((ise==1).and.(s==1)) then
    !  x1  = x
    !  call stderrss(n,x,ldat,dat,nboot, se)
    !endif
    !if (ise>1) then
    !  call stderrss(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)

  mahat       = mahat/float(nsimmle)
  sahat       = (sahat-float(nsimmle)*mahat*mahat)/float(nsimmle-1)
  sahat       = sqrt(sahat)
  mchat       = mchat/float(nsimmle)
  schat       = (schat-float(nsimmle)*mchat*mchat)/float(nsimmle-1)
  schat       = sqrt(schat)

  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 SS Model  '
  write(9,'(a41)') '========================================='
  write(9,'(a41)') '  Model:                                 '
  write(9,'(a41)') '    S[t+1] = A S[t] +B e[t],  Eee''=I    '
  write(9,'(a41)') '      Y[t] = C S[t]                      '
  write(9,'(a41)') '      S[t] = [log(khat[t]),s[t],s[t-1],1]'
  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 shocks (s):                   '
  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))')   '    A     = ',a(1,:)
  do i=2,ns
    write(9,'(a12,10(1x,f7.3))') '            ',a(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    B     = ',b(1,:)
  do i=2,ns
    write(9,'(a12,10(1x,f7.3))') '            ',b(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    C     = ',c(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',c(i,:)
  enddo
  write(9,*)
  write(9,'(a41)') '  ---------------------------------------'
  write(9,'(a41)') '         Statistics for True DGP         '
  write(9,'(a41)') '  ---------------------------------------'
  write(9,*)
  write(9,'(a43)') '    Matrix C*B 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 (Shock i in Column i)'
  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))')   '    A     = ',mahat(1,:)
  do i=2,ns
    write(9,'(a12,10(1x,f7.3))') '            ',mahat(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    C     = ',mchat(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',mchat(i,:)
  enddo
  write(9,*)
  write(9,'(a41)') '  Std. Dev. of Estimated Matrices:       '
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    A     = ',sahat(1,:)
  do i=2,ns
    write(9,'(a12,10(1x,f7.3))') '            ',sahat(i,:)
  enddo
  write(9,*)
  write(9,'(a12,10(1x,f7.3))')   '    C     = ',schat(1,:)
  do i=2,k
    write(9,'(a12,10(1x,f7.3))') '            ',schat(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 C*B 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 (Shock i in Column i)'
  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