*
* @KPSS( options ) series  start  end
*
* Calculates the KPSS unit root 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.
*
* Parameters:
*   series     series to analyze
*   start end  range of <<series>> to use. By default, the defined range of <<series>>
*
* 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
*
* Reference:
*   Kwiatowski, Phillips, Schmidt & Shin(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?", J. of Econometrics, vol 54, pp 159-178.
*
* Revision schedule:
*    01/1992 Written by Phil Rothman
*    08/1993 ETA(tau) part corrected by Phil Rothman
*    02/1994 Updated and revised by Tom Doan.
*    02/1996 On detrending regression, replaced "/" with "startl endl"
*    11/2006 Updated by Tom Doan to include PRINT, LAGS and DET options
*    03/2009 Switch to use of REPORT
*
Procedure KPSS series start end
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
local integer  row baserow tag
local vector   critvalues
*
fixed vect[vect] cv(2)(4)
  0.739 0.574   0.463  0.347
  0.216  0.176  0.146  0.119
*
* Check for invalid syntax
*
if .not.%defined(series) {
   disp "Syntax: @kpss(options)  SERIES start end"
   return
}

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

if %defined(lmax)
   compute lowlimit=0,highlimit=lmax
else
   compute lowlimit=lags,highlimit=lags
*
* Get deviations from deterministic components
*
if det==1 {
   diff(center) series startl endl res
   compute detlabel="about Level",testlabel="Eta(mu)",critvalues=cv(1)
}
else {
   set trend = t
   linreg(noprint) series startl endl res
   # constant trend
   compute detlabel="about Trend",testlabel="Eta(tau)",critvalues=cv(2)
}
*
*  Partial sums --> sums
*
acc res startl endl sums
sstats startl endl sums**2>>sum2
*
* Set up header of table
*
report(action=define)
report(atrow=1,atcol=1,span) "KPSS Test for Stationarity "+detlabel+", Series "+%l(series)
report(atrow=2,atcol=1,span) "From "+%datelabel(startl)+" to "+%datelabel(endl)
report(atrow=3,atcol=1,span) "Observations "+n
report(atrow=5,atcol=1) "Sig Level" "Crit Value"
report(atrow=6,atcol=1) "1%(**)"    critvalues(1)
report(atrow=7,atcol=1) "2.5%"      critvalues(2)
report(atrow=8,atcol=1) "5%(*)"     critvalues(3)
report(atrow=9,atcol=1) "10%"       critvalues(4)
report(atrow=11,atcol=1,align=center) "Lags" "TestStat"
*
*  Get consistent estimator of sigma^2, see KPSS,
*   for l = 1,2,...,lmax
*
compute row=baserow=12
do ikount = lowlimit,highlimit
   mcov(lwindow=newey,lags=ikount) startl endl res
   # constant
   compute %cdstat = (sum2/n)/%cmom(1,1)
   report(atrow=row,atcol=1) ikount
   compute tag=fix(1+(%cdstat>critvalues(1))+(%cdstat>critvalues(3)))
   report(atrow=row,atcol=2,special=tag) %cdstat
   compute row=row+1
end do ikount
compute %nobs=n
report(action=format,atrow=6,atcol=2,width=8,align=decimal)
if print
   report(action=show)
end



