*
* @PPUNIT  series  start  end
*
* Computes one of the Phillips-Perron modifications to the Dickey-
* Fuller unit root tests.
*
* Parameters:
*   series     series to analyze
*   start end  range of <<series>> to use. By default, the defined range of <<series>>
*
* Options:
*   [TTEST]/NOTTEST  Computes the regression t test, as opposed to the T(rho-1) test.
*   DET=[CONSTANT]/TREND
*   TREND/[NOTREND]
*     Either of these determine whether an intercept or intercept and trend is
*     included. With DET=TREND or TREND, estimated model includes 1 and t. By default,
*     only the 1 is included.
*   LAGS=number of lags in spectral estimation window [4]
*   TABLE/[NOTABLE] show sensitivity table (all lags 0 to LAGS)
*   [PRINT]/NOPRINT
*
* 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.
*
* References:
*   Phillips(1987), "Time Series Regressions with a Unit Root", Econometrica, pp.
*     277-301
*   Phillips and Perron(1988), "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
*
* Revision Schedule:
*   11/1989 Original writing
*   01/1992 Updated to version 4.0
*   06/2000 Updated to version 5.0
*   01/2003 Updated to improve output, change default to TTEST.
*   11/2003 Correction to reported value of "regressions run from"
*            output--should be "startl+1", not "starl+lags+1"
*   11/2006 Added DET and PRINT options for compatibility with other procedures.
*   03/2009 Updated to use REPORT, generate sensitivity table
*   01/2010 Declare report as local.
*
procedure ppunit series start end
type series series
type integer start end
*
option integer lags  4
option switch  trend 0
option switch  ttest 1
option choice  det   1  constant trend
option switch  print 1
option switch  table 0
*
local integer  startl endl regnobs nobs
local series   sdiff resids strend
local real     su st tstat dxx teststat vardiff
local string   descript
local vect     critvalues
local integer  dettrend
local integer  lwidth uwidth row baserow tlag tag
local report   ppunitr
*
* Check for invalid syntax
*
if .not.%defined(series) {
   disp "Syntax: @ppunit(options)  SERIES start end"
   return
}

inquire(series=series) startl<<start endl<<end
*
set sdiff  startl+1 endl = series-series{1}
set strend startl+1 endl = t

if trend.or.det==2
   compute dettrend=1
else
   compute dettrend=0
*
* Do the calculations that aren't dependent upon the window width
*
linreg(entries=5+dettrend,noprint) series startl+1 endl resids
# series{1} constant strend
compute nobs=%nobs+1
compute regnobs=%nobs
compute %seesq=%seesq*%ndf/%nobs
compute su=sqrt(%seesq)
compute tstat=(%beta(1)-1)/sqrt(%seesq*%xx(1,1))
if dettrend
   compute dxx=1./%det(%xx)
else {
   stats(noprint) series startl endl-1
   compute dxx=(%nobs-1)*%variance
}

if dettrend {
   if ttest {
      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 and trend, T(rho-1) variant"
      compute critvalues=||-29.5,-21.8,-18.3||
   }
}
else {
   if ttest {
      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)||
   }
   else {
      compute descript="With intercept, T(rho-1) variant"
      compute critvalues=||-20.7,-14.1,-11.3||
   }
}

if table
   compute lwidth=0,uwidth=lags
else
   compute lwidth=uwidth=lags

report(use=ppunitr,action=define)
report(use=ppunitr,atrow=1,atcol=1,span) "Phillips-Perron Test for a Unit Root for "+%label(series)
report(use=ppunitr,atrow=2,atcol=1,span) "Regression Run From "+%datelabel(startl+1)+" to "+%datelabel(endl)
report(use=ppunitr,atrow=3,atcol=1,span) "Observations "+%nobs
report(use=ppunitr,atrow=4,atcol=1,span) descript
report(use=ppunitr,atrow=6,atcol=1) "Sig Level" "Crit Value"
report(use=ppunitr,atrow=7,atcol=1) "1%(**)"    critvalues(1)
report(use=ppunitr,atrow=8,atcol=1) "5%(*)"     critvalues(2)
report(use=ppunitr,atrow=9,atcol=1) "10%"       critvalues(3)
report(use=ppunitr,atrow=11,atcol=1,align=center) "Lags" "Statistic"
compute row=baserow=12
do tlag=lwidth,uwidth
   mcov(lwindow=newey,lags=tlag) startl+1 endl resids
   # constant
   compute st=sqrt(%cmom(1,1)/%nobs)
   compute vardiff=(st**2-su**2)
   if .not.dettrend
      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))

   report(use=ppunitr,atrow=row,atcol=1) tlag
   compute tag=fix(1+(teststat<critvalues(1))+(teststat<critvalues(2)))
   report(use=ppunitr,atrow=row,atcol=2,special=tag) teststat
   compute row=row+1
end do tlag
report(use=ppunitr,action=format,atrow=baserow,atcol=2,width=8,align=decimal)
if print
   report(use=ppunitr,action=show)
*
compute %nobs=nobs,%cdstat=teststat,%rho=%beta(1)
end

