* rdiagrbst  version 1.0.0    Robert A. Yaffee    yaffee@nyu.edu      22 June 2012
* rdiagrbst  version 1.1.0       13 July 2012  fixed kd title
* primary objective   cross-sectional regression diagnostic analysis of standardized
*     residuals to test validity of final OLS regression models when the models are robust
*     This program overlooks the hettest and computes the scores based on all else
set more off
set linesize 140
cap drop tskew
cap scalar drop tskew
cap drop hat
cap drop resid2
cap drop myrstd
cap drop outout
cap drop badout2
cap drop badout
cap drop bdout
cap drop dfits
cap drop posout1
cap drop negout1
cap drop posoutliers
cap drop negoutliers
cap drop totaloutliers
cap drop outlier
cap drop badoutsum
cap drop totout
cap drop outlier
cap drop prboutliers
cap drop hatthreshld
cap drop rstd
cap drop rstd1
cap drop rstd2
cap drop r1
cap drop cd
cap drop cd1
cap program drop rdiag
cap drop dum
cap drop multc
cap drop strb
cap drop obsnum
cap drop MaXBF
cap drop strbrk
cap drop totout
cap scalar drop pcttot
cap scalar drop coltot
cap scalar drop multc
cap scalar drop reset
cap scalar drop het
cap scalar drop white
cap scalar drop Bpagan
cap scalar drop skrp
cap scalar drop sw
cap scalar drop condxNumb
cap scalar drop badout
cap scalar drop ckcrit
cap scalar drop strbrk
cap scalar drop MaxBoxF
cap scalar drop totout



cap program drop rdiagrbst
program rdiagrbst, rclass
version 12.1
set more off
MyPrepare
cap gen id = _n
capture confirm variable id
if _rc != 0 {
gen id = _n
}
local flag 0
capture confirm variable rstd1
if _rc != 0 {
local pos 0
local neg 0
local cksd 0
local skrp 0
local ht 0
local jr 0
predict `1'
local yvar e(depvar)
di "yvar = ",`yvar'
local cn7:colnames(e(b))
di "`cn7'"
local len7 = length("`cn7'")
di `len7'
local len7b = `len7'- 6
di `len7b'
local xvar = substr("`cn7'",1,`len7b')
di "xvar = ,""`xvar'"
}



matrix define Assumptns = J(10,1,0)
matrix rownames Assumptns =  POuts NOuts Badouts Collin SW KS Spec Strbrk Coltot Pcttot
matrix colnames Assumptns =  Problems
matlist Assumptns

* outlier diagnostic analysis
*  How many positive and negative outliers
*  How many have bad leverage


di "{hline}"
di as input  in red "Outlier dignostics"

di "{hline}"
summ `1', detail

capture confirm variable `1'std
if _rc !=0 {
gen `1'std = (`1'-r(mean))/r(sd)
}

local sampsize _N

tempvar q1
qui gen `q1'=0
replace `q1'= 1 if _n < `sampsize'/4


tempvar segmt1
qui gen `segmt1'=0
replace `segmt1'=1 if `q1'==1


tempvar q2
qui gen `q2' = 0
replace `q2' = 1 if _n >= .25*`sampsize' & _n  < .50*`sampsize'


tempvar segmt2
qui gen `segmt2'=0
replace `segmt2'=1 if `q2'==1



tempvar q3
qui gen `q3'= 0
replace `q3'=1 if _n >= .5*`sampsize' & _n < .75*`sampsize'


tempvar segmt3
qui gen `segmt3'=0
replace `segmt3'=1 if `q3'==1


tempvar q4
qui gen `q4' = 0
replace `q4' = 1 if _n >= .75*`sampsize' & _n != .


tempvar segmt4
qui gen `segmt4' = 0
replace `segmt4' = 1 if `q4'==1


