*
* @ERSTEST(OPTIONS) series start end
* performs Elliott, Rothenberg, Stock unit root tests
*
* Parameters:
*   series         Series to be tested
*
*   start, end     Range over which you want to perform the test.
*                  Defaults to defined range of the series.
*
* Options:
*   DET=NONE/[CONSTANT]/TREND, deterministic part (trend uses (1,t))
*   CBAR=(positive sign) value of c-bar to use [based upon choice for DET)
*   LAGS=number of augmenting lags [0]
*     LAGS indicates the degree of AR polynomial on the differences
*     to include in the unit root test on the residuals.  (Note - some
*     software will count the number of lags in the expanded autoregression,
*     which will be one more than the number of lags on the differences).
*   MAXLAGS=maximum number of additional lags to consider [number of observations**.25]
*   METHOD=[INPUT]/AIC/BIC/HQ/TTEST/GTOS/SBC/MAIC
*     Selects the method for deciding the number of additional lags. If
*     INPUT, the number of lags given by the MAXLAGS (or LAGS) option is used. If
*     AIC, the AIC-minimizing value between 0 and MAXLAGS is used; if BIC (or SBC),
*     it's the BIC-minimizing value, with MAIC the Ng-Perron Modified AIC minimizer,
*     and if TTEST (or GTOS), the number of lags for which the last included lag has
*     a marginal significance level less than the cutoff given by the SIGNIF option.
*   SIGNIF=cutoff significance level for METHOD=TTEST/GTOS[.10]
*
*   PRINT/[NOPRINT]
*   TITLE="title of report" ["DF-GLS Tests, Series xxxx"]
*
* Variables Defined:
*   %%pt        Value of PT test
*   %%dfgls     Value of DFGLS test
*   %%qt        Value of QT test
*   %%dfglsu    Value of DFGSLu test
*   %%autop     Lag length used in the DFGLS test
*
* A sub-procedure which can be used directly does the GLS detrending. The syntax
* for this is
*
* @ersdetrend series start end yd dy ut rhot
*
* Parameters:
*     series         Series to be tested
*
*     start, end     Range over which you want to perform the test.
*                    Defaults to defined range of the series.
*
*     yd             (Output) detrended series
*     dy             (Output) differenced series
*     ut             (Output) residuals from detrending regression
*     rhot           (Output) value of rho used in GLS
*
* Options:
*     DET=NONE/[INTERCEPT]/TREND, deterministic part (trend uses (1,t))
*     U/[NOU]  U indicates unconditional distribution is to be used
*     CBAR=value of c-bar to use [based upon choice for DET)
*
* The 5% critical values are supplied by the above papers except for
* the no detrending case with DFGLS, when standard DF critical values
* are used. The latter are taken from James Hamilton, _Time Series Analysis_,
* Princeton University Press, 1994, p. 763).
*
* References:
* Elliott, Rothenberg and Stock(1996), "Efficient Tests for an Autoregressive
* Unit Root", Econometrica, vol 64, no. 4, pp 813-836.
*
* Elliott(1999), "Efficient Tests for a Unit Root When the Initial Observation
* is Drawn from its Unconditional Distribution", International Economic Review,
* August 1999, vol 40, pp 767-783.
*
* Revision Schedule:
*  10/1998 RATS translation by Alan Taylor from Graham Elliott's GAUSS code
*          DEMO.GSS, ERS.PRC, POIU.PRC. I (Alan Taylor) thank Graham for his help.
*  07/2004 Syntax changed to be compatible with other UR tests. Updated to version
*          6 and generally cleaned up by Tom Doan, Estima.
*  03/2009 Correct typos in critical values. Switch to full use of REPORT.
*  09/2010 Add local report and title handling
*  06/2014 Add TITLE option
*  08/2019 Add options for lag length selection. Define %%AUTOP.
*
* Alan M. Taylor Department of Economics       tel:530-752-1572
* University of California                     fax:530-752-9382
* One Shields Ave
* Davis CA 95616-8578
* amtaylor@ucdavis.edu
*
*******************************************************************************
*
* @ERSTestCV( options ) critval
*
* Computes approximate critical values for ERS Unit Root test using response
* surface analysis.
*
* Parameters:
*   critval=[return] vector with 2 entries for the computed critical values at the
*     5% and 10% level (in this order). This doesn't have to be declared or
*     dimensioned.
*
* Options:
*   npts=number of observations [%NOBS]
*   lags=number of augmenting lags [0]
*   determ=noconst/[const]/trend
*
* DET=CONST is used for the DF-GLS mu test
* DET=TREND is used for the DF-GLS tau test
* DET=NOCONST is a simple DF test and returns the Mackinnov critical
*  values.
*
* Reference: Cheung and Lai(1995), "Lag Order and Critical Values of a
* Modified Dickey-Fuller test", Oxford Bulletin of Economics and
* Statistics, vol 53, no 3, pp 411-419.
*
procedure ERSTestCV critval
type vector *critval
*
option integer npts   %nobs
option integer lags   0
option choice  det    2      none constant trend
*
local integer i
local real    rn pn
*
fixed vect[rect] cl(3)(2,6)
 -1.6156   -0.1810    0.0000    0.0000    0.0000    0.0000
 -1.9393   -0.3980    0.0000    0.0000    0.0000    0.0000

 -1.9480  -17.8390  104.0860    0.8020    5.5580  -18.3320
 -1.6240  -19.8880  155.2310    0.7090    5.4800  -16.0550

 -2.8380  -20.3280  124.1910    1.2670   10.5300  -24.6000
 -2.5500  -20.1660  155.2150    1.1330    9.8080  -20.3130
