第4章：アルファーファクター研究、Kalman FilterとWavelet変換
今回の記事は時系列データを滑らかにする操作について紹介していきます。
英語を読める方はこちらも是非。
滑らかにする理由は、回帰モデルではよく残差の分布について注目されたりしますが、今回も同じで、統計的なモデルに落とし込むことによって、観測値と期待値の乖離に着目することはよく分析される手法であります。
カルマンフィルターとは
まずはWikipediaを引用させていただきます。
カルマンフィルター (英: Kalman filter) は、誤差のある観測値を用いて、ある動的システムの状態を推定あるいは制御するための、無限インパルス応答フィルターの一種である。
要するにギザギザしている株価を滑らかにするということなんですね。
ウェーブレット変換
Wikipedia引用
ウェーブレット変換（ウェーブレットへんかん、wavelet transformation）は、周波数解析の手法の一つ。基底関数として、ウェーブレット関数を用いる。フーリエ変換によって周波数特性を求める際に失われる時間領域の情報を、この変換においては残すことが可能である。フーリエ変換でも窓関数を用いる窓フーリエ変換で時間領域の情報は残せたが、窓幅を周波数に合わせて固定する必要があるため、広い周波数領域の解析には向かなかった。ウェーブレット変換では、基底関数の拡大縮小を行うので、広い周波数領域の解析が可能である。しかし、不確定性原理によって精度には限界がある。フーリエ変換では、N をデータのサイズとしたときに N logN のオーダーで計算量が増える(O(N logN))が、ウェーブレット変換では O(N) の計算量でできる利点がある。
非常に難しそうに書いてありますが、ウェーブレット関数という既に定義されている関数というものがございまして、それらを組み合わせて株価に近い形にして、周波数を解析するんですね。
インストール
pip install PyWavelets, pykalmanインポートと設定
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
from datetime import datetime
import itertools
import pandas as pd
import pandas_datareader.data as web
from pykalman import KalmanFilter
import pywt
import matplotlib.pyplot as plt
import seaborn as sns
warnings.filterwarnings('ignore')
sns.set_style('whitegrid')
idx = pd.IndexSlice
データ取得
assets.h5というファイルを読み込むのですが、ファイルの作成方法はこちらをご参照ください。
データの中身はQuandlから米国の株式市場2000-2018の調整済み終値のデータです。
DATA_STORE = '../data/assets.h5'
with pd.HDFStore(DATA_STORE) as store:
sp500 = store['sp500/stooq'].loc['2009': '2010', 'close']カルマンフィルター実装
#インスタンス初期化
kf = KalmanFilter(transition_matrices = [1],
observation_matrices = [1],
initial_state_mean = 0,
initial_state_covariance = 1,
observation_covariance=1,
transition_covariance=.01) #フィッティング
state_means, _ = kf.filter(sp500)プロット
sp500_smoothed = sp500.to_frame('close')
sp500_smoothed['Kalman Filter'] = state_means
for months in [1,2,3]:
sp500_smoothed[f'MA ({months}m)'] = sp500.rolling(window=months*21).mean()
ax = sp500_smoothed.plot(title='Kalman Filter vs Moving Average', figsize=(14,6), lw=1, rot=0)
ax.set_xlabel('')
ax.set_ylabel('S&P 500')
plt.tight_layout()
sns.despine();
オレンジ色の線がカルマンフィルターによって滑らかになったものです。MAと比較して短期MA（1か月分）に非常に近い印象があります。
ウェーブレット変換
wavelet = pywt.Wavelet('db6')
phi, psi, x = wavelet.wavefun(level=5)
df = pd.DataFrame({'$\phi$': phi, '$\psi$': psi}, index=x)
df.plot(title='Daubechies', subplots=True, layout=(1, 2), figsize=(14, 4), lw=2, rot=0)
plt.tight_layout()
sns.despine();
こちらがウェーブレット変換のサンプルになります。その他の波も紹介させていただきます。
plot_data = [('db', (4, 3)),
('sym', (4, 3)),
('coif', (3, 2))]
for family, (rows, cols) in plot_data:
fig = plt.figure(figsize=(24, 12))
fig.subplots_adjust(hspace=0.2, wspace=0.2, bottom=.02, left=.06,
right=.97, top=.94)
colors = itertools.cycle('bgrcmyk')
wnames = pywt.wavelist(family)
i = iter(wnames)
for col in range(cols):
for row in range(rows):
try:
wavelet = pywt.Wavelet(next(i))
except StopIteration:
break
phi, psi, x = wavelet.wavefun(level=5)
color = next(colors)
ax = fig.add_subplot(rows, 2 * cols, 1 + 2 * (col + row * cols))
ax.set_title(wavelet.name + " phi")
ax.plot(x, phi, color, lw=1)
ax.set_xlim(min(x), max(x))
ax = fig.add_subplot(rows, 2*cols, 1 + 2*(col + row*cols) + 1)
ax.set_title(wavelet.name + " psi")
ax.plot(x, psi, color, lw=1)
ax.set_xlim(min(x), max(x))
sns.despine()
for family, (rows, cols) in [('bior', (4, 3)), ('rbio', (4, 3))]:
fig = plt.figure(figsize=(24, 12))
fig.subplots_adjust(hspace=0.5, wspace=0.2, bottom=.02, left=.06,
right=.97, top=.94)
colors = itertools.cycle('bgrcmyk')
wnames = pywt.wavelist(family)
i = iter(wnames)
for col in range(cols):
for row in range(rows):
try:
wavelet = pywt.Wavelet(next(i))
except StopIteration:
break
phi, psi, phi_r, psi_r, x = wavelet.wavefun(level=5)
row *= 2
color = next(colors)
ax = fig.add_subplot(2*rows, 2*cols, 1 + 2*(col + row*cols))
ax.set_title(wavelet.name + " phi")
ax.plot(x, phi, color, lw=1)
ax.set_xlim(min(x), max(x))
ax = fig.add_subplot(2*rows, 2*cols, 2*(1 + col + row*cols))
ax.set_title(wavelet.name + " psi")
ax.plot(x, psi, color, lw=1)
ax.set_xlim(min(x), max(x))
row += 1
ax = fig.add_subplot(2*rows, 2*cols, 1 + 2*(col + row*cols))
ax.set_title(wavelet.name + " phi_r")
ax.plot(x, phi_r, color, lw=1)
ax.set_xlim(min(x), max(x))
ax = fig.add_subplot(2*rows, 2*cols, 1 + 2*(col + row*cols) + 1)
ax.set_title(wavelet.name + " psi_r")
ax.plot(x, psi_r, color, lw=1)
ax.set_xlim(min(x), max(x))
sns.despine()
plt.show()
これらの波形を使うために以下のコードを実行します。
pywt.families(short=False)
#['Haar', 'Daubechies', 'Symlets', 'Coiflets', 'Biorthogonal', 'Reverse biorthogonal', 'Discrete Meyer (FIR Approximation)', 'Gaussian', 'Mexican hat wavelet', 'Morlet wavelet', 'Complex Gaussian wavelets', 'Shannon wavelets', 'Frequency B-Spline wavelets', 'Complex Morlet wavelets']最後にS&P500にフィッティングさせます。
signal = (pd.read_hdf(DATA_STORE, 'sp500/stooq')
.loc['2008': '2009']
.close.pct_change()
.dropna())
fig, axes = plt.subplots(ncols=2, figsize=(14, 5))
wavelet = "db6"
for i, scale in enumerate([.1, .5]):
coefficients = pywt.wavedec(signal, wavelet, mode='per')
coefficients[1:] = [pywt.threshold(i, value=scale*signal.max(), mode='soft') for i in coefficients[1:]]
reconstructed_signal = pywt.waverec(coefficients, wavelet, mode='per')
signal.plot(color="b", alpha=0.5, label='original signal', lw=2,
title=f'Threshold Scale: {scale:.1f}', ax=axes[i])
pd.Series(reconstructed_signal, index=signal.index).plot(c='k', label='DWT smoothing}', linewidth=1, ax=axes[i])
axes[i].legend()
fig.tight_layout()
sns.despine();
まとめ
今回は二つの手法について紹介しました。MAは実装は簡単ですが、今後こちらの手法を使う機会があるかもしれないので、参考までに是非。
この記事が気に入ったらサポートをしてみませんか？