ころがる狸

ころがる狸のデータ解析ブログ

【python+株価+statsmodel】ARIMAモデルでダウ平均株価を解析してみた

こんにちは。ゴールデンウィーク始まりましたね。とりあえずこの5日間は、機械学習・統計解析の勉強をしつつどうぶつの森で島を開拓する予定でいます!
・・・・・・・・・
今日はpythonで時系列データの解析を行いたいと思います。時系列データというと株価や気温など多種多様なものがあるわけですが、ここでは経済指標の一つであるダウ平均株価を例にとって解析を進めていきたいと思います。csvファイルはyahoo!financeの以下のURLから取得しました。日、週、月毎のデータを取得できますが、1985年から2020年現在に至るまでの月ごとの株価情報を解析していきます。また、データ解析には主にstatsmodelというpythonのライブラリを用いました。これは様々な統計モデルを簡単に扱うことができる非常に便利なツールです。
finance.yahoo.com
www.statsmodels.org

また、本稿は2019年9月に出版された以下の参考書をもとに執筆させて頂きました。pythonでの時系列情報や状態空間モデルの取り扱いを実際のコードを参考にしながら学ぶことができます。コード付きの本だと自分の知らない書き方が必ず見つかるので、プログラミング技術向上という意味でも参考になると思います。

解析は以下の手順で行いました。まず、基本的なデータの可視化と平滑化。偏自己相関関数のプロットによる各時点のデータの相関の把握。続いて差分を可視化し時系列が定常過程か非定常過程かを推測し、単位根検定で定量的に評価。これらの結果を根拠にARIMAモデルによる予測です!

移動平均法による平滑化

まずは、pandasを使ってcsvを読み込み簡単なデータの分析を行いたいと思います。株の運用をしている方は見たことがあるかもしれませんが、株価のような変動の激しい時系列データでは移動平均法を取って平滑化しノイズを除去することがよく行われます。移動平均法とは、ある時点から過去の一定期間内のデータの平均値を取る操作のことを指します。ここでは12か月移動平均をとりましたが、確かに滑らかになっている様子が見て取れます。また一定期間の平均を取る必要があるので、移動平均線は実測のデータより12か月分未来側にシフトしています。そのため大暴騰や大暴落があったとき、それが移動平均線に反映されるのは遅れがでることになります。なおpandasでの移動平均の計算は以下のコードにもあるように非常に簡単です。

f:id:Dajiro:20200502190659p:plain
観測データと12か月移動平均線

import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
import statsmodels.api as sm
from datetime import datetime

df = pd.read_csv('data/ダウ平均株価(1985-2020).csv')
d = list(map(lambda x:datetime.strptime(x, '%Y/%m/%d'), df.Date))

plt.plot(d, df['Close'], label = '観測値')
plt.xlabel('年', fontsize = 16)
plt.ylabel('終値(USD)', fontsize = 16)
#ここで移動平均を計算すると同時にプロット
plt.plot(d, df.Close.rolling(48).mean(), label = '12か月移動平均')
plt.legend(fontsize = 12)
plt.savefig('fig1.png', bbox_inches='tight')

自己相関・ 偏自己相関の可視化

異なる時点に時系列データの相関関係を知るために、自己相関係数を計算して可視化することは非常に重要です。自分自身との自己相関係数は1となり、相関の強いデータほど自己相関係数は1に近づきます。1985年から株価は変動しつつもジワジワと上昇しており、前月より今月の方が株価は高いという相関関係があるように見えませんか?実際に1985年のデータ開始点から見たときの自己相関関数をプロットすると高い相関があることが分かります。この図のことをコレログラム(correlogram)と呼びます。青で塗られた領域はデータが無相関としたときの95%信頼区間で、この領域の外にある場合この仮定は棄却され、相関があると判断できます。

f:id:Dajiro:20200502190745p:plain
自己相関関数。ラグは時間差で、隣り合うデータ間はラグ1となる。
また、今月の値と先月の値が相関するとした場合、今月の値と先々月の値にも相関があると仮定したことになりますよね。この時先月の値の影響を排除し、今月と先々月の値の相関を直接比較できるようにしたものを偏自己相関関数と呼びます。偏自己相関関数に関しては以下の記事が非常に詳しかったので是非お読み下さい。
自己共分散と自己相関 | 時系列データにおけるある時刻の観測値とその過去のある時刻の観測値との相関
この偏自己相関関数をプロットすると、今月の値と翌々月以降の値には直接の相関がないことが分かります。つまり今月の値が分かったところで、株価の変動があるため翌々月の値を的確に予測することは困難である、と言えると思う(のですが、どうでしょうか…)。これらの自己相関関数の計算はstatsmodelを用いると非常に簡単に行うことができます。
f:id:Dajiro:20200502190933p:plain
偏自己相関関数。

close = df.Close
#自己相関係数
acf = sm.tsa.stattools.acf(close)
#偏自己相関係数
p_acf = sm.tsa.stattools.pacf(close,  method='ols')

sm.graphics.tsa.plot_acf(close, lags=30)
plt.xlabel('ラグ', fontsize=14)
plt.ylabel('自己相関係数', fontsize=14)

sm.graphics.tsa.plot_pacf(close, lags=30)
plt.xlabel('ラグ', fontsize=14)
plt.ylabel('偏自己相関係数', fontsize=14)

差分系列と単位根検定