*
dim critval(2)
*
compute rn=1.0/npts
compute pn=float(lags)/npts
*
do i=1,2
   compute critval(i)=cl(det)(i,1)+cl(det)(i,2)*rn+cl(det)(i,3)*rn^2+$
    cl(det)(i,4)*pn+cl(det)(i,5)*pn^2+cl(det)(i,6)*pn^3
end do i
end
*******************************************************************
*
* This does the lag selection if the user indicated a particular
* method. This is for a previously detrended series.
*
procedure ERSLagSelect s sdiff start end
type series  s sdiff
type integer start end
*
option integer maxlag
option choice  method      1 input aic bic hq ttest gtos sbc maic
option real    signif      .10
option integer *bestlag
*
local  real    icmult maicss ic icmin
local  integer lag

if method==1 {
   compute bestlag=maxlag
   return
}
cmom
# s{1} sdiff{0 to maxlag}

if method==2
   compute icmult=2.0/%nobs
else
if method==3.or.method==7
   compute icmult=log(%nobs)/%nobs
else
if method==4
   compute icmult=2.0*log(log(%nobs))/%nobs
else {
   *
   * This is needed for MAIC
   *
   compute maicss=%cmom(2,2)
}
linreg(cmom,noprint) sdiff
# s{1}
*
* Compute the IC's for lag=0
*
if method==8
   compute ic=log(%sigmasq)+2.0*(%beta(1)^2*maicss/%sigmasq)/%nobs
else
   compute ic=log(%sigmasq)+icmult*%nreg
compute icmin=ic,bestlag=0
*
* Walk backwards through the lags. The direction doesn't matter for
* the IC calculations, but does for GTOS.
*
do lag=maxlag,1,-1
   linreg(cmom,noprint) sdiff
   # s{1} sdiff{1 to lag}
   if method==5.or.method==6 {
      if %ttest(%tstats(%nreg),%ndf)<=signif {
         compute bestlag=lag
         break
      }
   }
   else {
      if method==8
         compute ic=log(%sigmasq)+2.0*(lag+%beta(1)^2*maicss/%sigmasq)/%nobs
      else
         compute ic=log(%sigmasq)+icmult*%nreg
      if ic<icmin
         compute icmin=ic,bestlag=lag
   }
