*
* @Breitung( options ) series start end
*
* Implements tests from Breitung(2000), "The local power of some unit
* root tests for panel data", from Baltagi, Fomby, Hill (eds),
* "Nonstationary Panels, Panel Cointegration, and Dynamic Panels
* (Advances in Econometrics, Volume 15)", Emerald Group Publishing,
* pp.161-177
*
* Parameters:
*
*   series     = series to analyze
*   start end  = range to use in performing tests [by default, maximum range of series]
*
* Options:
*
*   SMPL=Dummy series with 0's in entries to skip [include all]
*   DET=NONE/[CONSTANT]/TREND
*     Deterministic components
*
*   LAGS=(maximum) number of additional lags of the differenced series to
*    include [Schwert's suggestion]
*   CRIT=FIXED/[GTOS]/AIC/BIC/HQ
*     Criterion used in selecting the individual specific lags. FIXED uses
*     the value of the LAGS option for each. GTOS is General TO Specific,
*     starting at the value of the LAGS option and dropping lags as long as
*     the t-stat on the final one has a significance level below the SLSTAY
*     option. AIC, BIC and HQ the Akaike, Schwarz Bayesian and Hannan-Quinn,
*     respectively.
*   SLSTAY=significance level for keeping the marginal lag in CRIT=GTOS [.10]
*
*   ROBUST/[NOROBUST]
*     Uses robust (Eicker-White) standard errors in computing the test statistic.
*
*   [PRINT]/NOPRINT
*   TITLE=title for output ["Breitung Unit Root Test Series: ..."]
*
* Variables Defined:
*   %CDSTAT - the normalized t-statistic
*   %SIGNIF - one tailed significance level
*   %NGROUP - number of individuals
*
* Revision Schedule:
*  03/2012 Written by Tom Doan, Estima
*
procedure Breitung y start end
type series 	y
type integer 	start end
*
option choice  	det    		2  none constant trend
option integer 	lags
option choice  	crit   		2  fixed gtos aic bic hq
option real    	slstay 		.10
option series  	smpl
option switch  	print    	1
option string     title
option switch     robust      0
*
local integer 		startl endl
local integer 		i tnobs
local series 		dy y1 icount ehat vhat
local equation 	extras
local integer   	p pmax tmax tend starti endi lagsused slot
local series    	itrend
local vect[int] 	detvars xdetvars pchoice
local rect      	sweepxx
local real      	marginal_t ic bestic logl
local real        pt nt avglag lowt hight alpha beta
local vect      	sigmasq tsample
local report      breport
local string      ltitle descript lstyle

if .not.%defined(y) {
   disp "###Syntax: @BREITUNG( options ) >>>>>SERIES<<<<<"
   return
}
inquire(series=y) startl<<start endl<<end
*
* Compute first differenced data
*
clear dy
*
set dy = y-y{1}
*
* Figure out which individuals are covered by the range, and get the
* observation count for each.
*
compute starti=%indiv(startl),endi=%indiv(endl)
panel(icount=1,smpl=smpl,compress) y startl endl icount
sstats(max) starti endi icount>>pt
sstats      starti endi (icount>0)>>nt
compute tmax=fix(pt),tend=%panelobs()
*
dim pchoice(endi) sigmasq(endi) tsample(endi)
*
* Create the individual trend series in case we need it.
*
set itrend = %period(t)
*
* Create the set of regressors. The detvars are only used in choosing
* lag lengths.
*
if det==1 {
   dim detvars(0) xdetvars(0)
   compute descript="None"
}
else
if det==2 {
   compute detvars=||constant||
   dim xdetvars(0)
   compute descript="Constant"
}
else {
   compute detvars=||constant,itrend||
   compute xdetvars=||constant||
   compute descript="Constant+Trend"
}
*
* Figure out lag lengths for individual specific ADF regressions. Use
* Schwert's rule of thumb if the maximum lag length wasn't given.
*
if .not.%defined(lags)
   compute pmax=fix(%round(12*(.01*pt)^.25,0))
else
   compute pmax=lags
*
compute avglag=0.0
do i=starti,endi
   *
   * Skip an individual with no data
   *
   if icount(i)==0
      next
   if crit==1 {
      compute lagsused=pmax
      goto chosen
   }
   *
   * Set lagsused to 0, which will be the value if we don't find anything
   * significant at 1 or higher.
   *
   compute lagsused=0
   cmom(smpl=smpl) i//(pmax+2) i//tend
   # dy y{1} detvars dy{1 to pmax}
   *
   * Start by regressing out all the variables and work backwards.
   *
   compute sweepxx=%sweeplist(%cmom,%seq(2,%ncmom))
   compute slot=%ncmom
   do p=pmax,(crit==2),-1
      compute %ndf=%nobs-(slot-1)
      if crit==2 {
         *
         * Compute the (absolute) t-statistic from the information in
         * <<sweepxx>>. The pivot is at <<slot>>, which is the number of
         * regressors + 1 for the dependent variable.
         *
         compute marginal_t=abs(sweepxx(slot,1))/$
           sqrt(sweepxx(slot,slot)*sweepxx(1,1)/%ndf)
			*
         * If that's significant at the requested level, break the loop.
         *
         if marginal_t>%invttest(slstay,%ndf) {
            compute lagsused=p
            break
         }
      }
		else
		if crit==3.or.crit==4.or.crit==5 {
      	compute logl=-.5*%nobs*log(sweepxx(1,1)/%nobs)
      	if crit==3
         	compute ic=2.0*logl-2*p
      	else
      	if crit==4
         	compute ic=2.0*logl-log(%nobs)*p
      	else
         	compute ic=2.0*logl-2*log(log(%nobs))*p
			if p==pmax.or.ic>bestic
         	compute bestic=ic,lagsused=p
		}
		*
      * Sweep on the last lag to take it out of the regression.
   	*
   	compute sweepxx=%sweep(sweepxx,slot)
	   compute slot=slot-1
	end do p
   :chosen
   compute pchoice(i)=lagsused
   compute avglag=avglag+lagsused
end do i
compute avglag=avglag/nt
if crit==1
   compute lstyle="Using fixed lags "+pmax
else
   compute lstyle="With average lags "+%strval(avglag,"*.##")+" chosen from "+pmax
*
clear(zeros) ehat vhat
do i=starti,endi
   if icount(i)==0
      next
   equation extras *
   # xdetvars dy{1 to pchoice(i)}
   sstats(smpl=smpl,min) i//(pchoice(i)+1) i//tend %if(%valid(y),t,%na)>>lowT
   sstats(smpl=smpl,max) i//(pchoice(i)+1) i//tend %if(%valid(y),t,%na)>>highT
   *
   * Detrend/demean lagged regressor
   *
   compute alpha=y(fix(lowT))
   compute beta =(y(fix(highT))-y(fix(lowT)))/(highT-lowT)
   if det==1
      set y1 i//(pchoice(i)+2) i//tend = y{1}
   else
   if det==2
      set y1 i//(pchoice(i)+2) i//tend = y{1}-alpha
   else
      set y1 i//(pchoice(i)+2) i//tend = y{1}-alpha-beta*(t-lowT-1)
   *
   * Compute cross product with augmenting lags of original data.
   *
   cmom(smpl=smpl) i//(pchoice(i)+2) i//tend
   # dy y1 xdetvars dy{1 to pchoice(i)}
   compute sweepxx=%cmom
   compute tsample(i)=%nobs
   *
   * Sweep out the constant (if necessary) and augmenting lags
   *
   if %ncmom>=3
      compute sweepxx=%sweeplist(sweepxx,%seq(3,%ncmom))
   *
   * Get the residuals from projections on the augmenting variables.
   * (There might be none, hence the check). These drop another data point.
   *
   if %ncmom>=3 {
      set ehat i//(pchoice(i)+3) i//tend = dy + %dot(%xsubmat(sweepxx,1,1,3,%ncmom),%eqnxvector(extras,t))
      set vhat i//(pchoice(i)+3) i//tend = y1 + %dot(%xsubmat(sweepxx,2,2,3,%ncmom),%eqnxvector(extras,t))
   }
   else {
      set ehat i//(pchoice(i)+3) i//tend = dy
      set vhat i//(pchoice(i)+3) i//tend = y1
   }
   *
   * Breitung's Gauss code uses a variance estimator which doesn't
   * correct for degrees of freedom, so we'll do the same.
   *
   compute sigmasq(i)=sweepxx(1,1)/(%nobs-1)
   *
   * Forward de-mean the ehat if DET=TREND
   *
   if det==3
      panel(smpl=smpl,gls=forward,effects=indiv,vrandom=1.0,vindiv=%na) ehat i//(pchoice(i)+3) i//tend
   *
   * Standardize the projection residuals
   *
   set ehat i//(pchoice(i)+3) i//tend = ehat/sqrt(sigmasq(i))
   set vhat i//(pchoice(i)+3) i//tend = vhat/sqrt(sigmasq(i))
end do i
if robust {
   linreg(smpl=smpl,robust) ehat startl end
   # vhat
   compute %cdstat=%tstats(1)
}
else {
   cmom(smpl=smpl) startl endl
   # ehat vhat
   compute %cdstat=%cmom(1,2)/sqrt(%cmom(2,2))
}
compute tnobs=%nobs
compute %signif=%cdf(%cdstat)
compute %ngroup=fix(nt)
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Breitung Unit Root Test: Series "+%l(y)

report(use=breport,action=define,title=ltitle)
report(use=breport,atrow=1,atcol=1,span) ltitle
report(use=breport,atrow=2,atcol=1,span) "Test has large N, large T (sequential) asymptotics"
report(use=breport,atrow=3,atcol=1,span) "Null is Unit Root. Alternative is Common Stationary Root."
report(use=breport,atrow=4,atcol=1,span) "Individual Specific Components: "+descript
report(use=breport,atrow=5,atcol=1,span) lstyle
report(use=breport,atrow=7,atcol=1) "N"      fix(nt)
report(use=breport,atrow=8,atcol=1) "T"      fix(pt)
report(use=breport,atrow=9,atcol=1) "Test Statistic" %cdstat
report(use=breport,atrow=10,atcol=1) "Signif" %signif
report(use=breport,action=format,width=8,atcol=2)
if print
   report(use=breport,action=show)
end

