* rdiag  version 1.01    Robert A. Yaffee    yaffee@nyu.edu
* primary objective   cross-sectional regression diagnostic analysis of standardized
*     residuals



cap program drop rdiag
program rdiag, rclass

cap drop posout1
cap drop negout1
cap drop dum
cap drop hat1
cap drop cd1-prboutliers






version 12.1
set more off
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
local lhs neiw1 efradw3 BSIsoma
predict `1', rstandard
}
matrix define Assumptns = J(12,1,0)
matrix rownames Assumptns =  Sb POuts NOuts Badouts Collin SW KS Het Spec indep Coltot Pcttot
matrix colnames Assumptns =  Problems
matlist Assumptns



* structural break test
*  Is there a level shift in your data?


// a non parameteric structural break test for a level shift




di  as input in red "Testing for a structural break in the data"

centile `1', centile(33 67)
capture confirm variable resgroup
if _rc != 0 {
tempvar resgroup
qui gen `resgroup' = 1 if `1' < r(c_1)
qui replace `resgroup' = 2 if `1' >= r(c_1)
qui replace `resgroup' = 3 if `1' >= r(c_2)
qui kwallis `1', by(`resgroup')
local x2 `r(chi2_adj)'

loc Sb chi2tail(`r(df)',`r(chi2_adj)')



if `Sb'  < .05  {
scalar strucbk  = 1
return scalar strucbk = `Sb'
}
else {
scalar strucbk = 0
return scalar strucbk = `Sb'
}





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




di "{hline}"
di as input  in red "Outlier dignostics"
loc i 0
// First we count the outliers
local sampsize _N
predict cd1, cooksd
predict dfits1, dfits
predict hat1, lev
qui gen posout1 = 0
qui gen negout1 = 0
extremes `1'
qui replace posout1 =1 if `1' > 3.5 & `1' < .
qui count if posout1==1
noi di as input "There are `r(N)' positive outliers "
list id cd if posout1==1
qui count if posout1==1
qui gen posoutliers=r(N)
return scalar posoutliers = r(N)
qui replace negout1=1 if `1' < - 3.5
qui count if negout1==1
qui gen negoutliers=r(N)
noi di as input "There are `r(N)' negative outliers:  "
list id cd if negout1==1
return scalar negoutliers = `r(N)'
qui gen totout = r(negoutliers) + r(posoutliers)
count if posout1 > 0
scalar pos = r(N)
return scalar pos = `pos'
count if negout1 > 0
scalar neg = r(N)
return scalar neg = `neg'

local hlfss = (`sampsize'/2)
gen dum = 0
replace dum = 1 if id > `hlfss'

stem `1'`flg'
di "  "
di "  "
di as input "We now examine their leverage"
di "  "

// We determine whether these outliers have bad leverage (distort the predicted value much)?
// leverage analysis

tempvar cdcritval1
tempvar dfcritval
gen  `cdcritval1' = 3*e(df_m)/e(N)
scalar ckcrit = `cdcritval1'
gen `dfcritval' = 2*sqrt(e(df_m)/e(N))
return scalar cdcritval1 = `cdcritval1'
return scalar dfcritval = `dfcritval'
di "Cooks d critical value = " `cdcritval1'
di "  "
qui gen prboutliers=0
qui replace prboutliers = 1 if cd > `cdcritval1' & cd < . & totout > 0 & totout < .
qui count if prboutliers==1
scalar badout= `r(N)'
return scalar numbadout = `r(N)'
noi di as input "The number of problematic outliers is: ",  `badout'
noi di as input " and they are : "
list id if prboutliers==1


list id cd  if prboutliers==1


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

lvr2plot, mlabel(id)
gr save lvr2p.gph, replace



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

coldiag2
collin `lhs'
scalar condxnmbr = r(cn)
return scalar condxnmber = r(cn)

if condxnmbr > 4 {
scalar multc = 1
return scalar multc= multc
}
else {
scalar multc = 0
return scalar multc = multc
}



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'
if r(p) < 0.05 {
scalar sw = 1
return scalar sw = 1
}
else {
scalar sw = 0
return scalar sw=0
}

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

jb `1'

di "{hline}"
di as input in red "Residual heteroskedasticity tests:"
// Are the residual distributions homoskedastic?
// p >= .05 means H0 of homoskedasticity cannot be rejected for these tests


estat hettest
if r(p) < 0.05 {
scalar Bpagan = 1
return scalar Bpagan =1
}
else {
scalar Bpagan = 0
return scalar Bpagan = 0
}


estat imtest
if r(p) < .05 {
scalar white = 1
return scalar white=1
}
else {
scalar white = 0
return scalar white = 0
}


scalar het = (Bpagan + white)/2
return scalar het = (Bpagan + white)/2

// graphical review of residual distributions

qnorm `1', title(quantile normal plot of `1')
gr save qn`1'.gph, replace

gr box `1', title(Box plot of `1')
gr save bx`1'.gph, replace

hist `1', normal frequency title(histogram of `1')
gr save hist`1'.gph, replace


graph combine  qn`1'.gph  bx`1'.gph  lvr2p.gph  hist`1'.gph, col(2) title(Analysis of `1') ///
  subtitle(graphical residual analysis)
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 = 1
}
di "{hline}"
di as input in red  "Independence of Observations:"
cap drop obsnum

runtest `1'

if r(p) < .05 {
scalar indep = 1
return scalar indep = 1
}
else {
scalar indep = 0
return scalar indep = 0
}


scalar coltot = (strucbk + badout + multc + (sw + skrp)/2 + (het+ white)/2 + reset + indep)
return scalar coltot = (strucbk + badout + multc + (sw + skrp)/2 + (het+ white)/2 + reset + indep)
scalar pcttot = 100-(coltot/7)*100
return scalar pcttot = pcttot


di "{hline}"
matrix define Assumptns = (strucbk \ pos \ neg \ badout  \ multc  \ sw \ skrp \ het \ reset \ indep \coltot \ pcttot )
matrix rownames Assumptns = Structrl-Brk  Positive-outliers Negative-outliers Bad-outliers collinearity Shapiro-Wilk  ///
  Kolmogorov-Smirnov Heteroskedasticity Specification Independence_Obs Total_prob_type Percent_Valid
matrix colnames Assumptns =  Problems=1
di _skip(4)
di as input in red "Assessment of fulfillment of cross-sectional statistical regression assumptions: "
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 = SB + Bad-Outliers+collinearity + (SW + KS)/2 + het + Spec + independence"
di "Percent_Valid = 100 - (percent of violated assumptions)"
di _skip(3)
di "{hline}"
di "Legend:"
di "SB (Structural Break)  =  Significant diff of levels in residuals"
di "PosOutliers = positive outlier issue"
di "NegOutliers = negative outlier issue"
di "Badoutliers = Number of outliers with bad leverage"
di "collinearity =  Colinearity problem when cnd indx > 4"
di "SW = Significant Shapiro-Wilk test"
di "KS = Significant Kolmogorov-Smirnov test"
di "Het = (Significant Breusch-Pagan + significant white heteroskedasticity test)/2"
di "Spec = Significant Ramset23 reset test"
di "Independence of obs = runs test for independence"
di "Column total"
di "Percent Valid=  percent of regression model assumptions not violated"
di "{hline}"


erase qn`1'.gph
erase bx`1'.gph
erase hist`1'.gph
cap drop myrstd
cap drop cd
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
cap drop posout1
cap drop negout1
cap drop dum
cap drop posout1
cap drop negout1
cap drop dum
cap drop hat1
}
end