【Python】【機械学習】Prophetでの株価予測モデルの精度検証とハイパーパラメータ最適化

2023年6月3日 20:10

Prophetではチューニングできるパラメータ（ハイパーパラメータ）が多くあります。特に、各パラメータが大きい影響力を持っていて、値を少し変えただけでも予測値が大きく変化します。このハイパーパラメータを調整し最適な設定を探すタスクを、ハイパーパラメータチューニングと言います。

ここでは、ハイパーパラメータチューニングをベイズ最適化で行う代表的なライブラリ Oputuna を使って、より良いハイパーパラメータの組合せを探索し、予測精度を比較したいと思います。

なお、最低限のポイントのみの説明にするため、Pythonライブラリ、モジュール等のインストール方法については割愛させて頂きます。お使いのPC環境等に合わせてインストールしてもらえればと思います。

１．ライブラリをインポートして株価データを取得する

時系列機械学習ライブラリである fbprophet から Prophet をインポートします。oputuna、mean_absolute_error、mean_squared_errorをインポートします。

import pandas_datareader.data as web
import datetime

# visualization
import matplotlib.pyplot as plt
%matplotlib inline

from prophet import Prophet
import optuna
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error

２．株価データを取得してデータを確認

下記を参考にOHLCV（始値 / 高値 / 安値 / 終値 / 出来高）形式の日経平均株価(^NKX)データを取得します。データの取得期間は、2018年1月1日から2023年6月/1日までです。

start = '2018/01/01'
end = '2023/06/01'

df = web.DataReader('^NKX', 'stooq', start, end)

# 日付を昇順に並び替える
df.sort_index(inplace=True)

# インデックスをリセット
df = df.reset_index()

#　株価データを可視化
fig, ax = plt.subplots(figsize=(15,7))
ax.plot(df['Date'], df['Close'], lw=4)
ax.set_title("Stock price(Close)",fontsize=15)
plt.show()

３．基礎分析をする

まず、可視化して年、月、週の周期性がありそうか確認する

# 年、月、週の周期性がありそうか確認するため、年、月、週、日カラムを追加
df['year'] = df.Date.apply(lambda x: x.year)
df['month'] = df.Date.apply(lambda x: x.month)
df['weekday'] = df.Date.apply(lambda x: x.dayofweek)
df['day'] = df.Date.apply(lambda x: x.day)

年ごとに、月、週、日の終値の平均値を可視化する

# 年ごとに月、週、日の終値平均値を可視化
for i  in [2018,2019,2020,2021,2022,2023]:
    tmp_df = df[df.year==i].copy()
    tmp_df.sort_values('Date',inplace=True)
    
    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 4))
    
    # 月平均値のグラフ化
    tmp_df.groupby('month').mean()['Close'].plot(kind='line', title='average stock price(month)', lw=4, color='r', ax=ax1)
    
    # 週平均値のグラフ化
    tmp_df.groupby('weekday').mean()['Close'].plot(kind='line', title='average stock price(weekday)', lw=4, color='g', ax=ax2)
    
    # 日平均値のグラフ化
    tmp_df.groupby('day').mean()['Close'].plot(kind='line', title='average stock price(day)', lw=4, color='b', ax=ax3)
    
    fig.suptitle('Trends for year: {0}'.format(i), size=20, y=0.95)
    plt.tight_layout()
    plt.show()

2018年から2023年まで総合的に見た場合、明確な周期性はなさそうです。

４．Prophetで予測する

今回は、以下の条件で株価予測を行います。

対象：日経平均株価(^NKX)
学習期間：2018/1/1～2023/1/1
予測期間：2023/1/1以降、2023/6/1までの終値を予測

４. 1 デフォルトパラメータでの予測と精度検証

Prophetを用いるため、ds カラムと y カラムを追加し、train データと testデータに分割して確認のためグラフを描画します。trainデータセット（df_train）は青で、testデータセット（df_test）は橙で描かれており、2023年以降は、df_test になっておりデータが分割できていることが確認できます。

# 新たに ds列と y列を追加
df['ds'] = df['Date']
df['y'] = df['Close']

# trainデータとtestデータに分割
df_train = df[df['ds'] < '2023-01-01']
df_test = df[df['ds'] >= '2023-01-01']

#　株価データを可視化
fig, ax = plt.subplots(figsize=(15,7))
ax.plot(df_train.ds, df_train.y, label="actual(train dataset)", lw=4)
ax.plot(df_test.ds, df_test.y, label="actual(test dataset)", lw=4)
ax.set_title("Stock price(Close)",fontsize=15)
plt.legend()
plt.show()

予測モデルを構築し、test 期間を含めて予測します。そして予測値と実測値からProphetモデルを評価します。評価には、RMSE、MAE、MAPEという指標を使います。

・RMSE（Root Mean Squared Error）
予測値と実測値の誤差を二乗して平均し、平方根をとった値です。MSEに平方根をつけたバージョンです。RMSEは予測値と実測値の誤差の大きさを表す指標で、二乗を取るため、大きな誤差がある場合はその影響が大きく出ます。RMSEが小さいほど、予測精度が高いと言えます。MSEを損失関数として使って，RMSEを評価指標として使うことも多いです。

・MAE（Mean Absolute Error）
予測値と実測値の誤差の絶対値を平均した値です。MAEは予測値と実測値の誤差の大きさを表す指標で、二乗を取らないため、大きな誤差があっても影響が小さくなります。MAEが小さいほど、予測精度が高いと言えます。