tempvar minvar
summ `segmt1'
tempvar v1
gen `v1' = r(Var)
summ `segmt2'
tempvar v2
gen `v2' = r(Var)
summ `segmt3'
tempvar v3
gen `v3' = r(Var)
tempvar v4
summ `segmt4'
gen `v4'= r(Var)
gen `minvar' = min(`v1',`v2', `v3', `v4')
di `minvar'
tempvar maxvar
gen `maxvar'= max(`v1',`v2',`v3',`v4')
di `maxvar'

di as input in red "Split sample analysis"
di as input "There are 4 subsamples of equal sample size."
di _skip(2)
di as result "Variance 1 = ",`v1'
di as result "Variance 2 = ",`v2'
di as result "Variance 3 = ",`v3'
di as result "Variance 4 = ",`v4'
di _skip(1)
di as result "Minimum variance = ", `minvar'
di _skip(2)
di as result "Maximum variance = ", `maxvar'
di _skip(2)
di as result  "Computing Box's F test as indication of stuctural break:"
di _skip(1)
tempvar BoxF1
qui gen `BoxF1' = `v1'/`minvar'
tempvar BoxF2
qui gen `BoxF2' = `v2'/`minvar'
tempvar BoxF3
qui gen `BoxF3'= `v3'/`minvar'
tempvar BoxF4
qui gen `BoxF4' =  `v4'/`minvar'
di as result "for segment 1  Box's F = ",`BoxF1'
di as result "for segment 2  Box's F = ",`BoxF2'
di as result "for segment 3  Box's F = ",`BoxF3'
di as result "for segment 4  Box's F = ",`BoxF4'


local maxboxf = max(`BoxF1',`BoxF2',`BoxF3',`BoxF4')
local MAXBoxF = `maxboxf'


capture confirm variable strbrk
if _rc !=0 {
qui gen strbrk = 0
}
else {
drop strbrk
}
capture confirm variable strbrk
if _rc !=0 {
qui gen strbrk = 0
label var strbrk "structural break dummy"
}

qui replace strbrk = 1 if `MAXBoxF' >= 3.84
cap return scalar strbrk = strbrk



if strbrk ==1  {
di _skip(1)
di as input in red "There is a structural break and transformation may be nec.!"
}


else  {
di _skip(1)
di as input in red "There appears to be no structural break."
}

cap scalar drop tskew
tempvar posout1
tempvar negout1
local sampsize _N

qui gen `posout1'= 0
qui gen `negout1'= 0
extremes `1'std
qui replace `posout1' =1 if `1'std > 3.5 & `1'std < .
qui count if `posout1'==1

noi di as input "There are `r(N)' positive outliers "
list id `1'std if `posout1'==1
qui count if `posout1'==1
tempvar numpoutliers
qui gen `numpoutliers'=r(N)
return scalar numpouts = r(N)
qui replace `negout1'=1 if `1' < - 3.5
tempvar numnoutliers
qui count if `negout1'==1
qui gen `numnoutliers'=r(N)
noi di as input "There are `r(N)' negative outliers:  "
list id `1'std if `negout1'==1
return scalar numnegouts = `r(N)'
tempvar totouts
gen `totouts'=0
replace `totouts' = r(numpouts) + r(numnegouts)
count if `posout1' > 0 & `posout1' < .
local pos = r(N)
count if `negout1' > 0 & `negout1' < .
loc neg = r(N)

cap drop bdout
qui gen outlier= 0
replace outlier=1 if `1'std > 3.5 | `1'std < - 3.5
count if `1'std > 3.5 | `1' < -3.5
return scalar totoutliers = r(N)

capture confirm variable hatthreshld
if _rc != 0 {
gen hatthreshld = 3*e(df_m)/e(N)
}
else if _rc==0 {
cap drop badoutsum
}
summ outlier if `1'std > 3.5 | `1'std < 3.5 , detail
gen badoutsum = r(sum)

local sampsize _N
tempvar firsthlf
cap confirm variable `firsthlf'
if _rc ! = 0 {
local hlfss = `sampsize'/2
gen `firsthlf' = 0
replace `firsthlf'= 1 if _n > `hlfss'
}
tab `firsthlf'

stem `1'std`flg'
di "  "
di "  "
local N _N
di as input "Computing the leverage h"
qui gen tskew=0
// h_i = (x - xbar)^2/sum_(x - xbar)^2
// we use the standard error of skewness to determine whether an outlier is a problem
local ses sqrt((6*`N')*(`N'-1)/((`N'-2)*(`N'+1)*(`N'+2)))
summ `1'std, detail
di _skip(4)
di as input "standard error of snewness= ",`ses'
replace tskew = r(sknewness)/`ses'
return scalar tskew = r(skewness)/`ses'
di _skip(4)
di as input "t test for skewness = ", r(tskew)
if abs(r(tskew)) > 2  {
gen badout=1
}
else {
gen badout=0
}
scalar tskewed = tskew
return scalar badout=badout
scalar bdout1=badout
di as input "badout= ", badout

di as input "Positive outliers are : "

list id  if `posout1' == 1
count if `posout1'==1
di as input "Negative outliers are: "
list id  if `negout1'==1
count if `negout1'==1


di _skip(4)
di "{hline}"

di as input in red "multicollinearity tests:"


// Variance Inflation Factor = 1/(1 - R^2) =   1/Tolerance
// VIF levels above 20 are seriously problematic
// VIF levels above 10 can be problematic also
// if the condition number is above 30 there is significant collinearity
//  condition indices above 15 are problematic acc to R. Williams at Notre Dame
estat vif
coldiag2
local maxrow rowsof(r(cx))
matrix define cx = r(cx)
scalar condxNumb = el(cx,`maxrow',1)
return scalar condxNumb = condxNumb



di "{hline}"
di as input in red "Residual normality tests:"
// normality of the residual distribution p >= .05 means H0 of nonormality cannot be rejected
// for any of these tests


swilk `1'std
if r(p) < 0.05 {
scalar sw = 1
return scalar sw = 1
}
else {
scalar sw = 0
return scalar sw=0
}

sktest `1'std
if r(P_chi2) < 0.05 {
scalar skrp= 1
return scalar skrp = 1
}
else {
scalar skrp = 0
return scalar skrp = 0
}

jb `1'


qui replace strbrk = 1 if `MAXBoxF' >= 3.84

return scalar strbrk = strbrk



// graphical review of residual distributions

qnorm `1'std, title(quantile normal plot of `1') note("depvar = `e(depvar)' vs" "`xvar'")
gr save qn`1'std.gph, replace


gr box `1'std, title(Box plot of `1')  note("depvar = `e(depvar)' vs. " "`xvar'")
gr save bx`1'std.gph, replace

hist `1'std, normal frequency title(histogram of `1')  note(depvar = `e(depvar)' vs  "`xvar'")
gr save hist`1'std.gph, replace

kdensity `1'std, title("Kernel density distribution" "of `1'std")
gr save kd`1'std.gph, replace

graph combine  qn`1'std.gph  bx`1'std.gph kd`1'std.gph   hist`1'std.gph, col(2) title(Analysis of `1'std)
gr save combi4.gph, replace

