Pythonを用いて株式・長期債・コモディティの見通しを予測してみた

今回は、SARIMAモデルを用いた時系列分析を行い、未来予測をしてみようと思います。

★目次

０. 自己紹介・導入
１. データの読み込み
２. データの整理
３. データの可視化
４. データの周期の把握 (パラメーターsの決定)
５. s以外のパラメーターの決定
６. モデルの構築
７. データの予測とその可視化
８. 予測結果・考察・まとめ

０.自己紹介・導入

■自身のレベル

・株式会社Aidemyにて受講可能な「データ分析講座（３ヶ月）」を受講しました。
・データ分析の業務経験は無く、講座内容を参考に今回初めての実践として取り組みました。

■実行環境

・サービス
Google Colaboratory
・言語
Python 3.7.13
・ライブラリ
pandas-datareader 0.9.0
datetime （Python標準）
pandas 1.3.5
sklearn 0.0
matplotlib 3.2.2
numpy 1.21.6
statsmodels 0.10.2
itertools （Python標準）

■時系列分析の流れについて

・「時系列データ」は 時間の経過とともに変化するデータ のことをいいます。具体的には東京の毎時間の気温の記録、ある会社の毎月の売上高の記録、ある国のGDPの記録や日経平均株価の記録などは時系列データといえます。
・「時系列分析」は、会社の売り上げや商品の売り上げ予測、さらには来店者数の予測など、ビジネスにおいて非常に重要な分析技術となります。
・Pythonを用いた時系列分析にはさまざまな方法がありますが、今回はSARIMAモデルを用いたモデル構築・予測値の算出を行っていきます。
・SARIMAモデルを用いた時系列解析の手順は以下の７つで、それぞれを章立てして記載させていただきます。　

1. データの読み込み
2. データの整理
3. データの可視化
4. データの周期の把握 (パラメーターsの決定)
5. s以外のパラメーターの決定
6. モデルの構築
7. データの予測とその可視化
   　

■解析に用いた経済指標

（１）S&P500指数

「S&P500®」とは、S&Pダウ・ジョーンズ・インデックスLLCが公表している、アメリカの代表的な株価指数の1つです。
ニューヨーク証券取引所、NASDAQ等に上場している銘柄から代表的な500銘柄を時価総額で加重平均し、指数化しています。
★この指標を用いることでアメリカの代表的500社の株式（ひいてはアメリカ株式）に投資する場合の見通しを予測可能と考えています。

（２）アメリカ国債利回り（10年）

「アメリカ国債」とは、アメリカ合衆国財務省が発行するアメリカ政府の債務証券です。国債には満期設定があり、3年、5年、10年、20年などさまざまな期間があります。「利回り」とは、投資金額に対して、"利息による収益"と"償還差損益（投資した金額と額面金額との差額）"との合計がどのくらいの割合になるかを示すものです。
★この指標を用いることでアメリカ国債に投資する場合の見通しを予測可能と考えています。

（３）金先物価格

「金先物」取引とは、将来の期日に決められた量の金を現時点で買う（もしくは売る）取引のことです。「シカゴ・マーカンタイル（CME）取引所」グループ傘下の「ニューヨーク商品取引所（COMEX）」で取引される「NY金先物」が有名で、世界の金先物の指標となっています。
★この指標を用いることで金に投資する場合の見通しやコモディティ投資（商品先物市場で取引されているエネルギー、貴金属、穀物といったような商品投資すること）の価格変動の見通しを予測可能と考えています。

（４）原油先物価格

「原油先物」取引とは、将来の期日に決められた量の原油を現時点で買う（もしくは売る）取引のことです。ニューヨーク・マーカンタイル取引所（NYMEX）に上場している「WTI原油先物」は取引量、市場参加者ともに圧倒的に多く、流動性や価格決定における透明性が高いため、原油価格の主要指標のひとつにもなっています。
★この指標を用いることで原油に投資する場合の見通しやコモディティ投資の価格変動の見通しを予測可能と考えています。

１.データの読み込み