時系列データを扱う上で、定常性という概念が非常に重要になります。(私の理解の範囲で)一言で言うなら、時間経過に対してそう大きく変わらないデータのことを指します。例えばホワイトノイズのように平均ゼロの周辺で一定の分散をもってランダムに振動を続けるデータは定常性があると言いますし、今回のダウ平均株価は明らかな上昇傾向が存在し変化し続けているためこれは定常性があるとは言いません(非定常なデータですね)。
より細かい話をすると、定常性には弱定常性と強定常性という2種類の概念があります。弱定常性は常に平均と分散が一定である状態を指し、まさにホワイトノイズの性質と一致しますね。また、強定常性は時系列データのかたまり(yt, ..., yt+1)がある同時確率分布から生成されたとした場合、それらを任意の時刻だけシフトさせても同時確率分布は変わらないという性質のことを指します。
話が天下り的になってしまうのですが、定常性は解析を進める上で非常に重要な性質になります。そこで、見るからに非定常なこのダウ平均株価の推移を何とかして定常性のある時系列データ(定常過程と呼ぶことにします)に変換してみたいと思います。
ここで行われるのは、隣接するデータ間の差分を取るというかなり素朴な方法です。しかしこれによりデータを定常過程として扱うことができるようになります。差分化は簡単です。pandasで1行で処理できます。素晴らしきpandasですね。また対数を取った対数差分系列も可視化によく使われます。これも1行pandasです。

f:id:Dajiro:20200502191003p:plainf:id:Dajiro:20200502191017p:plain
左:差分系列 右:対数差分系列

close_d = close.diff().dropna()
close_log = (1 + close.pct_change()).apply(np.log)

さて、これが定常過程であることを検定するのに使われるのが、単位根検定になります。もともとのデータが非定常で、差分を取ったものが定常過程に従うときこれを単位根過程であると呼びます。そしてその時系列データが単位根過程に従うか否かを判定するのに、単位根検定が使われます。正確にはここで行うのはADF検定と呼ぶのですが、これをstatsmodelで使ってみましょう。単位根検定の詳細は冒頭で紹介した黄色い本を読んでみてください!
まずは、差分を取る前の株価情報を検定にかけてみましょう。ここでregressionという引数にアルファベットが与えられていますが、これはADF検定に使われる統計モデルの条件を指定しており、トレンド項を含めるか、定数項を含めるかという指定をしています。出力されるのは各検定で得られるp値で、これが大きいとき時系列は単位根過程であると判定できます。この場合pはどのモデルでも比較的大きい値をとることから、株価情報は単位根過程であると考えられます。

from statsmodels.tsa import stattools

#ADF検定
ctt = stattools.adfuller(close, regression='ctt')
ct = stattools.adfuller(close, regression='ct')
c = stattools.adfuller(close, regression='c')
nc = stattools.adfuller(close, regression='nc')

print(ctt[1], ct[1], c[1], nc[1])
#0.4375574529611433 0.3627093947736816 0.930086132887625 0.9826511637975629

なお、ここでは省略しますが差分を取った系列ではp値がほぼゼロとなり、単位根過程ではない、つまり定常過程であることが確認できます。

ARIMAモデル

ここまで長くなりましたが、ARIMAモデルを使ってフィッティングと予測をしてみましょう。これは自己回帰和分移動平均モデル(autoregressive integrated moving average model)の略称で、決して有馬教授という方が提案した手法というわけではないようです笑 これはデータの差分をモデル化する手法で、ARIMAモデルは上で紹介した弱定常性が備わっていることが分かっています。時系列データそのものは非定常でも問題ありませんが、差分値が定常過程である必要があります。そのため単位根検定を行い株価が非定常で、差分値が定常過程であることを確認する必要があったのです。それを確認しなければARIMAモデルを使う根拠がなくなってしまいますからね。また、この手法は機械学習のように膨大な学習データは必要でなく、予測以前の時系列データがあれば訓練と予測ができます。この手軽さが機械学習技術に勝る利点と言えるのではないでしょうか。
statsmodelを使うとモデルによるフィッティングと予測が簡単に行えます(冒頭で紹介した黄色い本ですとpredict関数で予測までしてますが手元のデータでは何故かエラーをはき未来のデータが予測できませんでした。forecast関数を使うと未来の値が予測できます)。全データのうち先頭の350データを持ちいてフィッティングを行い、残りのデータで予測をしてみましょう。上昇トレンドは確かに捉えられていることが分かりますね。また訓練データへのあてはめも良いのですが、近年の株価の上昇率の高さは予測できていません(まぁ、当たり前か)。ちなみに予測値に対する標準偏差と95%信頼区間も同時に確認できるので、予測値の信頼性評価も行えますね。参考までに図を張っておきましょう。

f:id:Dajiro:20200502191159p:plain
ARIMAモデルを用いた予測。傾向は掴めていると思われる。
f:id:Dajiro:20200502191247p:plain
ARIMAモデルを用いた予測。標準偏差を同時にプロットした。

from statsmodels.tsa.arima_model import ARIMA
import numpy as np
close_train = df.Close.dropna()[:350]
model_d1 = ARIMA(close_train, (3, 1, 1))
results_d1 = model_d1.fit()
res_d1 = results_d1.resid

tup = results_d1.forecast(steps=73)
plt.plot(close.values, label = '観測値')
plt.plot(results_d1.predict(1, 349, typ='levels').values, '--', label="訓練データへのあてはめ")
plt.plot(range(350, 423), tup[0], "k--", label = '予測値')
plt.legend(fontsize = 12)

まとめ

以上のように、ARIMAモデルを用いた株価の予測までを行うことができ、上昇トレンドを掴めていることを確認しました。次回は状態空間モデルを用いた予測や、機械学習技術の一つであるLSTMを用いた予測などを行ってみたいと考えています。