*
* @Zivot(options) series start end
*
* Implements the Zivot-Andrews unit root test, 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
*
*   CRIT=[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.
*
*   PI=fraction of data range to skip at either end when examining possible
*     break points [.15]
*
*   LAGS=number of additional lags (CRIT=INPUT)
*
*   MAXLAGS=number of additional lags (CRIT=INPUT) or the maximum number of
*     lags to consider (other crit's) [number of observations**.25]
*
*   SIGNIF=cutoff significance level for CRIT=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.
*
*   [PRINT]/NOPRINT
*
* Variables Defined:
*  %cdstat is the unit root test statistic
*  %%breakpoint is the entry at which %cdstat is achieved.
*
* Reference:
*   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/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
*   03/2009 Rename LAGMETHOD (to CRIT) and introduce LAGS option to be consistent
*           with similar instructions. Switch to use of REPORT.
*   10/2009 Fix * tags on test statistic
*   12/2009 Correct range on EXTREMUM.
*
procedure Zivot series start end
type series    series
type integer   start end
*
option choice  break     1 intercept trend both
option integer lags
option integer maxlags
option choice  crit      1 input aic bic ttest
option real    signif    .10
option series  smpl
option switch  graph     0
option switch  print     1
option real    pi     	 .15
*
local integer  maxlag bestlag startl endl lag znobs tstart tend time tag
local series   dy trend breakint breaktrend breakstats
local real     icmult ic icmin crit01 crit05
local string   lmethod
local vect[integer] others
*
* Check for invalid syntax
*
if .not.%defined(series) {
   disp "Syntax: @zivot(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
*
* 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)

if crit==1
   compute bestlag=maxlag
else {
   cmom(smpl=smpl) startl+maxlag endl
   # constant trend series{1} dy{0 to maxlag}
   compute icmult=%if(crit==2,2.0/%nobs,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 crit==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
   # series{1} others
   compute breakstats(time) = %tstats(1)
end do time
*
if graph {
   graph(header="Zivot-Andrews Unit Root Tests for "+%l(series))
   # breakstats tstart tend
}
ext(noprint) breakstats tstart tend
compute %cdstat=%minimum,%%breakpoint=%minent

report(action=define)
report(atrow=1,atcol=1,span) "Zivot-Andrews Unit Root Test, Series "+%l(series)
if break==1 {
   report(atrow=2,atcol=1,span) "Allowing for Break in Intercept Only"
   compute crit01=-5.34,crit05=-4.80
}
else
if break==2 {
   report(atrow=2,atcol=1,span) "Allowing for Break in Trend Only"
   compute crit01=-4.93,crit05=-4.42
}
else {
   report(atrow=2,atcol=1,span) "Allowing for Break in both Intercept and Trend"
   compute crit01=-5.57,crit05=-5.08
}
report(atrow=3,atcol=1,span) "Breaks Tested for "+%datelabel(tstart)+" to "+%datelabel(tend)

if crit==1
   compute lmethod="User"
else
if crit==2
   compute lmethod="AIC"
else
if crit==3
   compute lmethod="BIC"
else
if crit==4
   compute lmethod="Longest Lag with P-value<"+%strval(signif,".##")
report(atrow=4,atcol=1,span) "Including "+bestlag+" Lags of Difference"
report(atrow=5,atcol=1,span) "Selected by "+lmethod
report(atrow=7,atcol=1) "Sig Level" "Crit Value"
report(atrow=8,atcol=1) "1%(**)"    crit01
report(atrow=9,atcol=1) "5%(*)"     crit05
report(atrow=11,atcol=1) "Breakpoint" "TestStat"
report(atrow=12,atcol=1) %datelabel(%%breakpoint)
compute tag=fix(1+(%cdstat<crit01)+(%cdstat<crit05))
report(atrow=12,atcol=2,special=tag) %cdstat
report(action=format,atrow=7,atcol=2,tocol=2,width=8,align=decimal)
if print
   report(action=show)
end