・データの取得対象期間は2012/6/25～2022/6/24での10年間です。
・各指標の日足データの終値（取引終了時点の価格や指数）をその日の値とします。
・（１）（２）に関してはセントルイス連邦銀行が公表している「FRED®」より、pandas-datareaderライブラリを用いてデータを取得します。
・（３）（４）に関してはSBI証券提供のサービス「HYPER SBI Ⅱ」よりダウンロードしたcsvを、pandasライブラリを用いてデータを取得します。

#１.データの読み込み
#（共通）必要なモジュールのインポート
import pandas_datareader.data as web
import datetime as dt
import pandas as pd

#（１）S&P500指数（インターネット上より読み込み（再現可能））
start = dt.date(2012,6,25)
end = dt.date(2022,6,24)
code1 = 'SP500'
df1 = web.DataReader(code1, 'fred', start, end)

#（２）アメリカ国債利回り（10年）（インターネット上より読み込み（再現可能））
code2 = 'DGS10'
df2 = web.DataReader(code2, 'fred', start, end)

#（３）金先物価格（ローカル環境のドライブ上より読み込み（再現不可能））
data3 = pd.read_csv("/content/drive/MyDrive/成果物/金先物_day_20120625_20220624.csv", encoding='shift_jis')
# 不要な列削除
df3 = data3.drop(["始値", "高値", "安値", "5日平均", "25日平均", "75日平均", "VWAP", "出来高", "5日平均.1", "25日平均.1"], axis=1)
# 「日付」列の型変換
df3["日付"] = pd.to_datetime(df3["日付"])
# 「日付」列の昇順並べ替え
df3 = df3.sort_values(by="日付", ascending=True)
# 「日付」列をインデックスに指定
df3.set_index("日付", inplace=True)
# 「終値」列の型変換
#   列の要素の両端についている引用符の削除
df3["終値"] = df3["終値"].str.replace('"', '')
#   列の要素のカンマの削除
df3["終値"] = df3["終値"].str.replace(",", "")
#    列の要素の文字列型から数値型への変換
df3["終値"] = df3["終値"].astype("float")

#（４）原油先物価格（ローカル環境のドライブ上より読み込み（再現不可能））
data4 = pd.read_csv("/content/drive/MyDrive/成果物/原油_day_20120625_20220624.csv", encoding='shift_jis')
# 不要な列削除
df4 = data4.drop(["始値", "高値", "安値", "5日平均", "25日平均", "75日平均", "VWAP", "出来高", "5日平均.1", "25日平均.1"], axis=1)
# 「日付」列の型変換
df4["日付"] = pd.to_datetime(df4["日付"])
# 「日付」列の昇順並べ替え
df4 = df4.sort_values(by="日付", ascending=True)
# 「日付」列をインデックスに指定
df4.set_index("日付", inplace=True)

２. データの整理

・問題点として、金融商品取引市場の休日の行が取得したデータのままでは欠損しており、次章以降での「相関行列作成」や「予測値算出」にてエラーの原因となります。
・解決策として、データ取得対象期間の全ての日付をインデックスにもつダミーのDataframeを作成し、それと取得したデータを結合します。
・もともと休日だった行は欠損値となりますので、前営業日や翌営業日の値を取得するメソッドを用いて補完します。

#２.データの整理
#（共通）過去10年分の日付がインデックスとなっている空のDataframeを用意
df_date = pd.DataFrame(index=pd.date_range('2012-6-25', '2022-6-24'))

#（★）過去10年分の日付がインデックスとなっているDataframeへのマージにより、
#　　　・インデックスである日付の欠損の補完
#　　　・バリューである各指標値の欠損を前日値で補完
#（１）S&P500指数
df1_datefull = df1.merge(df_date, how="outer", left_index=True, right_index=True).fillna(method="bfill")
#（２）アメリカ国債利回り（10年）
df2_datefull = df2.merge(df_date, how="outer", left_index=True, right_index=True).fillna(method="ffill")
#（３）金先物価格
df3_datefull = df3.merge(df_date, how="outer", left_index=True, right_index=True).fillna(method="ffill")
#（４）原油先物価格
df4_datefull = df4.merge(df_date, how="outer", left_index=True, right_index=True).fillna(method="ffill")

