*! version 1.0.0 : 9-28-93  J. Hilbe
* GLM: Poisson based Extreme Value (Gumbel) regression 

program define glmgumb
  version 3.1
  local varlist "req ex"
  local options "`options' LTolera(real -1) ITerate(integer 0) Cen(string) Time Model(string)"
  local weight "fweight aweight"
  parse "`*'"
  parse "`varlist'",parse(" ")

  qui  {

  tempvar mu eta z dev oldev cens logt alpha a1 devw oldevw ty z0 at beta t_adj tt
  
* EXTRACT TIME VARIABLE; T
  local t "`1'"                /* survival time variable */
  mac shift
  gen `ty' = `t'  
  gen `tt' = `t'
  sum `ty'
  gen `t_adj' = _result(6)
  replace `ty'= `ty'/`t_adj' 
  

* DEFINE STATUS OF CENSOR VARIABLE
  if ("`cen'"!="") { gen byte `cens'=`cen' }   /* censor variable */
     else { gen byte `cens' = 1 }
  count if `cens'==1
  local ncens=_result(1)          /* number of non-censored cases */
  count  
  local nobs=_result(1)

  if `ltolera'<0  { local ltolera 0.0001  }
  if `iterate'<=0 { local iterate 50 }
  
*INITIALIZATIONS
  local i  = 1
  gen `z'  = 0
  gen `dev' = 0
  gen `oldev' = 1
  gen `alpha' = 1
  gen `beta' = 1
  local ddev = 1
  local j = 1
  gen `oldevw' = 1  
  local ddevw = 1
  gen `devw' = 0
  gen `a1'=0
  gen `at'=0
  
  gen `z0'=sum(`cens') 
  replace `z0'=`z0'[_N]
  replace `z0'=2*`z0' *log(`t_adj') 

  noi di in gr _col(11) "Deviance" _col(26) "Shape"
  
* MAIN LOOP
  while ((abs(`ddevw') > `ltolera' & `j'<`iterate') | (`j'==0))  { 
     cap drop `mu'
     cap drop `eta'
     egen `mu' = mean(`cens')           /* Poisson initial values */
     gen `eta' = ln(`mu')               /* Poisson link */
     local i = 1                        /* Reset */
     local ddev = 1
     replace `dev' = 0
     
* POISSON LOOP
     while ((abs(`ddev') > `ltolera' & `i'<`iterate') | (`i'==0))  {
         replace `z' = (`eta'+(`cens'-`mu')/`mu') - (`ty'*`alpha') 
         regress `z' `*' [iw=`mu'], mse1 dof(100000) dep(`t')
         cap drop `eta'
         predict `eta' 
         replace `eta' = `eta' + (`ty'*`alpha')
         replace `mu' = exp(`eta')
         replace `oldev' = `dev'
         replace `dev' = (`cens'*log(`cens'/`mu') - (`cens'-`mu'))
         replace `dev' = sum(`dev')
         replace `dev' = 2*`dev'[_N]
         local ddev = `dev' - `oldev'
         local i = `i'+1
     }
         
* CALCULATION OF EXTREME VALUE DEVIANCE
         replace `devw' = `cens'*log(`alpha'*`mu') - `mu'
         replace `devw' = sum(`devw')
         replace `devw' = `z0'- 2*`devw'[_N]
         replace `beta'=`alpha'/`t_adj' 
         local ddevw = `devw' - `oldevw'
         if "`model'"=="e" { noi di in gr _n "Generalized Deviance" }
         noi di in gr "Iter `j' : " in ye %9.4f `devw' " : "  %9.4f `beta'  
         if `j'==1 { local edev = `devw'}
         local j = `j'+1
         replace `oldevw' = `devw'

* RECALCULATION OF SHAPE PARAMETER
         replace `at' = `alpha'
         replace `a1' = sum((`mu'-`cens')*`ty'/`ncens')
         replace `a1' = `a1'[_N]                    
         replace `a1' = 1/`a1'
         replace `alpha'=(`alpha'+`a1')/2
 }

* LOG EXPECTED TIME FORM OPTION   
   if ("`time'"!="") {
    if (`cens'==1)  { replace `z' =(-(`eta'+(`cens'-`mu')/`mu')+(`tt'*`beta'))/`beta' }  
         else { replace `z' = 1-`eta'/log(`beta') }
    reg `z' `*' [iw=`mu'*`beta'^2], mse1 dep(`t') dof(100000)
  }
  
}
* DISPLAY
    di in gr _col(46)  "Observations         = " in ye %9.0f `nobs'
    di in gr "Shape  = " in ye %9.4f `beta'    in gr _col(46)  "Censored observ.     = " in ye %9.0f  `nobs'-`ncens' 
    di in gr "Scale  = " in ye %9.4f 1/`beta'  in gr _col(46)  "D(0) Deviance        = " in ye %9.4f `edev'
    di in gr _col(46)  "Scaled Deviance      = " in ye %9.4f `devw'
    local devdiff = `edev'- `devw'
    di in gr _col(46)  "Prob > Chi2          = " in ye %9.4f chiprob(1,`devdiff') _n
di in gr  "Extreme Value (Gumbel) regression" 
noi regress, noheader 
end