*
* @LPUnit(options) series start end
*
* Implements the Lumsdaine-Papell unit root test, allowing for two
* breaks in the intercept, the trend or both at unknown locations. The
* base model is
*
*  dy(t) = a + bt + r*y(t-1) + sum a(i) * dy(t-i)
*
* The break variables are c I(t>tb) or d max(0,t-tb) or both, where
* the break point tb is unknown and determined by finding the minimum
* value for the t-statistic for r=0. The number of lags of dy included
* can either be input, or can be selected automatically using AIC, BIC
* or general-to-specific pruning by t-tests. All of the automatic lag
* selections are applied to the base model
*
* Options:
*   BREAK=[INTERCEPT]/TREND/BOTH
*     Selects which variables are allowed to have a break
*
*   NBREAKS=0/1/[2] (actually can be any positive integer)
*
*   METHOD=[INPUT]/AIC/BIC/HQ/TTEST/GTOS
*     Selects the method for deciding the number of additional lags. If
*     INPUT, the number of lags given by the MAXLAGS option is used. If AIC,
*     the AIC-minimizing value between 0 and MAXLAGS is used; if BIC, it's
*     the BIC-minimizing value, and if TTEST or GTOS, the number of lags for
*     which the last included lag has a marginal significance level less than the
*     cutoff given by the SIGNIF option.
*   LAGS=number of additional lags (METHOD=INPUT)
*   MAXLAGS=number of additional lags (METHOD=INPUT) or the maximum number of lags
*     to consider (other METHOD's) [number of observations**.25]
*   SIGNIF=cutoff significance level for METHOD=TTEST or GTOS[.10]
*
*   SMPL=SMPL series with entries to skip
*
*   PI=fraction of entries on each end of data to exclude as break points, and minimum
*     gap between breakpoints [.10]
*
*   [PRINT]/NOPRINT
*   TITLE=title for report ["Lumsdaine-Papell Unit Root Test, Series ..."]
*   GRAPH/[NOGRAPH]
*     GRAPH does a graph of the series with grid lines at the break points
*
* 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
*
* References:
*   Lumsdaine and Papell(1997), "Multiple trend breaks and the unit root
*     hypothesis", Review of Economics and Statistics, vol 79, 212-218.
*
*   Zivot and Andrews(1992), "Further evidence on the Great Crash, the oil
*     price shock, and the unit-root hypothesis", JBES, vol. 10, 251-70.
*
* Revision Schedule:
*   08/2009 Written by Tom Doan, Estima
*   06/2011 Add asymptotic critical values, TITLE option.
*   10/2014 Switch to align=decimal on report.
*   05/2017 Correct output of maximum number of lags. Improve the description.
*   05/2017 Correct dummy trend (was 1 at BP, not BP+1)
*   05/2017 Correct t-stat on TREND in standard output
*
procedure LPUnit series start end
type series  series
type integer start end
*
option choice  break     1 intercept trend both
option integer nbreaks   2
option integer lags
option integer maxlags
option choice  method    1 input aic bic hq ttest gtos
option real    signif    .10
option series  smpl
option real    pi        .15
option switch  print     1
option switch  graph     0
option string  title
*
local integer maxlag bestlag startl endl lag znobs
local series  dy trend grid
local real    icmult ic icmin
local integer lower upper pinobs breakcomps done nobs tag pos
local integer lmethod
local real    mint
local vect    bestbeta besttstats critvalues
local string  descript ltitle llagreport llagmethod
local report  lpreport
local rect[series] dt
local vect[int] bps upperbound bestbreaks
*
declare integer %%autop
*
* First block is Zivot-Andrews (one break)
* Second is Lumsdaine-Papell (two breaks)
*
fixed rect[vect] cvtables(2,3)(3)

 -5.34 , -4.80, -4.58,
 -4.93 , -4.42, -4.11,
 -5.57 , -5.08, -4.82,

 -6.74 , -6.16 , -5.89,
 -7.19 , -6.62 , -6.37,
 -7.19 , -6.75 , -6.48
*
* Check for invalid syntax
*
if .not.%defined(series) {
   disp "Syntax: @lpunit( options )  SERIES start end"
   return
}
*
* Figure out the range available for the series and one lag
*
inquire(reglist) startl<<start endl<<end
# series{0 1}
*
set dy    startl endl = series-series{1}
set trend startl endl = t
*
compute lower=startl+maxlags+1
compute upper=endl
*
* Shortest interval between the breaks or between the ends and the breaks
*
compute pinobs=fix(pi*(upper-lower+1))
*
* Determine default value of maxlag
*
compute znobs=endl-startl+1

if %defined(lags)
   compute maxlag=lags
else
if %defined(maxlags)
   compute maxlag=maxlags
else
   compute maxlag=fix(znobs^.25)

compute lmethod=method
if lmethod==6
   compute lmethod=5

if lmethod==1
   compute bestlag=maxlag
else {
   cmom(smpl=smpl) startl+maxlag endl
   # constant trend series{1} dy{0 to maxlag}
   if lmethod==2
      compute icmult=2.0/%nobs
   else
   if lmethod==3
      compute icmult=log(%nobs)/%nobs
   else
      compute icmult=2.0*log(log(%nobs))/%nobs
   linreg(cmom,noprint) dy
   # constant trend series{1}
   compute ic=log(%rss/%nobs)+icmult*%nreg
   compute icmin=ic,bestlag=0
   do lag=maxlag,1,-1
      linreg(cmom,noprint) dy
      # constant trend series{1} dy{1 to lag}
      if lmethod==5 {
         if %ttest(%tstats(3+lag),%ndf)<=signif {
            compute bestlag=lag
            break
         }
      }
      else {
         compute ic=log(%rss/%nobs)+icmult*%nreg
         if ic<icmin
            compute icmin=ic,bestlag=lag
      }
   end do lag
}
if break<>3
   compute breakcomps=1
else
   compute breakcomps=2
*
* RECT[SERIES] for the dummies. By dimensioning this as breakcomps x
* nbreaks (rather than nbreaks x breakcomps), the regressors will be
* grouped by the break point number, since RECT's are expanded going
* down columns.
*
dim dt(breakcomps,nbreaks)
*
* lower and upper bounds of regressions, allowing for lags
*
compute lower=startl+bestlag+1
compute upper=endl
*
* These keep track of the break point information.
*
dim bps(nbreaks) upperbound(nbreaks) bestbreaks(nbreaks)
*
* Set up the starting break points (the leftmost legal values), and the
* upper bounds (rightmost legal values).
*
do i=1,nbreaks
   compute bps(i)=(lower-1)+pinobs*i
   compute upperbound(i)=endl+1-pinobs*(nbreaks+1-i)
end do i
*
* 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 required dummy components
   *
   do i=1,nbreaks
      if break<>2
         set dt(1,i) startl endl = t>bps(i)
      if break<>1
         set dt(1+(break==3),i) startl endl = %max(t-bps(i),0)
   end do i
   linreg(noprint,smpl=smpl) dy startl+bestlag endl
   # series{1} dt constant trend dy{1 to bestlag}
   *
   * 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.%tstats(1)<mint
      compute mint=%tstats(1),bestbreaks=bps,bestbeta=%beta,besttstats=%tstats
   *
   * Shift the break points if possible
   *
   compute done=1
   do i=nbreaks,1,-1
      compute bps(i)=bps(i)+1
      if bps(i)>=upperbound(i)
         next
      do j=i+1,nbreaks
         compute bps(j)=bps(j-1)+pinobs
      end do j
      compute done=0
      break
   end do i
}
compute %cdstat=mint
*
set grid startl endl = 0.0
do i=1,nbreaks
   set grid startl endl = grid.or.(t==bestbreaks(i))
