      subroutine tobit(nobs,nomax,nvar,nvmax,y,iy,xmat,beta,covmat,
     1 sigma,flogl,tol,ncall)
      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 Tobit estimation routine that uses the Hessian. Although the
c arguments are beta and sigma, the model is parametrized as
c gamma = beta/sigma and delta = 1/sigma, and the covariance
c matrix that is returned is for those parameters. Use routine
c tobcov to obtain a covariance matrix for beta and sigma.
c
c Routines called: olsqcm (chol2), toblf (ddnor)
c
c Arguments (input):
c
c  nobs:   number of observatons
c  nomax:  maximum allowable number of observations
c  nvar:   number of regressors
c  nvmax:  maximum allowable number of regressors
c  y:      vector of observations on dependent variable
c  iy:     integer vector. iy(i) = 1 if y(i) > 0, iy(i) = 0 otherwise.
c  xmat:   matrix of regessors
c  tol:    tolerance for convergence. 1.d-7 is a good choice.
c  ncall:  if ncall < 0, use beta and sigma as starting values
c  beta:   starting values for slope coefficients (optional)
c  sigma:  starting value for error standard deviation (optional)
c
c Arguments (output):
c
c  beta:   vector of slope coefficient estimates
c  covmat: estimated covariance matrix of gamma and delta
c  sigma:  estimate of standard deviation of error term
c  flogl:  maximized value of loglikelihood function
c  ncall:  number of times the loglikelihood function was evaluated
c
c As the program is currently written, there can be no more than
c 20 regressors and 2000 observations.
c
      real*8 y(nomax), xmat(nomax,nvmax), beta(nvmax), gold(20)
      real*8 xpy(20), gest(20), gest2(20), gamma(20)
      real*8 covmat(nvmax+1,nvmax+1)
      integer iy(nomax)
c
c If ncall => 0, regress y on xmat using ols to get starting values.
c
      if (ncall.lt.0) then
        sigold = sigma
        deltold = 1.d0/sigold
        do j=1,nvar
          gold(j) = beta(j)*deltold
        end do
      else
        call olsqcm(nobs,nomax,nvar,nvmax,beta,xpy,ssr,covmat,xmat,y)
        sigold = sqrt(ssr/nobs)
        deltold = 1.d0/sigold
        do j=1,nvar
          gold(j) = beta(j)*deltold
        end do
      end if
c
c Starting values are now in gold and deltold.
c Find logl for starting values.
c
      ncall = 0
  100 continue
      iopg = 1
      ncall = ncall + 1
      call toblf(nobs,nomax,nvar,nvmax,y,iy,xmat,gold,deltold,
     1 gest,delest,covmat,floglo,ess,iopg)
  101 continue
      if (ess.lt.tol) then
        sigma = 1.d0/deltold
        do 2 j=1,nvar
        beta(j) = gold(j)*sigma
    2   continue
        flogl = floglo
        return
      end if
c
c Now start searching in direction given by OPG regression
c Try step length 1, then smaller sizes
c
      step = 1.d0
  110 continue
      do 20 j=1,nvar
      gamma(j) = gold(j) + step*gest(j)
   20 continue
      delta = deltold + step*delest
      if (delta.le.0.d0) go to 150
      iopg = 1
      ncall = ncall + 1
      call toblf(nobs,nomax,nvar,nvmax,y,iy,xmat,gamma,delta,
     1 gest2,deles2,covmat,flogl,ess,iopg)
      if (flogl.gt.floglo) then
         deltold = delta
         delest = deles2
         do 21 j=1,nvar
         gold(j) = gamma(j)
         gest(j) = gest2(j)
   21    continue
         floglo = flogl
         go to 101
      end if
  150 continue
      step = .5d0*step
      if (step.lt.0.01d0) then
        sigma = 1.d0/delta
        do 3 j=1,nvar
        beta(j) = gold(j)*sigma
    3   continue
        flogl = floglo
        write(8,*) 'Search failed in tobit'
        return
      end if
      go to 110
      end
      subroutine toblf(nobs,nomax,nvar,nvmax,y,iy,xmat,gamma,delta,
     1 gest,delest,covmat,flogl,ess,iopg)
