pandas+matplotlibで世界時価総額トップ銘柄のアルファ値とベータ値をアニメ化する

#1. はじめに

　機械学習と金融工学の両方に興味があったので，一石二鳥を狙って[Machine learning for trading @ Udacity](https://www.udacity.com/course/machine-learning-for-trading--ud501)を受講しています．株価の評価指標である$\alpha$と$\beta$について学んだので，[pandas](http://pandas.pydata.org/)+[matplotlib](http://matplotlib.org/)で時変化をアニメ化して遊んでみました．対象は，2017年2月時点の[世界時価総額ランキング](http://www.180.co.jp/world_etf_adr/adr/ranking.htm)トップ10銘柄です．Githubは[こちら](https://github.com/haltaro/animating-alpha-and-beta)です．
　ご覧の通りド素人ですので，誤りや不備があればご指摘いただけると幸いです．

![ab_optimized.gif](https://qiita-image-store.s3.amazonaws.com/0/151195/97596e37-c42c-0356-2e24-f67d271cc94c.gif)



#2. 環境

##2-1. 株価データの入手

　日足の株価データを使いますので，事前に入手する必要があります．例えば，[YAHOO! FINANCE](https://finance.yahoo.com/)からは，次のような手順で入手できます．なお，**スクレイピングは規約違反ですのでご注意ください**．

1. 目的の銘柄のサマリーページにアクセス（例：[google（``GOOG.csv``）](https://finance.yahoo.com/quote/GOOG?p=GOOG)）
2. ページ上部の``Histrical Data``をクリック．
3. ``Time Period:``が所望の期間であること，また``Frequency:``が``Daily``であることを確認し，``Download Data``をクリック．
4. 必要に応じて，ファイル名を変更（例：``table.csv``→``GOOG.csv``）．

なお，本記事では，``GOOG.csv``等の株価データが``data/``ディレクトリに格納されていることを前提に話を進めます．

##2-2. python

　データ成形に[pandas](http://pandas.pydata.org/)，アニメ化に[matplotlib](http://matplotlib.org/)を使いますので，事前にインストールする必要があります．以下の環境で動作を確認しました．

* Ubuntu 14.04.5 LTS
* python 2.7.11
* NumPy 1.10.4
* matplotlib 1.5.1
* pandas 0.18.0


#3. 評価指標

## 3-1. daily-return

　daily-return $r$とは，簡単に言うと株価の前日比です．厳密には下式で定義されます．ここで，$x[t]$は日付$t$における調整後終値を表します．

```math
r = \frac{x[t]}{x[t-1]} - 1 
```
pandasを使うと，簡単にdaily-returnを計算することができます．詳細は[Machine learning for trading @ Udacity](https://www.udacity.com/course/machine-learning-for-trading--ud501)をご参照ください．ご参考までに，2006年12月1日から2016年12月1日の[`GOOG`](https://finance.yahoo.com/quote/GOOG/?p=GOOG)のdaily-returnを以下に示します．

![goog.png](https://qiita-image-store.s3.amazonaws.com/0/151195/f4f2b7c4-4811-e3f1-ff07-ca1ad4836a8a.png)


## 3-2. アルファ値とベータ値

　市場平均とdaily-returnを比較することで，その銘柄の特徴を評価することができます．ここでは，[`SPY`](https://finance.yahoo.com/quote/SPY?p=SPY)を市場平均と仮定し，[`SPY`](https://finance.yahoo.com/quote/SPY?p=SPY)と[`GOOG`](https://finance.yahoo.com/quote/GOOG/?p=GOOG)のdaily-returnの散布図を書き，回帰直線を引いてみます．期間は，2006年12月1日から2016年12月1日です．

![figure_3.png](https://qiita-image-store.s3.amazonaws.com/0/151195/b6edcc2e-b75f-38d8-0f8d-9c5752ce307a.png)

　この回帰直線$y=\beta x + \alpha$の切片を**アルファ値**，傾きを**ベータ値**と呼びます．銘柄が市場平均と一致する時，回帰直線は直線$y=x$と一致する，つまり$\alpha=0$かつ$\beta=1$となるはずです．$\alpha$が大きいほど，市場平均に対する超過リターンが大きく，$\beta$が大きいほど，市場平均に対する連動性が大きいことを示します．ちなみに上図の場合は，$\alpha=0.000309329688614 $，$\beta=0.917720842929$です．
　$\alpha$および$\beta$は，個別銘柄のアクティブリターンとパッシブリターンの評価に利用される指標だそうです．

#4. 評価

　matplotlibを使って，アルファ値とベータ値の推移をアニメ化します．対象は，以下の2017年2月時点の[世界時価総額ランキング](http://www.180.co.jp/world_etf_adr/adr/ranking.htm)トップ10銘柄です．業種はWikipediaを参考にしました．
　知識が足りないので，特に考察はありません．結果を見ていろいろ妄想するだけです．

|#|symbol|企業名|業種|
|:--|:--|:--|:--|
|1|`AAPL`|[Apple](http://www.apple.com/)|電気機器|
|2|`GOOG`|[Alphabet](https://abc.xyz/)|コングロマリット|
|3|`MSFT`|[Microsoft](https://www.microsoft.com)|情報・通信業|
|4|`BRK-A`|[Berkshire Hathaway](http://www.berkshirehathaway.com/)|保険業|
|5|`AMZN`|[Amazon.com](https://www.amazon.com/)|小売業|
|6|`FB`|[Facebook](https://www.facebook.com/)|インターネット|
|7|`XON`|[Exxon Mobil](http://corporate.exxonmobil.com/)|石油・石炭製品|
|8|`JNJ`|[Johnson & Johnson](https://www.jnj.com/)|サービス業|
|9|`JPM`|[JPMorgan Chase](https://www.jpmorganchase.com/)|その他金融業|
|10|`WFC`|[Wells Fargo](https://www.wellsfargo.com/)|金融業|

##4-1. 前処理

　csvファイルから株価情報を読みだして，daily-returnを計算します．本記記事では，[Udacity](https://www.udacity.com/course/machine-learning-for-trading--ud501)で紹介された以下の関数[^license]を使います．

* `get_data(symbol)`：`symbol`のCSVファイルから調整後終値を読みだして，`pandas.DataFrame`形式で出力する関数です．
* `fill_missing_values(stock_data)`：`stock_data`（`pandas.DataFrame`）の欠損を穴埋めする関数です．
* `compute_daily_returns(stock_data)`：`stock_data`（`pandas.DataFrame`）から，daily-returnを計算し，`pandas.DataFrame`で出力する関数です．

##4-2. アニメ化

　$\alpha$と$\beta$をアニメ化する関数を作りました．デフォルト設定では，2006年12月1日（`start_date`）から2016年12月1日（`end_date`）の上記10銘柄（`symbols`）について，直近2年間（`period`）のデータから値を計算します．10日（`interval`）ずつ計算範囲をずらしながら，`animate_polyfit()`で1フレームずつ描画し，`ab.gif`を出力します．凡例の大きさは`SPY`と当該銘柄の相関係数の大きさを表しています．グラフ右上の日付は，計算範囲の最終日を表します．つまり，右上の日付から過去`period`日分のデータをもとに計算した値を表示しています．
　なお，各銘柄のcsvファイルを`data`ディレクトリに置く必要があります．また，4.1節の関数を定義した`udacity.py`を同じディレクトリに置く必要がありますので，ご注意ください．

```py
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.animation as ani
import matplotlib.cm as cm
import udacity  # Functions defined in "Machine learning for trading"


def animate_a_b(
        symbols=["AAPL", "GOOG", "MSFT", "BRK-A", "AMZN",
                 "FB", "XOM", "JNJ", "JPM", "WFC"],
        start_date="2006-12-01", end_date="2016-12-01",
        period=252 * 2):

    """ --- Preprocess: You have to take Udacity! --- """
    # Read data
    dates = pd.date_range(start_date, end_date)  # date range as index
    stock_data = udacity.get_data(symbols, dates)  # get data for each symbol

    # Fill missing values
    udacity.fill_missing_values(stock_data)
    # Daily returns
    daily_returns = udacity.compute_daily_returns(stock_data)

    """ --- Make animation --- """
    interval = 10  # Nframe interval
    frames = (len(stock_data) - period) / interval  # Num of frames
    markers = ["o", "^", "s"]

    def animate_polyfit(nframe):

        plt.clf()
        daily_returns_p = daily_returns[
            nframe * interval: nframe * interval + period]
        corr = daily_returns_p.corr(method="pearson")

        xmin, xmax = -0.003, 0.003
        ymin, ymax = 0.0, 2.0

        plt.plot([0, 0], [ymin, ymax], '--', color='black')
        plt.plot([xmin, xmax], [1, 1], '--', color='black')
        for n, symbol in enumerate(symbols[1:]):
            beta, alpha = np.polyfit(daily_returns_p["SPY"],
                                     daily_returns_p[symbol], 1)
            plt.plot(alpha, beta, markers[n % len(markers)], alpha=0.7,
                     label=symbol, color=cm.jet(n * 1. / len(symbols)),
                     ms=np.absolute(corr.ix[0, n + 1]) * 25)
        plt.xlim([xmin, xmax])
        plt.ylim([ymin, ymax])
        plt.xlabel("Alpha")
        plt.ylabel("Beta")
        plt.text(xmax, ymax, str(daily_returns_p.index[-1]),
                 ha="right", va="bottom")
        plt.legend(loc="upper left")

    fig = plt.figure(figsize=(8, 8))
    anim = ani.FuncAnimation(fig, animate_polyfit, frames=frames)
    anim.save("ab.gif", writer="imagemagick", fps=18)

```

![ab_optimized.gif](https://qiita-image-store.s3.amazonaws.com/0/151195/574f65ed-211c-832f-cd51-c2435f2fe981.gif)

　2008年ごろのリーマン・ショック周辺で，金融業界の２銘柄（`JPM`，`WFC`）がエキサイトしていますね…．ちなみに，`FB`は比較的若い企業なので，アニメの途中から登場します．

#5. おわりに
　$\alpha$および$\beta$の推移をpandas+matplotlibでアニメ化してみました．使い道はさておき，二次元以上の指標の時変化を直感的に表現する方法として，アニメ化は有りかもしれません[^anime]．今後も受講を続け，記事にできるようなことがあれば，ご紹介したいと思います．
　最後まで読んで下さり，ありがとうございました！

[^license]: 権利上の問題を回避するため，ソースコードは掲載しません．無料講座ですので，ぜひ受講をば！

[^anime]: 実は昔TEDで，これの100倍格好いいアニメを使ったプレゼンを見たのを思い出したので，真似してみました．詳しくは覚えてませんが，GDPと何かを軸にとって，主要国の動きを解説してました．

#Appendix. Trendalyzer
　コメント欄でご指摘の通り，本記事の元ネタはHans Roslingさんの[Trendalyzer](https://www.gapminder.org/tag/trendalyzer/)です．Trendalyzerは，Googleスプレッドシートの[モーショングラフ](https://support.google.com/docs/answer/190718?hl=ja)として利用できますので，使ってみました．[こちら](https://docs.google.com/spreadsheets/d/187NFdwv4-KocUyrpSZRqgu4a2tdDrHtS1BdpHrsPqDc/edit?usp=sharing)でスプレッドシートを公開しておりますので，自由に使ってください．シート1にデータ本体，シート2にモーショングラフを置いてあります．

## A-1. 準備

　googleスプレッドシートで開けるよう，計算結果をエクセル形式で保存しました．モーショングラフ用の[フォーマット](https://support.google.com/docs/answer/190718?hl=ja)[^format]に合わせて，`pandas.DataFrame`を作成します．


```py:save_a_b()
def save_a_b(
        symbols=["AAPL", "GOOG", "MSFT", "BRK-A", "AMZN",
                 "FB", "XOM", "JNJ", "JPM", "WFC"],
        start_date="2006-12-01", end_date="2016-12-01",
        period=252 * 2):

    """ --- Preprocess: You have to take Udacity! --- """
    # Read data
    dates = pd.date_range(start_date, end_date)  # date range as index
    stock_data = udacity.get_data(symbols, dates)  # get data for each symbol

    # Fill missing values
    udacity.fill_missing_values(stock_data)
    # Daily returns
    daily_returns = udacity.compute_daily_returns(stock_data)

    """ --- Calculate spreadsheet of alpha and beta ---"""
    sheet = pd.DataFrame(columns=["Symbol", "Date", "Alpha", "Beta", "Color",
                                  "Size"])
    interval = 10  # Nframe interval
    frames = (len(stock_data) - period) / interval  # Num of frames

    for nframe in range(frames):

        daily_returns_p = daily_returns[
            nframe * interval: nframe * interval + period]
        corr = daily_returns_p.corr(method="pearson")

        for n, symbol in enumerate(symbols[1:]):
            beta, alpha = np.polyfit(daily_returns_p["SPY"],
                                     daily_returns_p[symbol], 1)

            new_row = pd.DataFrame(
                [[symbol, daily_returns_p.index[-1].strftime("%Y/%m/%d"),
                  alpha, beta, n, np.absolute(corr.ix[0, n + 1]) * 25]],
                columns=sheet.columns)
            sheet = sheet.append(new_row, ignore_index=True)

    sheet.to_excel("ab.xlsx")
```

[^format]: `%m/%d/%Y`形式で出力したのですが，google先生がモーショングラフに変換してくださらず，ハマりました．結局，`%Y/%m/%d`形式で出力し，googleスプレッドシートの`表示形式>数字`で`日付`を指定すると，解決しました．よくわかってません…:bow:

## A-2. 遊んでみた

　スタートボタンを押すと，アニメーションを開始します．アニメーションの速度を変更することも可能です．

![ab_resize.gif](https://qiita-image-store.s3.amazonaws.com/0/151195/df1418d2-4660-82b8-e602-a7537e632765.gif)

　タイムバーを使って，早送りや巻き戻しも思いのままです．特定の銘柄を選択すると，軌跡を表示することができます．すごい！

![select_resize.gif](https://qiita-image-store.s3.amazonaws.com/0/151195/94cad1bb-0483-4c97-1472-bd11d7d56030.gif)

#参考

　本記事の作成にあたっては，以下を参考にさせて頂きました．ありがとうございました！

* [Machine learning for trading @ Udacity](https://www.udacity.com/course/machine-learning-for-trading--ud501)：株価データの入手，加工について
* [投資信託のアルファとベータってなに？](https://allabout.co.jp/gm/gc/43752/)：アルファ値とベータ値について
* [Pythonを使った回帰分析の概念の解説 番外編１](http://qiita.com/kenmatsu4/items/c25b2003fb9faf7a783e)：アニメーションの作成について
* [Hans Rosling:
New insights on poverty](https://www.ted.com/talks/hans_rosling_reveals_new_insights_on_poverty)：本記事の元ネタ．