目次
・はじめに
・本記事の概要
・プログラム
・おわりに
はじめに
このブログは、通信講座の受講終了要件を満たすためのものです。
【自己紹介】
３７歳女。
営業事務や士業のアシスタントとして仕事をしてきましたが、
海外生活の経験からずっと、専門職・手に職というものに憧れていました。
子供が保育園に入るのを機に、最後の学習チャンスだと思い、
以前から興味があったPythonに挑戦しました。
本記事の概要
自分の学習の理解度の確認を主に、
Pythonの学習をしている初心者にも、わかりやすいことを心がけています。
SARIMAモデルを用いて日経平均の株価予測をしていきます。
実行環境
Google Colaboratory (Google Colab)
Python 3.10.12
プログラム作成開始！！
SARIMAモデルとは
とても簡単に言うと、季節性の影響を組み込んだ過去のデータに基づいて、未来の値を予測するコードのことです。
プログラムの流れ概要
①データを準備する
②データを可視化する
③データの前処理をする
④自己相関と偏自己相関関係を可視化し、データの系統を調べる
⑤データを更に分解する
⑥SARIMAモデルのパラメータを設定する
⑦statsmodelsライブラリを使ってSARIMAモデルをフィッティングする
⑧予測と評価
データの準備
まず、日経平均の株価データを取得します。
Yahooファイナンスから、1990年1月1日から2023年12月31日までの、
日経平均（ティッカーシンボル「^N225」）のデータを取得しました。
同時に、どういうデータなのかをまずは確認しました。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import yfinance as yf
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.statespace.sarimax import SARIMAX
import itertools
# 日経平均のティッカーシンボル
ticker = '^N225'
# データのダウンロード
data = yf.download(ticker, start='1990-01-01', end='2023-12-31', interval='1d')
#データの情報を出力
print(data.info())
出力結果
[*********************100%%**********************] 1 of 1 completed
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 8347 entries, 1990-01-04 to 2023-12-29
Data columns (total 6 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 Open 8347 non-null float64
1 High 8347 non-null float64
2 Low 8347 non-null float64
3 Close 8347 non-null float64
4 Adj Close 8347 non-null float64
5 Volume 8347 non-null int64
dtypes: float64(5), int64(1)
memory usage: 456.5 KB
None
出力結果自分用メモ
1. pandasのDataFrameオブジェクトである
2. 日付順に並んでいる、8347行ある、1990年1月4日から2023年12月29日までである
3. ６列ある。（Open, High, Low, Close, Adj Close, Volume）
4. 欠損値はない。
5. 浮動小数点数（float64）と整数（int64）がある。
データの前処理
取得したデータを前処理します。
終値の月平均を使用します。
【月平均にする理由】
データのノイズを減少させる効果があります。
日次の場合、短期的な変動が多く、予測に適していない場合があるためです。
また平均を取ることで、データが滑らかになるため、季節性やトレンドをより明確に捉えることができます。
# Close(終値)のみを抽出
close_data = data['Close']
# 月ごとの平均を計算
monthly_avg = close_data.resample('M').mean()
# データの情報を出力
print(monthly_avg.info())
出力結果
<class 'pandas.core.series.Series'>
DatetimeIndex: 408 entries, 1990-01-31 to 2023-12-31
Freq: M
Series name: Close
Non-Null Count Dtype
-------------- -----
408 non-null float64
dtypes: float64(1)
memory usage: 6.4 KB
None
出力結果自分用メモ
1.PandasのSeries型である
2.1990年1月31日から2023年12月31日まで408ヶ月分のデータである
3.月次（月ごと）の終値のデータである
4.Non-Null Count 欠損値無し
5.float64のデータである
プロットでデータの可視化
# プロットでデータの可視化
plt.figure(figsize=(12, 6))
plt.plot(monthly_avg.index, monthly_avg, label='Monthly Average of Closing Prices', color='blue')
plt.title('Monthly Average of Nikkei 225 Closing Prices')
plt.xlabel('Date')
plt.ylabel('Price')
plt.grid(True)
plt.legend()
plt.show()
出力結果
出力結果自分用メモ
1990年から2004年頃までは下降トレンド、2010年頃以降は概ね上昇トレンド。
このグラフからは、季節性はいまいちわからない。
なので、自己相関や偏自己相関で調べてみることにします。
自己相関と偏自己相関のプロットを作ってみる
# 自己相関プロット
fig, ax = plt.subplots(figsize=(12, 6))
plot_acf(monthly_avg, lags=40, ax=ax)
ax.set_title('Autocorrelation Function (ACF)')
ax.set_xlabel('Lag')
ax.set_ylabel('ACF')
plt.grid(True)
plt.show()
# 偏自己相関プロット
fig, ax = plt.subplots(figsize=(12, 6))
plot_pacf(monthly_avg, lags=40, ax=ax)
ax.set_title('Partial Autocorrelation Function (PACF)')
ax.set_xlabel('Lag')
ax.set_ylabel('PACF')
plt.grid(True)
plt.show()
自己相関プロット出力結果
自己相関係数（自己相関係数のプロット）
自己相関係数は、ある時点のデータが、別の時点のデータとどれくらい似ているかを示します。
上記の出力結果では、ACF（自己相関係数）が徐々に時間とともに下がっています。
これは、過去の値が現在の値に対する影響が時間とともに弱まっていることを示しています。
また、データにトレンドがあることを示唆している可能性があります。
日経平均の終値月平均のグラフで、下降トレンドと上昇トレンドが見られるので、
トレンドがある可能性というのは正しそうです。
周期性もなさそうですが、季節性についてはこの時点ではわかりません。
偏自己相関プロット出力結果
偏自己相関係数（偏自己相関係数のプロット）
あるラグの過去の値が現在の値にどれくらい影響を与えているかを示すとき、他のラグの影響を取り除いた後の影響を示します。つまり、「特定のラグが持つ純粋な影響」を測るものです。
上記のプロット結果では、最初の１ラグだけy軸の１に近く、それ以外は網掛けの中で０に近いプラスかマイナスの値を取っています。
これは、短期的なラグで自己相関が強く、長期になると自己相関が弱いことを示しています。
ラグ２以降すべて網掛けの中に入っているので、偶然の変動の範囲内にあることを示しています。
日経平均は、短期的には多少予測は可能だが、長期の予測は偶然でしか当たらないと示していると思われます。
もう少し詳しく分解してみる
自己相関プロットと、偏自己相関プロットで、トレンドがありそうなことと、周期性はなさそうなことがわかりましたが、季節性については明確ではありませんでした。
そこで、「sm.tsa.seasonal_decompose」という関数を使って、データを分解し、季節成分、トレンド成分、残差成分を個別に見てみます。
decomposition = sm.tsa.seasonal_decompose(monthly_avg, period=12)
decomposition.plot()
plt.show()
出力結果
上記結果を見ると、下降上昇のトレンドがあることと、季節性があることが明確になりました。
よって、SARIMAモデルを使用するのは問題なさそうです。
定常性について
予測精度の向上や、モデルの適用性を上げるため、
時系列のあるデータを分析する時、データに定常性を持たせる必要があります。
すなわち、データを加工するなどして、トレンドや季節性をデータから除去する必要があります。
しかしSARIMAモデルを使用する場合、定常性を持たせる必要がありません。
パラメータを設定することで、データが持っているトレンドや季節性を自動的に考慮してくれるので、定常性を持たせるのと同じ働きをしてくれます。
パラメータについては、以下で説明していきます。
SARIMAモデルのパラメータ設定
SARIMAモデルとは「未来予測をしてくれるコード」と先述しましたが、ただコードを書くだけではなく、自分で事前設定が必要です。
それは、p、d、q、P、D、Q、Sという７つの値（パラメータ）を決めるというものです。
この値を決めることで、データのトレンドや季節性を考慮できます。
各アルファベットには、下記の通りそれぞれ意味があります。
p：AR（自己回帰）成分の次数（通常０〜１０の間）
d：差分化の次数（通常０〜２の間）
q：MA（移動平均）成分の次数（通常０〜１０の間）
P：季節性AR成分の次数（通常０〜１０の間）
D：季節性差分化の次数（通常０〜１の間）
Q：季節性MA成分の次数（通常０〜１０の間）
s：季節性の周期（一般的に月次は１２、４半期は９０、週次は５２）
これらの数字を何にするか決めなければならないのですが、残念ながら自動でPythonが決めてくれないので、自分で設定します。
一つずつ組み合わせて検証するには、通常の数値を考慮しただけでも、２０万通り以上あるので、不可能に近いです。先に作成した自己相関プロットや偏自己相関プロットから予測を立てても良いですが、それも時間と労力を必要とします。
なのでここでは、学習の過程で教えてもらった「どの値が最も適切なのかを調べるコード」を使って、パラメータを設定していこうと思います。
下記がそのコードです。（簡易化のため、各パラメータの幅を０から１にしています。）
#SARIMAモデルパラメータの最適解を求める
def selectparameter(DATA, s):
# パラメータの範囲を拡張
p = d = q = range(0, 2) # 0から1まで
P = D = Q = range(0, 2) # 0から1まで
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], s) for x in list(itertools.product(P, D, Q))]
parameters = []
BICs = np.array([])
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
mod = SARIMAX(DATA,
order=param,
seasonal_order=param_seasonal)
results = mod.fit(disp=False) # disp=False でフィッティングの詳細表示を抑制
parameters.append([param, param_seasonal, results.bic])
BICs = np.append(BICs, results.bic)
except:
continue
best_params = parameters[np.argmin(BICs)]
print(f"Best parameters: {best_params}")
return best_params
一番適切なパラメータが決定！
Best parameters: [(0, 1, 1), (1, 1, 1, 12), 6641.980801367711]
モデルのフィッティングと実行
# SARIMAモデルのパラメータ設定
p, d, q = 0, 1, 1
P, D, Q, s = 1, 1, 1, 12
# SARIMAモデルのフィッティング
model = SARIMAX(monthly_avg, order=(p, d, q), seasonal_order=(P, D, Q, s))
results = model.fit()
# モデルの要約を表示
print(results.summary())
サマリー出力結果
上記サマリーを読み解いてみます。
２段落目の「P>|z|」の値で、係数の優位性を判断できます。
ma.L1、ar.S.L12、ma.S.L12のこれらがMA項（移動平均項）やAR項（自己回帰項）の係数です。
統計モデルにおいて、p値が小さい場合、その係数は統計的に有意であり、モデルにおいて重要な影響を持つとされます。これらのｐ値が全て「0.000」であることから、係数の有意性は高いことが分かります。
したがって、パラメータ設定はうまく機能していると言えそうです。
Jarque-Beraテストのｐ値（Prob(JB))が「0.00」は残差が正規分布していない可能性が高いことを示しています。
Heteroskedasticityのｐ値（Prob(H))が「0.00」は残差の分散が一定でない可能性があることを示しています。
これらは、残念ながらモデルの予測の信頼性や精度が低い可能性を示しています。
予測開始！
将来の予測が正しいか判断するために、予測と実績の範囲を少しばかり重ねることにしました。
偏自己相関プロットやサマリー結果から、予測が難しかったり、信頼性にかけると先述しましたが、の2020年から2023年の実績と、予測が近ければ、未来の予測も多少は信用できると思いたいからです。
#2023年1月〜2025年12月の終値を予測し、predに格納
pred = results.predict('2020-1', '2025-12')
#predと元の時系列データを可視化
#予測データは赤色で表示
fig=plt.figure(figsize=(10,6))
plt.title("Prediction of Close Price (monthly)")
plt.plot(monthly_avg)
plt.plot(pred, "r")
plt.show()
出力結果
上記の通り、赤色の予測が実績にとても近しい結果となりました！！
未来はわかりませんが、期待も込めて2025年は上昇トレンドとなることでしょう。
おわりに
講座学習で理解が及ばなかったことの復習になり、とても身になりました。
今回は、SARIMAモデルを使用しましたが、サマリーの結果、「モデルの予測の信頼性や精度が低い可能性」と出てしまったこともあり、遠い未来の予測は、これだけでは当然難しそうだなと感じました。
そこをもっと突き詰めてみたり、他のモデルや感情分析、更に講座では学習しなかったものなど、今後も学びを深めていきたいと思います。
ここまで読んでいただき、ありがとうございました。