end do i
*
if graph {
   graph(footer="Series with Breaks Highlighted",grid=grid)
   # series
}

if nbreaks>=1
   compute %minent=bestbreaks(1)
if nbreaks>=2
   compute %maxent=bestbreaks(2)
compute %%autop=bestlag


if print==0
   return
if nbreaks==0 {
   *
   * Critical values for zero breaks (Dickey-Fuller)
   *
   compute nobs=%nobs+1
   compute critvalues=||-3.9638-8.353/nobs-47.44/(nobs^2),$
                        -3.4126-4.039/nobs-17.83/(nobs^2),$
                        -3.1279-2.418/nobs-7.58/(nobs^2)||
}
else
if nbreaks<=2
   compute critvalues=cvtables(nbreaks,break)
else
   compute critvalues=||%na,%na,%na||

if %defined(title)
   compute ltitle=title
else
   compute ltitle="Lumsdaine-Papell Unit Root Test, Series "+%l(series)
report(use=lpreport,action=define,title=ltitle)
report(use=lpreport,atrow=1,atcol=1,span) ltitle
report(use=lpreport,atrow=2,atcol=1,span) "Regression Run From "+%datelabel(startl+bestlag+1)+" to "+%datelabel(endl)
report(use=lpreport,atrow=3,atcol=1) "Observations" %nobs
if nbreaks==0
   compute descript="None -- Dickey-Fuller Test"
