I incorporated a number of suggestions made by Junsoo Lee regarding the output, which includes output of all the coefficients other than the auxiliary lags.
Tom Doan
Estima
- Code: Select all
*
* @LSUnit( options ) series start end
*
* Implementation of
*
* (for 1 break)
* Lee and Strazicich, "Minimum LM Unit Root Test with One Structural Break", 2004
* Appalachian State University Working Paper,
* http://econ.appstate.edu/RePEc/pdf/wp0417.pdf
*
* (for 2 breaks)
* Lee and Strazicich, "Minimum LM Unit Root Test with Two Structural Breaks,"
* Review of Economics and Statistics 2003, vol 85, no. 4, pp 1082-1089.
*
* (for no breaks)
* Schmidt and Phillips, "LM Tests for a Unit Root in the Presence of
* Deterministic Trends." Oxford Bulletin of Economics and Statistics 1992, vol 54,
* pp 257-287.
* Parameters:
* series Series to test for unit roots
* start end Range of series to use
*
* Options:
* MODEL=[CRASH]/BREAK
* Model of structural break. MODEL=CRASH has an abrupt change in level, but no
* change in the trend rate. MODEL=BREAK allows for (simultaneous) changes in level
* and trend.
*
* BREAKS=Number of breaks [1]
* Number of breaks (allows for any positive number)
*
* LAGS=number of lags on the auxiliary regression [1]
* METHOD=[FIXED]/GTOS
* SLSTAY=significance level to keep lag in model with METHOD=GTOS [.10]
* If METHOD=FIXED, the number of lags on the LAGS option is included in all
* regressions. If METHOD=GTOS (general-to-specific), lags are dropped from the
* end so long as the significance level of the t-statistic is greater than the
* SLSTAY value.
*
* PI=fraction of entries on each end of data to exclude as break points, and minimum
* gap between breakpoints [.10]
* [PRINT]/NOPRINT
*
* Variables Defined:
* %CDSTAT=test statistic (minimum t-stat)
* %MINENTRY=entry of the (lower) break giving the minimum t
* %MAXENTRY=entry of the (higher) break for the 2 break model
*
* Revision Schedule:
* 02/2008 Written by Tom Doan, Estima
* 03/2008 Corrected calculation, added output of dummies, cleaned up handling of
* multiple breaks
*
*****************************************************************************************
* Procedure for pruning the number of lags. This is specific to this procedure,
* but can be adapted fairly easily to other situations where lags are being chosen
* using a general-to-specific method.
*
* sweepxx is a RECTANGULAR matrix formed from a cross-product matrix with the
* dependent variable in the first slot, any variables forced into all regressions
* following that, and finally, the lags. That needs to be "swept" on all the
* pivots (all explanatory variables) before calling the procedure.
*
procedure %%LSGToS sweepxx lagsused
type rect *sweepxx
type integer *lagsused
*
option integer lags 0
option real slstay .10
option integer ndet
*
local integer lagsused testlag slot
local real marginal_t
*
* General to specific pruning. Check backwards from <<lags>> to 1
*
* Set lagsused to 0, which will be the value if we don't find anything
* significant at 1 or higher.
*
compute lagsused=0
do testlag=lags,1,-1
*
* Number of regressors is <<testlag>> + <<ndet>>+1. That allows for the
* deterministics plus the the lagged stilde.
*
compute %ndf=%nobs-testlag-(ndet+1)
*
* 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 slot=ndet+2+testlag
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=testlag
break
}
*
* Otherwise, sweep on the rejected lag to take it out of the regression
*
compute sweepxx=%sweep(sweepxx,slot)
end do testlag
end
procedure LSUnit series start end
type series series
type integer start end
*
option choice model 1 crash break
option integer breaks 1
option integer lags 0
option choice method 1 fixed gtos
option real slstay .10
option real pi .10
option switch print 1
*
local integer startl endl pinobs lower upper
local integer i j
local integer ndet lagsused bestlag
local real tstat mint
local series dy ds stilde
local string descript
local rect sweepxx bestxx
local integer done
local rect[series] dt
local vect[int] bps upperbound bestbreaks
*
inquire(series=series) startl<<start endl<<end
*
set dy startl+1 endl = series-series{1}
*
* lower and upper bounds of regressions, allowing for lags
*
compute lower=startl+lags+1
compute upper=endl
*
* shortest interval between the breaks or between the ends and the breaks
*
compute pinobs=fix(pi*(upper-lower+1))
*
* RECT[SERIES] for the dummies. By dimensioning this as model x breaks (rather
* than breaks x models), the regressors will be grouped by the break point number,
* since RECT's are expanded going down columns.
*
dim dt(model,breaks)
*
* These keep track of the break point information.
*
dim bps(breaks) upperbound(breaks) bestbreaks(breaks)
*
* Set up the starting break points (the leftmost legal values), and the upper
* bounds (rightmost legal values).
*
do i=1,breaks
compute bps(i)=(lower-1)+pinobs*i
compute upperbound(i)=endl+1-pinobs*(breaks+1-i)
end do i
*
compute ndet=model*breaks+1
*
* The minimum t value is initialized as an NA, so we know to keep the first value
* we see.
*
compute mint=%na
*
compute done=0
while .not.done {
*
* Create the differences of the required dummy components
*
do i=1,breaks
set dt(1,i) startl endl = t==(bps(i)+1)
if model==2
set dt(2,i) startl endl = t>=(bps(i)+1)
end do i
*
* Run base regression on the deterministic components
*
linreg(noprint) dy startl+1 endl
# constant dt
*
* Generate stilde and its first difference
*
set(first=0.0) stilde startl endl = stilde{1}+%resids
set ds startl+1 endl = stilde-stilde{1}
*
* Compute the cross product matrix including any augmenting lags
*
cmom startl+lags+1 endl
# dy stilde{1} constant dt ds{1 to lags}
*
* Use %sweeplist to regress on the lags of stilde and the deterministics
*
compute sweepxx=%sweeplist(%cmom,%seq(2,%ncmom))
*
* If we're pruning lags, use %%LSGToS to figure out the number of lags
*
if method==2
@%%LSGToS(lags=lags,slstay=slstay,ndet=ndet) sweepxx lagsused
else
compute lagsused=lags
*
* Compute the t-stat (not absolute) on the first lag of stilde
*
compute %ndf=%nobs-lagsused-(ndet+1)
compute tstat=sweepxx(2,1)/$
sqrt(sweepxx(2,2)*sweepxx(1,1)/%ndf)
*
* If this is the smallest t-stat we've seen so far, save the number of
* lags, the break point array and the sweep array.
*
if .not.%valid(mint).or.tstat<mint
compute mint=tstat,bestlag=lagsused,$
bestbreaks=bps,bestxx=sweepxx
compute done=1
do i=breaks,1,-1
compute bps(i)=bps(i)+1
if bps(i)>=upperbound(i)
next
do j=i+1,breaks
compute bps(j)=bps(j-1)+pinobs
end do j
compute done=0
break
end do i
}
*
compute %ndf =%nobs-(ndet+1+bestlag)
compute %cdstat=mint
compute %minent=bestbreaks(1)
compute %maxent=bestbreaks(2)
dim %beta(1+ndet) %tstats(1+ndet)
ewise %beta(i)=bestxx(i+1,1)
ewise %tstats(i)=%beta(i)/sqrt(bestxx(1,1)*bestxx(i+1,i+1)/%ndf)
*
if print {
report(action=define)
report(atrow=1,atcol=1,span) "Lee-Strazicich Unit Root Test, Series "+%l(series)
report(atrow=2,atcol=1,span) "Regression Run From "+%datelabel(startl+lags+1)+" to "+%datelabel(endl)
report(atrow=3,atcol=1) "Observations" %nobs
if model==1
compute descript="Crash Model"
else
compute descript="Trend Break Model"
if breaks==0
report(atrow=4,atcol=1,span) "No Breaks - Schmidt-Phillips Test"
else
report(atrow=4,atcol=1,span) descript+" with "+breaks+" breaks"
if method==1
report(atrow=5,atcol=1,span) "Estimated with fixed lags "+lags
else
report(atrow=5,atcol=1,span) "With "+bestlag+" chosen from "+lags
report(atrow=7,atcol=1) "Variable" "Coefficient" "T-Stat"
report(atrow=8,atcol=1) "S{1}" %beta(1) %tstats(1)
report(atrow=9,atcol=1) "Constant" %beta(2) %tstats(2)
do i=1,breaks
report(atrow=model*(i-1)+10,atcol=1) "D("+%datelabel(bestbreaks(i))+")" $
%beta(model*(i-1)+3) %tstats(model*(i-1)+3)
if model==2
report(atrow=model*(i-1)+11,atcol=1) "DT("+%datelabel(bestbreaks(i))+")" $
%beta(model*(i-1)+4) %tstats(model*(i-1)+4)
end do i
report(action=format,picture="*.####",atrow=8)
report(action=show)
}
end LSUnit
The following runs a simple example
- Code: Select all
source lsunit.src
data(format=free,unit=input,org=cols) 1 38 year x
1960 11508
1961 11369
1962 11930
1963 12636
1964 13194
1965 13264
1966 13875
1967 14193
1968 15172
1969 15712
1970 16316
1971 16872
1972 17563
1973 18489
1974 18606
1975 18741
1976 18728
1977 18480
1978 19303
1979 19682
1980 19968
1981 20271
1982 19464
1983 20127
1984 20712
1985 21438
1986 21719
1987 22443
1988 23040
1989 22995
1990 22226
1991 21921
1992 22367
1993 23176
1994 24144
1995 24818
1996 25077
1997 25608
set logx = log(x)
@lsunit(breaks=2,pi=.10,lags=5,method=gtos,model=crash) logx
@lsunit(breaks=2,pi=.10,lags=5,method=gtos,model=break) logx