*
* @GregoryHansen(options) start end
* # endogenous variables (dependent variable listed first)
*
* NOTE --- MAJOR CHANGE OF SYNTAX, July 2012 (version 8.2) to match
* other cointegration procedures.
*
* Implements the Gregory-Hansen cointegration test. The cointegrating
* regression is allowed to have a trend or not, and can have either a
* break in the intercept only or a break in all coefficients. The break
* point t0 is unknown and determined by finding the minimum value for
* the ADF statistic on the residuals from the broken cointegration
* regression. The number of lags of the change in the residual used in
* computing the ADF statistic can either be input, or can be selected
* automatically using AIC, BIC or general-to-specific pruning by
* t-tests.
*
* Options:
*
* There are two ways to indicate the type of model: the combination of the
* DET and BREAK options, or the MODEL option.
*
*   DET=[CONSTANT]/TREND
*   BREAK=[INTERCEPT]/ALL
*     DET chooses the deterministic components in the cointegrating
*     regression. BREAK chooses what components are allowed to a break
*     at the unknown change point. BREAK=ALL has the intercept and slope
*     coefficients and the trend if present.
*
*   MODEL=[CONSTANT]/TREND/NOTREND/REGIMETREND/FULLBREAK
*     CONSTANT includes an intercept and a level shift dummy.
*     TREND includes an intercept, trend and a level shift dummy.
*     NOTREND includes no trend, but dummies for the intercept and the slopes
*     REGIMETREND includes breaks in intercept, trend and slopes
*    (FULLBREAK is an older synonym for NOTREND). We recommend using the
*     DET and BREAK options as it is clearer what your model is.
*
*   LAGMETHOD --- note, renamed as METHOD for compatibility with similar
*     procedures
*
*   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 LAGS option is used. If AIC,
*     the AIC-minimizing value between 0 and LAGS 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.
*
*   LAGS=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 [.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]
*
*   [PRINT]/NOPRINT
*   TITLE=title for report ["Gregory-Hansen Cointegration Test"]
*
* Variables Defined:
*  %cdstat 			unit root test statistic
*  %%breakpoint 	entry at which %cdstat is achieved.
*  %%autop			number of augmenting lags used
*
* References:
*   Gregory, Allan W. & Hansen, Bruce E., 1996. "Residual-based tests for
*    cointegration in models with regime shifts," Journal of Econometrics,
*    vol. 70, no 1, 99-126.
*
*   Gregory, Allan W & Hansen, Bruce E, 1996. "Tests for Cointegration in
*    Models with Regime and Trend Shifts," Oxford Bulletin of Economics and
*    Statistics, vol. 58, no 3, 555-560.
*
* Revision Schedule:
*  09/2003 Written by Tom Doan, Estima
*  10/2004 Correction made to independent variables in regression
*  11/2010 Switch to use of REPORT
*  01/2011 Add options to allow for trend with break from the OBES article.
*          Switch to use of FIXED for critical values. Change to select ADF
*          lag length for each tested break point, rather than once based on
*          the base regression.
*  07/2012 Change LAGMETHOD to METHOD, change MAXLAGS to LAGS and change
*          syntax to list all endogenous variables on the supplementary card
*          to match other cointegration procedures.
*  12/2015 Add .10 critical values and provide warning about lack of
*          critical values if number of variables is too large. Add description
*          of null, and sample range.
*  12/2018 Check for constant as lead dependent variable (probable old syntax).
*
procedure GregoryHansen start end
type integer start end
*
option choice  det      1 constant trend
option choice  break    1 intercept all
option choice  model    0 constant trend notrend regimetrend fullbreak
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  graph
option switch  print			1
option string  title
*
local vect[integer]  xlist datavec breaklist
local integer        depvar
local vect[series]   xdummies
local integer        maxlag bestlag startl endl lag znobs tstart tend m i time
local series         trend breakint breaktrend breakstats e de
local real           icmult ic icmin crit01 crit05 crit10 besttstat
local report         greport
local string         ltitle lmethod ldescript
local integer        lmodel
*
* These are the .01, .05 and .10 critical values. The rows are the
* number of x variables, and the columns are the model.
*
fixed rect c01(4,4)
 -5.13 -5.45 -5.47 -6.02
 -5.44 -5.80 -5.97 -6.45
 -5.77 -6.05 -6.51 -6.89
 -6.05 -6.36 -6.92 -7.31
fixed rect c05(4,4)
 -4.61 -4.99 -4.95 -5.50
 -4.92,-5.29,-5.50 -5.96
 -5.28,-5.57,-6.00 -6.32
 -5.56,-5.83,-6.41 -6.84
fixed rect c10(4,4)
 -4.34 -4.72 -4.68 -5.24
 -4.69 -5.03 -5.23 -5.72
 -5.02 -5.33 -5.75 -6.16
 -5.31 -5.59 -6.17 -6.58
*
* Get the list of series. Separate them into the dependent variable and
* the RHS variables.
*
enter(varying) datavec
compute m=%rows(datavec)-1
dim xlist(m)
compute depvar=datavec(1)
ewise xlist(i)=datavec(i+1)
if depvar==0 {
   disp "Syntax: @GregoryHansen(options)  start end"
   disp "        # list of endogenous variables, dependent first"
   disp "Deterministics are handled with DET option"
   return
}
*
* Figure out the range available for the series
*
inquire(reglist) startl<<start endl<<end
# datavec
*
* Figure out which model is wanted
*
if model>0 {
   compute lmodel=model
   if lmodel==5
      compute lmodel=3
}
else
   compute lmodel=(break-1)*2+det
*
* Prepare for the slope dummies, if needed
*
if lmodel==3.or.lmodel==4
   dim xdummies(m)
*
set trend startl endl = t
*
* Determine default value of maxlag
*
compute znobs=endl-startl+1
if .not.%defined(maxlags).and..not.%defined(lags)
   compute maxlag=fix(znobs^.25)
else
if .not.%defined(lags)
   compute maxlag=maxlags
else
   compute maxlag=lags
*
* Construct the list of regressors for the broken regression. All
* include the constant, the broken constant and the original slope
* variables.
*
compute breaklist=||constant,breakint||
do i=1,m
   compute breaklist=%rladdone(breaklist,xlist(i))
end do i
*
* Models 2 and 4 need the trend, and 4 needs the broken trend as well.
*
if lmodel==2.or.lmodel==4
   compute breaklist=%rladdone(breaklist,trend)
if lmodel==4
   compute breaklist=%rladdone(breaklist,breaktrend)
