統計,確率のお勉強

統計学を主に勉強しています。勉強したことをアウトプットしていきます。 (※数式はMathJaxにより描画されています。ロードに少し時間がかかることがあります。)

Study Probability & Statistics

確率統計の理論と実践

正規性の検定(Shapiro-Wilk検定)

正規性の検定

データを分析するにあたり、

  1. データが正規分布に従う
  2. データが独立な標本である

といった仮定を置くことは多い。そのような場合に分析をする際、これら二つの仮定が満たされているか確認する必要が出てくる。そのための手法として統計的仮説検定がある。今回はその中の Shapiro-Wilk検定Pythonで株価収益率に対してやってみようと思う。

二種類の誤り

仮説検定には2つの誤りがある。

第1種の誤り
帰無仮説が正しいにも関わらず、帰無仮説を棄却してしまう誤り。

第2種の誤り
帰無仮説が正しくないにも関わらず、帰無仮説を採択してしまう誤り。

通常、第一種の誤り確率をある水準(有意水準)以下に抑えた状態で、第二種の誤り確率をできるだけ小さくするようにする。そのため、帰無仮説が採択された場合は、積極的に帰無仮説を正しいとみているわけではなく、「帰無仮説が間違っているとは言えない」というような消極的な解釈となる。

Shapiro-Wilk検定

私が所持している本の中にはこの検定に関する記述のあるものが見つからなかったため、今回は
シャピロ–ウィルク検定 - Wikipedia
の内容でお茶を濁すことにする。

統計学における、シャピロ–ウィルク検定(シャピロ–ウィルクけんてい)とは、 標本  x_1, \ldots, x_nが正規母集団からサンプリングされたものであるという帰無仮説を検定する検定である。この検定方法は、サミュエル・シャピロとマーティン・ウィルクによって、1965年に発表された。検定統計量は、

$$
W = \frac{\left(\sum_{i=1}^n a_i x_{(i)} \right)^2}{\sum_{I=1}^n (x_i - \bar{x})^2}
$$

ただし

  • x_{(i)} は、i番目の順序統計量。
  • \bar{x} は標本平均。
  • 定数  a_i は次の式によって与えられる。

$$
(a_1, \ldots, a_n) = \frac{m' V^{-1} }{\left( m' V^{-1} V^{-1} m \right)^{1/2}}
$$

ただし、

$$
m = (m_1, \ldots, m_n)'
$$

m_1, \ldots, m_nは、標準正規分布からサンプリングされた独立同分布の確率変数の順序統計量の期待値であり、「V」は、この順序統計量の分散共分散行列である。帰無仮説は、「W」が小さすぎる場合に棄却される。

※一応Wikiの引用だが、少し書き換えてある。

兎にも角にもこれを用いて正規性を検定する。データは株式の収益率を使う。

Pythonでは次のようにすれば、Shapiro-Wilkの検定が可能である.

from spicy import stats
stats.shapiro(df['QUICK'])  # df['QUICK'] はデータフレームdfに入っているQUICKの収益率データを表す

結果は以下のようになった。

統計量W p値
QUICK 0.9260867834091187 1.0152270032202054e-22
日立 0.9507756233215332 1.2722628221227775e-18
武田薬品 0.9450637698173523 1.0990565268293669e-19

どれもp値がほとんど0に近いので、棄却される。つまり、どれも正規分布に従っているとは言えない事がわかる。

各業種と日経平均株価についてヒストグラムとガウシアンフィッティングで得られた曲線(もし、そのデータが正規分布に従うならばフィットするであろう正規分布な曲線)を一緒に描画した図をみてみよう。

f:id:doratai:20180207230849p:plain

これをみて見ると、正規分布よりも尖った分布になっているのがみて取れる。

参考文献の本では、4銘柄中3銘柄帰無仮説が採択されてたのに...

正規分布に従っているという仮定には無理がある事がわかった。

参考文献

[1] 横内大介, 青木義充(2014):『現場ですぐ使える時系列データ分析〜データサイエンティストのための基礎知識〜』技術評論社

収益率の相関を見る

相関係数

二つの銘柄X, Yの収益率データを \{x_1,\ ldots, x_T \}, \{y_1,\ldots, y_T\} とする。この時相関係数

$$
corr(X,Y) = \frac{ \sum_{t=1}^T (x_t - \bar{x})(y_t - \bar{y}) }{\sqrt{\sum_{t=1}^T (x_t - \bar{x})^2} \sqrt{\sum_{t=1}^T (y_t - \bar{y})^2} }
$$

で定義される。今回は散布図行列をプロットして見ることにする。
今回は

