      subroutine boottp(nobs,nomax,nvarb,nvmax,nvar,nboot,isboot,
     1 jsboot,xmat,btrue,sigtrue,tvect,pvect,nfail)
      implicit real*8 (a-h,o-z)
c
c Copyright (c) 1996, James G. MacKinnon
c This routine comes with absolutely no warranty!
c
c This routine computes bootstrap P values for LR and ES LM tests
c in tobit models. These are the tests that seem to perform best.
c
c Routines called: geny (ggnew2), tobit (toblf, ddnor),
c                  ormereg (olsqcm, chol2)
c Some of these routines are in the file tobit.f.
c
c Arguments (input):
c
c  nobs:    number of observatons (max 2000)
c  nomax:   maximum allowable number of observations (max 2000)
c  nvarb:   number of regressors under the alternative (max 20)
c  nvmax:   maximum allowable number of regressors (max 20)
c  nvar:    number of regressors under the null (max 20)
c  nboot:   number of bootstrap samples
c  isboot:  first seed for bootstrap sampling
c  jsboot:  second seed for bootstrap sampling
c  xmat:    matrix of regessors
c  btrue:   vector of parameter values to be used for the bootstrap
c  sigtrue: value of sigma to be used for the bootstrap
c  tvect:   2-vector of test statistics, LR then ES LM
c
c Arguments (output):
c
c  pvect:   2-vector of calculated P values
c  nfail:   number of bootstrap samples that had to be discarded
c
      real*8 xmat(nomax,nvmax), btrue(nvmax), y(2000)
      real*8 beta(20), covmat(21,21), best(21)
      real*8 tvect(2), pvect(2)
      integer iy(2000)
c
      nmin = nvarb + 1
      nvar1 = nvar + 1
      onboot = 1.d0/nboot
      lrgt = 0
      lmgt = 0
      do iboot=1,nboot
c
c generate the data under the null
c
  555   continue
        call geny(nobs,nomax,nvar,nvmax,isboot,jsboot,xmat,
     1   sigtrue,btrue,y,iy)
        tol = 1.d-7
c
c check number of ones
c
        isum = 0
        do i=1,nobs
          isum = isum + iy(i)
        end do
        if (isum.le.nmin) then
          nfail = nfail + 1
          go to 555
        end if
c
c estimate model under the null
c
        ncall = -1
        sigma = sigtrue
        do j=1,nvar
         beta(j) = btrue(j)
        end do
        call tobit(nobs,nomax,nvar,nvmax,y,iy,xmat,beta,covmat,
     1   sigma,floglr,tol,ncall)
c
c run LM test regression
c
        do j=nvar+1,nvarb
          beta(j) = 0.d0
        end do
        call ormereg(nobs,nomax,nvarb,nvmax,beta,sigma,best,
     1   sigest,ess,xmat,y)
        if (ess.gt.tvect(2)) lmgt = lmgt + 1
c
c estimate model under alternative
c
        do j=nvar1,nvarb
          beta(j) = 0.d0
        end do
        ncall = -1
        call tobit(nobs,nomax,nvarb,nvmax,y,iy,xmat,beta,covmat,
     1   sigma,floglu,tol,ncall)
        bootlr = 2.d0*(floglu - floglr)
        if (bootlr.gt.tvect(1)) lrgt = lrgt + 1
c
      end do
c
c end of bootstrap loop
c
      pvect(1) = lrgt*onboot
      pvect(2) = lmgt*onboot
      return
      end
      subroutine bttapp(nobs,nomax,nvarb,nvmax,nvar,nboot,isboot,
     1 jsboot,xmat,btrue,sigtrue,tvect,pvect,nround,nfail)
      implicit real*8 (a-h,o-z)
c
c Copyright (c) 1996, James G. MacKinnon
c This routine comes with absolutely no warranty!
c
c This routine computes approximate bootstrap P values for LR and
c ES LM tests c in tobit models. They are approximate because
c they are evaluated at estimates obtained from nround Newton
c steps. Note that nround could probably be lower if only the LR
c test were used.
c
c Routines called: geny (ggnew2), toblf (ddnor),
c                  ormereg (olsqcm, chol2)
c Some of these routines are in the file tobit.f.
c
c Arguments (input):
c
c  nobs:    number of observatons (max 2000)
c  nomax:   maximum allowable number of observations (max 2000)
c  nvarb:   number of regressors under the alternative (max 20)
c  nvmax:   maximum allowable number of regressors (max 20)
c  nvar:    number of regressors under the null (max 20)
c  nboot:   number of bootstrap samples
c  isboot:  first seed for bootstrap sampling
c  jsboot:  second seed for bootstrap sampling
c  xmat:    matrix of regessors
c  btrue:   vector of parameter values to be used for the bootstrap
c  sigtrue: value of sigma to be used for the bootstrap
c  tvect:   2-vector of test statistics, LR then ES LM
c
c Arguments (output):
c
c  pvect:   2-vector of calculated P values
c  nfail:   number of bootstrap samples that had to be discarded
c
      real*8 xmat(nomax,nvmax), btrue(nvmax), y(2000)
      real*8 best(21), tvect(2), gnextu(20), gnext(20)
      real*8 pvect(2), bnext(20), gest(20)
      real*8 covmat(21,21)
      integer iy(2000)
      nvar1 = nvar + 1
      nmin = nvarb + 1
      onboot = 1.d0/nboot
      lrgt = 0
      lmgt = 0
      do iboot=1,nboot
