*
* @PerronBreaks( options )  series start end
*
* Implementation of the general structure for analyzing breaks with unit
* roots in Perron(2006), "Dealing with Structural Breaks," in Palgrave
* Handbook of Econometrics, Vol. 1, pp 278-352, extended to allow more
* than one break.
*
* Parameters:
*   series      Series to test for unit roots
*   start end   Range of series to use
*
* Options:
*   AO=[NONE]/MEAN/CRASH/JOIN/BREAK
*   IO=[NONE]/MEAN/CRASH/BREAK
*
*   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/AIC/BIC/HQ
*   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
*   TITLE=title for report ["Unit Root Test, Series x"]
*
*  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
*    %%SHIFTS=RECT of estimated shifts (one column for each break).
*    %%AUTOP=number of lags selected
*
*  Revision Schedule:
*   06/2010 Written by Tom Doan, Estima
*   09/2011 Correct combination of AO and BREAKS>1. (Lag selection didn't
*           work correctly).
*   05/2017 Correct degrees of freedom calculations for AO models. (Didn't
*           allow for the fixed part removed in the detrending regression).
*   05/2017 Improve output of information regarding choice of lags.
*
*  References:
*    The two break case is used in David H. Papell & R. Prodan (2006)
*    "Additional Evidence of Long Run Purchasing Power Parity with
*    Restricted Structural Change", Journal of Money, Credit and Banking,
*    vol. 38(5), 1329-1349.
*
*************************************************************************
* 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 %%URGToS sweepxx lagsused
type rect    *sweepxx
type integer *lagsused
*
option integer lags    0
option choice  method  1  fixed gtos  aic  bic hq
option real    slstay  .10
option integer nbase
option integer ndrop   1		;* Number of regressors to drop when eliminating a lag
*
local integer lagsused testlag slot
local real    marginal_t bestic logl ic
*
if method==1 {
   compute lagsused=lags
   return
}
*
* Set lagsused to 0, which will be the value if we don't find anything
* significant at 1 or higher.
*
compute lagsused=0
*
* Stop at testlag=1 if we're doing GTOS; otherwise run down to testlag=0
*
do testlag=lags,(method==2),-1
   *
   * Number of regressors is <<testlag>>*<<ndrop>>+<<nbase>>+1. That
   * allows for the deterministics plus lagged dependent.
   *
   compute slot=nbase+2+testlag*ndrop
   compute %ndf=%nobs-(slot-1)
   if method==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=testlag
         break
      }
   }
   else
   if method==3.or.method==4.or.method==5 {
      compute logl=-.5*%nobs*log(sweepxx(1,1)/%nobs)
      if method==3
         compute ic=2.0*logl-2*(testlag*ndrop)
      else
      if method==4
         compute ic=2.0*logl-log(%nobs)*(testlag*ndrop)
      else
         compute ic=2.0*logl-2*log(log(%nobs))*(testlag*ndrop)
      if testlag==lags.or.ic>bestic
         compute bestic=ic,lagsused=testlag
   }
   *
   * Otherwise, sweep on the rejected lag to take it out of the
   * regression.
   *
   compute sweepxx=%sweeplist(sweepxx,%seq(slot-ndrop+1,slot))
end do testlag
end
************************************************************************
procedure PerronBreaks series start end
type series series
type integer start end
*
option choice  ao      1 	none  mean crash  join  break
option choice  io      1 	none  mean crash  break
option integer breaks  1
option integer lags    0
option choice  method  1  	fixed gtos  aic  bic  hq
option real    slstay  .10
option real    pi      .10
option switch  print   1
option string  title
*
local integer     startl endl pinobs lower upper
local integer     i j k next
local integer     nperbreak nbase lagsused bestlag bestndf
local real   		   tstat mint
local series 		   u du trend
local rect 			   sweepxx bestxx
local integer     done
local rect[series] dt
local vect[series] aoshift
local vect[int]   bps upperbound bestbreaks fixed aolags
local string      clabel
local string      llagmethod llagreport
local integer     spike meanbreak trendbreak trending
local string      ltitle
local report      breport
local integer     nperlag
local vector      thisShift bestshift
local rect[int]   aotable aochosenlags
dec   rect        %%shifts
*
inquire(series=series) startl<<start endl<<end
*
* 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))
*
* 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
*
if io>=2.and.ao>=2 {
   disp "###PerronBreaks-IO and AO options can't be used together"
   return
}

if io==2
   compute spike=1,meanbreak=0,trendbreak=0,trending=0,clabel="IO Model: Mean"
else
if io==3
   compute spike=1,meanbreak=1,trendbreak=0,trending=1,clabel="IO Model: Crash"
else
if io==4
   compute spike=1,meanbreak=1,trendbreak=1,trending=1,clabel="IO Model: Full Trend Break"
else
if ao==2
   compute spike=1,meanbreak=1,trendbreak=0,trending=0,clabel="AO Model: Mean"
else
if ao==3
   compute spike=1,meanbreak=1,trendbreak=0,trending=1,clabel="AO Model: Crash"
else
if ao==4
   compute spike=0,meanbreak=0,trendbreak=1,trending=1,clabel="AO Model: Joined Trend"
else
   compute spike=1,meanbreak=1,trendbreak=1,trending=1,clabel="AO Model: Full Trend Break"
*
* 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.
*
* For AO models, the "spike"'s are used in the final regression, not the
* detrending regression.
*
dim %%shifts(meanbreak+trendbreak,breaks)
if ao>=2 {
   compute nperbreak=meanbreak+trendbreak
   compute nbase    =0
   dim aoshift(breaks)
}
else {
   compute nperbreak=meanbreak+trendbreak+spike
   compute nbase    =1+(trending)+breaks*nperbreak
}
dim dt(nperbreak,breaks)
*
* Create the fixed regressors
*
set trend = t
if trending==0
   compute fixed=||constant||