の3銘柄で試してみた。

%matplotlib inline
import matplotlib.pyplot as plt
import quandl
import pandas as pd
import seaborn as sns
import numpy as np
from scipy.stats import norm

# データの取得

df_quick = quandl.get('TSE/4318')  # Quick
df_hitachi = quandl.get('TSE/6501') # 日立製作所
df_takeda = quandl.get('TSE/4502')  #武田薬品工業

# 対数差収益率に変換
Quick_RoR = np.diff(np.log(df_quick.loc["2013-07-16":, ['Close']]), axis=0)*100
Hitachi_RoR = np.diff(np.log(df_hitachi.loc["2013-07-16":, ['Close']]), axis=0)*100
Takeda_RoR = np.diff(np.log(df_takeda.loc["2013-07-16":, ['Close']]), axis=0)*100
Nikkei_RoR = np.diff(np.log(df_Nikkei.loc["2013-07-16":, ['Close Price']]), axis=0)*100

# データフレームに収益率をまとめる
df = pd.DataFrame(np.concatenate((Quick_RoR, Hitachi_RoR, Takeda_RoR), axis=1), columns=['QUICK', 'HITACHI', 'TAKEDA'])

# seaboard で散布図行列を出力
sns.pairplot(df, kind='reg',size=8)

得られた結果が以下

f:id:doratai:20180207003104p:plain

これを見る限り、武田薬品と日立の相関が比較的強いように見える。実際に相関を計算すると

df.corr()
QUICK HITACHI TAKEDA
QUICK 1.000000 0.324614 0.307317
HITACHI 0.324614 1.000000 0.511227
TAKEDA 0.307317 0.511227 1.000000

となる。日立と武田薬品は約0.5と相関があり、日立とQUICKと比べて大きい事がわかる。



少し短いけどここまでにします。

<時系列データ分析>ボラティリティを見る〜Pythonと株価データを使ってお勉強〜

ボラティリティ

ちゃんと本を読むまでボラティリティ のことを単純に標準偏差だと思っていたが、どうやら定義がいくつかあるようだ。詳しくはないのでWikipediaから引用する。

金融工学においてボラティリティ(volatility)とは、広義には資産価格の変動の激しさを表すパラメータ。

ボラティリティ - Wikipedia

細かいところはWikiをみてください。

とりあえずの理解としては株価のバラツキ具合だと思うことにする。株価のバラツキが大きいということは投資した場合に株価の値上がりで儲ける可能性も高いが、逆に失う可能性も高い。つまりはリスクが大きいということになる。ボラティリティは銘柄の特徴を捉えるのに重要な指標の一つのようだ。

ヒストリカル・ボラティリティ

ここではヒストリカル・ボラティリティの定義を確認する。ヒストリカル・ボラティリティは価格の対数差収益率の標準偏差で定義される。

t 時点の株価を P_t とした時の t 時点での対数差収益率を r_t とする。

$$
r_t = \log P_t - \log P_{t-1} = \log \frac{P_t}{P_{t-1}}
$$

そのため、ヒストリカル・ボラティリティn 個の収益率データ \{r_1, r_2, \ldots, r_n \} が得られているとき

$$
s = \sqrt{\frac{1}{n-1}\sum_{i=1}^n (r_i - \bar{r})^2 }
$$

で推定される。

対数差収益率をプロットする

各銘柄の変動のしやすさを視覚的に見るために対数差収益率をプロットしてみよう。
今回も前回同様銘柄は

の4つで行く。

import matplotlib.pyplot as plt
import quandl
import numpy as np

# データをquandlから取得
df_quick = quandl.get('TSE/4318')  # Quick
df_hitachi = quandl.get('TSE/6501') # 日立製作所
df_takeda = quandl.get('TSE/4502')  #武田薬品工業
df_Nikkei = quandl.get('NIKKEI/INDEX')  # 日経平均株価

# matplotlibの描画をかっこよく
plt.style.use('ggplot')

fig1 = plt.figure(figsize=(18, 10))
ax1_quick = fig1.add_subplot(221)
ax1_hitachi = fig1.add_subplot(222)
ax1_takeda = fig1.add_subplot(223)
ax1_Nikkei = fig1.add_subplot(224)

ax1_quick.plot(np.diff(np.log(df_quick.loc["2013-07-16":, ['Close']]), axis=0)*100, label='Quick', color='b')
ax1_hitachi.plot(np.diff(np.log(df_hitachi.loc["2013-07-16":, ['Close']]), axis=0)*100,label='HITACHI', color='r')
ax1_takeda.plot(np.diff(np.log(df_takeda.loc["2013-07-16":, ['Close']]), axis=0)*100 , label='Takeda', color='g')
ax1_Nikkei.plot(np.diff(np.log(df_Nikkei.loc["2013-07-16":, ['Close Price']])*100 , axis=0), label='Nikkei', color='purple')

ax1_quick.set_title('Quick')
ax1_quick.set_ylim=(-20,20)

ax1_hitachi.set_title('HITACHI')
ax1_hitachi.set_ylim(-20,20)

ax1_takeda.set_title('TAKEDA')
ax1_takeda.set_ylim(-20, 20)

ax1_Nikkei.set_title('NIKKEI')
ax1_Nikkei.set_ylim(-20, 20)

plt.show()

y軸の幅はQuickの変動幅が一番大きかったのを何回か実行する中で確認したので
それに合わせた。y軸の値は パーセント(%)

f:id:doratai:20180202170051p:plain

これを見るとボラティリティは Quick > 日立 > 武田薬品日経平均株価 となっている事が視覚的にわかる。

Quickは \pm 20 %と大きく変動しているが、日経平均武田薬品\pm 5 % の幅に収まっている。

つまり、この4つの中だとQuickが最もボラタイル(ボラティリティが大きい)である事がわかる。

データの要約

視覚的に見ることはできた。続いては、平均や分散のような統計量を求めて見ることにする。

収益率の平均

各銘柄の収益率の平均を求めてみよう。
次のようにすれば計算できる。

quick_mean = np.mean(np.diff(np.log(df_quick.loc["2013-07-16":, ['Close']]), axis=0)*100)
hitachi_mean = np.mean(np.diff(np.log(df_hitachi.loc["2013-07-16":, ['Close']]), axis=0)*100)
takeda_mean = np.mean(np.diff(np.log(df_takeda.loc["2013-07-16":, ['Close']]), axis=0)*100)
Nikkei_mean = np.mean(np.diff(np.log(df_Nikkei.loc["2013-07-16":, ['Close Price']]), axis=0)*100)

print("Quick : " + str(quick_mean))
print("HITACHI : " + str(hitachi_mean))
print("TAKEDA : " + str(takeda_mean))
print("Nikkei : " + str(Nikkei_mean))

結果は

銘柄名 QUICK 日立 武田薬品 日経平均
平均(%) 0.163 0.027 0.029 0.043

いづれの銘柄の平均値も正の値をとるが、どれも0に近く、小さな値になっている。

※計算結果(平均)
Quick : 0.162560953612
HITACHI : 0.0269943597075
TAKEDA : 0.0285518660969
Nikkei : 0.0426793473652

ボラティリティを計算する

続いては、各銘柄のボラティリティ標準偏差)を求める。
次のようにすれば計算できる。

quick_std = np.std(np.diff(np.log(df_quick.loc["2013-07-16":, ['Close']]), axis=0)*100)
hitachi_std = np.std(np.diff(np.log(df_hitachi.loc["2013-07-16":, ['Close']]), axis=0)*100)
takeda_std = np.std(np.diff(np.log(df_takeda.loc["2013-07-16":, ['Close']]), axis=0)*100)
Nikkei_std = np.std(np.diff(np.log(df_Nikkei.loc["2013-07-16":, ['Close Price']]), axis=0)*100)

print("Quick : " + str(quick_std))
print("HITACHI : " + str(hitachi_std))
print("TAKEDA : " + str(takeda_std))
print("Nikkei : " + str(Nikkei_std))

結果は

銘柄名 QUICK 日立 武田薬品 日経平均
ボラティリティ(%) 2.713 1.896 1.342 1.316

先ほどの図の印象と数字が一致している事がわかります。Quickのボラティリティは2.713% であり、武田薬品が 1.342% であることから大きさとしてはQuickは武田薬品の約2倍であり、より価格変動が大きい銘柄である事がわかる。


※計算結果(ボラティリティ)
Quick : 2.7130695301
HITACHI : 1.8958143491
TAKEDA : 1.34199583079
Nikkei : 1.31599308226

豆○ば〜

QUICKって、これを書く際に参考にしている「現場ですぐ使える時系列データ分析ーデータサイエンティストのための基礎知識ー」の著者が所属している会社なんだよ〜

<時系列データ解析> 株価をプロット〜Pythonと株価データを使って勉強する〜

主に参考にする本

「現場ですぐ使える時系列データ分析 データサイエンティストのための基礎知識」を中心に勉強して行きたいと思う。
この本にした理由は

  • 内容がそこまで難しくない
  • Rで実行し、手を動かしながら学べるスタイル
  • 扱うデータが株価と自分の興味と合致している

