*
* @LevinLin( options ) series start end
* Implements tests from Levin, Lin and Chu(2002), "Unit root tests in
* panel data: Asymptotic and finite-sample properties", Journal of
* Econometrics, vol 108, no 1, 1–24.
*
* Requires version 7.30 or later of RATS due to the TITLE option on
* REPORT(ACTION=DEFINE). That just affects the title of the window if
* you reload it with the Window--Restore Report operation. If you remove
* that, it will work with 7.1 or later.
*
* 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]
*
*   LWINDOW=[NEWEY]/BARTLETT/FLAT/PARZEN/QUADRATIC
*   BANDWIDTH=# of lags (or bandwidth for LWINDOW=QUADRATIC) [Schwert's]
*   RECOLOR/[NORECOLOR]
*     These choose the method used for computing the ratio of long- to short-
*     run variances. RECOLOR does this by using inverting the autoregressive
*     representation. By default, this is done using a lag window calculation
*     applied to the residuals from the differenced series on the deterministics
*     with lag window type and bandwidth determined by the LWINDOW and BANDWIDTH
*     options.
*
*   [PRINT]/NOPRINT
*   TITLE=title for output ["Levin-Lin Unit Root Test Series: ..."]
*
* Variables Defined:
*   %CDSTAT - the normalized t-statistic
*   %SIGNIF - one tailed significance level
*   %NGROUP - number of individuals
*   %%ILAGS - lags used in individual analysis
*
* Revision Schedule:
*  03/2012 Written by Tom Doan, Estima
*
procedure LevinLin 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 choice     lwindow  	2 newey bartlett flat parzen quadratic
option real       bandwidth
option switch     recolor   	0
option series  	smpl
option switch  	print    	1
option string     title
*
decl  vect[integer] %%ilags
*
local integer 		startl endl
local integer 		i tnobs tlookup
local series 		dy icount ehat vhat dyhat
local equation 	extras
local integer   	p pmax tmax tend starti endi lagsused slot
local series    	itrend
local vect[int] 	detvars pchoice
local rect      	sweepxx
local real      	snhat tavg xlookup mustar sigmastar
local real      	marginal_t ic bestic logl
local real        pt nt avglag klag
local vect      	sigmasq sigsqratio tsample
local report      lreport
local string      ltitle descript wstyle lstyle
local vector      lagcoeffs
*
* Labels for lag window estimation
*
fixed vect[string] lwindowtext(5)
Bartlett(Newey-West)
Bartlett(Newey-West)
Flat
Parzen
Quadratic Spectral
*
* Levin-Lin, Table 2, mean and standard deviation adjustments
*
* Individual record sample sizes
*
fixed vect[int] ss(13)
 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 250, 999999
*
* Mean (first row) and standard deviation (second row) adjustments for
* the three types of deterministics and various sample sizes.
*
fixed vect[rect[real]] msda(3)(2,13)

 0.004, 0.003, 0.002, 0.002 , 0.001, 0.001, 0.001, 0,         0,     0,     0,     0,     0
 1.049, 1.035, 1.027, 1.021,  1.017, 1.014, 1.011, 1.008, 1.007, 1.006, 1.005, 1.001, 1.000

 -.554,-.546,-.541,-.537,-.533,-.531,-.527,-.524,-.521,-.520,-.518,-.509,-.500
  .919, .889, .867, .850, .837, .826, .810, .798, .789, .782, .776, .742, .707

 -.703,-.674,-.653,-.637,-.624,-.614,-.598,-.587,-.578,-.571,-.566,-.533,-.500
 1.003, .949, .906, .871, .842, .818, .780, .751, .728, .710, .695, .603, .500

