Commit cdc097ff authored by Andreas Steinberg's avatar Andreas Steinberg

added kpss test

parent df3132b0
......@@ -2,9 +2,15 @@
import logging
from collections import defaultdict
import numpy as num
from pyrocko import util, pile, model, config, trace, snuffling, gui_util
from pyrocko.fdsn import enhanced_sacpz, station as fs
from pyrocko.guts import Object, Tuple, String, Float, dump_all, load_all
import statsmodels.api as sm
from scipy.interpolate import interp1d
import statsmodels.tsa.stattools
guts_prefix = 'grond'
......@@ -304,6 +310,96 @@ class Dataset(object):
raise NotFound('no response', (net, sta, loc, cha))
else:
raise NotFound('multiple responses', (net, sta, loc, cha))
def stationarity_test_adf(series):
#testing the stationarity of the data by checking if the null hypothesis of a unit root can be recjected
d_order0=statsmodels.tsa.stattools.adfuller(series)
#check if the p value is critcall (only esimation!) after MacKinnon, J.G. 2010.
if d_order0[0]> d_order0[4]['5%']:
#print 'Time Series is nonstationary'
sto= 'FALSE'
else:
#print 'Time Series is stationary'
sto= 'TRUE'
return sto
def kpssTest(x, regression="LEVEL",lshort = True):
"""
KPSS Test for Stationarity
"""
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],"likely non stationary")
if (v<x[-1]):
return (y[-1],"likley trend stationary (or level if adf test rejects null hypthosis)")
if y[0]> y[-1]:
print 'Time Series is nonstationary'
else:
print 'Time Series is stationary'
def get_waveforms_raw(self, obj, tmin=None, tmax=None, tpad=0.):
net, sta, loc = self.get_nsl(obj)
......@@ -352,6 +448,16 @@ class Dataset(object):
if toffset_noise_extract != 0.0:
for tr in trs:
tr.shift(-toffset_noise_extract)
noise_prop=True
if noise_prop != True:
extracted_obs = []
extracted_obs = tr.chop(tmin+(random.uniform(1.0, 3.0)),tmin+(random.uniform(5.0, 10.0)), inplace=False)
tshift = (random.uniform(tmin, tmax))
extracted_obs.shift(tshift)
tr.add(extracted_obs, interpolate=True)
if len(trs) == 1:
return trs[0]
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment