*! version 2.1.0 : 6-15-94 (c) J. Hilbe
* GLM: Log negative binomial without k estimation (Chi2 based convergence)
* expected IM initialization, observed thereafter: may be altered
* prior analytic/frequency weights and noconstant allowed; if/in allowed 
* defaults as geometric regression

program define negbink
  version 3.1
  local varlist "req ex"
  local options "`options' LTolera(real -1) ITerate(integer 0) Offset(string) EForm LEvel(integer $S_level) Resid K(real 1.0) NOLog noCONStant Expec(real 1)"
  local weight "fweight aweight"
  local in "opt"
  local if "opt"
  parse "`*'"
  parse "`varlist'",parse(" ")

  qui  {
  tempvar w wt mu dev n  z eta oldev ofs chi cdisp ocdisp alpha ofs llog touse
  local y "`1'"                      
  mac shift
  
  if `ltolera'<0  { local ltolera 1e-8  }
  if `iterate'<=0 { local iterate 50 }
   
  if `level'<10 | `level'>99  {
     local level 95
   }
  if "`constan'"!="" { local cons "nocons" }  
  
  if "`eform'"~="" { local eform "eform(IRR)" }
  
  if "`offset'"!="" { gen `ofs'=`offset'}
  else { gen `ofs' = 0 }
  
  gen `alpha' = `k'
  
  gen byte `touse'=1 `if' `in'
  replace `touse'=0 if `y'==. | `touse'==.
  if "`exp'"~="" {
     gen float `wt' `exp' if `touse'
     replace `wt'=. if `wt'<=0
     replace `touse'=0 if `wt'==.
     if ("`weight'"=="aweight") {
        sum `wt'
        replace `wt'=`wt'/_result(3)
     }
  }
  else gen float `wt'=`touse' if `touse'
  qui reg `y' `*' if `touse'
  mac list
  local nobs=_result(1)


* MAIN NEGATIVE BINOMIAL LOOP
  local ddisp = 1
  gen `ocdisp' = 1 if `touse'
  gen `dev'=1 if `touse' 
  gen `w' = 1 if `touse' 
      local i = 1
      local ddev = 1
      gen `z' = 0 if `touse' 
      replace `dev' = 0
      gen `oldev' = 1 if `touse' 
      egen `mu' = mean(`y') if `touse' 
      replace `mu' = (`y'+ `mu')/2
      gen `eta'= log(`mu') if `touse'         /* NB log link */
      while ((abs(`ddev') > `ltolera' & `i'<`iterate') | (`i'==0))  { 
          if `i'<=`expec' { 
             replace `w' = `mu'
             replace `z'=(`eta'+(`y'-`mu')/`w') -`ofs' 
             replace `w'=`w'^2/(`mu'+`alpha'*`mu'*`mu')
          }
          else {
            #delimit ;
            replace `w' =(`mu'/(1+`alpha'*`mu')) + (`y'-`mu')*
               ( (`alpha'*`mu')/(1+(2*`alpha'*`mu')+(`alpha'^2*`mu'^2)));
            #delimit cr
            replace `z'=(`eta'+(`y'-`mu')/( (1+`alpha'*`mu')*`w')) - `ofs'
          }
          regress `z' `*' [iw=`w'*`wt'] , mse1 dep(`y') `cons'  
          cap drop `eta' 
          predict `eta' if `touse' 
          replace `eta'=`eta'+`ofs'
          replace `mu'= exp(`eta')
          replace `oldev'=`dev'
          #delimit ;
          replace `dev' = `y'*log(`y'/`mu')-(1+`alpha'*`y')/`alpha' * 
                          log((1+`alpha'*`y')/(1+`alpha'*`mu')) if `y'>0;
          replace `dev' = log(1+`alpha'*`mu')/`alpha' if `y'==0;
          #delimit cr
          replace `dev' = sum(`dev'*`wt')
          replace `dev' = 2*`dev'[_N]
          if "`nolog'"=="" {
            noi di `i' in gr ":    Deviance: " in ye %9.4g  `dev'
          }
          local ddev = `dev'-`oldev'
          local i = `i'+1
       }
    gen `chi'=(`y'-`mu')^2/(`mu'+`alpha'*`mu'*`mu') if `touse' 
    replace `chi'=sum(`chi'*`wt')
    replace `chi'= `chi'[_N]
    local pred = _result(3)
    local df = `nobs'-`pred'-("`cons'"=="") /* degrees of freedom */
    if ("`weight'"=="fweight") {  
       tempvar fwt
       egen `fwt'= sum(`wt')
       local fobs =`fwt'
       local df = `fobs'-`pred'-1 
       if "`constan'"!="" { 
         local df = `fobs'-`pred'
       }
    } 
    gen `cdisp' = `chi'/`df'
  
