*! Atkinson residual normality plot by Stas Kolenikov. v.1.3
program define atkplot
   version 6
   syntax, [Level(int $S_level) RSTUdent DFIts Cooksd Reps(int 1000) SAving(passthru) Dots Symbol(str) COnnect(str) *]

   if "`e(cmd)'"~="regress" { error 301 }

   local pred `rstudent'`dfits'`cooksd'
   if "`pred'"=="" { local pred rstudent }
   if _N>2040 {
      di in red "too many observations; see help atkplot"
      exit 195
   }
   if "`symbol'"=="" { local symbol o }

   preserve

   tempfile topost tosave source
   tempvar yhat it index ysim itsim vobs vup vlow vexp
   tempname beta origin

   mat `beta'=e(b)
   local sigma=e(rmse)
   local reglist: colnames `beta'
   local iscons=index("`reglist'","_cons")
   if `iscons'==0 { local nocons nocons }
             else { local reglist=substr("`reglist'",1,`iscons'-1) }


   tokenize `reglist'
   local nreg : word count `reglist'
   if "``nreg''"=="_cons" {
       local reglist
       local i=1
       while `i'<`nreg' {
         local reglist `reglist' ``i''
         local i=`i'+1
       }
   }

   if "`dots'"~="" { di }

   * yhat    -- the predicted value
   * it      -- the actual values
   * ysim    -- simulated y
   * itsim   -- simulated it
   quietly {
     g long `index'=_n
     keep if e(sample)
     sort `index'
     save `source'
     local count=_N
     predict `yhat', xb
     predict `it', `pred'
     est hold `origin'
     save `tosave'
     keep `index' `it'
     replace `it'=abs(`it')
     sort `it'
     xpose, clear
     save `topost'
     * hence, the first obs is the obs no, and the second is the actual stat

     * simulation

     * cycle over the simulations
     local k=1
     while `k'<=`reps' {
        use `tosave', clear
        g `ysim'=`yhat'+invnorm(uniform())*`sigma'
        regress `ysim' `reglist', `nocons'
        predict `itsim', `pred'
        replace `itsim'=abs(`itsim')
        keep `itsim'
        sort `itsim'
        xpose, clear
        append using `topost'
        save `topost', replace
        if "`dots'"~="" & mod(`k',10)==0 { noi di in gre "." _c }
        local k=`k'+1
     }

     * simulation data are in memory; the last obs. is the actual data
     local lb=int((100-`level')*`reps'/200)
     local ub=`reps'-`lb'
     expand 3 in l
     * we shall have the last three observations as follows:
     * actual; lower bound; upper bound
     local i=1
     while `i'<=`count' {
       sort v`i' in 1/`reps'
       replace v`i'=v`i'[`lb'] in -2
       replace v`i'=v`i'[`ub'] in -1
       local i=`i'+1
     }
     keep in -4/l
     xpose, clear
     rename v1 `index'
     g `vexp'=invnorm( 0.5+(_n-0.375)/(2*_N+0.25) )
     lab var `vexp' "Expected value"
     rename v2 `vobs'
     rename v3 `vlow'
     rename v4 `vup'
     lab var `vobs' "Observed data"
     lab var `vlow' "Simulated lower bound"
     lab var `vup' "Simulated upper bound"
     sort `index'
     merge `index' using `source', nokeep
     local symstr `symbol'..
     sort `vexp'
     graph `vobs' `vlow' `vup' `vexp', s(`symstr') c(.ll) title("Atkinson simulation envelope") xlab ylab `saving' `options'
     est unhold `origin'
   }
   * end of quietly
   if "`dots'"~="" { di }
end

/*

   History:
    21.01.2001 -- more reliable detection of the constant
*/