日経ソフトウェア２018年9月号を読んでいたら恐怖指数がうんたらかんたらみたいな特集があって、サンプル通りに実行したら途中で、エラーになる謎コードがあり、前後関係も完全にイミフな内容だったので、書き直してみました
<開発環境>
- Python3.6
- JupyterLab
- Windows10
# ライブラリの読み込み
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pandas_datareader.data as web
import pandas as pd
# fredからVixを取得する
vix = web.DataReader("VIXCLS", "fred", "1949/1/1")
# fredからSP500を取得する
sp500 = web.DataReader("SP500", "fred", "1949/1/1")
FREDというのはアメリカにあるセントルイス中央銀行のデータベースのことです。ここには日経平均など主要な金融データが保存されています。
関連記事：【Python】pandasで日経平均の株価データをスクレイピングする
pandas-datareaderは仕様上時々死ぬことがあるので、コード実行に失敗した人は、こことここからcsvを落として以下のコードを実行して読み込んでください。(ディレクトリは任意で調整してください)
# pandasでcsvを読み込む
vix = pd.read_csv('vix.csv', index_col='DATE')
sp500 = pd.read_csv('sp500.csv', index_col='DATE')
ちなみに引数のindex_col='DATE'がないとインデックスがずれて後々エラーを吐くので、忘れないように注意してください。
# 対数変化率を計算する vol = (np.log(sp500).diff().dropna().rolling(20).std() * np.sqrt(365) * 100) # データを結合する ts = pd.concat([vix, vol], axis=1) # 欠損値を削除する ts = ts.dropna(how='any') # 列名を変更する ts.columns = [ 'vix', 'vol'] ts
<実行結果>
エルボー法で最適なクラスタ数を求める
K-Means法で分類する前に、最適なクラスタ数つまり分類数のあたりを付けていきます。
inertia =[]
x = range(1, 50)
# 1～50までで最適なクラスタ数を繰り返しで計算する
for i in x:
model = KMeans(n_clusters=i)
model.fit(tsvix)
inertia.append(model.inertia_)
plt.plot(x, inertia, marker='o')
<実行結果>
何となくだが、分類するクラスタ数は5くらいでいいと思う(小並感)。
k-Mean法でVIXをクラスタリングする
from sklearn.cluster import KMeans
n_clusters = 5
model = KMeans(n_clusters=n_clusters)
tsvix = ts.iloc[:, :1]
ts['pred'] = model.fit_predict(tsvix)
# クラスタリングしたデータを0~4まで一種類ずつ取り出して上書きプロットしていく
for i in range(n_clusters):
# i(0~4)に当てはまるデータをデータフレームから抽出する
ts2 = ts.loc[ts['pred'] == i]
# 散布図をプロットする
plt.scatter(ts2['vix'], ts2['vol']/100)
# ラベルを付ける
plt.xlabel("VIX")
plt.ylabel("Logarithmic rate of return")
<実行結果>
当たり前ではあるが、VIXが大きくなっていくと対数収益率(ボラティリティ)、つまり金融工学におけるリスクもそれにともなって大きくなっていくのが分かる
関連記事:【Python】pandasでビットコイン価格のローソク足と移動平均線を計算してプロットする
関連記事：【Python】機械学習で株価を予測する~Scikit-Learnの決定木アルゴリズムを使う
