subroutine stats1(k,b1,m,sig,alp,nwant,swant,lrand,randstat,nper,nsim,lstat, &
                  stat)
!
! Compute statistics for the estimated model:
!
!    Y[t] = (B1+M) Y[t-1] + e[t] - M e[t-1],  Eee'=Sig
!   (k x 1)  
! 
! Statistics reported:
!
!    (1) A0(ivar,ishk) from a SVAR of Y[t] regressed on infinite lags,
!
!        Y[t] = B1 Y[t-1] + M B1 Y[t-2] + M^2 B1 Y[t-3]+ ... + e[t]
!             = B1 Y[t-1] + B2 Y[t-2] + B3 Y[t-3]+ ... + e[t]  
!
!        A0   = inv(Cbar) * chol( Cbar* Sig *Cbar')'
!        Cbar = inv(I-B1-B2-...)
!             = inv(I-inv(I-M)*B1)
!
!        ==> stat(1:k*k) = [A0(1,:),A0(2,:),...,A0(k,:)]
!
!    (2) Population variance decomposition for Y using state space:
!
!        Z[t+1] = [B1+M,  I] Z[t] + [ I ] e[t+1],  
!                 [   0,  0]        [-M ]        
!
!          Y[t] = [   I,  0] Z[t]
!
!        done first with Eee'=A0*A0' and then with Eee'=L*L', where
!        L is Cholesky lower triangular
!
!        ==> stat(k*k+1:(3k+1)*k) = [var(Y[t]|all shocks)
!                                    var(Y[t]|shock 1 only, Sig=A0*A0')
!                                          :
!                                    var(Y[t]|shock k only, Sig=A0*A0')
!                                    var(Y[t]|shock 1 only, Sig=L*L')
!                                          :
!                                    var(Y[t]|shock k only, Sig=L*L') 
!                
!    (3) Second moments for variables of interest, W[t], where
!
!        W[t]   = S* X[t],   S is n x k
!        Y[i,t] = (1-alp(i)L) X[i,t] 
!               
!        and variables in W[t] HP-filtered.  The statistics reported
!        are standard deviations, autocorrelations, and contemporaneous 
!        correlations of HP-filtered W[t] for nsim simulations of the 
!        model for nper periods.
!     
!        ==> stat((3k+1)*k+1:lstat) = [std(HP W[t])
!                                      cor(HP W[t], HP W[t-1])
!                                      cor(HP W[:,t], HP W[1,t])
!                                            :
!                                      cor(HP W[:,t], HP W[nwant,t])]
!            where lstat = (3k+1)*k+(nwant+2)*nwant
!

  implicit none
  
  integer, parameter                    :: lwork=5000
  integer, intent(in)                   :: k,nwant,lrand,nper,nsim,lstat
  integer                               :: i,j,l,kk,info,n,s
  integer, dimension(k)                 :: ipiv
  integer, dimension(lwork)             :: iwork
  real, dimension(k,k), intent(in)      :: b1,m,sig
  real, dimension(lrand), intent(in)    :: randstat
  real, dimension(lstat), intent(out)   :: stat
  real, dimension(k,k)                  :: imat,cbar,cbari,a0,csig,csigc,ei, &
                                           phi,tem
  real, dimension(2*k,2*k)              :: a,at,v,vj,tem2,atem
  real, dimension(2*k,k)                :: b
  real, dimension(k,2*k)                :: tem1
  real, dimension(k,1)                  :: y,ylag,e,elag,x,alp
  real, dimension(nwant,1)              :: wtem
  real, dimension(nwant)                :: sum1,sum2,sum3,sum4,sum5,sum6,    &
                                           sum7,mean,meanf,meanl,std,var,    &
                                           varf,varl,acov,acor
  real, dimension(nwant*nwant)          :: sum8,cor
  real, dimension(nwant,nper)           :: w
  real, dimension(nwant,k)              :: swant
  real, dimension(nper)                 :: wn,wt,wd
  real, dimension((nwant+2)*nwant)      :: kp
  real, dimension(lwork)                :: work


  stat = 0.
  if (lstat<(3*k+1)*k+1+(nwant+2)*nwant) then
    write(*,*) 'ERROR: In call11.f90, fix ns in call to stats.f90'
  endif

  !========================
  ! First set of statistics
  !========================
  !
  ! Cbar^(-1) = I-inv(I-M)*B1
  !
  imat = 0.
  do i=1,k
    imat(i,i)  = 1.
  enddo
  cbari = imat-m
  call dgetrf(k,k,cbari,k,ipiv,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with dgetrf'
    write(*,*) 'info = ',info
    return
  endif
  call dgetri(k,cbari,k,ipiv,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with dgetri'
    write(*,*) 'info = ',info
    return
  endif
  cbari = imat-matmul(cbari,b1)
  !
  ! Invert Cbar^(-1)
  !
  cbar = cbari
  call dgetrf(k,k,cbar,k,ipiv,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with dgetrf'
    write(*,*) 'info = ',info
    return
  endif
  call dgetri(k,cbar,k,ipiv,work,lwork,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATIONARY: problem with dgetri'
    write(*,*) 'info = ',info
    return
  endif
  !
  ! Cholesky decomposition of Cbar* Sig *Cbar'
  !
  csigc = matmul(cbar,sig)
  csigc = matmul(csigc,transpose(cbar))
  call dpotrf('L',k,csigc,k,info)
  if (info.ne.0) then
    write(*,*) 'ERROR in STATS: problem with dpotrf'
    write(*,*) 'info = ',info
    return
  endif
  do j=2,k
    do i=1,j-1
      csigc(i,j) = 0.
    enddo
  enddo
  !
  ! A0   = inv(Cbar)*chol(Cbar*Sig*Cbar')
  !
  a0      = matmul(cbari,csigc)
  l       = 1
  do i=1,k
    do j=1,k
      stat(l) = a0(i,j)
      l       = l+1
    enddo
  enddo

  !=========================
  ! Second set of statistics
  !=========================
  !
  ! Compute V that satisfies V = A*V*A' + B*Sig*B'
  !
  kk            = 2*k
  a             = 0.
  a(1:k,1:k)    = b1+m
  a(1:k,k+1:kk) = imat
  b(1:k,:)      = imat
  b(k+1:kk,:)   = -m

  tem1    = matmul(sig,transpose(b))
  v       = matmul(b,tem1)
  at      = -transpose(a)
  atem    = a

  call sb04qd(kk,kk,atem,kk,at,kk,v,kk,tem2,kk,iwork,work,lwork,info)
  l       = k*k+1
  do i=1,k
    stat(l) = v(i,i)
    l       = l+1
  enddo
 
  !
  ! Variance decomposition (V(i) due to each e(i))
  !    V(i) = A*V(i)*A' + B*A0*ei*A0'*B', ei has 1 in i,i element
  !

  l         = k*k+k+1
  do j=1,k
    ei      = 0.
    ei(j,j) = 1.
    tem     = matmul(a0,ei) 
    tem     = matmul(tem,transpose(a0)) 
    tem1    = matmul(tem,transpose(b))
    vj      = matmul(b,tem1)
    atem    = a
    at      = -transpose(a)
    call sb04qd(kk,kk,atem,kk,at,kk,vj,kk,tem2,kk,iwork,work,lwork,info)
    do i=1,k
      stat(l)  = vj(i,i)/v(i,i)
      l        = l+1
    enddo
  enddo
  csig    = sig
  call dpotrf('L',k,csig,k,info)
  do j=2,k
    do i=1,j-1
      csig(i,j) = 0.
    enddo
  enddo
  do j=1,k
    ei      = 0.
    ei(j,j) = 1.
    tem     = matmul(csig,ei)
    tem     = matmul(tem,transpose(csig))
    tem1    = matmul(tem,transpose(b))
    vj      = matmul(b,tem1)
    atem    = a
    at      = -transpose(a)
    call sb04qd(kk,kk,a,kk,at,kk,vj,kk,tem2,kk,iwork,work,lwork,info)
    do i=1,k
      stat(l)  = vj(i,i)/v(i,i)
      l        = l+1
    enddo
  enddo
      
  !=========================
  ! Third set of statistics
  !=========================
  !
  ! KP Statistics for the variables of interest (W)
  !
  !    Example:  Suppse observed vector is
  !       Y[t] = [ log (y[t]/l[t]) - log(y[t-1]/l[t-1]) ]
  !              [ log l[t] - alpha log l[t-1]          ]
  !
  !       W[t] = [ log y[t] ]
  !              [ log l[t] ]
  !
  !    Then alp = [1,alpha]'
  !         S   = [1  1]
  !               [0  1]
  !       nwant = number of rows of W = 2

  phi  = b1+m
 
  !
  ! Set Y[-1]=e[-1]=0
  !
  l  = 1
  do i=1,k
    elag(i,1) = randstat(l)
    l         = l+1
  enddo
  ylag = elag 
  do i=1,k
    e(i,1)    = randstat(l)
    l         = l+1
  enddo
  e           = matmul(csig,e)
  ylag        = matmul(phi,ylag)+e-matmul(m,elag)
  elag        = e
  !
  ! Generate nsim time series for W[t] and HP filter each
  !
  kp          = 0.
  x           = 0.
  do s=1,nsim
    do j=1,nper
      do i=1,k
        e(i,1) = randstat(l)
        l      = l+1
      enddo
      e       = matmul(csig,e)
      y       = matmul(phi,ylag)+e-matmul(m,elag)
      x       = alp*x -y
      wtem    = matmul(swant,x)
      w(:,j)  = wtem(:,1)
      ylag    = y
      elag    = e
    enddo
    do j=1,nwant
      wn      = w(j,:)
      call hptrend(nper,wn,1600., wt,wd)
      w(j,:)  = wd
    enddo

    !
    ! With HP-filtered data, compute second moments (as in Kydland-Prescott)
    !
    sum1      = 0.
    sum2      = 0.
    sum3      = 0.
    sum4      = 0.
    sum5      = 0.
    sum6      = 0.
    sum7      = 0.
    sum8      = 0.
    do j=1,nper
      sum1    = sum1 + w(:,j)
      sum2    = sum2 + w(:,j)*w(:,j)
    enddo
    do j=2,nper
      sum3    = sum3 + w(:,j)
      sum4    = sum4 + w(:,j-1)
      sum5    = sum5 + w(:,j)*w(:,j-1)
      sum6    = sum6 + w(:,j)*w(:,j)
      sum7    = sum7 + w(:,j-1)*w(:,j-1)
    enddo
    mean      = sum1/float(nper)
    meanf     = sum3/float(nper-1)
    meanl     = sum4/float(nper-1)
    var       = (sum2-float(nper)*mean*mean)/float(nper-1)
    std       = sqrt(var)
    acov      = (sum5-float(nper-1)*meanf*meanl)/float(nper-2)
    varf      = (sum6-float(nper-1)*meanf*meanf)/float(nper-2)
    varl      = (sum7-float(nper-1)*meanl*meanl)/float(nper-2)
    acor      = acov/(sqrt(varf)*sqrt(varl)) 
    do i=1,nwant
      do n=1,nwant
        do j=1,nper
          sum8((i-1)*nwant+n) = sum8((i-1)*nwant+n)+w(i,j)*w(n,j)
        enddo
        cor((i-1)*nwant+n)    = (sum8((i-1)*nwant+n)-float(nper)*mean(i)*    &
                                 mean(n))/(float(nper-1)*std(i)*std(n))
      enddo
    enddo
                 
    n                   = (nwant+2)*nwant
    kp(1:nwant)         = kp(1:nwant)+std
    kp(nwant+1:2*nwant) = kp(nwant+1:2*nwant)+acor
    kp(2*nwant+1:n)     = kp(2*nwant+1:n)+cor

  enddo
  stat((3*k+1)*k+1:lstat)  = kp/float(nsim)

end subroutine stats1