di "{hline}"
di as input in red  "Specification error tests:"
estat ovtest

if r(p) < 0.05 {
scalar reset = 1
return scalar reset=1
}
else {
scalar reset = 0
return scalar reset = 0
}

if condxNumb > 20 {
scalar multc = 1
}
else {
scalar multc = 0
}

if `MAXBoxF' > 3.84 {
scalar strbrk = 1
return scalar strbrk = 1
}
else {
scalar strbrk = 0
return scalar strbrk = 0
}

capture confirm variable badout
if _rc != 0 {
cap gen badout=1 if r(badoutsum) > 0 & badoutsum < .
}
else if r(badoutsum)==0 {
cap gen badout = 0
}

scalar badout=badout

scalar coltot = (badout + multc + (sw + skrp)/2 +  reset + strbrk)
return scalar coltot = badout + multc + (sw + skrp)/2 +  reset + strbrk
scalar pctinvalid = (coltot/5)*100
scalar pctvalid = 100 - pctinvalid
return scalar pctinvalid = (coltot/5)*100
return scalar pctvalid = 100- pctinvalid






di "{hline}"
matrix define Assumptns = ( `pos' \ `neg' \ badout  \ multc  \ sw \ skrp \ reset \ strbrk \ coltot \ pctinvalid \ pctvalid)
matrix rownames Assumptns =  Positive-outliers Negative-outliers Bad-outliers collinearity Shapiro-Wilk  ///
  Kolmogorov-Smirnov Specification StructuralBreak TotalProbs Percent-Invalid Pct-Pseudo-Valid
matrix colnames Assumptns =  Problems=1
di _skip(4)
di as input in red "Assessment of fulfillment of cross-sectional statistical regression assumptions: "
di as input in green "Caveat: This score is a heuristicc only as these scores are not strictly additive"
matlist Assumptns, title(Regression model assumption test result matrix) left(18) twidth(18) lines(rows)  ///
   border format(%12.0g)
di _skip(3)

di "Problem type total score = Bad-Outliers+collinearity + (SW + KS)/2 + Spec + strbrk "
di "Percent_Valid = 100 - (percent of violated assumptions)"
di _skip(3)
di "{hline}"
di "Legend:"
di "PosOutliers = positive outlier issue"
di "NegOutliers = negative outlier issue"
di "Badoutliers = Number of outliers with bad leverage"
di "collinearity =  Collinearity problem when condition number > 20"
di "SW = Significant Shapiro-Wilk test"
di "KS = Significant Kolmogorov-Smirnov test"
di "Spec = Significant Ramset23 reset test"
di "strbrk = Existence of structural break in residual variance between 4 sample segments"
di "Column total = the total number of counted problems in the column vector"
di "Percent invalid: Percent of 6 regression assumptions violated "
di "Percent Pseudo-Validity =  percent of regression model assumptions not violated"
di "{hline}"


cap erase qn`1'.gph
cap erase bx`1'.gph
cap erase hist`1'.gph
cap drop myrstd
cap drop cd
cap drop cd1
cap drop posout
cap drop negout
cap drop numout
cap drop badout
cap drop totout
cap drop posoutliers
cap drop negoutliers
cap drop numoutliers
cap drop prboutliers
cap drop dfits
cap drop cd
cap drop numout
cap drop rstd-cd
cap drop hat
cap drop rstd

end