*
* @LSUnit( options )  series start end
*
* Implementation of
*
* (for 1 break)
* Lee and Strazicich(2013), "Minimum LM Unit Root Test with One Structural
* Break", Economics Bulletin, vol. 33, no 4, 2493-2492.(4), pp 2483-2492.
*
* (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. This allows for any positive number, but note that
*    tables of critical values exist only for 0, 1 and 2. The programming
*    for 3 or more breaks is effectively the same as for 2, which is why
*    we allow for it. Whether there is any actual value in that many breaks
*    is unclear.
*
*   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]
*   THIN=entry spacing on the grid of data points [1]
*     With the default of THIN=1, every possible set of breaks that are separated as
*     given by the PI option are examined. As described above, this can take a great
*     deal of time on a large data set if you're looking for 2 or more breaks. A
*     higher value for the THIN option (for example THIN=5) can reduce that time by
*     only looking at every 5th entry for each possible break.
*
*   [PRINT]/NOPRINT
*   TITLE=title for report ["Lee-Strazicich Unit Root Test, Series ..."
*
*  Variables Defined:
*    %CDSTAT=test statistic (minimum t-stat)
*    %MINENT=entry of the (lower) break giving the minimum t
*    %MAXENT=entry of the (higher) break for the 2 break model
*    %%AUTOP=number of lags chosen at the break values
*
*  Revision Schedule:
*   02/2008 Written by Tom Doan, Estima
*   03/2008 Corrected calculation, added output of dummies, cleaned up
*           handling of multiple breaks
*   07/2010 Add comment regarding lack of tables for three+ breaks
*   07/2012 Add progress box when breaks>0
*   05/2014 Define %%AUTOP
*   07/2016 Add SIGNIF as an alternative to SLSTAY.
*   10/2017 Add THIN option
*   10/2017 Add updated and expanded critical values
*
*************************************************************************
* 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 nbase
*
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>>+<<nbase>>+1.
   *
   compute %ndf=%nobs-testlag-(nbase+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=nbase+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    signif
option real    pi      .10
option switch  print   1
option string  title
option integer thin    1
*
declare integer   %%autop
*
local integer     startl endl pinobs lower upper
local integer     i j last1
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
local string      ltitle
local report      lreport
local real        gtossignif
local real        lambda1 lambda2 temp w1 w2 tw1 tw2 tnobs best1 best2
local integer     i1 i2 t1 t2
*
* For computing interpolated critical values
*
local vect        cv1(3) cv2(3) cv(3)
*
* Sample sizes for critical values
*
fixed vect nlookup(4)
 100 250 500 1000
*
* Critical values for model=break (cv's depend upon size of data set AND
* (approximate) location of break).
*
* Break points for one-break table
*
fixed vect bplookup1(9)
 .10 .15 .20 .25 .30 .35 .40 .45 .50
*
* Break points for two-break table. Combinations of lambda which sum to
* over 1.0 reflect back to 1-lambda2,1-lambda1.
*
fixed vect bplookup2(7)
 .20 .30 .40 .50 .60 .70 .80
*
* Zero break critical values
*
fixed vect[vect] cv0(4)(3)
 -3.597 -3.031 -2.745 ;* T=100
 -3.572 -3.023 -2.747 ;* T=250
 -3.570 -3.021 -2.748 ;* T=500
 -3.566 -3.023 -2.748 ;* T=1000
*
* One break MODEL=crash critical values
*
fixed vect[vect] cvc1(4)(3)
 -4.084 -3.487 -3.185 ;* T=100
 -3.987 -3.387 -3.076 ;* T=250
 -3.840 -3.277 -2.985 ;* T=500
 -3.798 -3.230 -2.925 ;* T=1000
*
* Two break MODEL=CRASH critical values
*
fixed vect[vect] cvc2(4)(3)
 -4.073 -3.563 -3.296 ;* T=100
 -4.101 -3.594 -3.345 ;* T=250
 -4.261 -3.647 -3.287 ;* T=500
 -3.828 -3.248 -3.005 ;* T=1000
*
* One break MODEL=BREAK critical values
*
fixed vect[rect] cvb1(4)(9,3)
-4.630	-4.064	-3.787
-4.704	-4.132	-3.843
-4.769	-4.199	-3.906
-4.820	-4.253	-3.963
-4.857	-4.293	-4.008
-4.891	-4.324	-4.042
-4.899	-4.338	-4.058
-4.910	-4.348	-4.071
-4.915	-4.351	-4.073 ;* T=100

-4.523	-3.974	-3.697
-4.536	-3.982	-3.702
-4.563	-4.010	-3.726
-4.609	-4.060	-3.777
-4.647	-4.098	-3.819
-4.667	-4.125	-3.846
-4.672	-4.130	-3.856
-4.681	-4.144	-3.870
-4.685	-4.145	-3.873 ;* T=250

-4.502	-3.959	-3.679
-4.489	-3.944	-3.664
-4.496	-3.950	-3.671
-4.534	-3.991	-3.714
-4.570	-4.031	-3.756
-4.587	-4.052	-3.777
-4.598	-4.065	-3.794
-4.612	-4.081	-3.809
-4.602	-4.074	-3.799 ;* T=500

-4.466	-3.928	-3.651
-4.455	-3.919	-3.640
-4.469	-3.933	-3.651
-4.501	-3.968	-3.687
-4.523	-3.993	-3.713
-4.546	-4.018	-3.740
-4.562	-4.033	-3.756
-4.576	-4.049	-3.773
-4.593	-4.065	-3.790 ;* T=1000
*
* Two break MODEL=BREAK critical values
*
fixed vect[rect] cvb2(4)(12,3)
-6.750 -6.108 -5.779
-7.196 -6.312 -5.893
-6.932 -6.175 -5.825
-7.004 -6.185 -5.828
-6.691 -6.152 -5.798
-6.821 -5.917 -5.541
-6.963 -6.201 -5.890
-6.821 -6.166 -5.832
-6.978 -6.288 -5.998
-7.032 -6.375 -6.011
-6.863 -6.268 -5.956
-7.014 -6.446 -6.072 ;* T=100

-5.667 -5.177 -4.921
-5.974 -5.342 -5.004
-5.869 -5.398 -5.112
-6.064 -5.462 -5.192
-5.957 -5.374 -5.104
-5.783 -5.272 -4.943
-5.934 -5.300 -4.988
-5.771 -5.357 -5.032
-5.958 -5.421 -5.179
-5.751 -5.349 -5.063
-5.971 -5.264 -5.018
-6.035 -5.484 -5.161 ;* T=250

-5.640 -4.855 -4.572
-5.392 -4.939 -4.691
-5.608 -4.947 -4.661
-5.602 -5.059 -4.782
-5.533 -4.989 -4.752
-5.619 -5.038 -4.723
-5.511 -4.898 -4.607
-5.494 -5.035 -4.742
-5.548 -5.043 -4.811
-5.716 -5.149 -4.912
-5.509 -4.963 -4.674
-5.655 -4.967 -4.705 ;* T=500

-5.116 -4.539 -4.195
-5.240 -4.578 -4.278
-5.437 -4.917 -4.533
-5.387 -4.884 -4.567
-5.099 -4.716 -4.481
-5.247 -4.776 -4.472
-5.085 -4.576 -4.301
-5.209 -4.697 -4.398
-5.212 -4.734 -4.534
-5.183 -4.773 -4.535
-5.100 -4.699 -4.364
-5.291 -4.803 -4.481 ;* T=1000
*
compute gtossignif=%if(%defined(signif),signif,slstay)
*
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
if breaks>0 {
   infobox(action=define,progress,lower=bps(1),upper=upperbound(1)) $
     "Lee-Strazicich Unit Root Test"
   compute last1=bps(1)
}
while .not.done {
   if breaks>0.and.bps(1)<>last1 {
      compute last1=bps(1)
      infobox(current=last1)
   }
   *
   * 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=gtossignif,nbase=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)+thin
      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
}
if breaks>0
   infobox(action=remove)
*
compute %ndf   =%nobs-(ndet+1+bestlag)
compute %cdstat=mint
compute %minent=bestbreaks(1)
compute %maxent=bestbreaks(2)
compute %%autop=bestlag
*
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 %defined(title)
   compute ltitle=title
else
   compute ltitle="Lee-Strazicich Unit Root Test, Series "+%l(series)

compute lambda1=float(bestbreaks(1)-lower)/(upper-lower+1)
compute lambda2=float(bestbreaks(2)-lower)/(upper-lower+1)

if breaks==1.and.model==2 {
   *
   * Critical values for MODEL=BREAK depend upon the position of the
   * break. The critical values will be linearly interpolated if they
   * fall between two of the tabled values.
   *
   * Reflect lambda1 through .5 if necessary
   *
   compute lambda1=%min(lambda1,1-lambda1)
   *
   * Figure out the weights on the consecutive values of lambda
   *
   * This handles a break outside the table.
   *
   if lambda1<bplookup1(1)
      compute w1=1.0,w2=0.0,i1=i,i2=i
   else
   if lambda1>bplookup1(%size(bplookup1))
      compute w1=1.0,w2=0.0,i1=%size(bplookup1),i2=i1
   else {
      do i=2,%size(bplookup1)
         if lambda1<bplookup1(i) {
            compute w2=(lambda1-bplookup1(i-1))/(bplookup1(i)-bplookup1(i-1))
            compute w1=1-w2,i1=i-1,i2=i
            break
         }
      end do i
   }
}
else
if breaks==2.and.model==2 {
   *
   * Critical values for a 2-break MODEL=BREAK depend upon the position
   * of both breaks. Unlike the 1-break model, the critical values
   * aren't monotonic, so instead of trying to do two-dimensional
   * interpolation, we simply take the closest match.
   *
   * Reflect lambda's through .5 if necessary
   *
   if lambda1+lambda2>=1.0 {
      compute temp=lambda1
      compute lambda1=1-lambda2
      compute lambda2=1-temp
   }

   compute best1=1.0,best2=1.0
   do i=1,%size(bplookup2)
      if abs(lambda1-bplookup2(i))<best1
         compute i1=i,best1=abs(lambda1-bplookup2(i))
      if abs(lambda2-bplookup2(i))<best2
         compute i2=i,best2=abs(lambda2-bplookup2(i))
   end do i
   *
   * The two can't be the same
   *
   if i1==i2 {
      if i1==1
         compute i2=i1+1
      else
      if i1==%size(bplookup2)
         compute i1=i2-1
      else
      if (lambda1-bplookup2(i1))+(lambda2-bplookup2(i2))>0
         compute i2=i1+1
      else
         compute i1=i2-1
   }
   *
   * Figure out the position in the lookup table
   *
   compute i1=(%size(bplookup2)-i1)*(i1-1)+i2-1
}

compute tnobs=%nobs
*
* Figure out how to linearly interpolate in the observations direction
*
compute tw1=1.0,tw2=0.0,t1=%size(nlookup),t2=%size(nlookup)
do i=1,%size(nlookup)
   if tnobs<nlookup(i) {
      if i==1
         compute tw1=1.0,tw2=0.0,t1=1,t2=1
      else {
         compute tw1=(tnobs-nlookup(i-1))/(nlookup(i)-nlookup(i-1)),t1=i-1
         compute tw2=1-tw1,t2=i
      }
      break
   }
end do i
*
if breaks==0 {
   ewise cv1(i)=cv0(t1)(i)
   ewise cv2(i)=cv0(t2)(i)
}
else
if breaks==1.and.model==1 {
   ewise cv1(i)=cvc1(t1)(i)
   ewise cv2(i)=cvc1(t2)(i)
}
else
if breaks==2.and.model==1 {
   ewise cv1(i)=cvc2(t1)(i)
   ewise cv2(i)=cvc2(t2)(i)
}
else
if breaks==1.and.model==2 {
   ewise cv1(i)=cvb1(t1)(i1,i)*w1+cvb1(t1)(i2,i)*w2
   ewise cv2(i)=cvb1(t2)(i1,i)*w1+cvb1(t2)(i2,i)*w2
}
else
if breaks==2.and.model==2 {
   ewise cv1(i)=cvb2(t1)(i1,i)
   ewise cv2(i)=cvb2(t2)(i1,i)
}
else
if breaks>=3 {
   ewise cv1(i)=%na
   ewise cv2(i)=%na
}
ewise cv(i)=tw1*cv1(i)+tw2*cv2(i)
*
report(use=lreport,action=define,title=ltitle)
report(use=lreport,atrow=1,atcol=1,span) ltitle
report(use=lreport,atrow=2,atcol=1,span) "Regression Run From "+%datelabel(startl+lags+1)+" to "+%datelabel(endl)
report(use=lreport,atrow=3,atcol=1) "Observations" %nobs
if model==1
   compute descript="Crash Model"
else
   compute descript="Trend Break Model"

if breaks==0
   report(use=lreport,atrow=4,atcol=1,span) "No Breaks - Schmidt-Phillips Test"
else
   report(use=lreport,atrow=4,atcol=1,span) descript+" with "+breaks+" breaks"

if method==1
   report(use=lreport,atrow=5,atcol=1,span) "Estimated with fixed lags "+lags
else
   report(use=lreport,atrow=5,atcol=1,span) "With "+bestlag+" lags chosen from "+lags

report(use=lreport,atrow=6,atcol=4,tocol=6,span,align=center) "Critical Values"
report(use=lreport,atrow=7,atcol=1) "Variable" "Coefficient" "T-Stat" ".01" ".05" ".10"
report(use=lreport,atrow=8,atcol=1) "S{1}" %beta(1) %tstats(1) cv(1) cv(2) cv(3)
report(use=lreport,atrow=9,atcol=1) "Constant" %beta(2) %tstats(2)
do i=1,breaks
   report(use=lreport,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(use=lreport,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(use=lreport,action=format,picture="*.####",atrow=8)
if print
   report(use=lreport,action=show)

end LSUnit
