統計,確率のお勉強

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

Study Probability & Statistics

確率統計の理論と実践

ベイズ統計:パラメータの点推定

久しぶりなので少し短めに・・・
大学卒業して少し統計離れてましたけど、土日やることもないので、前々から気になっていたベイズ統計の勉強でも少しづつ始めていこうと思います。

パラメータの点推定

点推定
未知であるパラメータの値をデータから推定すること

通常の統計学なら不偏推定量とか、最小分散不偏推定量とか、十分統計量とかありました。それのベイズ版です。

ベイズ分析におけるパラメータの点推定は、パラメータに真の値と点推定の乖離度をある尺度(損失関数)ではかり、この損失関数が出来る限り小さくなるようにすることを考えます。

ここからは、パラメータ  \pi の点推定を \delta 損失関数を  L(\pi, \delta) と書くことにする。ベイズ分析によく使われる損失関数としては次のような関数があります。

\begin{align}
2乗誤差損失({\rm quadratic \, loss}) &: L(\pi, \delta) = (\pi - \delta)^2 \\
絶対誤差損失({\rm absolute \, loss}) &: L(\pi, \delta) = |\pi - \delta|^2 \\
0 - 1 損失 ({\rm 0-1 \, loss}) &: L(\pi, \delta) = 1 - \boldsymbol{1}_{\pi} (\delta)
\end{align}

ここで、  \boldsymbol{1}_{\pi} (\delta) \delta = \pi の時 1 それ以外のとき 0となる指示関数。

パラメータの真の値はわかっていないので、ベイズ分析の点推定では、損失関数  L(\pi, \delta) の期待値を  \pi の事後分布で評価したもの

\begin{equation}
R(\delta | D) = E_{p (\pi | D)} [L(\pi, \delta)] = \int_0^1 L(\pi, \delta) p(\pi | D) d\pi
\end{equation}

を考え、これを出来るだけ小さくするように点推定を選択します。 R(\delta | D)期待損失 (expected loss) と呼びます。これを最小化問題として定式化したものが、未知のパラメータ  \pi の点推定  \delta^*

\begin{equation}
\delta^* = \arg \min_{0 < \delta < 1} R(\delta | D) = \arg \min_{0 < \delta < 1} \int_0^1 L(\pi, \delta) p (\pi | D) d\pi
\end{equation}

となります。"  \arg " は 「  \delta^* が \min_{0 < \delta < 1} R(\delta | D) という最小化問題の解である」という意味です。

参考文献

中妻照雄(2007):『ファイナンスライブラリー10 入門ベイズ統計学』, 朝倉書店

正規性の検定(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倍の差がある。

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




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