MCMC初心者がpymc3で株価の期待日次リターンを推定してみる

# 背景
 [Pythonで体験するベイズ推論](https://www.amazon.co.jp/Python%E3%81%A7%E4%BD%93%E9%A8%93%E3%81%99%E3%82%8B%E3%83%99%E3%82%A4%E3%82%BA%E6%8E%A8%E8%AB%96-PyMC%E3%81%AB%E3%82%88%E3%82%8BMCMC%E5%85%A5%E9%96%80-%E3%82%AD%E3%83%A3%E3%83%A1%E3%83%AD%E3%83%B3-%E3%83%87%E3%83%93%E3%83%83%E3%83%89%E3%82%BD%E3%83%B3-%E3%83%94%E3%83%AD%E3%83%B3/dp/4627077912)がMCMC(マルコフ連鎖モンテカルロ法)初心者の私にも分かりやすいうえに, 事例が豊富でかなりの良書だった.  だけど, 株価のボラティリティ(分散)の計算が一部動かないのと, 主観的な推定だったものを, より客観的な推定にしてみたかったので, やってみた. また, 書籍ではpymc2だったけど, APIが豊富なpymc3を使った.

# 問題設定
2017年のAAPL, GOOG, TSLA, AMZNの株価(終値)を使って, 期待日次リターンをMCMCで推定する. また, 共分散行列を計算して, 株価の相関を求める.

# 環境設定
+ python3.6
+ numpy==1.13.3
+ pandas==0.22.0
+ matplotlib==2.1.1
+ pymc3==3.2
+ pandas_datareader==0.5.0  

# 株価取得
pandas_datareaderを使って, Yahoo Financeから株価をもってくる.

```python3:
import datetime
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pandas_datareader.data as web


start = datetime.datetime(2017, 1, 1)
end = datetime.datetime(2017, 12, 31)
aapl = web.DataReader("AAPL", 'yahoo', start, end)
goog = web.DataReader("GOOG", 'yahoo', start, end)
tsla = web.DataReader("TSLA", 'yahoo', start, end)
amzn = web.DataReader("AMZN", 'yahoo', start, end)
```

# データ処理
今回使うのは終値だけ. それぞれの終値を別のDataFrameに格納する.

```python3:
close = pd.concat([goog['Close'], aapl['Close'], tsla['Close'], amzn['Close']], axis=1)
close.columns = ['goog', 'aapl', 'tsla', 'amzn']
close.plot(grid=True)
plt.ylabel('Close')
```
![close.png](https://qiita-image-store.s3.amazonaws.com/0/160903/7d52bb9a-ce0a-a21a-c2bd-876cf16072ca.png)

日次リターンを計算する.

```python3:
up_rate = close.diff(axis=0) / close
up_rate.dropna(inplace=True)
up_rate.plot(grid=True)
plt.ylabel('up_rate')
```

![up_rate.png](https://qiita-image-store.s3.amazonaws.com/0/160903/13c33cf4-6e29-a04c-1c31-51a6c65161d6.png)

わかりにくいので, 初期投資1ドルに対するリターンをプロットする.

```python3:
prod = (1 + up_rate).cumprod(axis=0) - 1
prod.plot(grid=True)
plt.ylabel('Return of $1 on first date, x100%')
```

![return.png](https://qiita-image-store.s3.amazonaws.com/0/160903/d73d1351-14ac-6da5-ab51-90f27e11f785.png)

日次リターンのヒストグラム.

```python3:
up_rate.hist(bins=30, sharex=True, sharey=True)
```
![hist.png](https://qiita-image-store.s3.amazonaws.com/0/160903/70b8b0d5-4437-6f13-71b8-9de85b6fc410.png)

# MCMCによる推定
いよいよ本題. MCMCについては[書籍](https://www.amazon.co.jp/Python%E3%81%A7%E4%BD%93%E9%A8%93%E3%81%99%E3%82%8B%E3%83%99%E3%82%A4%E3%82%BA%E6%8E%A8%E8%AB%96-PyMC%E3%81%AB%E3%82%88%E3%82%8BMCMC%E5%85%A5%E9%96%80-%E3%82%AD%E3%83%A3%E3%83%A1%E3%83%AD%E3%83%B3-%E3%83%87%E3%83%93%E3%83%83%E3%83%89%E3%82%BD%E3%83%B3-%E3%83%94%E3%83%AD%E3%83%B3/dp/4627077912)を見ていただきたい. 雑に説明すると, MCMCとは, 人工的に作り出したデータを利用して, ベイズの定理によって, 事後分布(今回の場合, 2017年の株価における期待日次リターンの分布)を得る手法である.

今回は, 期待日次リターンと株価のボラティリティを推定する. 書籍では, 決め打ちの事前分布を使っていたが, 株のプロではないので, もう少し客観的な事前分布を使いたい. そこで, 次のようにモデリングを行った.

+ 期待日次リターンは正規分布とする.
+ $\mu$は一様分布から推定する.
+ $\sigma^2$はLKJ相関分布から推定する.

初めは, Wishert分布を使おうとしたが, Wishert分布の分散に一様分布などを使うと, LKJ相関分布を使えと怒られた.　LKJ相関分布については, pymc3のドキュメント内の[このページ](http://docs.pymc.io/notebooks/LKJ.html)が分かりやすかった.

# 計算する

MCMCのパートは, 以下の通り.

```python3:
import pymc3 as pm

with pm.Model() as model:
    # LKJCholeskyCovは下三角行列の要素をリストで返す
    packed_L = pm.LKJCholeskyCov("packed_L", n=4, eta=1., sd_dist=pm.HalfCauchy.dist(beta=2.5))
    # 下三角行列に変換する
    L = pm.expand_packed_triangular(4, packed_L, lower=True)
    # 共分散行列にする
    sigma = pm.Deterministic('sigma', L.dot(L.T))
    
    prior_mu = pm.Uniform("prior_mu", -1, 1, shape=4)
    mu= pm.Normal("returns", mu=prior_mu, sd=1, shape=4)
    
    obs = pm.MvNormal("observed_returns", mu=mu, chol=L, observed=up_rate)
    step = pm.NUTS()
    trace = pm.sample(50000, step)
```

NUTSを使っているので, pymc2より遅いかな.
NUTSについては, [ドキュメント](http://docs.pymc.io/notebooks/getting_started.html)を見てください(あまり理解してない).

pymc3には強力な図化APIが用意されていて, 以下の通り.

```python3:
pm.traceplot(trace)
```
![traceplot.png](https://qiita-image-store.s3.amazonaws.com/0/160903/02b38063-bcf1-fe2c-6bda-53f0398f9ad0.png)


```python3:
pm.forestplot(trace)
```
![forestplot.png](https://qiita-image-store.s3.amazonaws.com/0/160903/8948fa74-f0bf-e132-0286-83c7981ba8f8.png)



pymc3では$\hat{R}$も簡単に出してくれる. pymc2でも計算してみたが, pymc3のほうが収束がすごくいい. NUTSのおかげ.
期待日次リターンの事後分布は以下のようになる.

```python3:
mu_samples = pd.DataFrame(trace['returns'],
                          columns=['goog', 'aapl', 'tsla', 'amzn'])
mu_samples.hist(bins=30, sharex=True, sharey=True)
```
![returns.png](https://qiita-image-store.s3.amazonaws.com/0/160903/301e471c-6a97-9ea8-9f33-e9fe9496b605.png)


それぞれの株価の相関や分散は次のようになる.

```python3:
sigma_post = trace["sigma"].mean(axis=0)

# 共分散行列から相関行列にする関数
def cov2corr(A):
    d = np.sqrt(A.diagonal())
    A = ((A.T / d).T) / d
    return A

plt.subplot(1, 2, 1)
plt.imshow(cov2corr(sigma_post), interpolation="none", cmap=plt.cm.hot)

plt.xticks(np.arange(4), up_rate.columns)
plt.yticks(np.arange(4), up_rate.columns)
plt.colorbar(orientation="vertical")
plt.title("(Mean posterior) correlation matrix")

plt.subplot(1, 2, 2)
plt.bar(np.arange(4), np.sqrt(np.diag(sigma_post)),
        color="#348ABD",
        alpha=0.7)
plt.xticks(np.arange(4), up_rate.columns)
plt.title("(Mean posterior) variances of daily stock returns")
plt.xlabel('Stock')
plt.ylabel('Variance')

plt.tight_layout()
```
![corr.png](https://qiita-image-store.s3.amazonaws.com/0/160903/0114e55e-51b4-e36a-3695-e0c8a0bee23d.png)


GOOGとAMZNが相関が高いのは興味深い. どちらもクラウドやってるからかな. 

# まとめ
このぐらいだったら, 頻度主義に走って計算したほうが早い気がする(普通に相関行列を計算してもほぼ同じ). だけど, 真骨頂は損失関数を使ったベイズ点推定や, さらにモデルを複雑にした階層ベイズとからしい. 期待値とかが分布で出てくるのも, アバウト人間の私としては精神衛生上よろしい.  
興味を持った方は, ぜひとも[書籍](https://www.amazon.co.jp/Python%E3%81%A7%E4%BD%93%E9%A8%93%E3%81%99%E3%82%8B%E3%83%99%E3%82%A4%E3%82%BA%E6%8E%A8%E8%AB%96-PyMC%E3%81%AB%E3%82%88%E3%82%8BMCMC%E5%85%A5%E9%96%80-%E3%82%AD%E3%83%A3%E3%83%A1%E3%83%AD%E3%83%B3-%E3%83%87%E3%83%93%E3%83%83%E3%83%89%E3%82%BD%E3%83%B3-%E3%83%94%E3%83%AD%E3%83%B3/dp/4627077912)を読んで欲しい. ダークマターのベイズ推定とかカッコイイ.