c
c generate the data
c
  555   continue
        call geny(nobs,nomax,nvar,nvmax,isboot,jsboot,xmat,
     1   sigtrue,btrue,y,iy)
c
c check number of ones
c
        isum = 0
        do i=1,nobs
          isum = isum + iy(i)
        end do
        if (isum.le.nmin) then
          nfail = nfail + 1
          go to 555
        end if
c
c run nround Hessian updates under the null
c
        dlnext = 1.d0/sigtrue
        do j=1,nvar
          gnext(j) = btrue(j)*dlnext
        end do
        signext = sigtrue
        do iround=1,nround
          iopg = 1
          call toblf(nobs,nomax,nvar,nvmax,y,iy,xmat,gnext,dlnext,
     1    gest,delest,covmat,floglr,essr,iopg)
          dlnext = dlnext + delest
          if (dlnext.lt.1.d-5) dlnext = 1.d-5
          do j=1,nvar
            gnext(j) = gnext(j) + gest(j)
          end do
        end do
c
c evaluate logl restricted
c
        iopg = 0
        call toblf(nobs,nomax,nvar,nvmax,y,iy,xmat,gnext,dlnext,
     1    gest,delest,covmat,floglr,essr,iopg)
        iopg = 1
c
c run nround Hessian updates under the alternative, starting at gnext
c
        do j=nvar1,nvarb
          gnext(j) = 0.d0
          bnext(j) = 0.d0
        end do
        do j=1,nvarb
          gnextu(j) = gnext(j)
        end do
        dlnextu = dlnext
        do iround=1,nround
          call toblf(nobs,nomax,nvarb,nvmax,y,iy,xmat,gnextu,dlnextu,
     1    gest,delest,covmat,floglu,essu,iopg)
          dlnextu = dlnextu + delest
          if (dlnextu.lt.1.d-5) dlnextu = 1.d-5
          do j=1,nvarb
            gnextu(j) = gnextu(j) + gest(j)
          end do
        end do
c
c evaluate logl unrestricted
c
        iopg = 1
        call toblf(nobs,nomax,nvarb,nvmax,y,iy,xmat,gnextu,dlnextu,
     1    gest,delest,covmat,floglu,essu,iopg)
        bootlr = 2.d0*(floglu - floglr)
        if (bootlr.gt.tvect(1)) lrgt = lrgt + 1
c
c run Orme regression to get LM test
c
        signext = 1.d0/dlnext
        do j=1,nvar
          bnext(j) = gnext(j)*signext
        end do
        call ormereg(nobs,nomax,nvarb,nvmax,bnext,signext,best,
     1     sigest,essu,xmat,y)
        if (essu.gt.tvect(2)) lmgt = lmgt + 1
c
c end of bootstrapping loop
c
      end do
c
c calculate bootstrap P values
c
      pvect(1) = lrgt*onboot
      pvect(2) = lmgt*onboot
      return
      end
      subroutine geny(nobs,nomax,nvar,nvmax,iseed,jseed,xmat,
     1 sigtrue,btrue,y,iy)
      implicit real*8 (a-h,o-z)
c
c Copyright (c) 1996, James G. MacKinnon
c This routine comes with absolutely no warranty!
c
c routine to generate data from a Tobit model
c
c Arguments (input):
c
c  nobs:    number of observatons (max 2000)
c  nomax:   maximum allowable number of observations (max 2000)
c  nvar:    number of regressors
c  nvmax:   maximum allowable number of regressors
c  iseed:   first seed
c  jseed:   second seed
c  xmat:    matrix of regessors
c  btrue:   vector of parameter values (slope coefficients)
c  sigtrue: value of sigma
c
c Argument (output):
c
c  y:       vector of observations on dependent variable
c  iy:      integer vector. iy(i) = 1 if y(i) > 0, iy(i) = 0 otherwise.
c
      real*8 xmat(nomax,nvmax), btrue(nvmax), y(nomax)
      real*8 uvect(2000)
      integer iy(nomax)
      call ggnew2(iseed,jseed,nobs,uvect)
c
      do i=1,nobs
        yz = uvect(i)*sigtrue
        do j=1,nvar
          yz = yz + btrue(j)*xmat(i,j)
        end do
        y(i) = yz
        iy(i) = 1
        if (yz.lt.0.d0) then
          y(i) = 0.d0
          iy(i) = 0
        end if
      end do
      return
      end
