*
*  PPUNIT  series  start  end
*     Computes one of the Phillips-Perron modifications to the Dickey-
*     Fuller unit root tests.  References:
*       Phillips, "Time Series Regressions with a Unit Root", Econometrica,
*         1987, pp. 277-301
*       Phillips and Perron, "Testing for a Unit Root in Time Series
*         Regressions", Biometrika, 1988.
*
*     The (approximate) critical values for t-test form are from
*       MacKinnon, "Critical Values for Cointegration Tests", Long-Run
*         Economic Relationships, R.F. Engle and C.W.J. Granger, eds,
*         London, Oxford, 1991, pp 267-276
*
*  Options:
*     [TTEST]/NOTTEST  Computes the regression t test, as opposed to the
*       T(rho-1) test.
*     TREND/[NOTREND]
*       This determine whether an intercept or intercept and trend is included.
*       With TREND, estimated model includes 1 and t. By default, only the 1 is included
*     LAGS=number of lags in spectral estimation window [4]
*
*   Defines:
*    %NOBS   = number of regression observations + 1 (tables are based upon this)
*    %RHO    = the lag coefficient
*    %CDSTAT = test statistic
*
*   Various other statistics (such as %DURBIN and %RSS) will be set to the values
*   from the regression run.
*
*   Revision Schedule:
*    Written November, 1989
*    Updated January, 1992 to Version 4.0
*    Updated June, 2000 for Version 5.0
*    Updated January 2003 to improve output, change default to TTEST.
*
procedure ppunit series start end
type series series
type integer start end
*
option integer lags  4
option switch  trend 0
option switch  ttest 1
*
local integer startl endl regnobs
local series sdiff resids strend
local real su st tstat dxx teststat vardiff
local string descript
local vect critvalues
*
inquire(series=series) startl>>start endl>>end
*
set sdiff  startl+1 endl = series-series{1}
set strend startl+1 endl = t
linreg(entries=5+(trend>0),noprint) series startl+1 endl resids
# series{1} constant strend
compute nobs=%nobs+1
if trend {
   compute descript='With intercept and trend'
   compute critvalues=||-3.9638-8.353/nobs-47.44/(nobs**2),-3.4126-4.039/nobs-17.83/(nobs**2),-3.1279-2.418/nobs-7.58/(nobs**2)||
}
else {
   compute descript='With intercept'
   compute critvalues=||-3.4335-5.999/nobs-29.25/(nobs**2),-2.8621-2.738/nobs-8.36/(nobs**2),-2.5671-1.438/nobs-4.48/(nobs**2)||
}
compute regnobs=%nobs
mcov(damp=1.0,lags=lags) startl+1 endl resids
# constant
compute %seesq=%seesq*%ndf/%nobs
compute su=sqrt(%seesq)
compute st=sqrt(%cmom(1,1)/%nobs)
compute tstat=(%beta(1)-1)/sqrt(%seesq*%xx(1,1))
if trend
  compute dxx=1./%det(%xx)
else {
  stats(noprint) series startl endl-1
  compute dxx=(%nobs-1)*%variance
  }
compute vardiff=(st**2-su**2)

display
display 'Phillips-Perron Unit Root Test, Series' %l(series)
display 'Regression Run From' %datelabel(startl+lags+1) 'to' %datelabel(endl)
display 'Observations' nobs
display descript ', with ' lags 'lags in the error process'
if .not.trend
  if .not.ttest
    compute teststat=regnobs*(%beta(1)-1)-.5*vardiff*regnobs**2/dxx
  else
    compute teststat=tstat*su/st-.5*vardiff*regnobs/(st*sqrt(dxx))
else
  if .not.ttest
    compute teststat=regnobs*(%beta(1)-1)-vardiff*float(regnobs)**6/(24*dxx)
  else
    compute teststat=tstat*su/st-vardiff*float(regnobs)**3/(4*st*sqrt(3*dxx))

if ttest {
   display 'T-test statistic' teststat
   display 'Critical values: 1%=' *.### critvalues(1) '5%=' critvalues(2) '10%=' critvalues(3)
}
else {
   display 'T(p-1) statistic' teststat
}
display
compute %nobs=nobs,%cdstat=teststat,%rho=%beta(1)
end


