*
* @BPPanelTests( options ) u
*
* Tests the residuals from a least squares regression for presence of
* random effects using variants of the Breusch-Pagan LM test. It
* produces the standard BP test (which has chi-squared asymptotics), the
* Honda one-sided form (which has N(0,1) asymptotics) and the SLM
* (standardized Lagrange Multiplier) which is similar to the Honda test
* but uses small-sample corrections for the mean and standard deviation.
*
* Parameters:
*   U = series of residuals from the OLS regression (not FE or any other
*       panel regression).
*
* Options:
*   EFFECTS=[INDIV]/TIME/BOTH
*   TITLE=title for report ["Breusch-Pagan Tests for Random Effects"]
*
* References:
*   Breusch & Pagan(1980), "The Lagrange Multiplier Test and Its
*     Applications to Model Specification in Econometrics," Review of
*     Economic Studies, vol. 47, no. 1, 239-53.
*
*   Baltagi(2008), "Econometric Analysis of Panel Data", 4th ed. for the
*     description of the Honda and SLM forms.
*
* Revision Schedule:
*   04/2011 Written by Tom Doan, Estima.
*
procedure BPPanelTests u
type series u
*
option choice effects 1 indiv time both
option string title "Breusch-Pagan Tests for Random Effects"
*
local real    tindiv   ttime
local real    bpindiv  bptime   bpjoint
local real    hindiv   htime    hjoint
local real    slmindiv slmtime
local integer start end nindiv ntime
local symm    xxols
local rect    xxadjust
local real    trdpz ed trdpzdpz vd
*
compute start =%regstart(),end=%regend()
compute nindiv=%indiv(end)-%indiv(start)+1
compute ntime =%panelobs()
*
if effects==1.or.effects==3 {
   mcov(cluster=%indiv(t))
   # u
   compute tindiv=%cmom(1,1)/%rss
}
if effects==2.or.effects==3 {
   mcov(cluster=%period(t)) start end
   # u
   compute ttime=%cmom(1,1)/%rss
}
*
* Breusch-Pagan LM tests
*
if effects==1.or.effects==3
   compute bpindiv=(tindiv-1)^2*%nobs/(2*(ntime-1))
if effects==2.or.effects==3
   compute bptime =(ttime-1)^2*%nobs/(2*(nindiv-1))
if effects==3
   compute bpjoint=bpindiv+bptime
*
* Honda tests
*
if effects==1.or.effects==3
   compute hindiv=(tindiv-1)*sqrt(%nobs/(2.0*(ntime-1)))
if effects==2.or.effects==3
   compute htime =(ttime-1)*sqrt(%nobs/(2.0*(nindiv-1)))
if effects==3
   compute hjoint=(hindiv+htime)/sqrt(2)
*
* SLM tests
*
compute xxols=%xx
if effects==1.or.effects==3 {
   mcov(cluster=%indiv(t),lastreg) start end
   compute xxadjust=xxols*%cmom
   compute trdpz=(%nobs-%trace(xxadjust))
   compute ed=trdpz/%ndf
   compute trdpzdpz=%nobs*ntime-2*ntime*%trace(xxadjust)+%trace(xxadjust*xxadjust)
   compute vd=2*(%ndf*trdpzdpz-trdpz^2)/(%ndf^2.0*(%ndf+2))
   compute slmindiv=(tindiv-ed)/sqrt(vd)
}
if effects==2.or.effects==3 {
   mcov(cluster=%period(t),lastreg) start end
   compute xxadjust=xxols*%cmom
   compute trdpz=(%nobs-%trace(xxadjust))
   compute ed=trdpz/%ndf
   compute trdpzdpz=%nobs*nindiv-2*nindiv*%trace(xxadjust)+%trace(xxadjust*xxadjust)
   compute vd=2*(%ndf*trdpzdpz-trdpz^2)/(%ndf^2.0*(%ndf+2))
   compute slmtime=(ttime-ed)/sqrt(vd)
}
*
report(action=define,title=title)
report(atrow=1,atcol=1,span) title
report(atrow=2,atcol=1,fillby=cols,align=left) "Breusch-Pagan" "" "Honda" "" "SLM" ""
if effects==1.or.effects==3
   report(atrow=1,atcol=2,fillby=cols) "Indiv" bpindiv %chisqr(bpindiv,1) $
      hindiv %cdf(-hindiv) slmindiv %cdf(-slmindiv)
if effects==2.or.effects==3
   report(atrow=1,atcol=%reportcol+1,fillby=cols) "Time" bptime %chisqr(bptime,1) $
     htime  %cdf(-htime)  slmtime %cdf(-slmtime)
if effects==3
   report(atrow=1,atcol=%reportcol+1,fillby=cols) "Joint" bpjoint %chisqr(bpjoint,1) $
     hjoint %cdf(-hjoint)
*
report(atrow=3,torow=3,atcol=2,tocol=4,action=format,special=parens)
report(atrow=5,torow=5,atcol=2,tocol=4,action=format,special=parens)
report(atrow=7,torow=7,atcol=2,tocol=4,action=format,special=parens)
report(action=format,atcol=2,atrow=2,picture="*.####",align=decimal)
report(action=show)
end