３. データの可視化

■主成分分析の手法を一部取り入れた可視化

・４つのデータを１つのグラフに集約することで、傾向を考えます。

#３.データの可視化
#３－１.データの正規化
#（共通）必要なモジュールのインポート
from sklearn.preprocessing import MinMaxScaler
# (共通)正規化で使用する最小値と最大値を定義し、インスタンスを用意
mmscaler = MinMaxScaler(feature_range=(0, 1), copy=True)

#（★）バリューの正規化により、絶対値の異なる各指標値の比較の準備
#（１）S&P500指数
df1_nom = mmscaler.fit_transform(df1_datefull)
#（２）アメリカ国債利回り（10年）
df2_nom = mmscaler.fit_transform(df2_datefull)
#（３）金先物価格
df3_nom = mmscaler.fit_transform(df3_datefull)
#（４）原油先物価格
df4_nom = mmscaler.fit_transform(df4_datefull)

#３－２.データのグラフ化と傾向把握
#（共通）必要なモジュールのインポート
import matplotlib.pyplot as plt
#（共通）各指標値が折れ線グラフにて表示されるよう設定
plt.title("Timeseries-Chart")
plt.xlabel("time(day)")
plt.ylabel("absolute value")
plt.plot(df1_nom, label='S&P500')
plt.plot(df2_nom, label='DGS10')
plt.plot(df3_nom, label='Gold')
plt.plot(df4_nom, label='Oil')
plt.legend(loc='lower right')
plt.show()

結果

・主成分分析というさまざまな特徴量（：その結果の原因となりえる構成要素。体力テストの結果に対する生徒の身長・体重・年齢などにあたる）を集約する手法の一部を用います。データの「相関行列作成」を行い、特徴量ごとの類似度を考えます。

#３－３.データの相関行列作成と関係性把握
#（共通）必要なモジュールのインポート
import numpy as np
#（共通）各指標値の相関行列が表示されるよう設定
X = np.concatenate([df1_nom, df2_nom, df3_nom, df4_nom], axis=1)
R = np.corrcoef(X.T)
print(R)

結果

	S&P500指数	米国債利回り（10年）	金先物価格	原油先物価格
S&P500指数	1	-0.3096479	0.63927007	-0.12074947
米国債利回り（10年）	-0.3096479	1	-0.63176752	0.43031394
金先物価格	0.63927007	-0.63176752	1	0.16702285
原油先物価格	-0.12074947	0.43031394	0.16702285	1

■傾向　

1. S&P500指数と金先物価格間には正の比例関係が考えられる
2. アメリカ国債利回り（10年）と金先物価格間には負の比例関係が考えられる
3. アメリカ国債利回り（10年）と原油先物価格間には負の比例関係が考えられる
   　

４. データの周期の把握 (パラメーターsの決定)

・今回用いるSARIMAモデルでは最適なパラメーターセット(p,d,q)(sp,sd,sq)を求める必要があります。
・pは自己相関度といい、モデルが直前p個の値を用いて予測されるのかを表します。
・dは誘導といい、時系列データを定常にするためにd次の階差が必要だったことを表します。
・qは移動平均といい、モデルが直前q個の値に影響を受けることを表します。
・sは季節性といい、季節単位の周期性がある場合にその周期を表します。
・sを決め打ちするために、自己相関係数の算出、偏自己相関係数の算出、季節調整の利用を行います。

#４.データの周期の把握 (パラメーターsの決定)
#（共通）必要なモジュールのインポート
import statsmodels.api as sm
import itertools
#（１）S&P500指数
fig = sm.graphics.tsa.plot_acf(df1_datefull, lags=200)
plt.title("Autocorrelation-S&P500")
plt.show()
fig = sm.graphics.tsa.plot_pacf(df1_datefull, lags=25)
plt.title("Partial-Autocorrelation-S&P500")
plt.show()
sm.tsa.seasonal_decompose(df1_datefull, freq=100).plot()
plt.title("Seasonal-Index-S&P500")
plt.show()

結果

５. s以外のパラメーターの決定

