* 0-TRUNCATED NEGATIVE BINOMIAL REGRESSION : Joseph Hilbe 17Feb1998
* OPTIONS: tolerance,log,iteration,offset, robust, cluster, score
program define trnbin0
 version 5.0
 local options "Level(integer $S_level) IRr"
 if "`*'"=="" | substr("`1'",1,1)=="," {
   if "$S_E_cmd"=="trnbin0" {
     error 301
   }
   parse "`*'"
   if `level'<10 | `level'> 99 {
     di in red "level() must be between 10 and 99"
     exit 198
   }
 }
 else {
   local varlist "req ex"
   #delimit ;
   local options "`options' LTolerance(real -1) noLOg ITerate
   OFfset(string) SCore(string) SCORE2 Robust CLuster(string) *";
   #delimit cr
   local in "opt"
   local if "opt"
   local weight "fweight aweight"
   parse "`*'"
   parse "`varlist'",parse(" ")
   if "`log'"!="" { local log "quietly" }
   global S_mloff "`offset'"

 tempvar touse mysamp
 tempname b f V
 mark `touse' `if' `in'
 markout `touse' `varlist' `offset'

  if "`weight'" != "" {
        tempvar wvar
        gen double `wvar' `exp'
  }
  else  local wvar   1

  if ("`weight'"=="aweight")  {
        qui sum `wvar'
        qui replace `wvar' = `wvar'/_result(3)
  }

  if "`offset'" != "" {
        local ovar "`offset'"
  }
  else  local ovar   0
  global S_cen "`censor'"
* ESTIMATION OF LL0
 quietly {
   qui poisson `1'
   mat coef=get(_b)
   eq `1': `1'
   eq lnalpha:
   ml begin
   ml function nbin0lf
   ml method lf
   ml model `b'= `1' lnalpha, from(coefs) depv(10)
   ml sample `mysamp' [`weight' `exp'] if `touse'
   `log' ml maximize `f' `V', `options'
   local lf0 = `f'
   drop `mysamp'
 }

* ESTIMATION OF FULL MODEL
   qui poisson `varlist'
   local pll = $S_E_ll
   mat coef=get(_b)
   eq `1': `varlist'
   eq lnalpha:
   ml begin
   ml function nbin0lf
   ml method lf
   ml model `b' = `1' lnalpha, from(coefs) depv(10)
   ml sample `mysamp' [`weight' `exp'] if `touse'
   `log' ml maximize `f' `V', `options'
   #delimit ;
   ml post nbin0, title("0-Truncated Negative Binomial Estimates")
   lf0(`lf0') pr2;
   #delimit cr
   local nbll = $S_E_ll
   global S_1 = 2*(`nbll' - `pll')

* ROBUST CALCULATIONS
   local y "`1'"
   mac shift
   global rhs `*'
   tempname b V
   tempvar xb e1 e2 muvar alpha amuvar digam1 digam2

   mat `b' = get(_b)
   mat `V' = get(VCE)
   predict double `xb', index

* SCORE VARIABLES
   qui gen `muvar' = exp(`xb' + `ovar')
   qui gen `alpha' = exp([lnalpha]_b[_con])
   qui gen `amuvar' = `muvar'*`alpha'
   #delimit ;
   qui gen `digam1' = ( lngamma( (`y'+1/`alpha') +.0001)-
   lngamma((`y'+1/`alpha')-.0001) )/.0002;
   qui gen `digam2' = ( lngamma(1/`alpha' +.0001)-
   lngamma(1/`alpha'-.0001) )/.0002;
   #delimit cr

*SCORES
   tempvar tr0 tra ialpha
   qui gen double `ialpha' = -1/`alpha' if `touse'
   #delimit ;
   qui gen double `tr0' = ( `muvar' * -(1+`amuvar')^(-(1+`alpha')/
   `alpha') ) / (-1+(1+`amuvar')^(`ialpha')) if `touse';
   qui gen double `tra' = - (1 +`amuvar')^ (`ialpha')  *
   ( 1/(`alpha'^2) * ln(1+`amuvar') - (`muvar'/ (`alpha'*
   (1+`amuvar')))) / ( 1-(1+`amuvar')^(`ialpha')) if `touse';
   qui gen double `e1' = (`y' - `muvar')/(1+ `amuvar') - `tr0'
   if `touse';
   qui gen double `e2' = `ialpha'* ( (`alpha' * (`muvar'-`y'))/
   (1+`amuvar')) - ln(1+`amuvar') + `digam1' - `digam2' -
   `tra' if `touse';
#delimit cr

*  SCORE FUNCTIONS TO FILE
   if "`score'"!="" {
     gen `score' = `e1'
   }

* CLUSTER OPTIONS
   if "`robust'"!="" & "`cluster'"=="" {
      _robust `e1' `e2' if `touse'
      qui test $rhs, min
      local ch2 = _result(6)
   }
   if "`cluster'"!="" {
      _robust `e1' `e2' if `touse', cluster(`cluster')
      qui test $rhs, min
      local ch2 = _result(6)
   }
}
* OUTPUT
if "`irr'"!="" {
      local eform "eform(IRR)"
   }
ml mlout nbin0, level(`level') `eform'
di in gr "   alpha    " in ye %9.0g exp([lnalpha]_b[_cons]) /*
   */ in gr "   [lnalpha]_cons = ln(alpha)"
di _col(21) in gr "(LR test against Poisson, chi2(1) = " in ye %9.0g /*
   */ $S_1 /*
   */ in gr " P = ", in ye %6.4f chiprob(1,$S_1) in gr ")"

end