subroutine stderr11(n,x,ldat,dat,nboot, se)
!
! STDERR11  Computes standard errors of Phi, Theta, and Sig
! for the VARMA(1,1)
!
!     Y[t] = Phi(x) Y[t-1] + e[t] -Theta(x) e[t-1] 
!
! where Y[t] is a vector of length k, e[t] ~ N(0,Sig), and the
! sample size is for Y[t] is T, e.g., t=1,...T.
!
! Inputs:
!
!     n     = length of parameter vector (x)
!     x     = parameter vector
!     ldat  = length of dat
!     dat   = [T,k,Y(:)], where Y is T x k
!
! The vector x is
!
!     x     = A11,A12,A21,A22, B11,B12,B21,B22, C11,C21,C22
!
! where Mij is element (i,j) of matrix M and
!
!     Phi   = inverse[chol(I+AA')'] * A
!     Theta = inverse[chol(I+BB')'] * B
!     Sig   = CC'

!
! Ellen McGrattan, 5-25-06
!

  implicit none
  integer, parameter                 :: lwork=1000,tmax=1000
  integer, intent(in)                :: n,ldat,nboot
  integer, dimension(lwork)          :: iwork,ipiv
  integer                            :: i,j,k,l,m,t,iu,info
  real, dimension(n), intent(in)     :: x
  real, dimension(n), intent(out)    :: se
  real, dimension(ldat), intent(in)  :: dat
  real, dimension(ldat)              :: newdat
  real, dimension(lwork)             :: work
  real, dimension(n)                 :: xb,sumx,sumx2
  real, dimension(:,:), allocatable  :: ahat,bhat,cs,sige,phi,thet,y,ylag,e,  &
                                        elag,ehat
  real                               :: lnl,u
                                                                                
  external loglike11

  t         = dat(1)
  k         = dat(2)
  if (t>tmax) then
    write(*,*) 'ERROR in stderr11: increase tmax in stderr11.f90'
    se      = 0.
  else
    allocate(ahat(k,k))
    allocate(bhat(k,k))
    allocate(cs(k,k))
    allocate(sige(k,k))
    allocate(phi(k,k))
    allocate(thet(k,k))
    allocate(y(k,t))
    allocate(ylag(k,1))
    allocate(elag(k,1))
    allocate(e(k,1))
    allocate(ehat(k,t))
    do j=1,k
      y(:,j) = dat((j-1)*t+3:j*t+2)
    enddo
    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
    l    = 2*k*k+1
    do j=1,k
      do i=j,k
        cs(i,j) = x(l)
        l       = l+1
      enddo
    enddo
    sige    = matmul(cs,transpose(cs))
                                                                                                                    
    call stationary(k,ahat,sige, phi)
    call stationary(k,bhat,sige, thet)
    se          = 0.
    sumx        = 0.
    sumx2       = 0.
    ylag        = 0.
    elag        = 0.
    do j=1,t
      e(:,1)    = y(:,j)
      e         = e-matmul(phi,ylag)+matmul(thet,elag)
      ehat(:,j) = e(:,1)
      ylag(:,1) = y(:,j)
      elag      = e
    enddo
    do m=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(phi,ylag)+e-matmul(thet,elag)
        y(:,j)   = ylag(:,1)
        elag     = e
      enddo  
      do i=1,k
        newdat((i-1)*t+5:i*t+4) = y(i,:)
      enddo
      call uncmnd(n,x,ldat,newdat,loglike11,xb,lnl,info,work,lwork)

      l = 1
      do i=1,k
        do j=1,k
          ahat(i,j) = xb(l)
          l         = l+1
        enddo
      enddo
      do i=1,k
        do j=1,k
          bhat(i,j) = xb(l)
          l         = l+1
        enddo
      enddo
      l    = 2*k*k+1
      do j=1,k
        do i=j,k
          cs(i,j) = xb(l)
          l       = l+1
        enddo
      enddo
      sige    = matmul(cs,transpose(cs))

      call stationary(k,ahat,sige, phi)
      call stationary(k,bhat,sige, thet)
      l = 1
      do i=1,k
        do j=1,k
          xb(l)  = phi(i,j)
          l      = l+1
        enddo
      enddo
      do i=1,k
        do j=1,k
          xb(l)  = thet(i,j)
          l      = l+1
        enddo
      enddo
      sumx       = sumx+xb
      sumx2      = sumx2+xb*xb
    enddo
    se   = sqrt((sumx2-sumx*sumx/float(nboot))/float(nboot-1))
    deallocate(ahat)
    deallocate(bhat)
    deallocate(cs)
    deallocate(sige)
    deallocate(phi)
    deallocate(thet)
    deallocate(y)
    deallocate(ylag)
    deallocate(elag)
    deallocate(e)
    deallocate(ehat)
  endif

end subroutine stderr11