end do lag
end ERSLagSelect
***********************************************************************
*
* This is a modified DF test with no deterministics. We also need an
* estimate of spectral density at frequency zero = what. It's assumed
* that hard values for start and end are included.
*
procedure ersdf s sdiff startl endl
type series  s sdiff
type integer startl endl
*
option integer lags  0
option real    *dftstat
option real    *what
*
local real    a1
*
linreg(noprint) sdiff startl+lags+1 endl
# s{1} sdiff{1 to lags}
if lags>0 {
   summarize(noprint)
   # sdiff{1 to lags}
   compute a1 = %sumlc ; * sum of lag diff coefficients
}
else
   compute a1 = 0.0
*
compute dftstat = %tstats(1)
compute what    = sqrt(%seesq/(1-a1)^2) ; * spectral density

end ersdf

*******************************************************************************
procedure ersdetrend series start end yd dy ut
type series    series *yd *dy *ut
type integer   start end
*
* Local to unity detrending routine
*
option choice  det      2  none constant trend
option switch  u        0
option real    cbar
option real    *rhot
*
local integer  startl endl nobs
local series   y zt zt2 qy
*
* Check parameters
*

if .not.%defined(series).or..not.%defined(yd).or..not.%defined(dy).or..not.%defined(ut) {
   display "##@ERSDETREND Syntax: series  start  end  yd  dy  ut"
   return
}

*
* Adjust the range and determine the number of observations
*

inquire(series=series) startl<<start endl<<end
stats(noprint) series startl endl ; compute nobs=%nobs
set y startl endl = series

*
* do GLS detrending
* Set the local to unit autoregressive parameter for GLS
*

if %defined(cbar)
   compute rhot = 1 - cbar/nobs

*
* deterministic=()
*

if det==1
{
   if u {
      disp "##@ERSDETREND: Can't Use U option with DET=NONE"
   }
   if .not.%defined(cbar) ; compute rhot=1-10.0/nobs
   set  yd  startl   endl   = y
   set  qy  startl   endl   = %if(t==startl,y,y-rhot*y{1})
   set  dy  startl   endl   = %if(t==startl,y,y-y{1}) ; * OLS detrending under null, keep y(1)
   set  ut  startl   endl   = qy
}

*
* deterministic = (1)
*

if det==2
{
   if .not.%defined(cbar) ; compute rhot=%if(u,1-10.0/nobs,1-7.0/nobs)
   set  qy  startl endl = y-rhot*y{1}
   set  zt  startl endl = 1-rhot
   if u {
      set qy startl startl = y*sqrt(1-rhot^2)
      set zt startl startl = sqrt(1-rhot^2)
   }
   else {
      set qy startl startl = y
      set zt startl startl = 1
   }
   linreg(noprint) qy startl endl ut; # zt
   set yd startl endl = y - %beta(1)
   set dy startl endl = %if(t==startl,0,y-y{1}); * OLS detrending under null
}

*
* deterministic = (1,t)
*

if det==3
{
   if .not.%defined(cbar) ; compute rhot=%if(u,1-10.0/nobs,1-13.5/nobs)
   set  qy  startl   endl   = y-rhot*y{1}
   set  zt  startl   endl   = 1-rhot
   set  zt2 startl   endl   = (1-rhot)*(t-startl+1)+rhot
   if u {
      set qy  startl startl = y*sqrt(1-rhot^2)
      set zt  startl startl =   sqrt(1-rhot^2)
      set zt2 startl startl =   sqrt(1-rhot^2)
   }
   else {
      set qy  startl startl = y
      set zt  startl startl = 1
      set zt2 startl startl = 1
   }
   linreg(noprint) qy startl endl ut; # zt zt2
   set yd startl endl = y - %beta(1) - %beta(2)*(t-startl+1)
   set dy startl endl = %if(t==startl,0,y-y{1})
   diff(center) dy startl endl dy               ; * OLS detrending under null
}