の3つある。私は普段はRをほとんど使わないで言語はPython(Anaconda)、環境はJupyter notebookを使っているので、この本に載っている内容をPythonでやって行きたいと思う。

株価をプロットする

何を分析するにしても、まずはグラフを描画して見てみることはデータの特性を見るのにはいい手段だ。今回は株価データを単純に時系列プロットすることにする。

株価データを取得する

少し前だとyahoo financeやgoogle finance、あとは2017年で終わってしまったk-dbなどからデータを取ってきていたが今は現実的ではない。今だとQuandlというサイトからデータを取得するのがいいだろう。英語のサイトだがpythonから簡単なコードでデータを取得できるので使いやすい。


まずはライブラリのインポート

%matplotlib inline
import matplotlib.pyplot as plt
import quandl  #pipで別途インストールが必要
import pandas as pd

quandlがインストールされていない場合はターミナルで

>pip install quandl 

とやればインストールできる(おそらく)

株価データを取得する。今回は

の4つのデータを使うことにした。

データは次のコードで取得できる

df_quick = quandl.get('TSE/4318')  # Quick
df_hitachi = quandl.get('TSE/6501') # 日立製作所
df_takeda = quandl.get('TSE/4502')  #武田薬品工業
df_Nikkei = quandl.get('NIKKEI/INDEX')  # 日経平均株価

試しにQuickのデータを見ると

df_quick.head()
Date Open High Low Close Volume
2013-07-16 314.0 314.0 305.0 306.0 22100.0
2013-07-17 304.0 312.0 300.0 304.0 25700.0
2013-07-18 305.0 312.0 300.0 311.0 23300.0
2013-07-19 308.0 310.0 304.0 310.0 3600.0
2013-07-22 310.0 315.0 300.0 310.0 29100.0

株価をプロットする

のようになっている。今回は終値、つまり'Close'の値をプロットする。

# matplotlibの描画をかっこよく
plt.style.use('ggplot')

# 以下描画
fig = plt.figure(figsize=(18,8))
ax_quick = fig.add_subplot(221)
ax_hitachi = fig.add_subplot(222)
ax_takeda = fig.add_subplot(223)
ax_Nikkei = fig.add_subplot(224)

ax_quick.plot(df_quick.loc["2013-07-16": ,['Close']], label='Quick', color='b')
ax_hitachi.plot(df_hitachi.loc["2013-07-16": ,['Close']], label='HITACHI', color='r')
ax_takeda.plot(df_takeda.loc["2013-07-16": ,['Close']], label='Takeda', color='g')
ax_Nikkei.plot(df_Nikkei.loc["2013-07-16": ,['Close Price']], label='Nikkei', color='purple')

ax_quick.set_title('Quick')
ax_hitachi.set_title('HITACHI')
ax_takeda.set_title('TAKEDA')
ax_Nikkei.set_title('NIKKEI')

plt.show()

そしてこれが今回プロットした結果である。
f:id:doratai:20180202001816p:plain

どの銘柄も2016年の6月あたりから今に至るまで上がり続けている事が見て取れる。

ただ、縦軸のスケールが違うので、見た目の上がり方は同じでも、変動の幅は大きく違うことに注意したい。例えば、日立は2016年6月頃は400円ほどであり、2018年に入ることには900円程度になっているためその変動幅は約500円であるのに対し、武田薬品は2016年6月頃は約4000円であり、2018年に入る頃は6500円と、その幅は2500円ある。日立とは5倍の差がある。

銘柄ごとに株価の水準が異なるので、これをそのまま比較することはできない。そのためそのことを考慮した変動の割合を考える必要がある。




ただプロットしただけだけど今回はここまでにします。

学生時代にやっておくべきこと1位はなぜかいつも「海外旅行」

マ○ナビからメールが来た

社会人に聞いた!内定者が入社までにやっておいた方がいい遊びTOP10!
という記事がメールマガジンで紹介されてたので気になって見てみました。
そこには次のようにランキング形式でやっておいた方がいい遊びが紹介されていました。

1位 海外旅行をする 25人(9.6%)
2位 自分の好きな趣味に没頭する 17人(6.5%)
3位 国内旅行をする 13人(5%)
3位 合コン 13人(5%)
3位 飲み会 13人(5%)
6位 美術館・博物館めぐり 12人(4.6%)
7位 カラオケ 6人(2.3%)
7位 DVD鑑賞 6人(2.3%)
9位 映画鑑賞 5人(1.9%)
9位 アウトドア・スポーツ 5人(1.9%)
9位 キャンプ 5人(1.9%)

インドア派で休日は家から出ないのが基本な私にはなかなか難しい遊びが多い...

