前編の続きで時系列データを用いて分析し、モデルにフィッティングしました。
PythonのStatsModelsを用いて業種別株価指数の分析をしました。
今回は業種別株価指数(サービス業)を例に分析を行います。
以下コード作成の流れです。
///////////////////////////////////////////////////////////////////////////////////////////////
① 時系列データの取得
② 自己相関係数と偏自己相関係数の調査
③ 時系列データをトレンド、季節性、残差に分解
④ ADF検定で調査
⑤ 時系列データの相関
⑥ モデルにフィッティング
///////////////////////////////////////////////////////////////////////////////////////////////
① 時系列データの取得
前回と同様の方法で日経500種平均株価と業種別株価指数(サービス業)の時系列データを取得しました。
2018年度から2021年7月までの株価指数になります。
import pandas as pd
from io import StringIO
import urllib.request
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from datetime import datetime
import warnings
import itertools
%matplotlib inline
#日経500種平均株価と業種別株価指数(サービス業)を取得
url = "https://indexes.nikkei.co.jp/nkave/historical/nikkei_500_stock_average_daily_jp.csv"
def read_csv(url):
res = urllib.request.urlopen(url)
res = res.read().decode('shift_jis')
df = pd.read_csv(StringIO(res))
df = df.drop(df.shape[0]-1)
return df
df = read_csv(url)
df["データ日付"] = pd.to_datetime(df["データ日付"], format='%Y/%m/%d')
df = df.set_index('データ日付')
df2 = df.iloc[:,[0, -1]]
df2 = df2.sort_index(ascending=True)
#名前の変更
df2.rename(columns={'終値':'Nikkei500', '業種別（サービス）終値':'Services'}, inplace=True)
② 自己相関係数と偏自己相関係数の調査
次にPythonを用いて自己相関係数と偏自己相関係数を求めました。
自己相関係数は、過去のデータとどれほど相関あるのかを示す係数となります。
例えばt日間差の自己相関係数を求めたときに
t日前 → t-1日前 → t-2日前 … 1日前 → 今日と時間をずらして相関を求めていきます。
偏自己相関係数は、例えばt日前と今日の相関関係のように、
t-1日前から1日前までの間の影響を取り除いて相関を求めたものになります。
業種別株価指数(サービス業)データの自己相関係数と偏自己相関係数を可視化しました。
そのためのコードは以下となります。
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot()
#自己相関係数
sm.graphics.tsa.plot_acf(df2['Services'], lags=100,ax = ax)
plt.show()
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot()
#偏自己相関係数
fig = sm.graphics.tsa.plot_pacf(df2['Services'], lags=100, ax=ax)
plt.show()
可視化したグラフはコレログラムと呼ばれます。
自己相関(Autocorrelation Function)のコレログラムを以下に示します。
横軸はラグで、ラグが1離れるときは、1つ前のデータが該当データにどのくらい影響するかを示しています。
青で塗られる領域は95%信頼区間になります。
領域内にプロットされると相関があるとは言い切れない一方で、領域外だと相関があると判断できます。
このグラフを見ると、ラグが離れると自己相関係数は少しずつ下がる傾向がありますが、長期にわたり相関関係があることを示唆しています。
偏自己相関(Partial Autocorrelation Function)のグラフでは、ラグ1で大きな相関がみられます。
ラグが2以上離れると相関が弱いことを示唆しています。
③ 時系列データをトレンド、季節性、残差に分解
時系列データでは以下のコードを用いて、トレンド、季節性、残差に分解できます。
業種別株価指数(サービス業)の時系列データを分解してみました。
今回扱う指数は右肩上がりのトレンドがあることが分かります。
res = sm.tsa.seasonal_decompose(df2['Services'],freq=51)
fig = res.plot()
plt.xlabel("Date")
plt.show()
④ ADF検定で調査
時系列データの解析を行うときに「定常性」という概念は重要になります。
定常性とは「時間の経過によらずに、ある一定の値を中心に同程度の幅で振れて変化する」ことを意味します。
多くの時系列データはトレンドや季節変動といった要因が加わり非定常になっているため、
複雑で解析がしにくくなります。
扱う時系列データには定常性があるかを検定する方法として、
単位根検定(代表的なものとしてADF検定)があります。
今回扱う時系列データは定常性を有するかどうかADF検定にかけて調べます。
adf_result = sm.tsa.stattools.adfuller(df2['Services'],autolag='AIC')
adf2 = pd.Series(adf_result[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
print(adf)
結果は
Test Statistic -0.508646
p-value 0.890318
#Lags Used 0.000000
Number of Observations Used 851.000000
dtype: float64
p値が大きいことから、今回の時系列データは定常性がないと言えます。
③でトレンドがあることが確認できたので、定常性がないと言えばそうでしょう。
時系列データが定常性を持たないとき、データに定常性を持たせるように加工していくことでデータが解析しやすくなります。
具体的な方法の1つとして、差分をとり、トレンド・季節変動を除去させる方法があります。
1日前と当日の差分を以下のコードでとり、ADF検定で確認しました。
df2_diff = df2.diff().dropna()
adf2_result = sm.tsa.stattools.adfuller(df2_diff['Services'],autolag='AIC')
adf2 = pd.Series(adf2_result[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
print(adf2)
値は省略しますが、p値は0に近い値になり、定常過程であることが確認できます。
元々の時系列データ(原系列と呼びます)が非定常で、差分を取ったときに定常過程である時は単位根過程と呼ばれます。
今回の時系列データは単位根過程であると考えられます。
単位根過程では以下に述べるARIMAモデルや季節変動を考慮したSARIMAモデルでの解析ができます。
⑤ 時系列データの相関
③で時系列データを用いてトレンド、季節性、残差に分解しました。
日経500種平均株価と業種別株価指数(サービス業)の残差から互いの相関関係を分析しました。
日経500種平均株価のデータを1として、
対応するデータ(今回はサービス業の株価指数)の相関係数は-1から1の範囲で値を持ちます。
1は正の相関があり、1は負の相関(逆向きな関係性)があると言われます。
また、0に近いと相関が無いといわれます。
日経500種平均株価と業種別株価指数(サービス業)の互いの相関係数を出すために、
corr()を呼び出します。
時系列データの残差を出して互いの相関係数を求めました。
df2_diff_corr = df2_diff.corr()
print(df2_diff_corr)
日経500を1として0.88という値が得られました。
日経500の値と高い相関があることを示唆しています。
⑥ モデルにフィッティング
ARIMAモデルを用いて業種別株価指数(サービス業)の時系列データを分析しました。
ARIMAモデルは、"自己回帰和分移動平均モデル"と呼ばれるものです。
自己回帰モデル(ARモデル)と移動平均モデル(MAモデル)を合成したモデル(ARMAモデル)に
和分過程を導入したモデルになります。
ARIMAモデルは非定常過程にも適応可能ですが、差分値が定常過程である必要があります。
④で差分が定常過程であることを確認しました。
直近約1年間の株価指数を用い、1か月前までの値をフィッティングさせ、
直近1か月の値を予測値に用いました。
from statsmodels.tsa.arima_model import ARIMA
import numpy as np
N=853
stock_train = df2['Services'][600:N]
ARIMA_stock_price = ARIMA(stock_train,order=(1, 1, 1)).fit()
pred = ARIMA_stock_price.predict(1, 230, typ="levels")
pred2 = ARIMA_stock_price.forecast(steps=23)
plt.plot(df2['Services'].values[600:N], label = 'Observations')
plt.plot(pred.values, linestyle='--', color='y', label="Applying to training data")
plt.plot(range(230, 253), pred2[0], linestyle='--', color='k', label="Predicted value")
plt.legend(fontsize = 10)
plt.xlabel('Days', fontsize = 12)
plt.ylabel('Price', fontsize = 12)
青が学習データ、黄緑が学習させた結果、そして黒い点線が予測値です。
ARIMAモデルの引数になるorderのp,qの値は、
statesmodelのsm.tsa.arma_order_select_ic()のパラメーター推定関数を参考に決定しました。
予測値と観測値の値の差分を出し、値が解離していないか求めました。
コードは以下になります。
df3_diff = pred2[0] - df2['Services'][830:N]
print(df3_diff)
データ日付
2021-06-03 20.731258
2021-06-04 34.520736
2021-06-07 10.841220
2021-06-08 4.674962
2021-06-09 11.623379
2021-06-10 15.567811
2021-06-11 1.396759
2021-06-14 9.561169
2021-06-15 1.054556
2021-06-16 25.679921
2021-06-17 54.131306
2021-06-18 45.412091
2021-06-21 88.314864
2021-06-22 17.305176
2021-06-23 26.381403
2021-06-24 33.700101
2021-06-25 36.372292
2021-06-28 38.533688
2021-06-29 53.572140
2021-06-30 43.620196
2021-07-01 74.537581
2021-07-02 47.222283
2021-07-05 40.876516
グラフを見るに傾向は捉えていると思われますが、
予測開始日から現在に近づくにつれて、予測値と観測値の値が解離しており、
まだまだ工夫が必要なようです。
今回は手軽に行えるARIMAモデルを用いましたが、
深層学習の1つのモデルであるLSTMなどを用いた他モデルの実装にも取り組んでみたいです。
以上業種別株価指数を試しにpythonを用いて分析しました。
ここでは1業種に限定しましたが、他業種と比較/分析することで、
指数の値動きの関連性や、
株価に影響を与えるイベントが起きたときに与える業種の値動きの違い、
がより具体的に見えるかもしれません。
今回は業種別日経平均株価の時系列データを用いて分析を行いました。
他の業種別株価指数である東証業種別株価指数や米国市場の業種別株価など、
興味があるものはまだまだあるので今後も解析していきたいと思います。