program call
!
! CALLFIML   is the driver for a laboratory testing the model in 
!            `Business Cycle Accounting.'
!
! ishk = 1 x k listing the shocks to include
! ishk = [1,3] implies a 2 shock model with shocks to 
!              technology and investment taxes
! iobs = 1 x k listing the observables to include from:
!              1: difference in log labor productivity
!              2: quasi-difference in labor (alp)
!              3: labor
!              4: investment-output
!              5: consumption-output
!              6: difference in log output
! iobs = [3,6] implies the observables are labor and the change in output
!
!           Ellen McGrattan, 6-2-06
!           Revised, ERM, 10-2-09

  use random
  implicit none
  integer, parameter                  :: t=200,k=3,ns=2*k+2,                 &
                                         nsimmle=1000,nboot=500,nperkp=200,  &
                                         nsimkp=500,lwork=5000,              &
                                         lrandstat=k*(nperkp*nsimkp),        &
                                         lparam=42,ldat=(t+2)*k+2+4*lparam 
  integer                             :: i,j,l,iseed,s,nseed,info,ise,istat, &
                                         nwant,nstat,n
  integer, dimension(lparam)          :: ind,icon
  integer, dimension(k)               :: ishk,iobs
  integer, dimension(ns)              :: istate
  integer, dimension(:), allocatable  :: seed
  real, dimension(ns,ns)              :: a,ahat,mahat,sahat
  real, dimension(ns,k)               :: b,bhat,mbhat,sbhat
  real, dimension(k,ns)               :: c,chat,mchat,schat
  real, dimension(ns,1)               :: state,s0,s0hat,ms0hat,ss0hat
  real, dimension(k,k)                :: mkbyk,skbyk
  real, dimension(k,1)                :: eps,alp
  real, dimension(k,t)                :: y
  real, dimension(ldat)               :: dat
  real, dimension(lwork)              :: work
  real, dimension(nsimmle)            :: logl
  real, dimension(lrandstat)          :: randstat
  real, dimension(lparam)             :: param,par0,lb,ub
  real, dimension(6)                  :: gamma,ss,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, dimension(:), allocatable     :: stat0,stat,mstat,sstat,x,x0,x1,     &
                                         xhld,se,mx,sx,mse,sse
  real, dimension(:,:), allocatable   :: dwant,est
  real                                :: lnl,mlnl,slnl,alpha,tem
  external loglikefiml

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

  !
  ! Read in model parameters (irate,delta,gn,psi,sigma,theta,sbar,p,q) 
  ! and options:
  !
  n  = 0
  do i=1,lparam
    read(5,*) param(i),ind(i),icon(i),lb(i),ub(i)
    if (ind(i)==1) n = n+1
  enddo
  read(5,*) ishk
  read(5,*) iobs
  read(5,*) ise
  read(5,*) istat
  read(5,*) alp
  read(5,*) nwant
  istate   = (/ 1, 1+ishk, 5+ishk, 10 /)
  alpha    = .99

  !
  ! Allocate space to parameter vectors
  !
  allocate(x(n))
  allocate(x0(n))
  allocate(x1(n))
  allocate(xhld(n))
  allocate(se(n))
  allocate(mx(n))
  allocate(sx(n))
  allocate(mse(n))
  allocate(sse(n))
  allocate(est(n,nsimmle))
  
  !
  ! Compute equilibrium for initial parameters
  !

  call model(lparam,param, gamma,ss,psil,psiy,psix,psic,info)

  !
  ! Set up state space for initial parameters
  !

  sbar(:,1)      = param(7:10)
  l              = 11
  do i=1,4
    imat(i,i)    = 1.
    do j=1,4
      p(i,j)     = param(l)
      q(i,j)     = param(l+16)
      l          = l+1
    enddo
  enddo
  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)
  par0           = param

  !
  ! Allocate space to matrices used when computing statistics
  !
  if (istat==0) nwant = 1
  nstat  = (2*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
  !

  l  = 1
  do i=1,lparam
    if (ind(i)==1) then
      if (icon(i)==0) then
        x0(l)  = param(i)
      else
        tem    = (2.*param(i)-ub(i)-lb(i))/(ub(i)-lb(i))
        x0(l)  = tem/sqrt(1.-tem*tem)
      endif 
      l      = l+1
    endif
  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 statfiml(ns,k,a,b,c,s0,alp,nwant,dwant,lrandstat,randstat,nperkp,   &
                  nsimkp,nstat,  stat0)
  endif

  !
  ! Simulate data, estimate parameters, and compute statistics 
  !
  !l          = 2+2*k+4*lparam
  !dat(1:l)   = (/ float(t), float(k), float(ishk), float(iobs), float(ind), &
  !                float(icon), lb, ub /)
  l           = 2+2*k+2*lparam
  dat(1:l)    = (/ float(t), float(k), float(ishk), float(iobs), float(ind), &
                   float(icon) /)
  dat(l+1:l+lparam) = lb
  l           = 2+2*k+3*lparam
  dat(l+1:l+lparam) = ub
 
  mlnl        = 0.
  slnl        = 0.
  mx          = 0.
  sx          = 0.
  mse         = 0.
  sse         = 0.
  mstat       = 0.
  sstat       = 0.
  mahat       = 0.
  sahat       = 0.
  mbhat       = 0.
  sbhat       = 0.
  mchat       = 0.
  schat       = 0.
  ms0hat      = 0.
  ss0hat      = 0.
  do s=1,nsimmle
    !
    ! Generate a Txk matrix of observables y
    !
    state     = s0
    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+2*k+4*lparam
    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,loglikefiml,x,lnl,info,work,lwork)
    mx        = mx + x
    sx        = sx + x*x
    mlnl      = mlnl + lnl
    slnl      = slnl + lnl*lnl

    logl(s)   = lnl

    !
    ! Fill in new estimates
    !

    l  = 1
    do i=1,lparam
      if (ind(i)==1) then
        if (icon(i)==0) then
          param(i)    = x(l)
        else
          param(i)    = .5*(ub(i)+lb(i))+.5*(ub(i)-lb(i))*x(l)/              &
                         sqrt(1.+x(l)*x(l))
        endif
        est(l,s)      = param(i)
        l             = l+1
      endif
    enddo
   

    call model(lparam,param, gamma,ss,psil,psiy,psix,psic,info)

    sbar(:,1)         = param(7:10)
    l                 = 11
    do i=1,4
      do j=1,4
        p(i,j)        = param(l)
        q(i,j)        = param(l+16)
        l             = l+1
      enddo
    enddo
    ahat              = 0.
    bhat              = 0.
    ahat(1,1)         = gamma(1)
    do i=2,k+1
      l               = ishk(i-1)
      ahat(1,i)       = gamma(l+1)
      ahat(i+k,i)     = 1.
    enddo
    ahat(1,ns)        = gamma(6)
    p0                = matmul(imat-p,sbar)
    ahat(2:k+1,ns)    = p0(ishk,1)
    ahat(2:k+1,2:k+1) = p(ishk,ishk)
    ahat(ns,ns)       = 1.
    bhat(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)
    chat              = allc(iobs,istate)
    s0hat(:,1)        = alls0(istate)
    mahat             = mahat + ahat
    sahat             = sahat + ahat*ahat
    mbhat             = mbhat + bhat
    sbhat             = sbhat + bhat*bhat
    mchat             = mchat + chat
    schat             = schat + chat*chat
    ms0hat            = ms0hat+ s0hat
    ss0hat            = ss0hat+ s0hat*s0hat

    !
    ! Compute statistics
    ! 
    if (istat==1) then
      call statfiml(ns,k,ahat,bhat,chat,s0hat,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 stderr(n,x,ldat,dat,nboot, se)
    !endif
    !if (ise>1) then
    !  call stderr(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)
  mbhat       = mbhat/float(nsimmle)
  sbhat       = (sbhat-float(nsimmle)*mbhat*mbhat)/float(nsimmle-1)
  sbhat       = sqrt(sbhat)
  mchat       = mchat/float(nsimmle)
  schat       = (schat-float(nsimmle)*mchat*mchat)/float(nsimmle-1)
  schat       = sqrt(schat)
  ms0hat      = ms0hat/float(nsimmle)
  ss0hat      = (ss0hat-float(nsimmle)*ms0hat*ms0hat)/float(nsimmle-1)
  ss0hat      = sqrt(ss0hat)

  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 DSGE 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,6(1x,i2))')   '    Observables: ',iobs
  write(9,'(a17,6(1x,i2))')   '    Shocks:      ',ishk
  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 utility/technology parameters: '
  write(9,'(a32,f7.3)') '    interest rate             = ', par0(1)
  write(9,'(a32,f7.3)') '    capital depreciation      = ', par0(2)
  write(9,'(a32,f7.3)') '    growth in population      = ', par0(3)
  write(9,'(a32,f7.3)') '    utility weight on leisure = ', par0(4)
  write(9,'(a32,f7.3)') '    utility curvature         = ', par0(5)
  write(9,'(a32,f7.3)') '    capital share             = ', par0(6)
  write(9,*)
  write(9,'(a58)') '  Initial Argument of Log(L), Bounds, and Constrain Option'
  write(9,*) 
  j  = 1
  do i=1,lparam
    if (ind(i)==1) then
      write(9,'(a12,3(1x,f7.3),3x,i1)') '            ',x0(j),lb(i),          &
                                                       ub(i),icon(i)
      j = j+1
    endif
  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,'(a12,10(1x,f7.3))')   '    S0    = ',s0(:,1)
  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))')   '    B     = ',mbhat(1,:)
  do i=2,ns
    write(9,'(a12,10(1x,f7.3))') '            ',mbhat(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,'(a12,10(1x,f7.3))')   '    S0    = ',ms0hat(:,1)
  write(9,*) 
  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))')   '    B     = ',sbhat(1,:)
  do i=2,ns
    write(9,'(a12,10(1x,f7.3))') '            ',sbhat(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,'(a12,10(1x,f7.3))')   '    S0    = ',ss0hat(:,1)
  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,'(20(1x,f11.6)') est(:,s)
    write(8,*) logl(s)
  enddo

  deallocate(x)
  deallocate(x0)
  deallocate(x1)
  deallocate(xhld)
  deallocate(se)
  deallocate(mx)
  deallocate(sx)
  deallocate(mse)
  deallocate(sse)
  deallocate(est)
  deallocate(seed)
  deallocate(dwant)
  deallocate(stat)
  deallocate(mstat)
  deallocate(sstat)

end program call