時系列の適合度検定
疑問:(離散)時系列データがあった時に、Wiener過程や幾何Brown運動で記述するのが適切か検定する方法
Mandelbrotが、綿花の価格の変化率がパレート分布に従うことを示した云々はよく聞くけど、実際、Brown運動・幾何Brown運動であるかどうかの検定をしてみたくなった
とりあえず、時系列は置いておいて、標本データの正規性検定。標本データが、正規分布(や指定した確率分布)に従っているか検定する手段として
(1)Shapiro-Wilk検定(正規性検定。標本数が小さい時向き？)
(2)Anderson-Darling検定(任意の分布に適用可能)
(3)一標本Kolmogorov-Smirnov検定(任意の分布に適用可能)
(4)Lilliefors検定(正規性に特化してKS検定を改良。標本数が大きい正規性検定向き？)
(5)D'AgostinoのK2乗検定(歪度と尖度による正規性検定。実用性は低い？)
(6)Jarque-Bera検定(歪度と尖度による正規性検定。実用性は低い？)
(7)カイ二乗適合度検定(正規性のみ。区間の切り方の任意性が気持ち悪い)
などがあるらしい。多分、大体Rに実装されてるし、scipyにもあるっぽい。意味はないけど、Shapiro-Wilk検定をpythonで実装してみた。そこらに落ちてる実装をみると、サンプル数が5000を超えるときは、やめとけとコメントに書いてあることが多い
from math import *
#正規累積分布関数
def normcdf(u):
d1=0.0498673470
d2=0.0211410061
d3=0.0032776263
d4=0.0000380036
d5=0.0000488906
d6=0.0000053830
if u > 0:
temp=(((((d6*u+d5)*u+d4)*u+d3)*u+d2)*u+d1)*u+1.0
return 1.0 - 0.5/(temp**16)
else:
u=-u
temp=(((((d6*u+d5)*u+d4)*u+d3)*u+d2)*u+d1)*u+1.0
return 0.5/(temp**16)
#指定した下側確率に対応する正規分布のZ値(qnorm(normcdf(x)) = x)
def qnorm(p):
A0_p,A1_p,A2_p,A3_p = 3.3871327179, 5.0434271938E+01, 1.5929113202E+02, 5.9109374720E+01
B1_p , B2_p , B3_p = 1.7895169469E+01, 7.8757757664E+01, 6.7187563600E+01
C0_p,C1_p,C2_p,C3_p = 1.4234372777E+00, 2.7568153900E+00, 1.3067284816E+00, 1.7023821103E-01
D1_p , D2_p = 7.3700164250E-01, 1.2021132975E-01
E0_p,E1_p,E2_p,E3_p = 6.6579051150E+00, 3.0812263860E+00, 4.2868294337E-01,1.7337203997E-02
F1_p , F2_p = 2.4197894225E-01, 1.2258202635E-02
q = p -0.5
if fabs(q) <= 0.425:
r = 0.180625 - q*q
return q*(((A3_p*r + A2_p)*r+A1_p)*r+A0_p)/(((B3_p*r+B2_p)*r+B1_p)*r+1.0)
else:
if q<0.0:
r=p
else:
r = 1.0 -p
if r<=0.0:return 0.0
r = sqrt(-log(r))
if r <= 5.0:
r = r - 1.6
normal_dev = (((C3_p * r + C2_p) * r + C1_p) * r + C0_p) / ((D2_p * r + D1_p) * r + 1.0)
else:
r = r-5.0
normal_dev = (((E3_p * r + E2_p) * r + E1_p) * r + E0_p) / ((F2_p * r + F1_p) * r + 1.0)
if q < 0.0:
return -normal_dev
else:
return normal_dev
#Calculates Shapiro-Wilk normality test and P-value for sample sizes 3 <= n <= 5000.
"""
Ported from FORTRAN77 code
ALGORITHM AS R94 APPL. STATIST. (1995) VOL.44, NO.4
"""
def swtest(dat):
def poly(cs , x):
if len(cs)==1:return c[0]
p = x *cs[-1]
if len(cs)!=2:
for i in xrange(len(cs)-2):
p = (p+cs[-(i+2)])*x
return cs[0]+p
N = len(dat)
assert(N > 2)
#-- compute normalize coeffcients
weights = [] #coeffients for test
if N==3:
weights = [sqrt(2)/2]
else:
weights = [qnorm((n-0.375)/(N+0.25)) for n in xrange(1,1+int(N/2))]
summ2 = 2.0*sum([x*x for x in weights])
ssumm2 = sqrt(summ2)
rsn = 1.0/sqrt(N)
a1 = poly([0.0, 0.221157, -0.147981, -2.07119, 4.434685, -2.706056] , rsn) - weights[0]/ssumm2
if N>5:
i1 = 2
a2 = -weights[1]/ssumm2 + poly([0.0, 0.042981, -0.293762, -1.752461, 5.682633, -3.582633],rsn)
fac = sqrt((summ2 - 2.0*weights[0]*weights[0] - 2.0*weights[1]*weights[1])/(1.0 - 2.0*a1*a1 - 2.0*a2*a2))
weights[0] = a1
weights[1] = a2
else:
i1 = 1
fac = sqrt((summ2 - 2.0*weights[0]*weights[0])/(1.0 - 2.0*a1*a1))
weights[0] = a1
for i in xrange(i1 , N/2):
weights[i] = -weights[i]/fac
#-- compute W statistics
tmp = sorted(dat)
mean = sum(dat)/N
try:
W = (sum([weights[i]*(tmp[N-i-1] -tmp[i]) for i in xrange(N/2)])**2)/sum([(x-mean)*(x-mean) for x in dat])
except ZeroDivisionError:
return 0.0
#-- compute significance level
if N==3:
return 6*( asin(sqrt(W)) - asin(sqrt(0.75)) )/pi
elif N<=11:
y = log(1-W)
gamma = poly([-0.2273E1, 0.459E0] , N)
if y>=gamma:return 1E-19
y = -log(gamma - y)
m = poly([0.5440E0, -0.39978E0, 0.25054E-1, -0.6714E-3], N)
s = exp(poly([0.13822E1, -0.77857E0, 0.62767E-1, -0.20322E-2], N))
else:
y,xx = log(1-W),log(N)
m = poly([-0.15861E1, -0.31082E0, -0.83751E-1, 0.38915E-2],xx);
s = exp(poly([-0.4803E0, -0.82676E-1, 0.30302E-2], xx))
return 1.0-normcdf((y-m)/s)
"""
import random
if __name__=="__main__":
alpha = 0.05
fail_cnt = 0
for _ in xrange(1000):
dat = [random.gauss(0.0 , 1.0) for _ in xrange(50)]
if swtest(dat) <= alpha:
fail_cnt += 1
print "%2.2f(%%) failed" % (fail_cnt*0.1)
"""
α=0.05とすると、当然ながら、5%くらいの確率で、正規乱数を、正規でないと判定する。人為的に外れ値を入れてやると、急激にP値が下がるので、変化点検出とかに使えそうでもある
次に、時系列が、Brown運動や幾何Brown運動であるかどうかの判定。まず、(弱)定常性の検定をやらないといけない気がする。これも色々方法があるようだけど、今は置いておく。で、定常性が成立するなら、増分ないし対数値の増分が正規分布に従っていることを言えばいい気がする。11月くらいから現在までの日経平均は、割と定常的な気がするので、日経平均(の終値)100日分を使って実験して見る。以下は、コードで、P値は0.844827691028らしいので、対数価格の変化率は正規分布に従っている可能性を棄却できない
from datetime import *
import urllib2,traceback,time
from bs4 import BeautifulSoup
#-- Yahoo!ファイナンスからデータ取得
def getYF(kcode , npage=1):
def S(s):
return int(s.replace(',' , ''))
def F(s):
return float(s.replace(',' , ''))
ret = []
try:
today = datetime.today()
params = [datetime.strftime(today,'%Y') , datetime.strftime(today,'%m') , datetime.strftime(today,'%d')]
url_format = 'http://info.finance.yahoo.co.jp/history/?code=%s&sy=2003&sm=01&sd=01&ey=%s&em=%s&ed=%s&tm=d&p=%d'
for p in xrange(1,npage+1):
url = url_format % (str(kcode) , params[0] , params[1] , params[2] , p)
res = urllib2.urlopen(url)
soup = BeautifulSoup( res.read() )
for t in soup.find_all('table'):
if 'boardFin' in t['class']:
for c,line in enumerate( t.find_all('tr') ):
if c==0:continue
dat = [e.text for e in line.find_all('td')]
if len(dat)==7: #-- 一般銘柄
t = time.strptime(dat[0].encode('utf-8') , '%Y年%m月%d日')
d = date(t.tm_year , t.tm_mon , t.tm_mday).toordinal()
ret.append( (d ,S(dat[1]) , S(dat[2]) , S(dat[3]), S(dat[4]) , S(dat[5]) ,S(dat[6])) )
elif len(dat)==5: #-- 日経平均など
t = time.strptime(dat[0].encode('utf-8') , '%Y年%m月%d日')
d = date(t.tm_year , t.tm_mon , t.tm_mday).toordinal()
ret.append( (d ,F(dat[1]) , F(dat[2]) , F(dat[3]), F(dat[4]) , None ,F(dat[4])) )
except:
print "error@%s" % str(kcode)
traceback.print_exc()
return []
return ret
from math import *
def test(kcode):
tmp = getYF(kcode , npage=2) #-- 100日分
return swtest([log(cur[-1])-log(prev[-1]) for (cur,prev) in zip(tmp , tmp[1:])])けど、対数をとらない、普通の増分も正規分布に従うことを棄却できない(P値0.640581714301)ので、これだと、Brown運動か幾何Brown運動か区別できない。個別銘柄も、いくつかやってみると、こっちは結構微妙な結果になる
日経平均:0.8448276910278005
(株)三菱ＵＦＪフィナンシャル・グループ:0.03019741434584633
日本マクドナルドホールディングス(株):7.70282535533795e-08
トヨタ自動車(株):0.013263024199066442
ガンホーオンラインエンターテイメント(株):0.05917450794749601
まぁ、そもそもこんなやり方でいいのか分からないけど、飽きたので、終わり
