*     breg.ado performs quick-and-dirty bounded-influence robust
* regression.  See REGRESSION WITH GRAPHICS (1991, Brooks/Cole),
* by Lawrence C. Hamilton.
*
program define breg
version 2.1
     if "%_*"=="" {
          #delimit ;
               di in gr "Usage:  " in wh  "breg"
               in gr " lhsvar rhsvar [rhsvar ..] ["
               in wh "if"
               in gr "] [" in wh "," in gr " ["
               in wh "graph" in gr "] ["
               in wh "tolerance("
               in gr "#" in wh ")" in gr "]" ;
               #delimit cr
               exit
           }
     mac def _varlist "req ex min(2)"
     mac def _if "opt"
     mac def _options "Graph Tolerance(real .01)"
     parse "%_*"
     parse "%_varlist", parse(" ")
     mac def _lhs "%_1"
     mac shift
     mac def _rhs "%_*"
     capture drop _weight
     capture drop _hat
     capture drop wH
     capture drop MAXD
     capture drop TUNE
     capture drop SCALE
     capture drop RES
     capture drop NN
     capture drop ABSDEV
                        generate TUNE=4.685
     quietly {
          regress %_lhs %_rhs %_if
          generate NN=_result(1)
          predict RES %_if, residuals
          summarize RES %_if, detail
          gen SCALE=. %_if
          gen ABSDEV=abs(RES-_result(10)) %_if
          gen MAXD=1 %_if
          gen _weight=1 %_if
          replace _weight=. if RES==.
          predict _hat %_if, hat
          summ _hat %_if, detail
          generate wH=1
          replace wH=(_result(12)/_hat)^2 if _hat>_result(12)
          }
     mac def _it 1
     quietly summarize MAXD %_if
     while _result(6)>%_toleran & %_it<51 {
          capture drop OLDW
          rename _weight OLDW
          quietly {
               summarize ABSDEV %_if, detail
               gen _weight=1
               #delimit ;
               replace _weight=wH*2*_result(10)/abs(RES)
                           if abs(RES)>2*_result(10);
               #delimit cr
               replace _weight=. if RES==.
               regress %_lhs %_rhs =_weight %_if
               if "%_graph"~="" {
                    #delimit ;
                    graph _weight OLDW OLDW %_if,
                         sy(oi) sort c(.l)
                         l1("New Weight")
                         b2("Old Weight") ;
                    #delimit cr
               }
               drop RES
               predict RES %_if, residual
               summarize RES %_if, detail
               replace ABSDEV=abs(RES-_result(10)) %_if
               replace MAXD=abs(_weight-OLDW) %_if
               summarize MAXD %_if
          }
          #delimit ;
          di in gr "Huber iteration %_it:  "
               "maximum difference in weights = "
               in ye _result(6) ;
          #delimit cr
          macro define _it = %_it+1
     }
     quietly {
          replace MAXD=1
          summ MAXD
          }
     capture drop OLDW
     rename _weight OLDW
     quietly {
          summarize ABSDEV %_if, detail
          replace SCALE=_result(10)/.6745 %_if
          gen _weight= wH*max(1-(RES/(TUNE*SCALE))^2,0)^2 %_if
          replace _weight=. if RES==.
          regress %_lhs %_rhs =_weight %_if
          if "%_graph"~="" {
               #delimit ;
               graph _weight OLDW OLDW %_if,
               sy(oi) sort c(.l)
               l1("New Weight")
               b2("Old Weight") ;
               #delimit cr
               }
          drop RES
          predict RES %_if, residual
          summarize RES %_if, detail
          replace ABSDEV=abs(RES-_result(10)) %_if
          replace MAXD=abs(_weight-OLDW) %_if
          summarize MAXD %_if
          }
          #delimit ;
          di in gr "Biweight iteration %_it:  "
               "maximum difference in weights = "
               in ye _result(6) ;
          #delimit cr
     regress
     quietly {
     label variable _weight "Bounded-influence weights"
     label variable _hat "OLS leverage (hat diagonals)"
     drop wH RES SCALE MAXD NN OLDW TUNE ABSDEV
     }
di in gr "NOTE:  Standard errors and tests are not correct"
end