*
* @EndersGranger( options)  y  start  end
* does various types of the unit root regressions with threshold
* breaks
*
* The regression run is
*
*  dy = (rho1 D1 x (y{1} - attractor) + rho2 D2 x (y{1} - attractor)) + lags of dy
*
* where D1 is the dummy for threshold series < tau and D2 is
* 1-D1.
*
* It produces three test statistics:
*
*   t-max is the larger of the two t-statistics on rho1=0 and rho2=0.
*   phi is the joint F-test for rho1=rho2=0
*   equal is the test for rho1=rho2.
*
* Parameters:
*   y   = series to be tested.
*   start end  = range to use [defined range of y]
*
* You can either input the threshold value or request a search.
*
* Options:
*   LAGS=# of lags on the differences [1]
*   ATTRACTOR=ZERO/[CONSTANT]/TREND
*
*   MODEL=[TAR]/MTAR/CMTAR
*   With MODEL=TAR, the threshold series is the lag of y. With MODEL=MTAR,
*   it's the lagged difference of y.
*
*   TAU=threshold value [not used]
*   PI=fraction of high and low empirical values for the threshold series
*    which are omitted in the search.
*
*   [PRINT]/NOPRINT
*   TITLE=title for test report ["Enders-Granger Test"]
*
* If you provide a value of TAU, the regressions are done with that
* fixed value for tau. If you don't, the values of the threshold series
* (excluding the PI fraction at each end) are searched for the one that
* minimizes the sum of squared residuals.
*
* Variables Defined:
*  %%BREAKVALUE    Threshold value chosen or input
*  %%EQUALSTAT     F-test for rho1=rho2
*  %%PHISTAT       joint F-test for rho1=0 and rho2=0
*  %%TMAX          maximum t-statistic between rho1=0 and rho2=0
*
* Reference:
* Enders and Granger(1998), "Unit-Root Tests and Asymmetric Adjustment
* with an Example Using the Term Structure of Interest Rates", JBES, vol
* 16, pp 304-11.
*
* Revision Schedule:
*  10/2011 Written by Tom Doan, Estima
*
procedure EndersGranger yp start end
type series  yp
type integer start end
*
option integer lags 1
option choice  attractor 2 zero constant trend
option choice  model 1 tar mtar cmtar
option real    tau
option real    pi   	.15
option switch  print  1
option string  title
*
local series 	y dy trend a0test
local vec[series] xsplit
local integer  i ttau
local vector   testvalues
local integer  startl endl pistart piend
local real     rssbest ltau taubest a1best
local report   esreport
local string   ltitle
local real     a1
local vect     tgrid
*
* Make a copy of the input series in case it's %resids and would get
* overwritten.
*
set y  = yp
set dy = y-y{1}
*
inquire(reglist) startl<<start endl<<end
# dy{0 to lags}
*
dim xsplit(2)
*
if attractor==1 {
   compute ltau=0.0
   if model==1 {
      set xsplit(1) = %max(y{1},0.0)
      set xsplit(2) = %min(y{1},0.0)
   }
   else {
      set xsplit(1) = y{1}*(dy{1}>=0.0)
      set xsplit(2) = y{1}*(dy{1}<0.0)
   }
}
else
if attractor==2 {
   if %defined(tau)
      compute ltau=tau
   else {
      @UniqueValues(values=testvalues) y startl endl
      compute pistart=fix(pi*%rows(testvalues)),piend=%rows(testvalues)-pistart+1
      do ttau=pistart,piend
         compute ltau=testvalues(ttau)
         if model==1 {
            set xsplit(1) = %max(y{1}-ltau,0.0)
            set xsplit(2) = %min(y{1}-ltau,0.0)
         }
         else {
            set xsplit(1) = (y{1}-ltau)*(dy{1}>=0.0)
            set xsplit(2) = (y{1}-ltau)*(dy{1}<0.0)
         }
         *
         linreg(noprint) dy startl endl
         # xsplit dy{1 to lags}
         if ttau==pistart.or.%rss<rssbest
            compute taubest=ltau,rssbest=%rss
      end do tau
      compute ltau=taubest
   }
   if model==1 {
      set xsplit(1) = %max(y{1}-ltau,0.0)
      set xsplit(2) = %min(y{1}-ltau,0.0)
   }
   else {
      set xsplit(1) = (y{1}-ltau)*(dy{1}>=0.0)
      set xsplit(2) = (y{1}-ltau)*(dy{1}<0.0)
   }
}
else
if attractor==3 {
   set trend = t
   linreg(noprint) y startl endl
   # constant trend
   compute tgrid=%seqa(%beta(2)-4.0*%stderrs(2),8.0*%stderrs(2)/50.0,50)
   compute rssbest=%na
   dofor a1 = tgrid
      set a0test startl endl = y-a1*t
      @UniqueValues(values=testvalues) a0test startl endl
      do ttau=1,%rows(testvalues)
         compute ltau=testvalues(ttau)
         if model==1 {
            set xsplit(1) = %max(y{1}-ltau-a1*(t-1),0.0)
            set xsplit(2) = %min(y{1}-ltau-a1*(t-1),0.0)
         }
         else {
            set xsplit(1) = (y{1}-ltau-a1*(t-1))*(dy{1}>=0.0)
            set xsplit(2) = (y{1}-ltau-a1*(t-1))*(dy{1}<0.0)
         }
         linreg(noprint) dy startl endl
         # xsplit dy{1 to lags}
         if .not.%valid(rssbest).or.%rss<rssbest
            compute taubest=ltau,a1best=a1,rssbest=%rss
      end do ttau
   end dofor a1
   compute ltau=taubest
   if model==1 {
      set xsplit(1) = %max(y{1}-ltau-a1best*(t-1),0.0)
      set xsplit(2) = %min(y{1}-ltau-a1best*(t-1),0.0)
   }
   else {
      set xsplit(1) = (y{1}-ltau-a1best*(t-1))*(dy{1}>=0.0)
      set xsplit(2) = (y{1}-ltau-a1best*(t-1))*(dy{1}<0.0)
   }
}
*
linreg(noprint) dy startl endl
# xsplit dy{1 to lags}
*
compute %%breakvalue=ltau
compute %%tmax=%max(%tstats(1),%tstats(2))
test(zeros,noprint)
# 1 2
compute %%phistat=%cdstat
restrict(noprint) 1
# 1 2
# 1.0 -1.0 0.0
compute %%equalstat=%cdstat
*
if .not.print
   return