海外旅行       →そもそも日本から一回も出たこと無い。海外怖い。
趣味に没頭      →いつも思うんだけど趣味って何?
国内旅行       →お金無い。
合コン        →経験0。一度くらいは行ってみたいが呼ばれる可能性0。
飲み会        →どんだけ飲みたがりだよ...
美術館・博物館巡り  →これは面白そう。行ってみたい(博物館)。
カラオケ       →苦痛
DVD鑑賞       →寝る
アウトドア・スポーツ →インドアだっつてんだろ。
キャンプ       →あり得ない。

むしろ社会人になったほうがお金あるだろうし、旅行楽しそうだけどなあ...なんで学生時代に限定したがるんだ。
ふと、なぜこういうランキングでは海外旅行がいつも1位に来るんだろう?というどうでもいい疑問を抱いた。だってそこまでして行きたい理由分からないモン。必ず学生時代にやっておいたほうがいい事ランキングって必ず海外旅行1位だよね。もういいよ。そのランキングいらないよ。

理由1:「学生時代にやるべきこと=海外旅行」のイメージがつきすぎてる

いろいろなメディアで学生時代に海外旅行は必ず行くべき!みたいなことを言っている。これは先程のようなランキング形式のやつだったら絶対と言っていいほど出てくる。そういうものを普段から継続的に(この手の内容は毎年定期的に書かれてるよね?)目にすることで、多くの人に学生時代=海外旅行!みたいなイメージが既についてしまっている。

さて、アンケートの回答者は次のような質問を受ける。

「あなたが思う、学生時代にやっておいたほうがいいことはなんですか?」

ここで自分の周りの友達を思い出してほしい。そんなに皆が皆海外旅行に行ってる?少なくとも私の周りには片手で数えられるくらいしかいない。実際、学生の内に海外旅行に行く人達なんて一部だろう。その場合、アンケートの回答者は海外旅行に行ったことは無いが、「学生時代に海外旅行に行くべき!」という意識を刷り込まれた人達が必然的に多くなる。では先程の質問の回答はどうなるか?そう、

「海外旅行」

となるだろう。アンケートを回答する際、普通そんな真剣に考えない。そしたら当然こうなる用に思える。

理由2:社会人になって忙しくなる人達のほうが多い

当然だろ?と思うかもしれないがそうでないこともある。私は対して忙しく無いが、理系の大学に在籍しているため、特に実験などが多い研究室は週に休みがないのが当たり前みたいな教授の奴隷生活を送っている学生の話耳にすることもある。まあ、あくまで学生なので責任的な意味での大変さは無いだろうが時間的な意味での忙しさは、「社会人のほうが楽そう...」という言葉からなんとなく察する。まあ、ほとんど家から出ない私は低みの見物をするのみだが(紙とペンとパソコンと本あればなんとかなるのでそもそも研究室に行く意味ないし)。そんな私でもゼミのために勉強することは多いわ卒論書かなきゃだわと思ったより時間が無い(決して昔読んでた漫画やライトノベル等を読み返したり、アニメ見たりしているせいではない。そう、決して)。そんな中「時間のある学生時代の内に~」という言葉には違和感がある。学生時代はあくまで大学で「勉強」するための時間である以上海外「旅行」ではなく「留学」しろよ。どうせ海外いくなら。

大多数の人は社会人になってからのほうが時間的余裕がなくなるので、学生時代に海外旅行行っておけばよかったと後悔の念からそのような結果になるのだろう。でも思うんだよね。学生時代に海外行くようなアクティブな人なら、社会人になって忙しくても行くよね?

結局、学生時代も社会人になっても海外旅行行ってない人が答えている可能性...

理由3:1,2ときたら3がなきゃいけない気がする。

特に思いつかない。1、2ときたから一応様式美の3。

皆が学生時代に海外旅行に行くべきって言ってたから学生時代に海外旅行にいくべき。

つまりは、みんな(ソースは不明、ネットのランキングとか)が学生時代は海外旅行に行くべきって言ってたから学生時代は海外旅行にいくべきだと社会人が答え、それを聞いて社会に出た社会人がまた学生時代は海外旅行に行くべきって言ってたから学生時代は海外旅行にいくべきだと答え(以下無限ループ)により学生時代に海外旅行に行かなきゃならない感じになってる。

とりあえず学生時代にやるべきことランキングのまとめで

学生時代は海外旅行がおすすめ!!
海外旅行は学生時代のうちに!!

みたいなこと書くのはなんかなあ。旅行代理店か何かですか?それなら旅行プランへのリンク貼っとけよ。

注意

以上全て偏見でした。