Procedure KPSS series start end
*
*  Revision schedule:
*    01/92 Written by Phil Rothman
*    08/93 ETA(tau) part corrected by Phil Rothman
*    02/94 Updated and revised by Tom Doan.
*    02/96 On detrending regression, replaced "/" with "startl endl"
*    11/06 Updated by Tom Doan to include PRINT, LAGS and DET options
*
*  Procedure calculates Kwiatowski, Phillips, Schmidt & Shin ETA(mu) and
*  ETA(tau) statistics.  With ETA(mu), test H0: {X(t)} is stationary around
*  a level.  With ETA(tau), test H0: {X(t)} is trend stationary.  Details
*  found in KPSS (1992), "Testing the Null Hypothesis of Stationarity Against
*  the Alternative of a Unit Root:  How Sure Are We That Economic Time Series
*  Have a Unit Root?", JOURNAL OF ECONOMETRICS, 54, 159-178.
*
*  Syntax:
*    @KPSS( options ) series  start  end
*
*  Options:
*     LMAX=maximum number of lags in window [not used]
*     LAGS=specific number of lags to use [4]
*     [PRINT]/NOPRINT
*     DET=[CONSTANT]/TREND
*
*  Defines:
*    %NOBS   = number of observations
*    %CDSTAT = (final) test statistic
*
type series series
type integer start end
*
option integer lmax
option integer lags 4
option switch  print 1
option choice  det   1  constant trend
*
local series sums trend res
local real sum2
local integer n startl endl
local string detlabel testlabel
local integer lowlimit highlimit ikount
*
inquire(series=series) startl>>start endl>>end
compute n = endl-startl+1

*
* Get deviations from deterministic components
*
if det==1 {
   diff(center) series startl endl res
   compute detlabel="about level",testlabel="Eta(mu)"
}
else {
   set trend = t
   linreg(noprint) series startl endl res
   # constant trend
   compute detlabel="about trend",testlabel="Eta(tau)"
}
*
*  Partial sums --> sums
*
acc res startl endl sums
sstats startl endl sums**2>>sum2
*
*  Get consistent estimator of sigma^2, see KPSS,
*   for l = 1,2,...,lmax
*
if print {
   display
   display "KPSS Test for Stationarity "+detlabel+", Series" %l(series)
   display "From " %datelabel(startl) "to" %datelabel(endl)
   display "Observations" n
}

if %defined(lmax)
   compute lowlimit=0,highlimit=lmax
else
   compute lowlimit=lags,highlimit=lags

do ikount = lowlimit,highlimit
   mcov(lwindow=newey,lags=ikount) startl endl res
   # constant
   compute %cdstat = (sum2 / n) / %cmom(1,1)
   if print
      display "Lags=" ikount @12 testlabel+"=" ###.### %cdstat
end do ikount
compute %nobs=n
if .not.print
   return
*
*  First: set up test for ETA(mu)
*
display
if det==1 {
   display "Critical Level:  0.10   0.05   0.025  0.01"
   display "Critical Value:  0.347  0.463  0.574  0.739"
}
else {
   display "Critical Level:  0.10   0.05   0.025  0.01"
   display "Critical Value:  0.119  0.146  0.176  0.216"
}
display
end


