Stationarity Test with KPSS (Kwiatkowski–Phillips–Schmidt–Shin) in Python

In addition to augmented Dickey–Fuller test (ADF), KPSS (Kwiatkowski–Phillips–Schmidt–Shin) is widely used to verify stationarity of a signal. Python statsmodels contains ADF test, but I could not find any implementation of KPSS in python. Therefore, I converted R implementation of kpss.test in tseries library to python as follows (Click here to download the source code) :

"""
Created on Sat Jan-03-2015
@author: Deniz Turan (http://denizstij.blogspot.co.uk/)

"""
import numpy as np
import statsmodels.api as sm
from scipy.interpolate import interp1d

def kpssTest(x, regression="LEVEL",lshort = True):
    """
    KPSS Test for Stationarity

    Computes the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test for the null hypothesis that x is level or trend stationary.

    Parameters
    ----------
    x : array_like, 1d
        data series
    regression : str {'LEVEL','TREND'} 
        Indicates the null hypothesis and must be one of "Level" (default) or "Trend". 
    lshort : bool
        a logical indicating whether the short or long version of the truncation lag parameter is used.

    Returns
    -------
    stat : float
        Test statistic
    pvalue : float
        the p-value of the test.
    usedlag : int
        Number of lags used.

    Notes
    -----
    Based on kpss.test function of tseries libraries in R.
    
    To estimate sigma^2 the Newey-West estimator is used. If lshort is TRUE, then the truncation lag parameter is set to trunc(3*sqrt(n)/13), otherwise trunc(10*sqrt(n)/14) is used. The p-values are interpolated from Table 1 of Kwiatkowski et al. (1992). If the computed statistic is outside the table of critical values, then a warning message is generated.

    Missing values are not handled.

    References
    ----------
    D. Kwiatkowski, P. C. B. Phillips, P. Schmidt, and Y. Shin (1992): Testing the Null Hypothesis of Stationarity against the Alternative of a Unit Root. Journal of Econometrics 54, 159--178.

    Examples
    --------
    x=numpy.random.randn(1000)  #   is level stationary    
    kpssTest(x)

    y=numpy.cumsum(x)           # has unit root    
    kpssTest(y)

    z=x+0.3*arange(1,len(x)+1)   # is trend stationary
    kpssTest(z,"TREND")

    """
    x = np.asarray(x,float)
    if len(x.shape)>1:
        raise ValueError("x is not an array or univariate time series")
    if regression not in ["LEVEL", "TREND"]:
        raise ValueError(("regression option %s not understood") % regression)

    n = x.shape[0]
    if regression=="TREND":
        t=range(1,n+1)
        t=sm.add_constant(t)
        res=sm.OLS(x,t).fit()
        e=res.resid
        table=[0.216, 0.176, 0.146, 0.119]
    else:
        t=np.ones(n)
        res=sm.OLS(x,t).fit()
        e=res.resid
        table=[0.739, 0.574, 0.463, 0.347]

    tablep=[0.01, 0.025, 0.05, 0.10]
    s=np.cumsum(e)
    eta=np.sum(np.power(s,2))/(np.power(n,2))
    s2 = np.sum(np.power(e,2))/n
    if lshort:
        l=np.trunc(3*np.sqrt(n)/13)
    else:
        l=np.trunc(10*np.sqrt(n)/14)
    usedlag =int(l)
    s2=R_pp_sum(e,len(e),usedlag ,s2)
    
    stat=eta/s2

    pvalue , msg=approx(table, tablep, stat)
    
    print "KPSS Test for ",regression," Stationarity\n"
    print ("KPSS %s=%f" % (regression, stat))
    print ("Truncation lag parameter=%d"% usedlag )
    print ("p-value=%f"%pvalue )

    if msg is not None:
        print "\nWarning:",msg 
    
    return ( stat,pvalue , usedlag )


def R_pp_sum (u, n, l, s):
    tmp1 = 0.0
    for i in range(1,l+1):
        tmp2 = 0.0
        for j in range(i,n):
            tmp2 += u[j]*u[j-i]
        tmp2 = tmp2*(1.0-(float(i)/((float(l)+1.0))))
        tmp1 = tmp1+tmp2
    
    tmp1 = tmp1/float(n)
    tmp1 = tmp1*2.0
    return s + tmp1

def approx(x,y,v):
    if (v>x[0]):
        return (y[0],"p-value smaller than printed p-value")
    if (v

Stationary Tests : Augmented Dickey–Fuller (ADF), Hurst Exponent, Variance Ratio (VRTest) of Time Series within R

Several trading strategies (momentum, mean reverting, ...) are based on if data is stationary or not. In this text, i demonstrate how to test it statistically.

library("quantmod") # for downloading fx data
library("pracma") # for hurst exponent 
library("vrtest") # variance ratio test
library("tseries") # for adf test
library("fUnitRoots")  # for adf test

## first lets fetch USD/CAD data for last 5 years
getFX("UDSCAD")
usdCad=unclass(USDCAD) # unwrap price column
# estimate log return
n=length(usdCad)
usdcadLog=log(usdCad[1:n])

## First use Augmented Dickey–Fuller Test (adf.test) to test USD/CAD is statationary
>adfTest(usdCad, lag=1)

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 1
  STATISTIC:
    Dickey-Fuller: 0.2556
  P VALUE:
    0.6978 

Description:
 Sun Nov 03 16:47:27 2013 by user: deniz

## As you see above, null hypothesis (unit root) can not be rejected with p-value ~70 %

## So we demonstrated it is not stationary. So if there is a trend or mean reverting. Hurst exponent (H) can be used for this purpose (Note Hursy exponent relies on that random walk diffuse in proportion to square root of time.). 
#Value of H can be interpreted such as: 
#H=0.5:Brownian motion (Random walk)
#H<0.5:Mean reverting  
#H>0.5:Trending 
> hurst(usdcadLog) Hurst exponent
[1] 0.9976377

#So, USDCAD is in trending phase. 

## Another way to test stationary is to use V:
> vrtest::Auto.VR(usdcadLog)
[1] 83.37723
> vrtest::Lo.Mac(usdcadLog,c(2,4,10))
$Stats
           M1       M2
k=2  22.10668 15.98633
k=4  35.03888 25.46031
k=10 56.21660 41.58861

## Another way to analyse stationary condition is via linear regression in which we will try to establish if there is a link between data and diff(data(t-1))

>deltaUsdcadLog=c(0,usdcadLog[2:n]-usdcadLog[1:n-1])
> r=lm(deltaUsdcadLog ~ usdcadLog)
> summary(r)

Call:
lm(formula = deltaUsdcadLog ~ usdcadLog)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0121267 -0.0013094 -0.0000982  0.0012327  0.0103982 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -7.133e-05  1.306e-04  -0.546   0.5853  
usdcadLog    8.772e-03  5.234e-03   1.676   0.0944 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.002485 on 498 degrees of freedom
Multiple R-squared:  0.005608, Adjusted R-squared:  0.003611 
F-statistic: 2.808 on 1 and 498 DF,  p-value: 0.0944

> r$coefficients[2]
  usdcadLog 
0.008771754 

## Coefficient (Beta) gives clue about if there is mean reverting. If it is negative, there is a mean reverting.  As you see above, it is positive, therefore as we already concluded before, it is trending. If it was negative, we would use following to find out half life of mean revertion:

>-log(2)/r$coefficients[2]