c
c Copyright (c) 1996, James G. MacKinnon
c This routine comes with absolutely no warranty!
c
c Subroutine to evaluate tobit loglikelihood function and, optionally,
c take Newton step. Uses gamma-delta parametrization.
c
c Routines called: ddnor, chol2
c
c Arguments (input):
c
c  nobs:   number of observatons
c  nomax:  maximum allowable number of observations
c  nvar:   number of regressors
c  nvmax:  maximum allowable number of regressors
c  y:      vector of observations on dependent variable
c  iy:     integer vector. iy(i) = 1 if y(i) > 0, iy(i) = 0 otherwise.
c  xmat:   matrix of regessors
c  gamma:  vector of slope coefficients, = beta/sigma
c  delta:  1/sigma
c  iopg:   iopg = 0 ===> only loglikelihood evaluated
c          iopg = any other value ===> Newton step taken
c
c Arguments (output):
c
c  gest:   vector of slope coefficient estimates
c  delest: estimate of delta
c  covmat: estimated covariance matrix of gamma and delta
c  flogl:  maximized value of loglikelihood function
c  ess:    g'(H^{-1})g, where g = gradient and H = Hessian.
c
c Must have nobs <= 2000, nvar <= 20 because that is how arrays
c defined internally are dimensioned.
c
      implicit real*8 (a-h,o-z)
      real*8 y(nomax), xmat(nomax,nvmax), gamma(nvmax), gest(nvmax)
      real*8 fits(2000), bigphi(2000), resid(2000), smphi(2000)
      real*8 covmat(nvmax+1,nvmax+1), grad(21)
      integer iy(nomax)
      data const/.398942280401432678d0/, big/7.d0/
c
      nvar1 = nvar + 1
      flogl = 0.d0
      do 1 i=1,nobs
      fits(i) = 0.d0
      do 2 j=1,nvar
      fits(i) = fits(i) + xmat(i,j)*gamma(j)
    2 continue
      if(iy(i).eq.0) then
        arg = - fits(i)
        call ddnor(arg,bigphi(i))
        contr = log(bigphi(i))
      else
        resid(i) = delta*y(i) - fits(i)
        arg = resid(i)
        aarg = abs(arg)
        if (aarg.gt.big) arg = (arg/aarg)*big
        z = .5d0*arg**2
        smphi(i) = const*dexp(-z)
        contr = log(delta*smphi(i))
      end if
      flogl = flogl + contr
c
c flogl is value of loglikelihood function
c
    1 continue
