*
* @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-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
*   TITLE=title for report ["Zivot-Andrews Unit Root Test, Series ..."]
*
* Variables Defined:
*  %cdstat is the unit root test statistic
*  %%breakpoint is the entry at which %cdstat is achieved.
*  %%autop is the chosen number of lags at the breakpoint
*
* 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.
*   04/2010 Add display of selected regression
*   06/2011 Add TITLE option
*   05/2014 Add %%AUTOP variable
*   11/2014 Add (more standard) METHOD option as alternative to older CRIT option.
*   05/2017 Fix reported break data (was off by one)
*   05/2017 Add more standard GTOS as synonym for TTEST for METHOD
*   05/2017 Improve description of chosen lag length
*
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      0 input aic bic gtos ttest
option choice  method    1 input aic bic gtos ttest
option real    signif    .10
option series  smpl
option switch  graph     0
option switch  print     1
option real    pi     	 .15
option string  title
option switch  trace     0
*
declare integer %%breakpoint %%autop
*
local integer  lmethod
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   llagmethod llagreport
local string   ltitle
local report   zreport
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)
*
* Figure out whether to obey the CRIT or METHOD option
*
if crit>0
   compute lmethod=crit
else
   compute lmethod=method

*
* Remap GTOS
*

if lmethod==5
   compute lmethod=4

if lmethod==1
   compute bestlag=maxlag
else {
   cmom(smpl=smpl) startl+maxlag endl
   # constant trend series{1} dy{0 to maxlag}
   compute icmult=%if(lmethod==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
   if trace.and.lmethod<>4
      disp "Lag" 0 "IC" ic
   do lag=maxlag,1,-1
      linreg(cmom,noprint) dy
      # constant trend series{1} dy{1 to lag}
      if lmethod==4 {
         if trace
            disp "Lag" lag "TStat" %tstats(3+lag)
         if %ttest(%tstats(3+lag),%ndf)<=signif {
            compute bestlag=lag
            break
         }
      }
      else {
         compute ic=log(%rss/%nobs)+icmult*%nreg
         if trace
            disp "Lag" lag "IC" ic
         if ic<icmin
            compute icmin=ic,bestlag=lag
      }
   end do lag
}
enter(varying) others
# constant trend
*
compute %%autop=bestlag
*
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-time,0)
   linreg(noprint,smpl=smpl) dy startl+bestlag endl
   # series{1} others
   compute breakstats(time) = %tstats(1)
end do time
*
if trace
   print / breakstats
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

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

report(use=zreport,action=define,title=ltitle)
report(use=zreport,atrow=1,atcol=1,span) ltitle
if break==1 {
   report(use=zreport,atrow=2,atcol=1,span) "Allowing for Break in Intercept Only"
   compute crit01=-5.34,crit05=-4.80
}
else
if break==2 {
   report(use=zreport,atrow=2,atcol=1,span) "Allowing for Break in Trend Only"
   compute crit01=-4.93,crit05=-4.42
}
else {
   report(use=zreport,atrow=2,atcol=1,span) "Allowing for Break in both Intercept and Trend"
   compute crit01=-5.57,crit05=-5.08
}
report(use=zreport,atrow=3,atcol=1,span) "Breaks Tested for "+%datelabel(tstart)+" to "+%datelabel(tend)

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="GTOS/t-tests("+%strval(signif,"#.##")+")"

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

report(use=zreport,atrow=4,atcol=1,span) llagreport
report(use=zreport,atrow=5,atcol=1,span) "Selected by "+llagmethod
report(use=zreport,atrow=7,atcol=1) "Sig Level" "Crit Value"
report(use=zreport,atrow=8,atcol=1) "1%(**)"    crit01
report(use=zreport,atrow=9,atcol=1) "5%(*)"     crit05
report(use=zreport,atrow=11,atcol=1) "Breakpoint" "TestStat"
report(use=zreport,atrow=12,atcol=1) %datelabel(%%breakpoint)
compute tag=fix(1+(%cdstat<crit01)+(%cdstat<crit05))
report(use=zreport,atrow=12,atcol=2,special=tag) %cdstat
report(use=zreport,action=format,atrow=7,atcol=2,tocol=2,width=8,align=decimal)
if print {
   report(use=zreport,action=show)
   set breakint   startl endl = t>%%breakpoint
   set breaktrend startl endl = %max(t-%%breakpoint,0)
   labels breakint breaktrend
   # "D("+%datelabel(%%breakpoint)+")" "DT("+%datelabel(%%breakpoint)+")"
   linreg(print,smpl=smpl,title="Zivot-Andrews---Selected Regression") dy startl+bestlag endl
   # series{1} others
}
end