・MAPE（Mean Absolute Percentage Error）
予測値と実測値の差を実測値で割った値（パーセント誤差）の絶対値を平均した値です。MAPEは、予測値と実測値の相対的な誤差を表す指標で、予測値が小さい場合や0の場合には使用できません。100％の確率値にするため、一般的には最後に100をかけます。MAPEが小さいほど、予測精度が高いと言えます。

# 予測モデル構築
m = Prophet()
m.fit(df_train)

# 予測の実施（学習期間＋テスト期間）
df_future = m.make_future_dataframe(periods=len(df_test))
df_pred = m.predict(df_future)

# 元のデータセットに予測値を結合
df['Predict'] = df_pred['yhat']

# 予測精度（テストデータ期間）
## 予測値と実測値
preds = df.iloc[len(df_train):].loc[:, 'Predict']
y = df.iloc[len(df_train):].loc[:, 'y']

## 指標出力
print('RMSE:')
print(np.sqrt(mean_squared_error(y, preds)))
print('-----------------')

print('MAE:')
print(mean_absolute_error(y, preds))
print('-----------------')

print('MAPE(%):')
print(np.mean(abs(y - preds)/y)*100)

最後にグラフを描画して予測時と実測値を比較します。

#　株価データを可視化
fig, ax = plt.subplots(figsize=(15,7))

# テスト期間の実測値と予測値を描画
ax.plot(df.iloc[len(df_train):]['Date'], df.iloc[len(df_train):]['y'], label="Actual")
ax.plot(df.iloc[len(df_train):]['Date'], df.iloc[len(df_train):]['Predict'], label="Forecast")
ax.set_title("Forecast vs Actual",fontsize=15)

plt.legend()
plt.show()

４. ２ Optunaによるパラメータチューニング

Optunaでハイパーパラメータをチューニングします。まずは、以下の手順で最適パラメータの探索します。

目的関数の設定
ハイパーパラメータの探索の実施
最適パラメータの出力

・changepoint_prior_scale
Prophetのハイパーパラメータの1つで最も影響力のあるパラメータです。時系列データの変化点を検出する際の柔軟性を調整するためのスケーリングパラメータです。このパラメータの値が大きいほど、モデルはより柔軟になり、小さい値の場合はより制限されます。[0.001, 0.5]の範囲がほぼ正しいでしょう。

・seasonality_prior_scale
Prophetのハイパーパラメータの1つで、季節性の柔軟性を制御します。値が大きいと季節性が大きな変動に適合し、値が小さいと季節性の大きさが小さくなります。デフォルトは 10 で、基本的に正則化は適用されません。調整の妥当な範囲はおそらく [0.01, 10] です。

・seasonality_mode
Prophetのハイパーパラメータの1つで、季節性のモードを指定するパラメータです。オプションは [ ‘additive’、’multiplicative’] です。デフォルト‘additive’です。additiveは季節性が加算され、multiplicativeは乗算されます。

・changepoint_range
Prophetのハイパーパラメータの1つで、変化点を検出する期間を制限するパラメータです。このパラメータの値が小さいほど、モデルは変化点をより早く検出するようになります。このパラメーターは、おそらく多数の時系列を除いて、調整しないほうがよいでしょう。設定では、[0.8, 0.95]が妥当な範囲かもしれません。

・n_changepoints
Prophetのハイパーパラメータの1つで、変化点の数を指定するパラメータです。このパラメータの値が大きいほど、モデルはより多くの変化点を検出するようになります。

まず目的関数を設定し、ハイパーパラメータの探索を実施して最適なパラメータを出力します。

# 目的関数の設定
def objective(trial):

    params = {
        'changepoint_prior_scale' : trial.suggest_uniform('changepoint_prior_scale', 0.001,0.5),
        'seasonality_prior_scale' : trial.suggest_uniform('seasonality_prior_scale', 0.01,10),
        'seasonality_mode' : trial.suggest_categorical('seasonality_mode', ['additive', 'multiplicative']),
        'changepoint_range' : trial.suggest_discrete_uniform('changepoint_range', 0.8, 0.95, 0.001),
        'n_changepoints' : trial.suggest_int('n_changepoints', 20, 35),
    }
    
    m = Prophet(**params)
    m.fit(df_train)
    
    df_future = m.make_future_dataframe(periods=len(df_test), freq='D')
    df_pred = m.predict(df_future) 
    preds = df_pred.tail(len(df_test))
    val_rmse = np.sqrt(mean_squared_error(df_test.y, preds.yhat))
    return val_rmse


# 目的関数の最適化を実行する
study = optuna.create_study(direction="minimize")
study.optimize(objective, n_trials=100)


# 最適パラメータの出力
print(f"The best value is : \n {study.best_value}")
print(f"The best parameters are : \n {study.best_params}")

次に、上記の最適化したパラメータを設定して予測と構築したProphetモデルを評価します。評価には、RMSE、MAE、MAPEという指標を使います。

# 最適パラメータで予測モデル構築
m = Prophet(**study.best_params)
m.fit(df_train)

# 予測の実施（学習期間＋テスト期間）
df_future = m.make_future_dataframe(periods=len(df_test))
df_pred = m.predict(df_future)

# 元のデータセットに予測値を結合
df['Predict'] = df_pred['yhat']

# 予測精度（テストデータ期間）
## 予測値と実測値
preds = df.iloc[len(df_train):].loc[:, 'Predict']
y = df.iloc[len(df_train):].loc[:, 'y']

## 指標出力
print('RMSE:')
print(np.sqrt(mean_squared_error(y, preds)))
print('-----------------')

print('MAE:')
print(mean_absolute_error(y, preds))
print('-----------------')

print('MAPE(%):')
print(np.mean(abs(y - preds)/y)*100)