c
      if (iopg.eq.0) return
      nvmax1 = nvmax + 1
      sigma = 1.d0/delta
      do 49 j=1,nvar1
      grad(j) = 0.d0
      do 49 k=j,nvar1
      covmat(j,k) = 0.d0
   49 continue
      do 50 i=1,nobs
      if (iy(i).eq.1) then
        grad(nvar1) = grad(nvar1) + sigma - resid(i)*y(i)
        do 51 j=1,nvar
        grad(j) = grad(j) + resid(i)*xmat(i,j)
        covmat(j,nvar1) = covmat(j,nvar1) + xmat(i,j)*y(i)
        do 51 k=j,nvar
        covmat(j,k) = covmat(j,k) -  xmat(i,j)*xmat(i,k)
   51   continue
        covmat(nvar1,nvar1) = covmat(nvar1,nvar1) - y(i)**2 - sigma**2
      else
        arg = -fits(i)
        aarg = abs(arg)
        if (aarg.gt.big) arg = (arg/aarg)*big
        z = .5d0*arg**2
        smphi(i) = const*dexp(-z)
        obphi = 1.d0/bigphi(i)
        rmills = smphi(i)*obphi
        do 52 j=1,nvar
        grad(j) = grad(j) - rmills*xmat(i,j)
        do 52 k=j,nvar
        crap = (smphi(i)*fits(i)*bigphi(i)
     1   - smphi(i)**2)*xmat(i,j)*xmat(i,k)
        covmat(j,k) = covmat(j,k) + crap*(obphi**2)
   52   continue
      end if
   50 continue
      do 60 j=1,nvar1
      do 60 k=j,nvar1
      covmat(j,k) = - covmat(j,k)
      covmat(k,j) = covmat(j,k)
   60 continue
      call chol2(covmat,nvmax1,nvar1,kf2)
      do 61 j=1,nvar
      gest(j) = 0.d0
      do 61 k=1,nvar1
      gest(j) = gest(j) + grad(k)*covmat(k,j)
   61 continue
      delest = 0.d0
      do 62 k=1,nvar1
      delest = delest + grad(k)*covmat(k,nvar1)
   62 continue
      ess = delest*grad(nvar1)
      do 63 j=1,nvar
      ess = ess + gest(j)*grad(j)
   63 continue
      return
      end
      subroutine ddnor(ystar,gauss)
c
c Copyright (c) 1993, James G. MacKinnon
c This routine come with absolutely no warranty!
c
c This routine uses Cody's method to evaluate the cumulative normal
c distribution. It should be accurate to almost full double precision,
c up to 19 or 20 significant digits. It was originally written by
c James MacKinnon late in 1977, based on the Cody article referred to
c in the documentation for IMSL subroutine mdnor.
c
c Modified 1993 to avoid changing the argument.
c
c Arguments (input):
c
c  ystar:  argument of the normal c.d.f.
c
c Arguments (output);
c
c  gauss:  value of the normal c.d.f. evaluated at ystar.
c
      implicit real*8(a-h,o-z)
      real*8 p(6), q(5), a(9), b(8), c(5), d(4)
      data p(1)/-6.58749161529837803157d-04/,
     1     p(2)/-1.60837851487422766278d-02/,
     2     p(3)/-1.25781726111229246204d-01/,
     3     p(4)/-3.60344899949804439429d-01/,
     4     p(5)/-3.05326634961232344035d-01/,
     5     p(6)/-1.63153871373020978498d-02/
      data q(1)/2.33520497626869185443d-03/,
     1     q(2)/6.05183413124413191178d-02/,
     2     q(3)/5.27905102951428412248d-01/,
     3     q(4)/1.87295284992346047209d00/,
     4     q(5)/2.56852019228982242072d00/
      data a(1)/1.23033935479799725272d03/,
     1     a(2)/2.05107837782607146532d03/,
     2     a(3)/1.71204761263407058314d03/,
     3     a(4)/8.81952221241769090411d02/,
     4     a(5)/2.98635138197400131132d02/,
     5     a(6)/6.61191906371416294775d01/,
     6     a(7)/8.88314979438837594118d00/,
     7     a(8)/5.64188496988670089180d-01/,
     8     a(9)/2.15311535474403846343d-08/
      data b(1)/1.23033935480374942043d03/,
     1     b(2)/3.43936767414372163696d03/,
     2     b(3)/4.36261909014324715820d03/,
     3     b(4)/3.29079923573345962678d03/,
     4     b(5)/1.62138957456669018874d03/,
     5     b(6)/5.37181101862009857509d02/,
     6     b(7)/1.17693950891312499305d02/,
     7     b(8)/1.57449261107098347253d01/
      data c(1)/3.209377589138469472562d03/,
     1     c(2)/3.774852376853020208137d02/,
     2     c(3)/1.138641541510501556495d02/,
     3     c(4)/3.161123743870565596947d00/,
     4     c(5)/1.857777061846031526730d-01/
      data d(1)/2.844236833439170622273d03/,
     1     d(2)/1.282616526077372275645d03/,
     2     d(3)/2.440246379344441733056d02/,
     3     d(4)/2.360129095234412093499d01/
      data orpi/.5641895835477562869483d0/,
     1   root2/.70710678118654752440083d0/
      isw = 1
      y = ystar
      if (ystar.lt.-16.d0) y = -16.d0
      if (ystar.gt.16.d0) y = 16.d0
      x = -y*root2
      if(x.gt.0.d0) go to 1
      if(x.lt.0.d0) go to 2
      gauss = .5d0
      return
    2 continue
      x = - x
      isw = -1
    1 continue
      if(x.lt..477d0) go to 10
      if(x.le.4.d0) go to 20