・sを決め打ちしたら、(p,d,q)を求めます。
・SARIMAモデルのパラメーターを自動で最も適切にしてくれる機能はありません。そのため 情報量規準 (今回の場合は BIC(ベイズ情報量基準) )によってどの値が最も適切なのか調べます。BICの場合は この値が低ければ低いほどパラメーターの値は適切です。

#５.s以外のパラメーターの決定
# (共通)パラメーターセット算出のための関数を定義（情報量規準 (今回の場合は BIC(ベイズ情報量基準) )によってどの値が最も適切なのか調べるプログラム）
def selectparameter(DATA,s):
    p = d = q = range(0, 2)
    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 = sm.tsa.statespace.SARIMAX(DATA,
                                            order=param,
                                            seasonal_order=param_seasonal)
                results = mod.fit()
                parameters.append([param, param_seasonal, results.bic])
                BICs = np.append(BICs,results.bic)
            except:
                continue
    return parameters[np.argmin(BICs)]

#（１）S&P500指数
selectparameter(df1_datefull, 24)

結果
[(1, 1, 0), (0, 1, 1, 24), 33951.97624536817]

6. モデルの構築

・前章までで求められたパラメーターを用いてSARIMAモデルの構築を行います。

#６.モデルの構築
#（１）S&P500指数
SARIMA_df1 = sm.tsa.statespace.SARIMAX(df1_datefull, order=(1, 1, 0),seasonal_order=(0, 1, 1, 24)).fit()

７. データの予測とその可視化

・構築したモデルを用いて未来の予測を行います。
・予測期間は2022/6/25～2023/6/24での1年間です。

#７.データの予測とその可視化
#（１）S&P500指数
pred_df1 = SARIMA_df1.predict("2022-6-1", "2023-6-24")
plt.title("Predicted-Timeseries-Chart-S&P500")
plt.xlabel("time(day)")
plt.ylabel("value")
plt.plot(df1_datefull)
plt.plot(pred_df1, color="r")
plt.show()

★４～７までを（２）～（４）の指標に対しても同様に実施

#（２）アメリカ国債利回り（10年）
fig = sm.graphics.tsa.plot_acf(df2_datefull, lags=200)
plt.title("Autocorrelation-DGS10")
plt.show()
fig = sm.graphics.tsa.plot_pacf(df2_datefull, lags=25)
plt.title("Partial-Autocorrelation-DGS10")
plt.show()
sm.tsa.seasonal_decompose(df2_datefull, freq=50).plot()
plt.title("Seasonal-Index-DGS10")
plt.show()
#selectparameter(df2_datefull, 12)
SARIMA_df2 = sm.tsa.statespace.SARIMAX(df2_datefull, order=(0, 1, 0),seasonal_order=(0, 0, 0, 12)).fit()
pred_df2 = SARIMA_df2.predict("2022-6-1", "2023-6-24")
plt.title("Predicted-Timeseries-Chart-DGS10")
plt.xlabel("time(day)")
plt.ylabel("value")
plt.plot(df2_datefull)
plt.plot(pred_df2, color="r")
plt.show()

#（３）金先物価格
fig = sm.graphics.tsa.plot_acf(df3_datefull, lags=200)
plt.title("Autocorrelation-Gold")
plt.show()
fig = sm.graphics.tsa.plot_pacf(df3_datefull, lags=25)
plt.title("Partial-Autocorrelation-Gold")
plt.show()
sm.tsa.seasonal_decompose(df3_datefull, freq=100).plot()
plt.title("Seasonal-Index-Gold")
plt.show()
#selectparameter(df3_datefull, 12)
SARIMA_df3 = sm.tsa.statespace.SARIMAX(df3_datefull, order=(0, 1, 1),seasonal_order=(1, 1, 1, 24)).fit()
pred_df3 = SARIMA_df3.predict("2022-6-1", "2023-6-24")
plt.title("Predicted-Timeseries-Chart-Gold")
plt.xlabel("time(day)")
plt.ylabel("value")
plt.plot(df3_datefull)
plt.plot(pred_df3, color="r")
plt.show()

