*
*  @SPUNIT(options) series start end
*
*  Parameters:
*   Series      The series to be analyzed.
*   Start end   The range of the series to use (defaults to the entire range)
*
*  Options:
*   P=maximum polynomial trend order considered (1 to 4) [1].
*   LAGS=(maximum) # of lags in the lag windows for the non-parametric
*     estimates of the long-run variance [Schwert's fix(12*(T/100)^(1/4))]
*
*   [PRINT]/NOPRINT
*   TITLE="title for report" ["Schmidt-Phillips Test for a Unit Root in xxxx"]
*
*   SIGNIF=(no longer used - table of values is provided)
*
* @SPUNIT computes the "Schmidt-Phillips" test (TAU) for a unit root
* which is constructed considering the following parameterization -DGP-:
*
*       yt = psi + epsilon * trend + X(t), X(t) = beta * X(t-1) + e(t)
*
* Under the assumption that 'e(t)' are iid N(0,sigma^2e) and the initial
*  X(0) is taken as fixed.
*
* The MLE's of epsilon and psi(x) = psi + X(0) are:
*
*       eps(estimated) = mean delta(y) = (y(T)-y(1))/(T-1)
*       psi(x)(estimated) = y(1) - epsilon(estimated)*
*
* The 'residuals' are defined as follows:
*
*       Sbar(t) = y(t) - psi(x)(estimated) - epsilon(estimated), t=1,2,...,T.
*
* The TAU statistic corresponds to the t-statistic for phi=0 in the next
* regression:
*
*       delta(y(t)) = intercept + phi * Sbar(t-1) + error(t=2,...,T).
*
* This procedure also computes the semiparametric corrections to the TAU
* statistic Z(TAU) and a consistent estimate of the variance proposed by
* Schmidt and Phillips (1992) following the lines of Phillips (1987) and
*  Phillips and Perron(1988) that are asymptotically robust to error
* autocorrelation and heteroskedasticity so:
*
*       Z(TAU) = Tau/omega,
*       omega = s^2/s^2(l),
*
* Expressions for S^2 and s^2(l) are in Schmidt and Phillips (1992) pp
* 266. Z(TAU) statistics are calculated until (l), the truncation
* parameter of the Bartlett's Window, equals to the value of
* [fix(12*(T/100)^(1/4))] proposed by Schwert(1989).
*
* Furthermore the procedure computes (based upon the choice of P) the
* corresponding statistics TAU and Z(TAU) in the case of extended
* parametrizations to higher order polynomial trends considered by
* Schmidt and Phillips (1992) like:
*
*   delta(y(t)) = intercept + c(1)*trend + Phi*Sbar(t-1) + error
*   delta(y(t)) = intercept + c(1)*trend + c(2)* trend^(2) + Phi*Sbar(t-1) + error
*   delta(y(t)) = intercept + c(1)*trend + c(2)* trend^(2) + c(3)* trend^(3)
*                                                          + Phi*Sbar(t-1) + error
*  The maximum polynomial trend order considered is p = 4.
*
* The procedure provides the critical values tabulated by Schmidt and
* Phillips for TAU, at 0.01, 0.025, 0.05, and 0.1  significance levels
* for different sample sizes and for p=1,2,3 and 4. In addition the
* Ljung-Box Q statistic with significance level for residual serial
* correlation at the last lag calculated form (T/4) is displayed.
*
* This test allows for trend under both the null and the alternative,
* without introducing any parameters that are irrelevant under either,
* as happens in the Dickey-Fuller test.
*
* References:
*
* Schmidt, Peter. and Phillips P.C. "LM Test for a Unit Root in the
* Presence of Deterministic Trends". Oxford Bulletin of Economics and
* Statistics, 54,3.(1992) pp 257-287.
*
* Schwert, William G. "Test for Unit Roots: a Montecarlo Investigation".
* Journal of Business and Economic Statistics, April 1989, vol 7, No.2.
*
* Revision Schedule:
*   01/1998 Written by Diego Mauricio Vasquez, Santafe de Bogota, Colombia.
*   10/2008 Updated to use newer instructions. Correction of error in serial
*           correlation calculations.
*   03/2009 Updated for version 7.20 wizard. LAG option added.
*   08/2019 Add information to output
**********************************************************************************************
Procedure SPUNIT series start end
type series    series
type integer   start end
*
option real    signif 0.05
option integer p      1
option switch  print  1
option integer lags
option string  title
*
local series  Sbar
local series  dseries resids
local real    sum2 omega2 tau ztau
local integer n nrange startl endl l4 l12 tcount i row lagsToUse
local integer tag
local vect[series] trends
local report  sreport
local string  ltitle

fixed vect[int] nlimits(6)
 28 53 150 250 550 1500
*
* Blocked by p
*
fixed vect[rect] cv(4)(7,4)

-3.90 -3.50 -3.18 -2.85
-3.73 -3.39 -3.11 -2.80
-3.63 -3.32 -3.06 -2.77
-3.61 -3.30 -3.04 -2.76
-3.59 -3.29 -3.04 -2.76
-3.58 -3.28 -3.02 -2.75
-3.56 -3.27 -3.02 -2.75

-4.52 -4.09 -3.77 -3.41
-4.28 -3.93 -3.65 -3.34
-4.16 -3.84 -3.60 -3.31
-4.12 -3.81 -3.55 -3.28
-4.08 -3.80 -3.55 -3.26
-4.06 -3.78 -3.53 -3.26
-4.06 -3.78 -3.52 -3.26

-5.07 -4.62 -4.26 -3.89
-4.73 -4.37 -4.08 -3.77
-4.59 -4.29 -4.03 -3.72
-4.53 -4.24 -3.99 -3.69
-4.50 -4.20 -3.96 -3.68
-4.49 -4.19 -3.95 -3.68
-4.44 -4.16 -3.93 -3.67

-5.57 -5.09 -4.70 -4.30
-5.13 -4.78 -4.47 -4.15
-4.99 -4.66 -4.39 -4.10
-4.90 -4.60 -4.33 -4.06
-4.85 -4.56 -4.31 -4.03
-4.83 -4.54 -4.31 -4.03
-4.81 -4.52 -4.29 -4.01

if .not.%defined(series) {
   disp "##Syntax: @spunit(options)  SERIES start end"
   return
}

if p<1.or.p>4 {
   disp "##@SPUNIT. P option must be >=1 and <= 4"
   return
}

inquire(series=series) startl<<start endl<<end
compute n = endl-startl+1

compute l4  = fix(4*(N/100.0)^0.25)
compute l12 = fix(12*(N/100.0)^0.25)
if %defined(lags)
   compute lagsToUse = lags
else
   compute lagsToUse = l12

set dseries startl endl = series-series{1}

dim trends(p-1)
do i=1,p-1
   set trends(i) = t^i
end do i
*
* Compute residuals under null of unit root
*
linreg(noprint) dseries startl+1 endl resids
# constant trends
*
* Build the sbar variable, and redo the regression including its lag
*
set(first=0.0) sbar startl endl = sbar{1}+%resids
linreg(noprint) dseries startl+1 endl resids
# sbar{1} constant trends
compute tau=%tstats(1)
*
* Figure out which row to use
*
compute nrange=7
do i=1,6
   if n<=nlimits(i) {
      compute nrange=i
      break
   }
end do i
compute %cdstat=tau

if %defined(title)
   compute ltitle=title
else
   compute ltitle="Schmidt-Phillips Test (TAU) for a Unit Root in "+%LABEL(SERIES)

report(use=sreport,action=define,title=ltitle)
report(use=sreport,atrow=1,atcol=1,span) ltitle
report(use=sreport,atrow=2,atcol=1,span) "Regression Run From "+%datelabel(startl+1)+" to "+%datelabel(endl)
report(use=sreport,atrow=3,atcol=1) "Observations" %nobs
report(use=sreport,atrow=4,atcol=1,span) "Null of unit root. Reject in left tail."
report(use=sreport,atrow=6,atcol=1) "Signif. Level" "" "Critical Value"
report(use=sreport,atrow=7,atcol=1) "1%(**)" "" cv(p)(nrange,1)
report(use=sreport,atrow=8,atcol=1) "2.5%"   "" cv(p)(nrange,2)
report(use=sreport,atrow=9,atcol=1) "5%(*)"  "" cv(p)(nrange,3)
report(use=sreport,atrow=10,atcol=1) "10%"   "" cv(p)(nrange,4)
report(use=sreport,atrow=12,atcol=1,span) "LM Test Auxiliary Regression"
report(use=sreport,atrow=13,atcol=1) "Variable" "Coefficient" "T-Stat"
report(use=sreport,atrow=14,atcol=1) "SBAR{1}" %beta(1)
compute tag=fix(1+(%tstats(1)<cv(p)(nrange,3))+(%tstats(1)<cv(p)(nrange,1)))
report(use=sreport,atrow=14,atcol=3,tag=tag) %tstats(1)
report(use=sreport,atrow=15,atcol=1) "Constant" %beta(2) %tstats(2)
if p>=2
   report(use=sreport,atrow=16,atcol=1) "Trend" %beta(3) %tstats(3)
if p>=3
   report(use=sreport,atrow=17,atcol=1) "Trend^2" %beta(4) %tstats(4)
if p>=4
   report(use=sreport,atrow=18,atcol=1) "Trend^3" %beta(5) %tstats(5)

compute row=16+p
report(use=sreport,atrow=row,atcol=1,span) "Semiparametric Corrections for TAU"
report(use=sreport,atrow=row+1,atcol=1,span) "Bartlett Window Estimates of sigma^2"
report(use=sreport,atrow=row+2,atcol=1,span) "Schwert value of L4 ="+l4+", L12="+l12
report(use=sreport,atrow=row+4,atcol=1,align=right) "Lags" "Sigma^2" "tau"
sstats startl+1 endl resids^2>>sum2

do tcount = 0,lagsToUse
   mcov(lwindow=newey,lags=tcount) startl+1 endl resids
   # constant
   compute omega2=sum2/%cmom(1,1)
   compute ztau  =tau/sqrt(omega2)
   compute tag=fix(1+(ztau<cv(p)(nrange,3))+(ztau<cv(p)(nrange,1)))
   report(use=sreport,atrow=row+5+tcount,atcol=1) tcount %cmom(1,1)/%nobs
   report(use=sreport,atrow=row+5+tcount,atcol=3,special=tag) ztau
end do tcount
*
* Format the coefficients
*
report(use=sreport,action=format,atrow=14,torow=14+p,atcol=2,tocol=2,width=8)
*
* Variance estimates
*
report(use=sreport,action=format,atrow=row+5,atcol=2,tocol=2,width=8)
*
* t-statistics/critical values for t-stats
*
report(use=sreport,action=format,atcol=3,tocol=3,picture="*.###",align=decimal)
*
* Fix the alignment on the labels in column 3
*
report(use=sreport,action=format,atcol=3,tocol=3,atrow=13,torow=13,align=right)
report(use=sreport,action=format,atcol=3,tocol=3,atrow=row+4,torow=row+4,align=right)

if print==0
   return

report(use=sreport,action=show)
end spunit
