*!version 1.2.1 : 2-3-95  Joseph Hilbe -  revised:  6-19-96
* Maximum-likelihood probit regression, binomial/bernoulli : OIM-based
program define eprobit
  version 4.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 Group RESidual Expec(real 2) DISP /*
     */ SCale(string)";
  #delimit cr
  local weight "fweight aweight"
  local in "opt"
  local if "opt"
  parse "`*'"
  parse "`varlist'",parse(" ")
  qui  {
    tempvar w wt mu eta m ym touse u ll dev oldev n z ofs cdisp chi g1 g2
    tempvar v v1 v2 w1 ws g2 ndf wo ll0
  local y "`1'"
  mac shift
  if "`group'" != "" {
    local m = "`1'"
    mac shift
  }
  else { local m = 1 }
  global S_E_cmd
  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'
  qui reg `y' `*' if `touse'
  mac list
  local nobs=_result(1)
 
  replace `ofs'=`ofs' if `touse'
  gen `ym'=`y'/`m' if `touse'
  gen `w' = 1 if `touse'
  gen `ll' = 0 if `touse'
  gen `ll0'= 0 if `touse'
 
* BINOMIAL INITIALIZATION
  local i = 1
  local ddev = 1
  gen `z'    = 0 if `touse'
  gen `dev'  = 0 if `touse'
  gen `oldev'= 1 if `touse'
  gen `u'    = 0 if `touse'
  gen `g1'   = 0 if `touse'
  gen `g2'   = 0 if `touse'
  gen `v'    = 0 if `touse'
  gen `v1'   = 0 if `touse'
  gen `v2'   = 0 if `touse'
  gen `w1'   = 0 if `touse'
  gen `ws'   = 0 if `touse'
  gen `wo'   = 0 if `touse'
  gen `ndf'  = 0 if `touse'
 
 
* MODULE TO CALCULATE LL0
  gen `mu' = (`y'+0.5)/(`m'+1) if `touse'
  gen `eta' = invnorm(`mu') if `touse'
  while ((abs(`ddev') > `ltolera' & `i'<`iterate') | (`i'==0))  {
         if `i'<=`expec'  {
             replace `w'= 1/sqrt(2*_pi)*exp(-(`eta'^2/2))
             replace `z' = (`eta' + (`ym'-`mu')/`w') -`ofs'
             replace `w'=`w'^2/(`mu'*(1-`mu'))
             regress `z' [iw=`w'*`m'*`wt'], mse1  dep(`y') `cons'
         }
         else {
           replace `ndf'= 1/sqrt(2*_pi)*exp(-(`eta'^2/2))
           replace `v' = `mu'*(1-`mu')
           replace `v1' = 1-2*`mu'
           replace `v2'=(`mu'-`mu'^2)^2
           replace `w' = `ndf'^2/`v'
           replace `u' = (`ym'-`mu')/`ndf'
           replace `g1' = 1/`ndf'
           replace  `g2' = `eta'/`ndf'^2
           replace `ws' = `w'+(`ym'-`mu')*(`v'*`g2'+`v1'*`g1')/(`v2'*`g1'^3)
           replace `z' = (`eta'+1/`ws'*`w'*`u') - `ofs'
           regress `z'  [iw=`ws'*`m'*`wt'], mse1  dep(`y') `cons'
 
         }
         cap drop `eta'
         predict `eta' if `touse'
         replace `eta' = `eta' + `ofs'
         replace `mu' = normprob(`eta')
         replace `oldev' = `dev'
         #delimit ;
         replace `dev' = (`m'-`y')*log((`m'-`y')/(`m'-(`m'*`mu')))
                         if (`y'==0);
         replace `dev' = `y'*log(`y'/(`m'*`mu')) if (`y'==`m');
         replace `dev' = `y'*log(`y'/(`m'*`mu')) +
                          (`m'-`y')*log((`m'-`y')/(`m'-(`m'*`mu')))
                          if (`y'~=0 & `y'~=`m');
         #delimit cr
         replace `dev' = sum(`dev'*`wt')
         replace `dev' = 2*`dev'[_N]
         replace `ll0'=`y'*ln(`mu'/(1-`mu'))+`m'*ln(1-`mu') if `touse'
         replace `ll0'=sum(`ll0'*`wt')
         replace `ll0'=`ll0'[_N]
         local ddev = `dev' - `oldev'
         local i = `i'+1
  }
 
* MAIN MODULE
  replace `mu' = (`y'+0.5)/(`m'+1) if `touse'
  local i=1
  local ddev=1
  replace `eta' = invnorm(`mu') if `touse'
  while ((abs(`ddev') > `ltolera' & `i'<`iterate') | (`i'==0))  {
         if `i'<=`expec'  {
             replace `w'= 1/sqrt(2*_pi)*exp(-(`eta'^2/2))
             replace `z' = (`eta' + (`ym'-`mu')/`w') -`ofs'
             replace `w'=`w'^2/(`mu'*(1-`mu'))
             regress `z' `*' [iw=`w'*`m'*`wt'], mse1  dep(`y') `cons'
         }
         else {
           replace `ndf'= 1/sqrt(2*_pi)*exp(-(`eta'^2/2))
           replace `v' = `mu'*(1-`mu')
           replace `v1' = 1-2*`mu'
           replace `v2'=(`mu'-`mu'^2)^2
           replace `w' = `ndf'^2/`v'
           replace `u' = (`ym'-`mu')/`ndf'
           replace `g1' = 1/`ndf'
           replace  `g2' = `eta'/`ndf'^2
           replace `ws' = `w'+(`ym'-`mu')*(`v'*`g2'+`v1'*`g1')/(`v2'*`g1'^3)
           replace `z' = (`eta'+1/`ws'*`w'*`u') - `ofs'
           regress `z' `*' [iw=`ws'*`m'*`wt'], mse1  dep(`y') `cons'
 
         }
         cap drop `eta'
         predict `eta' if `touse'
         replace `eta' = `eta' + `ofs'
         replace `mu' = normprob(`eta')
         replace `oldev' = `dev'
         #delimit ;
         replace `dev' = (`m'-`y')*log((`m'-`y')/(`m'-(`m'*`mu')))
                         if (`y'==0);
         replace `dev' = `y'*log(`y'/(`m'*`mu')) if (`y'==`m');
         replace `dev' = `y'*log(`y'/(`m'*`mu')) +
                          (`m'-`y')*log((`m'-`y')/(`m'-(`m'*`mu')))
                          if (`y'~=0 & `y'~=`m');
         #delimit cr
         replace `dev' = sum(`dev'*`wt')
         replace `dev' = 2*`dev'[_N]
         replace `ll'=`y'*ln(`mu'/(1-`mu'))+`m'*ln(1-`mu') if `touse'
         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 ddev = `dev' - `oldev'
         local i = `i'+1
  }
 
local pdev = `dev'
gen `chi' = (`y'-`mu'*`m')^2 / (`mu'*`m'*(1-(`mu'*`m')/`m')) 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 pdisp = `cdisp'
local ll = `ll'[_N]
local ll0 = `ll0'[_N]
local pr2 = 1-(`ll'/`ll0')
local ch2 = 2*(`ll'-`ll0')
local cdisp = `cdisp'
}
 
qui regress,  noheader `eform' level(`level')
tempname b V
mat `b' = get(_b)
mat `V' = get(VCE)
mat post `b' `V', depname(`y') obs(`nobs')
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_cdis `cdisp'
global S_E_ddis `dev'/`df'
global S_E_cmd "eprobit"
 
 
}
 
else  {
     if "$S_E_cmd"!="eprobit" {
          error 301
     }
     parse "`*'"
}
if `level'<10 | `level'>99 {
    local level 95
}
if "`eform'" !=""  { local eform "eform(EF)" }
 
noi di in gr _n "Probit 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'
if "`disp'"!="" {
    noi di in gr "chi2 dispersion     = " in ye %9.4f $S_E_cdis
    noi di in gr "deviance dispersion = " in ye %9.4f $S_E_ddis
}
 
if"`residual'"!="" {
  set type float
  qui {
  tempvar yhat stdp Hat Presid Dresid Aresid Spear Sdev pp
  gen `Dresid'=1
  gen `yhat' = `mu'
  gen `pp' = `yhat'/`m'
  predict `stdp',stdp
  gen `Hat'=(`stdp'^2 * `yhat'*(1-`pp'))
  gen `Presid' = (`y'-`yhat')/sqrt(`yhat'*(1-`pp'))
  replace `Dresid' = `m'*log((`m'-`y')/(`m'-`yhat')) if (`y'==0)
  replace `Dresid' = `y'*log(`y'/`yhat') if (`y'==`m')
  #delimit ;
  replace `Dresid' = `y'*log(`y'/`yhat') + (`m'-`y') *
          log((`m'-`y')/(`m'-`yhat')) if (`y'~=0 & `y'~=`m');
  #delimit cr
  replace `Dresid' = sqrt(2*`Dresid')
  replace `Dresid' = -`Dresid' if (`y'-`yhat')<0
  gen `Aresid' = (2.0533902*ibeta(2/3,2/3,(`y'/`m')))- /*
     */ (2.0533902*ibeta(2/3,2/3,(`pp')))
  replace `Aresid' = `Aresid'/(((`pp'*(1-`pp'))^(1/6))*sqrt((1-`Hat')/`m'))
  gen `Spear' = `Presid'/sqrt(1-`Hat')
  gen `Sdev'  = `Dresid'/sqrt(1-`Hat')
  gen _mu = `mu'
  gen _lp = `eta'
  gen _Pear = `Presid'
  gen _Dev  = `Dresid'
  gen _Hat  = `Hat'
  gen _Anscom = `Aresid'
  gen _StandP = `Spear'
  gen _StandD = `Sdev'
  noi di in gr _n "Variables created:  "
  noi di in ye " _mu (fit),                  _lp (linear predictor)"
  noi di in ye " _Pear (Pearson resid),      _Dev (deviance resid)"
  noi di in ye " _Anscom (Anscombe resid),   _Hat (hat matrix diag.)"
  noi di in ye " _StandP (standard Pearson), _StandD (standard deviance)"
  set type double
  }
}
set type float
end