*!version 1.0 : 12-16-96  Joseph Hilbe    STB-41 sg53.2
* Complementary log-log regression; like Stata logit/probit
program define cloglog
  version 5.0
  set type double                              
  local options "Level(integer $S_level) EForm" 
  if "`1'"!="" & substr("`1'",1,1)!="," {
  local varlist "req ex"
  #delimit ;
  local options "`options' LTolera(real -1) ITerate(integer 0) Offset(string)

     noCONStant NOLog CLuster(string) Robust SCore(string) *";
  #delimit cr
  local weight "fweight aweight"
  local in "opt"
  local if "opt"
  parse "`*'"
  parse "`varlist'",parse(" ")
  qui  {
    tempvar w wt mu eta touse u ll n z ofs cdisp chi g1 g11  
    tempvar v wo ll0 dll oll count e
    local y "`1'"
  mac shift
  global rhs `*'
  count
  local nobs = _result(1)
  if `ltolera'<0  { local ltolera 1e-6 }
  if `iterate'<=0 { local iterate 50 }
  if `level'<10 | `level'>99  { local level 95 }
  if "`constan'"  !=""  { local cons "nocons" }
  if "`eform'"    !=""  { local eform "eform(EF)" }
  if "`offset'"   !=""  { gen `ofs' = `offset' }
     else { gen `ofs' = 0 }
  gen byte `touse'=1 `in' `if'
  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'
  markout `touse' `varlist'
  if "`cluster'" ~= "" {
    confirm variable `cluster'
    markout `touse' `cluster', strok
  }
  qui reg `y' `*' if `touse'
  mac list
  local nobs=_result(1)

  replace `ofs'=`ofs' if `touse'
  gen `w' = 1 if `touse'
  gen `ll' = 0 if `touse'
  gen `ll0'= 0 if `touse'

* BINOMIAL INITIALIZATION
  local i = 1
  local dll = 1
  gen `z'    = 0 if `touse'
  gen `u'    = 0 if `touse'
  gen `g1'   = 0 if `touse'
  gen `g11'  = 0 if `touse'
  gen `v'    = 0 if `touse'
  gen `wo'   = 0 if `touse'
  gen `oll'  = 1 if `touse'  
  
* MODULE TO CALCULATE LL0
  gen `mu' =(`y'+0.5)/2 if `touse'    
  gen `eta' = ln(-ln(1-`mu')) if `touse'  
  while ((abs(`dll')>`ltolera' & `i'<`iterate') | (`i'==0))  {
    replace `u' = -(`y'-`mu')/((1-`mu')*ln((1-`mu')))
    replace `w'= ((1-`mu')*(ln((1-`mu'))))^2 /(`mu'*(1-`mu'))
    replace `g1' = -1/((1-`mu')*ln((1-`mu')))
    replace `g11' = (-1-ln((1-`mu')))/((`mu'-1)^2* /* 
       */ (ln((1-`mu')))^2)
    replace `v' = `mu'*((1-`mu'))
    replace `wo' = `w'+(`y'-`mu')*(`v'*`g11'+((1-2*`mu'))*`g1')/ /* 
       */ ( (`v'^2 * `g1'^3))
    replace `z' = `eta' + ((1/`wo')*`w'*`u') -`ofs'
    regress `z' [iw=`wo'*`wt'], mse1 dep(`y') `cons' dof(100000) 
    cap drop `eta'
    predict `eta' 
    replace `eta' = `eta' + `ofs'
    replace `mu' = 1-exp(-exp(`eta')) 
    replace `oll' = `ll0'
    replace `ll0'  = `y'*ln(`mu')+(1-`y')*ln(1-(`mu'))
    replace `ll0'  = sum(`ll0'*`wt')
    replace `ll0'  = `ll0'[_N]
    local dll = `ll0' - `oll'
    local i = `i' + 1
}

* MAIN MODULE
  replace `mu' =(`y'+0.5)/2 if `touse'    
  replace `eta' = ln(-ln(1-`mu')) if `touse'  
  local i=1
  local dll=1
  while ((abs(`dll')>`ltolera' & `i'<`iterate') | (`i'==0))  {
    replace `u' = -(`y'-`mu')/((1-`mu')*ln((1-`mu')))
    replace `w'= ((1-`mu')*(ln((1-`mu'))))^2 /(`mu'*(1-`mu'))
    replace `g1' = -1/((1-`mu')*ln((1-`mu')))
    replace `g11' = (-1-ln((1-`mu')))/((`mu'-1)^2* /* 
       */ (ln((1-`mu')))^2)
    replace `v' = `mu'*((1-`mu'))
    replace `wo' = `w'+(`y'-`mu')*(`v'*`g11'+((1-2*`mu'))*`g1')/ /* 
        */ ( (`v'^2 * `g1'^3))
    replace `z' = `eta' + ((1/`wo')*`w'*`u') -`ofs'
    regress `z' `*' [iw=`wo'*`wt'], mse1 dep(`y') `cons' dof(100000) 
    cap drop `eta'
    predict `eta' 
    replace `eta' = `eta' + `ofs'
    replace `mu' = 1-exp(-exp(`eta')) 
    replace `oll' = `ll'
    replace `ll'  = `y'*ln(`mu')+(1-`y')*ln(1-(`mu'))
    replace `ll'  = sum(`ll'*`wt')
    replace `ll'  = `ll'[_N]
    if "`nolog'"==""  {
       noi di in gr "Iteration  " `i'-1  in gr ": Log Likelihood = " /* 
         */ in ye %12.6f `ll'
    }
    local dll = `ll' - `oll'
    local i = `i' + 1
}
gen `chi' = (`y'-`mu')^2 / (`mu'*(1-`mu')) if `touse'
replace `chi' = sum(`chi'*`wt')
replace `chi' = `chi'[_N]

local pred = _result(3)
local df = `nobs'-`pred'-("`cons'"=="")
if ("`weight'" == "fweight")  {
   tempvar fwt
   egen `fwt' = sum(`wt')
   local fobs = `fwt'
   local df = `fobs' - `pred'-1
   if "`constan'"!=""  {
      local df = `fobs' - `pred'
   }
 }
local df = `nobs'-`pred'-("`cons'"=="")
gen `cdisp' = `chi'/`df' if `touse'
local pr2 = 1-(`ll'[_N]/`ll0'[_N])
local ch2 = -2*(`ll0'[_N]-`ll'[_N])
local ll = `ll'[_N]
local ll0 = `ll0'[_N]
}

qui regress,  noheader `eform' level(`level')  

* GET MAT AND SCORE (e) STATS
tempname b V
mat `b' = get(_b)
mat `V' = get(VCE)
local n = _result(1)
local k = colsof(`V')
tempvar xb
predict double `xb', index
qui replace `xb'=`xb'
qui gen double `e' =( exp(`xb') * exp(-exp(`xb')) / (1-exp(-exp(`xb'))) ) /*
     */  if `y'!=0 & `touse'
qui replace `e' =  -exp(`xb') if `y'==0 & `touse'
if "`score'"!="" {
   gen `score' = `e'
}

* CLUSTER
if "`cluster'"!="" {
  qui {
    sort `touse' `cluster'
    by `touse' `cluster': gen byte `count' = 1 if _n==1 & `touse'
    summ `count', meanonly
    local nclust=_result(18)
    local dof = `nclust'-1
    local clopt "cluster(`cluster')"
    }
}
else local dof = `n' - `k'

* CLUSTER
if "`cluster'"!="" {
   mat post `b' `V', depn(`y')    /* dof(`dof') */
   _robust `e' if `touse', minus(1) `clopt'
   qui test $rhs, min
   local ch2 = _result(6)
}

* ROBUST
else if "`cluster'"=="" & "`robust'"!=""  {
   mat post `b' `V', depn(`y')    /*  dof(`dof')  */
   _robust `e' if `touse', minus(1)
   qui test $rhs, min
   local ch2 = _result(6)
}

* STANDARD
else {
    mat post `b' `V', depname(`y') obs(`nobs') 
}
* END ROBUST
global S_E_depv "`y'"
global S_E_ll `ll'
global S_E_ll0 `ll0'
global S_E_pr2 `pr2'
global S_E_chi2 `ch2'
global S_E_nobs `nobs'
global S_E_pred `pred'
global S_E_mdf  `dof'
global S_E_cmd "cloglog"
}

else  {
     if "$S_E_cmd"!="cloglog" {
          error 301
     }      
     parse "`*'"
}
if `level'<10 | `level'>99 {
    local level 95
}
if "`eform'" !=""  { local eform "eform(EF)" } 

noi di in gr _n "Complementary log-log Estimates"  /*
   */ _col(56) "Number of obs = " in ye %7.0g $S_E_nobs
noi di in gr _col(56) "Chi2(" in ye $S_E_pred in gr ") " _col(69) " = " /* 
   */ in ye %7.2f $S_E_chi2
noi di in gr _col(56) "Prob>Chi2     = " /* 
   */ in ye %7.4f chiprob($S_E_pred,$S_E_chi2) 
noi di in gr "Log Likelihood = " %9.3g /* 
   */ in ye $S_E_ll in gr _col(56) "Pseudo R2     = " in ye %7.4f $S_E_pr2 _n


mat mlout, level(`level') `eform'

set type float
end