c
c  evaluate erfc for x.gt.4.0
c
      x2 = x*x
      xm2 = 1.d0/x2
      xm4 = xm2*xm2
      xm6 = xm4*xm2
      xm8 = xm4*xm4
      xm10 = xm6*xm4
      top = p(1) + p(2)*xm2 + p(3)*xm4 + p(4)*xm6 + p(5)*xm8 + p(6)*xm10
      bot = q(1) + q(2)*xm2 + q(3)*xm4 + q(4)*xm6 + q(5)*xm8 + xm10
      crap = orpi + top/(bot*x2)
      erfc = dexp(-x2)*crap/x
c
      if(isw.eq.-1) erfc = 2.d0 - erfc
      gauss = erfc*.5d0
      return
   20 continue
c
c  evaluate erfc for .477.lt.x.le.4.0
c
      x2 = x*x
      x3 = x2*x
      x4 = x2*x2
      x5 = x3*x2
      x6 = x3*x3
      x7 = x3*x4
      x8 = x4*x4
      top = a(1) + a(2)*x + a(3)*x2 + a(4)*x3 + a(5)*x4 + a(6)*x5 +
     1 a(7)*x6 + a(8)*x7 + a(9)*x8
      bot = b(1) + b(2)*x + b(3)*x2 + b(4)*x3 + b(5)*x4 + b(6)*x5 +
     1 b(7)*x6 + b(8)*x7 + x8
      erfc = dexp(-x2)*top/bot
c
      if(isw.eq.-1) erfc = 2.d0 - erfc
      gauss = erfc*.5d0
      return
   10 continue
c
c  evaluate erf for x.lt..477
c
      x2 = x*x
      x4 = x2*x2
      x6 = x4*x2
      x8 = x4*x4
      top = c(1) + c(2)*x2 + c(3)*x4 + c(4)*x6 + c(5)*x8
      bot = d(1) + d(2)*x2 + d(3)*x4 + d(4)*x6 + x8
      erf = x*top/bot
c
      erf = erf*isw
      erfc = 1.d0 - erf
      gauss = erfc*.5d0
      return
      end
      subroutine olsqcm(nobs,nomax,nvar,nvmax,beta,xpy,ssr,
     1 xpxi,xmat,yvect)
      implicit real*8 (a-h,o-z)
