*
* @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:
*   DET=[CONSTANT]/TREND
*     The component around which the null of stationarity is assumed.
*
*   LMAX=maximum number of lags in window [not used]
*   LAGS=specific number of lags to use [4]
*     The KPSS test requires an estimate of the long-run variance, which is
*     computed using a Bartlett (Newey-West) window. The value of the test
*     statistic will change based upon the number used. If you use the LMAX
*     option, the procedure will show the test statistic for the range of
*     lags up to LMAX. Otherwise, it will do only the test with the number
*     of lags chosen.
*
*   [PRINT]/NOPRINT
*   TITLE=title for report ["KPSS Test for Stationarity"]
*
* Defines:
*   %NOBS   = number of observations
*   %CDSTAT = test statistic (using highest lag)
*
* 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
*    01/2010 Declare report as local.
*    07/2010 Switch to use FILTER(REMOVE=...)
*    07/2010 Add TITLE option, TITLE to REPORT
*    08/2019 Add additional information to output
*
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
option string  title
*
local series   sums 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
local report   kpssr
local string   ltitle ldesc
*
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
*
clear res sums
if det==1 {
   filter(remove=mean) series startl endl res
   compute detlabel="about Level",testlabel="Eta(mu)",critvalues=cv(1)
}
else {
   filter(remove=trend) series startl endl res
   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
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="KPSS Test for Stationarity, Series "+%l(series)

compute ldesc="Null is Stationary "+detlabel+". Reject for large values."

report(use=kpssr,action=define,title=ltitle)
report(use=kpssr,atrow=1,atcol=1,span) ltitle
report(use=kpssr,atrow=2,atcol=1,span) "From "+%datelabel(startl)+" to "+%datelabel(endl)
report(use=kpssr,atrow=3,atcol=1,span) "Observations "+n
report(use=kpssr,atrow=4,atcol=1,span) ldesc
report(use=kpssr,atrow=6,atcol=1,align=center) "Sig Level" "Crit Value"
report(use=kpssr,atrow=7,atcol=1) "1%(**)"    critvalues(1)
report(use=kpssr,atrow=8,atcol=1) "2.5%"      critvalues(2)
report(use=kpssr,atrow=9,atcol=1) "5%(*)"     critvalues(3)
report(use=kpssr,atrow=10,atcol=1) "10%"       critvalues(4)
report(use=kpssr,atrow=12,atcol=1,align=center) "Lags" "TestStat"
*
*  Get consistent estimator of sigma^2, see KPSS,
*   for l = 1,2,...,lmax
*
compute row=baserow=13
do ikount = lowlimit,highlimit
   mcov(lwindow=newey,lags=ikount) startl endl res
   # constant
   compute %cdstat = (sum2/n)/%cmom(1,1)
   report(use=kpssr,atrow=row,atcol=1) ikount
   compute tag=fix(1+(%cdstat>critvalues(1))+(%cdstat>critvalues(3)))
   report(use=kpssr,atrow=row,atcol=2,special=tag) %cdstat
   compute row=row+1
end do ikount
compute %nobs=n
report(use=kpssr,action=format,atrow=6,atcol=2,picture="*.###",align=decimal)
if print
   report(use=kpssr,action=show)
end
