PKϨMLp
codemeta.json{"@context": ["https://doi.org/doi:10.5063/schema/codemeta-2.0", "http://schema.org"], "@type": "SoftwareSourceCode", "provider": {"url": "https://www.comses.net", "@type": "Organization", "name": "CoMSES Network (CoMSES)", "@id": "https://www.comses.net"}}PKϨMLV
docs/readme_stationarity.txtStationarity test
###########################Readme########################################
# #
# Jakob Grazzini #
# #
# jakob.grazzini(at)unito.it #
# #
#########################################################################
Loading the statonaity_test.py module you can use the wwRunTest Function.
The module imports the following python modules which you need to install to run the test:
random
rpy --> to be upgraded to rpy2
math
operator
numpy
In the module there is a function "process" used to create the experimental processes, AR(1). The arguments are the length of the process (len_p) and teta.
to run the test:
from stationarity_test import *
wwRunTest(time_series, len_sub)
the wwRunTest function requires a time series and a parameters that indicate how long shall the sub time series be (it has to be an natural number)
The function
1) divide the time series in w sub time series of the specified length.
2) creates a list with the mean of each sub time series
3) compute the mean of the overall time series
If the moment we are testing is constant then the sample averages should be well explained by the overall mean, that is the sample averages are randomply distributed around the overall mean.
The function creates a new list. Then it checks in order, from the first to the last, whether the sample averages are above or below the overall mean. If the sample average it adds a 1 to the new list, if below it adds a -1 to the list (the numbers are symbols, so it is indifferent what we use, we could also use a and b). Given the new list it computes the Runs and the U statistic. Given the mean, the variance and the asymptotic normality, the function computes the p_value of U, that means the probability that U comes from the distribution that is true under the null. If the p_value < alfa, we reject the null of stationarity.
To test different moments we have to transform the time series. If for example we want to test the secon non centered moment we can take the time series and compute a new time series consisting in the squared elements of the first time series. The we can use the test on the new time series.
The test automatically standardizes the number of runs. The outcome is directly Z, the standardized value of the number of Runs. Therefore you will reject the two tail test if abs(Z)>1.96, ans will reject the one tail (left) test if Z < -1.645. Both the values refers to an alfa=0.05
The experiments visualized in the paper can be reproduced using the expe.py
PKϨMLa,code/stationarity_test.py_.txt# -*- coding: utf-8 -*-
### I run the non parametric tests:
### median test on raw data for the strict stationarity (Gibbons 1985)
### median test on moments
### runs test (fit-test) on moments
import random
#import rpy
#from statlib import stats
import math
import numpy
from operator import itemgetter
def process(len_p,teta):
proc=[]
y0=0
for i in range(len_p):
y= teta*y0 + random.uniform(-1,1)
proc.append(y)
y0=y
return proc
def wwRunTest(time_serie,sub_len):
'run test, checks whether the ensemble moments explains the sample moments. For second third.. moments pre transform the time serie'
r = []
rr = []
for i in range(len(time_serie)):
r.append(time_serie[i])
if (i+1)%sub_len==0:
rr.append(r)
r = []
sample_means=[]
for i in rr:
sample_means.append(sum(i)*len(i)**-1)
#rpy.r.plot(sample_means,ylab='')
overallMean = sum(time_serie)*len(time_serie)**-1
#overallMean = rpy.r.median(time_serie)
runList=[]
runs = 0
state = 0
above=0
below=0
for i in sample_means:
if i < overallMean:
below+=1
runList.append(-1)
if i < overallMean and state != -1:
state = -1
runs+=1
if i > overallMean:
above+=1
runList.append(+1)
if i > overallMean and state != 1:
state = 1
runs += 1
#mu = ((2*above*below)*(len(sample_means))**-1)+1
#print mu
alfa= above*below**-1
mu = 2*above*(1+alfa)**-1 + 1
#print mu1
#print runs
#print mu, 'mean'
## using exact mean and variance instead of asymptotic because the test is more precise with few windows
s = ((mu-1)*(mu-2))*(len(sample_means)-1)**-1
sigma = math.sqrt(s)
#print s
#s1 = 4*alfa*above*(1+alfa)**-3
#print s1
#sigma = math.sqrt(s1)
z = (runs-mu)*sigma**-1
#p_value=1-rpy.r.pnorm(abs(z),0,1)
#rpy.r.plot(sample_means, ylab='', type='l')
mm=[]
for i in range(len(sample_means)):
mm.append(overallMean)
#rpy.r.lines(mm,col = 'red')
return z
PKϨMLp
codemeta.jsonPKϨMLV
,docs/readme_stationarity.txtPKϨMLa,gcode/stationarity_test.py_.txtPK