*
* Models 3 and 4 need the broken x variables
*
if lmodel==3.or.lmodel==4 {
   do i=1,m
      compute breaklist=%rladdone(breaklist,xdummies(i))
   end do i
}
compute tstart=startl+fix(pi*znobs),tend=endl-fix(pi*znobs)
set breakstats tstart tend = 0.0
compute besttstat=%na
do time=tstart,tend
   *
   * Create the dummied out variables
   *
   set breakint   startl endl = t>time
   set breaktrend startl endl = %max(0,t-time)
   if lmodel==3.or.lmodel==4 {
      do i=1,m
         set xdummies(i) startl endl = xlist(i){0}*(t>time)
      end do i
   }
   *
   * Do the broken regression and save the residuals
   *
   linreg(noprint,smpl=smpl) depvar startl endl e
   # breaklist
   *
   * Do an ADF test on the residuals
   *
   set de = e-e{1}
   if method==1
      compute bestlag=maxlag
   else {
      cmom(smpl=smpl) startl+maxlag+1 endl
      # constant e{1} de{0 to maxlag}
      if method==2
         compute icmult=2.0/%nobs
      else
      if method==3
         compute icmult=log(%nobs)/%nobs
      else
         compute icmult=2.0*log(log(%nobs))/%nobs
      linreg(cmom,noprint) de
      # constant e{1}
      compute ic=log(%rss/%nobs)+icmult*%nreg
      compute icmin=ic,bestlag=0
      do lag=maxlag,1,-1
         linreg(cmom,noprint) de
         # constant e{1} de{1 to lag}
         if method==5.or.method==6 {
            if %ttest(%tstats(2+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
   }
   *
   * Re-run the regression over the range permitted by the choice
   * for <<bestlag>>.
   *
   if bestlag>0 {
      linreg(noprint,smpl=smpl) de startl+bestlag+1 endl
      # e{1} de{1 to bestlag}
   }
   else {
      linreg(noprint,smpl=smpl) de startl+1 endl
      # e{1}
   }
   *
   * Save the break statistic
   *
   compute breakstats(time) = %tstats(1)
   *
   * See if this is the smallest t-stat seen so far. If so, save it, the
   * time period and the number of lags used.
   *
   if .not.%valid(besttstat).or.besttstat>breakstats(time)
      compute besttstat=breakstats(time),%%breakpoint=time,%%autop=bestlag
end do time
*
compute %cdstat=besttstat
if graph {
   graph(header="Gregory-Hansen Cointegration Test",grid=(t==%%breakpoint))
   # breakstats tstart tend
}
*
* Handle description of lag choice
*
if method==2
   compute lmethod="AIC"
else
if method==3
   compute lmethod="BIC/SBC"
else
if method==4
   compute lmethod="Hannan-Quinn"
else
   compute lmethod="GTOS/t-tests("+%strval(signif,"#.###")+")"

if method==1
   compute lmethod="Using fixed lags "+%%autop
else
   compute lmethod="With "+%%autop+" lags chosen from "+maxlag+" by "+lmethod
*
* Handle description of model
*
if lmodel==1
   compute ldescript="Break in Intercept. No Trend."
else
if lmodel==2
   compute ldescript="Break in Intercept. Trend Included."
else
if lmodel==3
   compute ldescript="Full Structural Break. No Trend."
else
   compute ldescript="Full Structural Break. Trend Included."

if %defined(title)
   compute ltitle=title
else
   compute ltitle="Gregory-Hansen Cointegration Test"

report(use=greport,action=define,title=ltitle)
report(use=greport,atrow=1,atcol=1,span) ltitle
report(use=greport,atrow=2,atcol=1,span) "Null is no cointegration (residual has unit root)"
report(use=greport,atrow=3,atcol=1,span) "Regression Run From "+%datelabel(%regstart())+" to "+%datelabel(%regend())
report(use=greport,atrow=4,atcol=1) "Observations" %nobs
report(use=greport,atrow=5,atcol=1,span) ldescript
report(use=greport,atrow=6,atcol=1,span) lmethod

report(use=greport,atrow=8,atcol=1) "Variables"
report(use=greport,atrow=9,atcol=1) %l(depvar)
do i=1,m
   report(use=greport,atrow=9+i,atcol=1) %l(xlist(i))
end do i
*
compute crit01=c01(m,lmodel),crit05=c05(m,lmodel),crit10=c10(m,lmodel)
*
report(use=greport,atrow=11+m,atcol=1) "Minimum T-Statistic" %cdstat
report(use=greport,atrow=12+m,atcol=1) "Achieved At" %datelabel(%%breakpoint)
if m>4
   report(use=greport,atrow=13+m,atcol=1,span) "Critical values not available when m>4"
else {
   report(use=greport,atrow=13+m,atcol=1) "1% Critical Value" crit01
   report(use=greport,atrow=14+m,atcol=1) "5% Critical Value" crit05
   report(use=greport,atrow=15+m,atcol=1) "10% Critical Value" crit10
}
*
report(use=greport,action=format,atrow=11+m,atcol=2,tocol=2,picture="*.###")
if print
   report(use=greport,action=show)
end