end ersdetrend
*******************************************************************************
*
* ersdfproc should be called only after doing detrending. It uses the
* output series from ersdetrend.
*
procedure ersdfproc yd dy ut start end
type series    yd dy ut
type integer   start end
*
option integer lags    0
option real    rhot
option real    *dfgls
option real    *pt
*
local real     what sdy sqy
local integer  nobs
*
@ersdf(lags=lags,dftstat=dfgls,what=what) yd dy start end
sstats start end dy^2>>sdy ut^2>>sqy
compute nobs=%nobs,sdy=sdy/nobs,sqy=sqy/nobs
compute pt = nobs*(sqy-rhot*sdy)/(what^2)
end
*******************************************************************************
*
* Main entry point for erstest
*
procedure erstest series start end
type series    series
type integer   start end
*
option choice  det    2   none constant trend

option integer lags
option integer maxlags
option choice  method      1 input aic bic hq ttest gtos sbc maic
option real    signif      .10

option choice  range       1 data estimation

option switch  print  1
option string  title
option real    cbar
*
local integer  startl endl maxlag nobs rstart bestlag bestlagu znobs tag
local series   s dy yd ut
local real     rhot
local string   ltitle lmethod ltrend
local report   ereport
*
* Global Variables Defined by Procedure
*
declare real   %%pt %%dfgls %%qt %%dfglsu
*
* set up arrays with 1%,2.5%,5%,10% asymptotic cvs for each test
* cvs from GE 10/98 and 6/99
* store the cvs -- one row per test, detrending method in i, sig lev in j (1%,2.5%,5%,10%)
*

fixed rect cvpt(3,4)
  .    .    .    .
 1.99 2.55 3.26 4.48
 3.96 4.78 5.62 6.89
fixed rect cvdfgls(3,4)
 -2.58 -2.23 -1.95 -1.62
 -2.58 -2.23 -1.95 -1.62
 -3.48 -3.15 -2.89 -2.57
fixed rect cvqt(3,4)
   .     .     .     .
  3.06  3.80  4.65  5.94
  2.05  2.44  3.15  3.44
fixed rect cvdfglsu(3,4)
   .     .     .     .
 -3.28 -2.98 -2.73 -2.46
 -3.71 -3.41 -3.17 -2.91

if .not.%defined(series)
{
   display "##Syntax: @ERSTEST( options ) <<series>>  start  end"
   return
}

inquire(series=series) startl<<start endl<<end
compute znobs=endl-startl+1
*
if %defined(lags)
   compute maxlag=lags
else
if %defined(maxlags)
   compute maxlag=maxlags
else
if method==1
   compute maxlag=0
else
   compute maxlag=fix(znobs^.25)
*
* Create the string describing the method of lag selection
*
if method==1
   compute lmethod="User Input"
else
if method==2
   compute lmethod="AIC"
else
if method==3.or.method==7
   compute lmethod="BIC/SBC"
else
if method==4
   compute lmethod="Hannan-Quinn"
else
if method==8
   compute lmethod="Ng-Perron MAIC"
else
   compute lmethod="GTOS/t-tests("+%strval(signif,"#.###")+")"

*
* Adjust the range
*

if range==2.and.%defined(start) {
   compute rstart=start
   compute startl=startl-(maxlag+1)
}
else
   compute rstart=-1
*
* Make a copy of the input series in case it gets overwritten (if, for
* instance, series=%resids).
*
clear s
set s startl endl = series
statistics(noprint) s startl endl
compute nobs = %nobs
if %nmiss>0 {
   display "##@ERSTEST on "+%l(series)+" over "+%datelabel(startl)+" to "+%datelabel(endl)
   display "##has missing observations"
   return
}
*
* DFGLS and PT test
*
@ersdetrend(det=det,cbar=cbar,rhot=rhot) s startl endl yd dy ut
@erslagselect(maxlags=maxlag,method=method,signif=signif,bestlag=bestlag) yd dy startl endl
@ersdfproc(lags=bestlag,rhot=rhot,dfgls=%%dfgls,pt=%%pt) yd dy ut startl endl
*
* DFGLSu and QT test
*
@ersdetrend(det=det,cbar=cbar,u,rhot=rhot) s startl endl yd dy ut
@erslagselect(maxlags=maxlag,method=method,signif=signif,bestlagu=bestlagu) yd dy startl endl
@ersdfproc(lags=bestlagu,rhot=rhot,dfgls=%%dfglsu,pt=%%qt) yd dy ut startl endl
*
compute %%autop=bestlag
*
* display results
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="DF-GLS Tests, Dependent Variable "+%l(series)

