*
* @LSUnit( options )  series start end
*
* Implementation of
*
* (for 1 break)
* Lee and Strazicich(2004), "Minimum LM Unit Root Test with One Structural
* Break", Appalachian State University Working Paper,
* http://econ.appstate.edu/RePEc/pdf/wp0417.pdf
*
* (for 2 breaks)
* Lee and Strazicich(2003), "Minimum LM Unit Root Test with Two Structural
* Breaks," Review of Economics and Statistics, vol 85, no. 4, pp 1082-1089.
*
* (for no breaks)
* Schmidt and Phillips(1992), "LM Tests for a Unit Root in the Presence
* of Deterministic Trends." Oxford Bulletin of Economics and Statistics,
* 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,bestbreaks=bps,bestlag=lagsused,bestxx=sweepxx
   *
   * Update the break points. Add one to the final slot, until we hit
   * its upper bound. When we hit that, add one to the second to last
   * slot, and reset the final slot to its smallest possible value. When
   * the second to last slot is exhausted, add one to the one earlier,
   * etc. If <<done>> is zero at the end, we still have a valid
   * combination that we haven't checked.
   *
   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