c
c Copyright (c) 1996, James G. MacKinnon
c
c This is a fast, but not always accurate, least squares routine
c based on the Cholesky decomposition. It should not be used if the
c data are ill-conditioned. Residuals are not returned.
c
c Arguments (input):
c
c  nobs:   number of observatons
c  nomax:  maximum allowable number of observations
c  nvar:   number of regressors
c  nvmax:  maximum allowable number of regressors
c  yvect:  vector of observations on dependent variable
c  xmat:   matrix of regessors
c
c Arguments (output):
c
c  beta:  vector of coefficient estimates
c  xpy:   X'y vector
c  xpxi:  (X'X)^{-1} matrix
c  ssr:   sum of squared residuals
c
      real*8 xmat(nomax,nvmax), yvect(nomax), beta(nvmax)
      real*8 xpy(nvmax), xpxi(nvmax,nvmax)
c
c  form and invert xpx matrix
c
      do i=1,nvar
        do j=i,nvar
          xpxi(i,j) = 0.d0
        end do
      end do
c
c only one triangle of xpxi computed.
c
      do i=1,nvar
        do j=i,nvar
          do k=1,nobs
            xpxi(i,j) = xpxi(i,j) + xmat(k,i)*xmat(k,j)
          end do
          xpxi(j,i) = xpxi(i,j)
        end do
      end do
c
c  now invert xpx
c
      call chol2(xpxi,nvmax,nvar,kf2)
c
c preliminary work has now been done.
c
c form xpy
c
      do i=1,nvar
        xpy(i) = 0.d0
      end do
      do i=1,nvar
        do k=1,nobs
          xpy(i) = xpy(i) + xmat(k,i)*yvect(k)
        end do
      end do
c
c  now form estimates of beta.
c
      do i=1,nvar
        beta(i) = 0.d0
        do j=1,nvar
          beta(i) = beta(i) + xpxi(i,j)*xpy(j)
        end do
      end do
c
c Compute ssr.
c
      ssr = 0.d0
      do k=1,nobs
        resid = yvect(k)
        do i=1,nvar
          resid = resid - xmat(k,i)*beta(i)
        end do
        ssr = ssr + resid**2
      end do
      return
      end
      subroutine chol2(a,m,n,kf2)
c
c This routine uses the cholesky decomposition to invert a real
c symmetric matrix.
c
c Modified 1996 to reduce the number of divides.
c
      implicit real*8 (a-h,o-z)
      real*8 a(m,m)
      kf2 = 0
      do 8 i=1,n
      kl = i - 1
      do 7 j=i,n
      if (i.gt.1) then
        do 3 k=1,kl
    3   a(i,j) = a(i,j) - a(k,i)*a(k,j)
      else
        if (a(i,i).le.0.d0) then
        kf2 = i
        go to 20
      end if
      end if
      if (i.eq.j) then
        a(i,i) = dsqrt(a(i,i))
      else
        if (j.eq.i+1) ooa = 1.d0/a(i,i)
        a(i,j) = a(i,j)*ooa
      end if
    7 continue
    8 continue
      do 13 i=1,n
      do 12 j=i,n
      ooa = 1.d0/a(j,j)
      if (i.ge.j) then
        t = 1.d0
        go to 12
      end if
      kl = j - 1
      t = 0.d0
      do 11 k=i,kl
   11 t = t - a(i,k)*a(k,j)
   12 a(i,j) = t*ooa
   13 continue
      do 16 i=1,n
      do 15 j=i,n
      t = 0.d0
      do 14 k=j,n
   14 t = t + a(i,k)*a(j,k)
      a(i,j) = t
   19 a(j,i) = t
   15 continue
   16 continue
   20 return
      end
      subroutine tobcov(nvar,nvmax1,covmat,bhat,sigest,gmat,altcov)
      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 transforms a covariance matrix in the alternative
c (gamma-delta) parametrization into a covarinnce matrix in the natural
c (beta-sigma) one.
c
c Arguments (input):
c
c  nvar:    number of betas
c  nvmax1:  maximum nunber of betas, plus 1 (21 if nvmax = 20)
c  covmat:  covariance matrix in gamma-delta parametrization.
c  bhat:    vector of coefficient estimates
c  sigest:  estimate of sigma
c
c Arguments (output):
c
c  gmat:    matrix of derivatives used by routine
c  altcov:  covariance matrix in beta-sigma parametrization
c
      real*8 covmat(nvmax1,nvmax1), altcov(nvmax1,nvmax1)
      real*8 bhat(nvmax1), gmat(nvmax1,nvmax1)
      nvar1 = nvar + 1
c
c Form matrix of derivatives
c
      do i=1,nvar
        do j=1,nvar
          gmat(i,j) = 0.d0
        end do
        gmat(i,i) = sigest
      end do
      do j=1,nvar
        gmat(nvar1,j) = 0.d0
        gmat(j,nvar1) = -bhat(j)*sigest
      end do
      gmat(nvar1,nvar1) = - sigest**2
c
c Multiply things out to get new covariance matrix
c
      do i=1,nvar1
      do l=1,nvar1
        sum = 0.d0
        do j=1,nvar1
        do k=1,nvar1
          sum = sum + gmat(i,j)*gmat(l,k)*covmat(j,k)
        end do
        end do
        altcov(i,l) = sum
      end do
      end do
      return
      end
      subroutine ormereg(nobs,nomax,nvar,nvmax,bstar,sigstar,best,
     1 sigest,ess,xmat,y)
      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 run Orme's artificial regression for the Tobit model
c
c Routines called: olsqcm (chol2), ddnor
c
c Arguments (input):
c
c  nobs:    number of observatons
c  nomax:   maximum allowable number of observations
c  nvar:    number of regressors
c  nvmax:   maximum allowable number of regressors
c  xmat:    matrix of regessors
c  y:       vector of observations on dependent variable
c  bstar:   vector of coefficients (beta) to use for regression
c  sigstar: value of sigma to use for regression
c
c Arguments (output):
c
c  best:   vector of estimates from artificial regression
c  sigest: best(nvar1)
c  ess:    ESS from artificial regression
c
c Must have nobs <= 2000, nvar <= 20.
c
      real*8 xmat(nomax,nvmax), y(nomax), bstar(nvmax), best(21)
      real*8 rdep(4000), bigr(4000,21), xpxi(21,21), xpy(21)
      data const/.398942280401432678d0/
      osigma = 1.d0/sigstar
      osigsq = osigma**2
      nvar1 = nvar + 1
      nobs2 = nobs*2
c
c Generate fitted values, normalized by sigstar.
c Then generate normal density and distribution, and ratio, and
c everything else needed for artificial regression
c
      tss = 0.d0
      do i=1,nobs
        sum = 0.d0
        do j=1,nvar
          sum = sum + bstar(j)*xmat(i,j)
        end do
        fits = osigma*sum
        if (fits.lt.-7.d0) fits = -7.d0
        arg = -fits
        call ddnor(arg,bmphi)
        sphi = const * exp(-.5d0*fits**2)
        ratio = sphi/bmphi
        omee = osigsq*(1.d0 - bmphi*(1.d0 - ratio*(ratio - fits)))
        omgg = osigsq*(2.d0 - bmphi*(2.d0 + fits*ratio*
     1   (1.d0 - fits*(ratio - fits))))
        omeg = osigsq*sphi*(1.d0 - fits*(ratio - fits))
        ac = sqrt(omgg)
        ab = omeg/ac
        aa = (1.d0/ac)*sqrt(omee*omgg - omeg**2)
        if (y(i).gt.0.d0) then
           cc = 0.d0
        else
           cc = 1.d0
        end if
        uu = osigma*y(i) - fits
        eps = osigma*((1.d0 - cc)*uu - cc*ratio)
        gamma = osigma*((1.d0 - cc)*
     1   (uu**2 - 1.d0) + cc*fits*ratio)
c
c generate regressand
c
        rdep(i) = eps/aa - ab*gamma/(aa*ac)
        rdep(i+nobs) = gamma/ac
        tss = tss + rdep(i)**2 + rdep(i+nobs)**2
c
c generate regressors
c
        do j=1,nvar
          bigr(i,j) = xmat(i,j)*aa
          bigr(i+nobs,j) = xmat(i,j)*ab
        end do
        bigr(i,nvar1) = 0.d0
        bigr(i+nobs,nvar1) = ac
      end do
c
c run artificial regression
c
      call olsqcm(nobs2,4000,nvar1,21,best,xpy,ssr,xpxi,bigr,rdep)
      sigest = best(nvar1)
      ess = tss - ssr
      return
      end