else
   compute fixed=||constant,trend||
*
* Create the lag list for additive outliers. This needs to include the
* spike dummies for all models except the joined trend.
*
compute aolags=%rlempty()
if ao>=2 {

   if meanbreak {
      do i=1,breaks
         compute aolags=%rladdone(aolags,aoshift(i))
      end do i
      compute nperlag=breaks+1
   }
   else
      compute nperlag=1

   do k=1,lags
      if meanbreak {
         do i=1,breaks
            compute aolags=%rladdlag(aolags,aoshift(i),k)
         end do i
      }
      compute aolags=%rladdlag(aolags,du,k)
   end do k
}
*
* 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,breaks
      if ao>=2 {
         set aoshift(i) startl endl = (t==bps(i)+1)
         compute next=1
      }
      else {
         set dt(1,i) startl endl    = (t==bps(i)+1)
         compute next=2
      }
      if meanbreak {
         set dt(next,i) startl endl = (t>=bps(i)+1)
         compute next=next+1
      }
      if trendbreak
         set dt(next,i) startl endl = %max(0,t-bps(i))
   end do i
   *
   * For additive outliers, detrend first. We save into NBASE the number
   * of degrees of freedom lost in this for correction of degrees of
   * freedom calculations later.
   *
   if ao>=2 {
      linreg(noprint) series startl endl u
      # fixed dt
      compute nbase=%nreg
      compute thisShift=%xsubvec(%beta,%rows(fixed)+1,%nreg)
   }
   else
      set u startl endl = series

   set du = u-u{1}

   if io>=2 {
      *
      * IO models
      *
      cmom
      # du u{1} fixed dt du{1 to lags}
      *
      * Use %sweeplist to regress on the lags of stilde and the
      * deterministics.
      *
      compute sweepxx=%sweeplist(%cmom,%seq(2,%ncmom))
      @%%URGToS(lags=lags,slstay=slstay,method=method,nbase=nbase) sweepxx lagsused
      *
      * Redo the regression with the chosen number of lags and maximal range
      *
      cmom
      # du u{1} fixed dt du{1 to lagsused}
      compute sweepxx=%sweeplist(%cmom,%seq(2,1+nbase+lagsused))
      linreg(noprint) du
      # u{1} fixed dt du{1 to lagsused}
      compute %ndf=%nobs-lagsused-(nbase+1)
      compute thisshift=%xsubvec(%beta,2+%size(fixed),1+%size(fixed)+%size(dt))
   }
   else {
      cmom
      # du u{1} aolags
      compute sweepxx=%sweeplist(%cmom,%seq(2,%ncmom))
      @%%URGToS(lags=lags,slstay=slstay,nbase=nbase,method=method,ndrop=nperlag) sweepxx lagsused
      *
      * Redo the regression with the chosen number of lags and maximal range
      *
      compute aotable=%tablefromrl(aolags)
      dim aochosenlags(2,breaks+nperlag*lagsused)
      ewise aochosenlags(i,j)=aotable(i,j)
      cmom
      # du u{1} %rlfromtable(aochosenlags)
      compute sweepxx=%sweeplist(%cmom,%seq(2,%ncmom))
      *
      * Correct the degrees of freedom for the deterministic components removed in
      * detrending, plus the various lags and dummies, plus the lagged dep var.
      *
      compute %ndf=%nobs-nbase-(breaks+nperlag*lagsused)-1
   }
   *
   * Compute the t-stat (not absolute) on the first lag of u
   *
   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
      compute bestbreaks=bps
      compute bestlag=lagsused
      compute bestxx=sweepxx
      compute bestshift=thisshift
      compute bestndf=%ndf
   }
   *
   * 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   =bestndf
compute %cdstat=mint
compute %minent=bestbreaks(1)
compute %maxent=bestbreaks(2)
compute %%shifts=%reshape(bestshift,%size(bestshift)/breaks,breaks)
compute %%autop=bestlag
*
dim %beta(1+nbase) %tstats(1+nbase)
ewise %beta(i)=bestxx(i+1,1)
ewise %tstats(i)=%beta(i)/sqrt(bestxx(1,1)*bestxx(i+1,i+1)/%ndf)
*
if .not.print
   return

if %defined(title)
   compute ltitle=title
else
   compute ltitle="Unit Root Test, Series "+%l(series)

report(use=breport,action=define,title=ltitle)
report(use=breport,atrow=1,atcol=1,span) ltitle
report(use=breport,atrow=2,atcol=1,span) "Regression Run From "+%datelabel(startl+lags+1)+" to "+%datelabel(endl)
report(use=breport,atrow=3,atcol=1) "Observations" %nobs

report(use=breport,atrow=4,atcol=1,span) clabel+" with "+breaks+" breaks"

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

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

report(use=breport,atrow=5,atcol=1,span) llagreport
report(use=breport,atrow=6,atcol=1,span) "Selection by "+llagmethod
report(use=breport,atrow=8,atcol=1) "Variable" "Coefficient" "T-Stat"
report(use=breport,atrow=9,atcol=1) "Y{1}" %beta(1) %tstats(1)
if breaks>0 {
   report(use=breport,atrow=11,atcol=1) "Break point(s)"
   do i=1,breaks
      report(use=breport,row=new,atcol=1) %datelabel(bestbreaks(i))
   end do i
}
report(use=breport,action=show)
end PerronBreaks
