*
* @Zivot(options) depvar start end
*
* Implements the Zivot-Andrews unit root test ("Further evidence on the Great Crash, the
* oil price shock, and the unit-root hypothesis." JBES 10 (1992), 251–70), allowing for a
* single break in the intercept, the trend or both. 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>=t0) or d max(0,t+1-t0) or both, where the break
* point t0 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
*
*   LAGMETHOD=[INPUT]/AIC/BIC/TTEST
*     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 t-test, the number of lags for which the last included lag has a marginal
*     significance level less than the cutoff given by the SIGNIF option.
*
*   MAXLAGS=number of additional lags (LAGMETHOD=INPUT) or the maximum number of lags
*     to consider (other LAGMETHOD's) [number of observations**.25]
*
*   SIGNIF=cutoff significance level for LAGMETHOD=TTEST [.10]
*
*   SMPL=SMPL series with entries to skip
*
*   GRAPH/[NOGRAPH]
*     GRAPH requests a graph of the unit-root statistics with the different break points
*
*   PI=fraction of data range to skip at either end when examining possible break points [.15]
*
* Variables Defined:
*  %cdstat is the unit root test statistic
*  %%breakpoint is the entry at which %cdstat is achieved.
*
*
* Revision Schedule:
*  08/2003 Written by Tom Doan, Estima
*  11/2003 Correction to BREAK=BOTH which didn't include the trend break. Thanks to Kit Baum
*          for pointing this out.
*  07/2007 Suppress LINREG output on the lag selection regressions
*
procedure Zivot depvar start end
type series depvar
type integer start end
*
option choice  break     1 intercept trend both
option integer maxlags
option choice  lagmethod 1 input aic bic ttest
option real    signif    .10
option series  smpl
option switch  graph
option real    pi     .15
*
local integer maxlag bestlag startl endl lag znobs tstart tend time
local series dy trend breakint breaktrend breakstats
local real icmult ic icmin crit01 crit05
local vect[integer] others
*
* Figure out the range available for the series and one lag
*
inquire(reglist) startl>>start endl>>end
# depvar{0 1}
*
set dy    startl endl = depvar-depvar{1}
set trend startl endl = t
*
* Determine default value of maxlag
*
compute znobs=endl-startl+1
if .not.%defined(maxlags)
   compute maxlag=fix(znobs**.25)
else
   compute maxlag=maxlags

if lagmethod==1
   compute bestlag=maxlag
else {
   cmom(smpl=smpl) startl+maxlag endl
   # constant trend depvar{1} dy{0 to maxlag}
   compute icmult=%if(lagmethod==2,2.0/%nobs,log(%nobs)/%nobs)
   linreg(cmom,noprint) dy
   # constant trend depvar{1}
   compute ic=log(%rss/%nobs)+icmult*%nreg
   compute icmin=ic,bestlag=0
   do lag=maxlag,1,-1
      linreg(cmom,noprint) dy
      # constant trend depvar{1} dy{1 to lag}
      if lagmethod==4 {
         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
}
enter(varying) others
# constant trend
if bestlag>0 {
   enter(varying) others
   # others dy{1 to bestlag}
}
if break<>2 {
   enter(varying) others
   # others breakint
}
if break<>1 {
   enter(varying) others
   # others breaktrend
}
compute znobs=endl-startl+1-bestlag
compute tstart=startl+bestlag+fix(pi*znobs),tend=endl-fix(pi*znobs)
set breakstats tstart tend = 0.0
do time=tstart,tend
   set breakint   startl endl = t>=time
   set breaktrend startl endl = %max(t+1-time,0)
   linreg(noprint,smpl=smpl) dy startl+bestlag endl
   # depvar{1} others
   compute breakstats(time) = %tstats(1)
end do time
*
if graph {
   graph(header="Zivot-Andrews Unit Root Tests for "+%l(depvar))
   # breakstats tstart tend
}
ext(noprint) breakstats startl+bestlag endl
compute %cdstat=%minimum,%%breakpoint=%minent
disp
disp "Zivot-Andrews Unit Root Test"
disp "Dependent Variable" %l(depvar)
if break==1 {
   disp "Allowing for Break in Intercept Only"
   compute crit01=-5.34,crit05=-4.80
}
else
if break==2 {
   disp "Allowing for Break in Trend Only"
   compute crit01=-4.93,crit05=-4.42
}
else {
   disp "Allowing for Break in both Intercept and Trend"
   compute crit01=-5.57,crit05=-5.08
}
disp(hold) "Including" bestlag "Lags of Difference Selected by "
if lagmethod==1
   disp "User"
else
if lagmethod==2
   disp "AIC"
else
if lagmethod==3
   disp "BIC"
else
if lagmethod==4
   disp "Truncation at Longest Lag with P-value<" signif

disp "Minimum T-Statistic" %cdstat "at" %datelabel(%%breakpoint)
disp "Critical Values are 1%" *.## crit01 "and 5%" crit05
end