gen `llog' = `y'*log(`alpha'*`mu')-(`y'+1/`alpha')*log(1+`alpha'*`mu')+lngamma(`y'+1/`alpha')-lngamma(`y'+1) - lngamma(1/`alpha') if `touse' 
replace `llog'=sum(`llog'*`wt')
replace `llog'=`llog'[_N]

noi di in gr _n "Alpha (k)  = " in ye %9.4g `alpha' in gr _col(55) "No obs.     = " in ye %9.0g `nobs' _n  
noi di in gr "Chi2       = " in ye %9.0g `cdisp'*`df' in gr _col(55) "Deviance    = " in ye %9.0g `dev'
noi di in gr "Prob>chi2  = " in ye %9.0g chiprob(`df',`cdisp'*`df') in gr _col(55) "Prob>chi2   = " in ye %9.0g chiprob(`df',`dev')
noi di in gr "Dispersion = " in ye %9.0g `cdisp' in gr _col(55) "Dispersion  = " in ye %9.0g `dev'/`df' _n
scalar S_al = `alpha'
noi di in gr  "Log-link negative binomial regression"
noi regress,  noheader `eform' level(`level')
noi di in gr "Loglikelihood = " `llog'

if "`resid'"!="" {
   tempvar Presid Dresid stdp nbv Hat Spear Sdev delta
   gen `delta'=1
   cap drop _mu _eta _Presid _Dresid _hat _Spear _Sdev
   gen `Presid' = (`y'-`mu')/sqrt(`mu'+`alpha'*`mu'* `mu') * sqrt(`wt') if `touse' 
   gen `Dresid' =  sqrt(2*(`y'*log(`y'/`mu')-(1+`alpha'*`y')/`alpha' * log((1+`alpha'*`y')/(1+`alpha'*`mu')))) if `y'>0 if `touse' 
   replace `Dresid' =  sqrt(2*(log(1+`alpha'*`mu')/`alpha')) if `y'==0
   replace `Dresid' = -`Dresid' if (`y'-`mu')<0 
   predict `stdp',stdp if `touse' 
   gen `nbv' = `mu'+`alpha'*`mu'*`mu' if `touse' 
   gen `Hat' = (`stdp'^2*`nbv')/`delta' if `touse' 
   gen `Spear' = (`y'-`mu')/sqrt(`nbv' * `delta'* (1-`Hat')) if `touse' 
   gen `Sdev' =  sqrt(2*(`y'*log(`y'/`mu')-(1+`alpha'*`y')/`alpha' * log((1+`alpha'*`y')/(1+`alpha'*`mu')))/(`delta'*(1-`Hat'))) if `y'>0 if `touse' 
   replace `Sdev' =  sqrt(2*(log(1+`alpha'*`mu')/`alpha')/(`delta'*(1-`Hat'))) if `y'==0
   replace `Sdev' = -`Sdev' if (`y'-`mu')<0 
   gen _mu = `mu'
   gen _eta = `eta'
   gen _Presid = `Presid'
   gen _Dresid = `Dresid'
   gen _hat    = `Hat'
   gen _Spear  = `Spear'
   gen _Sdev   = `Sdev'
   noi di in gr _n "Created variables: " in ye "_mu, _eta, _Presid, _Dresid, _hat, _Spear, _Sdev"
   }

}
end