if .not.%defined(y) {
   disp "###Syntax: @LEVINLIN( options ) >>>>>SERIES<<<<<"
   return
}
inquire(series=y) startl<<start endl<<end
*
* Compute first differenced data
*
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) sigsqratio(endi) tsample(endi)
*
* Create the individual trend series in case we need it.
*
set itrend = %period(t)
*
* Create the set of regressors.
*
if det==1 {
   dim detvars(0)
   compute descript="None"
}
else
if det==2 {
   compute detvars=||constant||
   compute descript="Constant"
}
else {
   compute detvars=||constant,itrend||
   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
*
* Determine choice for calculation of long-run variance
*
if %defined(bandwidth)
   compute klag=bandwidth
else
   compute klag=fix(%round(4.0*(.01*pt)^(2.0/9.0),0))
if recolor
   compute wstyle="Long-run variances by recolored AR"
else
   compute wstyle="Long-run variances by "+lwindowtext(lwindow)+"("+fix(klag)+")"
*
clear ehat vhat
do i=starti,endi
   if icount(i)==0
      next
   equation extras *
   # detvars dy{1 to pchoice(i)}
   cmom(smpl=smpl) i//(pchoice(i)+2) i//tend
   # dy y{1} detvars dy{1 to pchoice(i)}
   compute tsample(i)=%nobs
   *
   * First sweep out the deterministics to get the dyhat for computing
   * the long-run variance.
   *
   if %rows(detvars)==0 {
      set dyhat i//2 i//tend = dy
      compute sweepxx=%cmom
   }
   else {
      compute sweepxx=%sweeplist(%cmom,%seq(3,2+%rows(detvars)))
      set dyhat i//2 i//tend = dy + $
         %dot(%xsubmat(sweepxx,1,1,3,2+%rows(detvars)),%xsubvec(%eqnxvector(extras,t),1,%rows(detvars)))
   }
   *
   * Next, sweep out the augmenting lags
   *
   if pchoice(i)>0
      compute sweepxx=%sweeplist(sweepxx,%seq(3+%rows(detvars),%ncmom))
   *
   * Save (minus) the coefficients on the augmenting lags
   *
   if pchoice(i)>0
      compute lagcoeffs=%xsubmat(sweepxx,1,1,%ncmom-pchoice(i)+1,%ncmom)
   *
   * Get the residuals from projections on the augmenting variables.
   * (There might be none, hence the check).
   *
   if %ncmom>=3 {
      set ehat i//(pchoice(i)+2) i//tend = dy  +%dot(%xsubmat(sweepxx,1,1,3,%ncmom),%eqnxvector(extras,t))
      set vhat i//(pchoice(i)+2) i//tend = y{1}+%dot(%xsubmat(sweepxx,2,2,3,%ncmom),%eqnxvector(extras,t))
   }
   else {
      set ehat i//(pchoice(i)+2) i//tend = dy
      set vhat i//(pchoice(i)+2) i//tend = y{1}
   }
   *
   * Then sweep out the lag so we can get the ADF variance estimate.
   *
   compute sweepxx=%sweep(sweepxx,2)
   compute sigmasq(i)=sweepxx(1,1)/(%nobs-%eqnsize(extras)-1)
   *
   * And standardize the projection residuals
   *
   set ehat i//(pchoice(i)+2) i//tend = ehat/sqrt(sigmasq(i))
   set vhat i//(pchoice(i)+2) i//tend = vhat/sqrt(sigmasq(i))
   *
   if recolor {
      *
      * Use the coefficients from the AR regression. (LAGCOEFFS
      * has - the augmenting lag coefficients).
      *
      if pchoice(i)==0
         compute sigsqratio(i)=1.0
      else
         compute sigsqratio(i)=1.0/(1.0+%sum(lagcoeffs))^2
   }
   else {
      *
      * Use a lag window applied to the DYHAT series
      *
      mcov(lwindow=lwindow,lags=klag,smpl=smpl) i//2 i//tend
      # dyhat
      compute sigsqratio(i)=(%cmom(1,1)/%nobs)/sigmasq(i)
   }
end do i
*
* Compute the regression on the standardized series.
*
linreg(smpl=smpl,noprint) ehat startl endl
# vhat
compute tnobs=%nobs
*
sstats(mean,smpl=(icount>0)) starti endi sqrt(sigsqratio(t))>>snhat tsample(t)>>tavg
sstats(min)  1 %cols(ss) %if(tavg<ss(t),t,%na)>>xlookup
compute tlookup=fix(xlookup)+1
compute mustar=msda(det)(1,tlookup)
compute sigmastar=msda(det)(2,tlookup)
*
compute %cdstat=(%tstats(1)-tnobs*snhat/%seesq*%stderrs(1)*mustar)/sigmastar
compute %signif=%cdf(%cdstat)
compute %ngroup=fix(nt)
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Levin-Lin Unit Root Test: Series "+%l(y)

report(use=lreport,action=define,title=ltitle)
report(use=lreport,atrow=1,atcol=1,span) ltitle
report(use=lreport,atrow=2,atcol=1,span) "Test has large N, N/T-->0"
report(use=lreport,atrow=3,atcol=1,span) "Null is Unit Root. Alternative is Common Stationary Root."
report(use=lreport,atrow=4,atcol=1,span) "Individual Specific Components: "+descript
report(use=lreport,atrow=5,atcol=1,span) lstyle
report(use=lreport,atrow=6,atcol=1,span) wstyle
report(use=lreport,atrow=8,atcol=1) "N"      fix(nt)
report(use=lreport,atrow=9,atcol=1) "T"      fix(pt)
report(use=lreport,atrow=10,atcol=1) "t-unadjusted" %tstats(1)
report(use=lreport,atrow=11,atcol=1) "t-adjusted" %cdstat
report(use=lreport,atrow=12,atcol=1) "Signif" %signif
report(use=lreport,action=format,width=8,atcol=2)
if print
   report(use=lreport,action=show)

dim %%ilags(%ngroup)
compute slot=0
do i=starti,endi
   if icount(i)==0
      next
   compute slot=slot+1,%%ilags(slot)=pchoice(i)
end do i

end