else
if break==1
   compute descript="Intercept Only"
else
if break==2
   compute descript="Trend Only"
else
   compute descript="Intercept and Trend"
report(use=lpreport,atrow=4,atcol=1,span) "Breaks in "+descript

compute descript="Breaks at "
do i=1,nbreaks
   compute descript=descript+%datelabel(bestbreaks(i))+" "
end do i
*
if nbreaks>0
   report(use=lpreport,atrow=5,atcol=1,span) descript

if lmethod==1
   compute llagmethod="User"
else
if lmethod==2
   compute llagmethod="AIC"
else
if lmethod==3
   compute llagmethod="BIC"
else
if lmethod==4
   compute llagmethod="HQ"
else
if lmethod==5
   compute llagmethod="GTOS/t-tests("+%strval(signif,"#.##")+")"

if lmethod==1
   compute llagreport="Using fixed lags "+bestlag
else
   compute llagreport="With "+bestlag+" lags chosen from "+maxlag

report(use=lpreport,atrow=6,atcol=1,span) llagreport
report(use=lpreport,atrow=7,atcol=1,span) "Selected by "+llagmethod

report(use=lpreport,atrow=8,atcol=1) "Sig Level" "Crit Value"
report(use=lpreport,atrow=9,atcol=1) "1%(**)"    critvalues(1)
report(use=lpreport,atrow=10,atcol=1) "5%(*)"     critvalues(2)
report(use=lpreport,atrow=11,atcol=1) "10%"      critvalues(3)

report(use=lpreport,atrow=13,atcol=1) "Variable"
report(use=lpreport,atrow=13,atcol=2,align=center) "Coefficient" "T-Stat"
report(use=lpreport,atrow=14,atcol=1) "Y{1}" bestbeta(1) besttstats(1)
compute tag=fix(1+(%cdstat<critvalues(1))+(%cdstat<critvalues(2)))
report(use=lpreport,action=format,atrow=15,atcol=3,special=tag)
compute pos=2
do i=1,nbreaks
   if break<>2 {
   	report(use=lpreport,row=new,atcol=1) "D("+%datelabel(bestbreaks(i))+")" $
        bestbeta(pos) besttstats(pos)
      compute pos=pos+1
   }
   if break<>1 {
      report(use=lpreport,row=new,atcol=1) "DT("+%datelabel(bestbreaks(i))+")" $
        bestbeta(pos) besttstats(pos)
      compute pos=pos+1
   }
end do i
report(use=lpreport,row=new,atcol=1) "Constant" bestbeta(pos)   besttstats(pos)
report(use=lpreport,row=new,atcol=1) "Trend"    bestbeta(pos+1) besttstats(pos+1)
report(use=lpreport,action=format,picture="*.####",atrow=8)
report(use=lpreport,action=format,align=decimal,atcol=3,atrow=8)
report(use=lpreport,action=show)
end
