program define grubbs
* this is a revised version of the original command, it no longer deletes missing obs or outlier
* instead, it sets "tag_grubbs" =1 if it believes the obs to be an outlier
* usage: "grubbs myvar .05 10"
version 8.0
*arguments:
* 1= Name of variable
* 2= Confidence interval (0.05 or 0.01)
* 3= Max number of iterations
args xvar conf maxit
tempvar dev missx
gen byte tag_grubbs = 0
local i = 1
di "deleting missing values"
gen byte missx =  `xvar'==.
* initial guess for critical value
scalar Gcrit = 10
* start with G > Gcrit (otherwise loop will not begin)
scalar G = Gcrit +1
di "maxit = " `maxit'
di "G= " G
di "Gcrit = " Gcrit
         while G > Gcrit & `i'<= `maxit' {
                   sum `xvar' if tag_grubbs == 0
                   local nobs = r(N)
                   gen `dev' = (abs(`xvar' -r(mean)))/r(sd)
                   gsort -`missx' -tag_grubbs `dev'
                   scalar G = `dev'[_N]
                   local ct = `conf'/(2*`nobs')
                   local ts = invttail(`nobs'-2,`ct')
                   scalar Gcrit = (`nobs'-1)*sqrt(`ts'^2/(`nobs'*(`nobs'-2+`ts'^2)))
                   di "Iteration = " `i' "   Critical G = " Gcrit " Current G = " G
                   if (G > Gcrit) di `xvar'[_N] " is an outlier, so tag_grubbs = 1"
                   replace tag_grubbs = 1 if  `dev' == G & G > Gcrit
                   local i = `i'+1
                   drop `dev'
}

if (`i'<=`maxit') di "Grubbs procedure terminated: no more outliers"
else di "Maximum iterations exceeded: Use larger maxit"
end