if det==1
   compute ltrend="Detrend = none, z(t)=0"
else
if det==2
   compute ltrend="Detrend = constant, z(t)=(1)"
else
   compute ltrend="Detrend = constant and linear time trend, z(t)=(1,t)"

report(use=ereport,action=define,title=ltitle)
report(use=ereport,atrow=1,atcol=1,span) ltitle
report(use=ereport,atrow=2,atcol=1,span) "From "+%datelabel(startl)+" to "+%datelabel(endl)

report(use=ereport,atrow=3,atcol=1,span) "Lag Length Chosen from "+maxlag+" by "+lmethod
report(use=ereport,atrow=4,atcol=1,span) ltrend
report(use=ereport,atrow=5,atcol=1,span) "Tests for a unit root null. All tests reject null in lower tail"
report(use=ereport,atrow=6,atcol=4,tocol=7,span) "Critical values (asymptotic)"

report(use=ereport,atrow=7,atcol=1,span) "Elliott et al (1996 Econometrica)"
report(use=ereport,atrow=8,atcol=2,align=center) "Stat" "Lags" "1%(**)" "2.5%" "5%(*)" "10%"
compute tag=fix(1+(%%pt<cvpt(det,1))+(%%pt<cvpt(det,3)))
report(use=ereport,atrow=9,atcol=1) "PT" "" "" cvpt(det,1) cvpt(det,2) cvpt(det,3) cvpt(det,4)
report(use=ereport,atrow=9,atcol=2,special=tag) %%pt
compute tag=fix(1+(%%dfgls<cvdfgls(det,1))+(%%dfgls<cvdfgls(det,3)))
report(use=ereport,atrow=10,atcol=1) "DFGLS" "" "" cvdfgls(det,1) cvdfgls(det,2) cvdfgls(det,3) cvdfgls(det,4)
report(use=ereport,atrow=10,atcol=2,special=tag) %%dfgls
report(use=ereport,atrow=10,atcol=3) bestlag

report(use=ereport,atrow=11,atcol=1,span) "Elliott (IER 1999)"
if det==1 {
   report(use=ereport,atrow=12,atcol=1) "QT" %na
   report(use=ereport,atrow=13,atcol=1) "DFGLSu" %na
}
else {
   compute tag=fix(1+(%%qt<cvqt(det,1))+(%%qt<cvqt(det,3)))
   report(use=ereport,atrow=12,atcol=1) "QT" "" "" cvqt(det,1) cvqt(det,2) cvqt(det,3) cvqt(det,4)
   report(use=ereport,atrow=12,atcol=2,special=tag) %%qt
   compute tag=fix(1+(%%dfglsu<cvdfglsu(det,1))+(%%dfglsu<cvdfglsu(det,3)))
   report(use=ereport,atrow=13,atcol=1) "DFGLSu" "" ""  cvdfglsu(det,1) cvdfglsu(det,2) cvdfglsu(det,3) cvdfglsu(det,4)
   report(use=ereport,atrow=13,atcol=2,special=tag) %%dfglsu
   report(use=ereport,atrow=13,atcol=3) bestlagu
}
report(use=ereport,action=format,atrow=9,atcol=4,tocol=7,picture="*.##",width=8,align=decimal)
report(use=ereport,action=format,atrow=9,atcol=2,tocol=2,picture="*.###",width=10,align=decimal)

if print==0
   return

report(use=ereport,action=show)

end