#（４）原油先物価格
fig = sm.graphics.tsa.plot_acf(df4_datefull, lags=200)
plt.title("Autocorrelation-Oil")
plt.show()
fig = sm.graphics.tsa.plot_pacf(df4_datefull, lags=25)
plt.title("Partial-Autocorrelation-Oil")
plt.show()
sm.tsa.seasonal_decompose(df4_datefull, freq=100).plot()
plt.title("Seasonal-Index-Oil")
plt.show()
#selectparameter(df4_datefull, 24)
SARIMA_df4 = sm.tsa.statespace.SARIMAX(df4_datefull, order=(0, 1, 1),seasonal_order=(1, 0, 0, 24)).fit()
pred_df4 = SARIMA_df4.predict("2022-6-1", "2023-6-24")
plt.title("Predicted-Timeseries-Chart-Oil")
plt.xlabel("time(day)")
plt.ylabel("value")
plt.plot(df4_datefull)
plt.plot(pred_df4, color="r")
plt.show()

結果

８．予測結果・考察・まとめ

■予測結果

（１）S&P500指数：2022年夏頃に現在の下落が底打ちして、以降は高い水準圏でゆるやかに上昇
（２）アメリカ国債利回り（10年）：2022年夏頃に現在の上昇が天井をつけ、以降は高い水準圏である3%前後で推移
（３）金先物価格：2022年夏頃まで小幅に上下し、以降は高い水準圏でゆるやかに上昇
（４）原油先物価格：2022年夏頃まで上下を繰り返しながら、以降は高い水準で推移

■考察

・アメリカ株式は、2022年夏以降ゆるやかだが高確率で上昇トレンドに入ると予測されるため、1年間の運用期間であれば投資すべき。【傾向1+予測結果(1)+予測結果(3)】
・アメリカ国債利回りは、2022年夏の水準で推移されると予測されるものの、負の比例関係にあるコモディティ価格の上昇幅によっては変化する可能性が高いため、様子見すべき。【傾向2+傾向3+予測結果(2)+予測結果(3)+予測結果(4)】
・金先物は、2022年夏以降ゆるやかだが高確率で上昇トレンドに入ると予測されるため、1年間の運用期間であれば投資すべき。【傾向1+傾向2+予測結果(1)+予測結果(2)+予測結果(3)】
・原油先物は、2022年夏の水準で推移されると予測されるものの、価格変動の激しいエネルギー分野という特性もあり、様子見すべき。【傾向3+予測結果(2)+予測結果(4)】

■まとめ

以前から関心のあったPythonについて講座を通してデータ分析に用いる基本知識を学べたと感じています。コードを書く中でやりたいと思ったデータの操作について、講座内容で見つからなくても調べることで対応するメソッドがほとんど見つかったことで、ライブラリの可用性の高さに驚きました。他の分析手法の利用の検討としましては、LSTMなどを候補と考えていまして、今後より精度の高い分析の実施を行っていきたいと思います。

参考

SBI証券　ＳＢＩ－ＳＢＩ・Ｖ・Ｓ＆Ｐ５００インデックス・ファンド　目論見書
(https://search.sbisec.co.jp/v2/popwin/info/connect/fund/8931119900000011.pdf)

SBI証券　債権の利率と利回り
(https://site3.sbisec.co.jp/ETGate/WPLETmgR001Control?OutSide=on&getFlg=on&burl=search_bond&cat1=bond&cat2=none&dir=info&file=bond_lecture_03.html)

SBI証券　金・銀・プラチナについてもっと知ろう！「金・銀・プラチナへの投資方法」
(https://site3.sbisec.co.jp/ETGate/WPLETmgR001Control?OutSide=on&getFlg=on&burl=search_gold&cat1=gold&cat2=guide&dir=guide&file=gold_guide_03.html)

みんかぶ　原油先物とは
(https://fu.minkabu.jp/beginner/column/11?utm_source=ca&utm_medium=beginner&utm_campaign=20220316)

SMBC日興証券　初めてでもわかりやすい用語集　コモディティ投資
(https://www.smbcnikko.co.jp/terms/japan/ko/J0244.html)