はじめに
前回の投稿でScipyを使用してpingの応答時間の分布を確認しました。
結果として、指数分布だと思っていた分布が対数正規分布によりフィットする結果となりました。どの分布に従うかは、見た目だけで判断できないことがわかります。
Scipyで処理時間の分布を確かめる
本来は理論から「これこれこうだからこの分布に従うはず」と仮説を立て、実測値と比較して検証できれば良いのですが、そのような知識を残念ながら持ち合わせていません。どうすればよいのか調べていたところ、下記の記事を見つけました。これなら理論を導き出す知識がなくても、フィットする分布を見つけられそうです。
とりあえず経験分布になんらかのモデルをフィットさせたい Python/Scipy
本記事では、日経平均株価の変化率の分布を調べてみようと思います。
日経平均株価の取得
まずは日経平均株価の取得です。
株価の取得はpandas_datareaderを使用して米国のYahoo!Financeから取得しています。
2011年から2020年の10年分の株価を取得してみます。
import pandas_datareader
import datetime
ticker = "^N225"
start_date = datetime.date(2011, 1, 1)
end_date = datetime.date(2020, 12, 31)
df = pandas_datareader.data.DataReader(
ticker, 'yahoo', start_date, end_date)
調整済み終値を使って10年間の株価の推移をグラフにしてみます。
import matplotlib.pyplot as plt
df["Adj Close"].plot()
plt.show()
2011年は10000円くらいだった株価が、2020年には27500円くらいまであがってます。
変化率の分布
では、変化率のヒストグラムを作成してみます。
df["pct_change"] = df["Adj Close"].pct_change()
df["pct_change"].plot.hist(bins=50)
plt.show()
左右対称の釣り鐘型のグラフになりました。正規分布よりも尖度が大きそうです。
さて、これは何分布でしょうか？
分布のフィッティング
では、下記記事を真似して分布の種類を一つずつフィッティングしてきます。
とりあえず経験分布になんらかのモデルをフィットさせたい Python/Scipy
まずはbest_fit_distribution関数です。オリジナルから少し修正しています。
import numpy as np
from scipy import stats
import warnings
def best_fit_distribution(data, bins=50, ax=None):
"""Model data by finding best fit distribution to data"""
# データからヒストグラムを作成する
y, x = np.histogram(data, bins=bins, density=True)
x = (x + np.roll(x, -1))[:-1] / 2.0
# 以下の分布でフィッティングする
DISTRIBUTIONS = [
stats.alpha,
stats.anglit,
stats.arcsine,
stats.beta,
stats.betaprime,
stats.bradford,
stats.burr,
stats.cauchy,
stats.chi,
stats.chi2,
stats.cosine,
stats.dgamma,
stats.dweibull,
stats.erlang,
stats.expon,
stats.exponnorm,
stats.exponweib,
stats.exponpow,
stats.f,
stats.fatiguelife,
stats.fisk,
stats.foldcauchy,
stats.foldnorm,
stats.frechet_r,
stats.frechet_l,
stats.genlogistic,
stats.genpareto,
stats.gennorm,
stats.genexpon,
stats.genextreme,
stats.gausshyper,
stats.gamma,
stats.gengamma,
stats.genhalflogistic,
stats.gilbrat,
stats.gompertz,
stats.gumbel_r,
stats.gumbel_l,
stats.halfcauchy,
stats.halflogistic,
stats.halfnorm,
stats.halfgennorm,
stats.hypsecant,
stats.invgamma,
stats.invgauss,
stats.invweibull,
stats.johnsonsb,
stats.johnsonsu,
stats.ksone,
stats.kstwobign,
stats.laplace,
stats.levy,
stats.levy_l,
# stats.levy_stable, ← 計算が終わらないのでコメントアウト
stats.logistic,
stats.loggamma,
stats.loglaplace,
stats.lognorm,
stats.lomax,
stats.maxwell,
stats.mielke,
stats.nakagami,
stats.ncx2,
stats.ncf,
stats.nct,
stats.norm,
stats.pareto,
stats.pearson3,
stats.powerlaw,
stats.powerlognorm,
stats.powernorm,
stats.rdist,
stats.reciprocal,
stats.rayleigh,
stats.rice,
stats.recipinvgauss,
stats.semicircular,
stats.t,
stats.triang,
stats.truncexpon,
stats.truncnorm,
stats.tukeylambda,
stats.uniform,
stats.vonmises,
stats.vonmises_line,
stats.wald,
stats.weibull_min,
stats.weibull_max,
stats.wrapcauchy,
]
# [分布名,分布パラメータ,誤差]を格納するためのリスト
l = []
# 分布毎にそれぞれパラメータを推測してみる
for i, distribution in enumerate(DISTRIBUTIONS):
print(f"{i+1}/{len(DISTRIBUTIONS)} {distribution.name: <20}", end="")
# 分布によってはフィットできないこともあるので、
# フィットできなければ、passする
try:
with warnings.catch_warnings():
warnings.filterwarnings('ignore')
# 分布をフィットさせる
params = distribution.fit(data)
# わかりやすい用にパラメータをばらばらにする
arg = params[:-2]
loc = params[-2]
scale = params[-1]
# フィットさせた確率密度関数を計算して、分布にフィットさせる
pdf = distribution.pdf(x, loc=loc, scale=scale, *arg)
# 残差平方和を誤差として計算する
sse = np.sum(np.power(y - pdf, 2.0))
print(f"done sse={sse}")
# l
l.append([
distribution,
params,
sse,
])
except Exception:
print("not fit")
l.sort(key=lambda x: x[2])
return l
フィッティングします。
x = df["pct_change"].dropna().values
fit_list = best_fit_distribution(x)
では、順位の発表です。
for i, fit in enumerate(fit_list):
print(f"{i+1} {fit[0].name}")
結果
1 dweibull
2 dgamma
3 gennorm
4 laplace
5 johnsonsu
6 t
7 nct
8 hypsecant
9 tukeylambda
10 burr
11 genlogistic
12 logistic
13 foldcauchy
14 cauchy
15 exponnorm
16 powerlognorm
17 powernorm
18 exponweib
19 loggamma
20 johnsonsb
21 pearson3
22 beta
23 lognorm
24 foldnorm
25 norm
26 gengamma
27 rice
28 fatiguelife
29 chi
30 nakagami
31 invgamma
32 alpha
33 recipinvgauss
34 invgauss
35 chi2
36 frechet_r
37 weibull_min
38 frechet_l
39 weibull_max
40 gumbel_l
41 gompertz
42 invweibull
43 gumbel_r
44 exponpow
45 betaprime
46 fisk
47 loglaplace
48 kstwobign
49 triang
50 maxwell
51 gausshyper
52 rayleigh
53 halfgennorm
54 genhalflogistic
55 genexpon
56 powerlaw
57 mielke
58 ksone
59 uniform
60 bradford
61 wald
62 truncexpon
63 gilbrat
64 levy_l
65 halfnorm
66 halflogistic
67 lomax
68 halfcauchy
69 levy
70 pareto
71 expon
72 ncf
73 genpareto
74 anglit
75 cosine
76 rdist
77 semicircular
78 vonmises_line
79 genextreme
80 ncx2
81 arcsine
82 erlang
83 gamma
84 f
85 reciprocal
86 truncnorm
87 vonmises
88 wrapcauchy
ベスト5は下記となりました。
- dweibull
- dgamma
- gennorm
- laplace
- johnsonsu
pdf(確率密度)を計算する関数を準備して、56位までをグラフ化してみます。
def make_pdf(distribution, params, start, end, size=500):
arg = params[:-2]
loc = params[-2]
scale = params[-1]
x = np.linspace(start, end, size)
y = distribution.pdf(x, loc=loc, scale=scale, *arg)
return x, y
COL_NUM = 4
fl = fit_list[0:56]
fig, axes = plt.subplots(math.ceil(len(fl) / COL_NUM), COL_NUM,
figsize=(10, 30), sharex=True, sharey=True, dpi=300)
plt.subplots_adjust(top=0.95, bottom=0.05, hspace=0.5)
for i, fit in enumerate(fl):
distribution, params, sse = fit
x_fit, y_fit = make_pdf(distribution, params, min(x), max(x), bins)
ax, ay = divmod(i, COL_NUM)
axes[ax, ay].hist(x, bins=bins, density=True)
axes[ax, ay].plot(x_fit, y_fit)
axes[ax, ay].set_title(f"No.{i+1} {distribution.name}")
plt.savefig("stock_distribution_result.png")
最後に
今回の結果では、ベスト5は下記のようになりました。
- dweibull
- dgamma
- gennorm
- laplace
- johnsonsu
ここでは詳細を書きませんが、私の好きなAppleで同じ事を行ったところ、結果は
- tukeylambda
- johnsonsu
- t
- nct
- hypsecant
となりました。
銘柄によってフィットする分布は違ってくるようです。
ソースは GitHub に置いています。
こちらのソースでは柄名を引数で指定できるようにしています。