if %defined(title)
   compute ltitle=title
else
   compute ltitle="Enders-Granger Test"

report(use=esreport,action=define,title=ltitle)
report(use=esreport,atrow=1,atcol=1,span) ltitle
if model==1
   report(use=esreport,row=new,atcol=1) "TAR Model"
else
   report(use=esreport,row=new,atcol=1) "M-TAR Model"
report(use=esreport,row=new,atcol=1) "Lags" lags
report(use=esreport,row=new,atcol=1) "Attractor" ltau
if attractor==3
   report(use=esreport,row=current,atcol=3) "+" a1best "t"
report(use=esreport,row=new,atcol=1) "T-Max" %%tmax
report(use=esreport,row=new,atcol=1) "Phi"   %%phistat
report(use=esreport,row=new,atcol=1) "Equality" %%equalstat
report(use=esreport,row=new,atcol=2,align=center) "Coeff" "Std Error"
report(use=esreport,row=new,atcol=1) "Above" %beta(1) %stderrs(1)
report(use=esreport,row=new,atcol=1) "Below" %beta(2) %stderrs(2)
do i=1,lags
   report(use=esreport,row=new,atcol=1) "DY{"+(i)+"}" %beta(2+i) %stderrs(2+i)
end do i
report(use=esreport,action=format,picture="*.####")
report(use=esreport,action=show)
end
