【SCM分析6】AICを指標とするSARIMAモデルの自動構築で出荷数量を需要予測

2023年10月16日

前回は重回帰分析を使って、物流センターからの出荷数量の需要予測を行いました。

今回は時系列データの予測手法として有名なSARIMAモデルを使って、前回と同じ6か月分の出荷データから需要予測を行ってみました。

最初の5か月分のデータを使ってSARIMAモデルを構築し、残り1か月分のデータで予測値の精度を検証しました。

 

SARIMAモデルの原理についてはネット上でいろいろな記事がアップされていますが、こちらの記事がわかりやすかったです。

Pythonで時系列解析・超入門(その3)ARIMA系モデルで予測する方法

 

出荷データをPythonに読み込んでグラフにする

560アイテムについて取得した6ヶ月間の日々の出荷データをデータフレームraw_dfに読み込みます。

raw_df = pd.read_csv('demand.csv')
raw_df['Date'] = pd.to_datetime(raw_df['Date']) # 日付データを日付型に変換
raw_df

今回も前回に続いて、この中のS00239のアイテムについて分析することにします。

data_df = raw_df[['Date', 'S00239']]
data_df = data_df.set_index('Date')
data_df

グラフにすると次のようになります。

plt.figure(figsize=(10, 5))
plt.plot(data_df)
plt.show()

 

 

手動で決めたパラメータで需要予測

SARIMAモデルにはARIMAモデルのパラメータpdqと、季節性を考慮するためのパラメータPDQsの合計7つがあります。

これらのパラメータは元データである出荷データを分析することによって予め当たりを付けることができます。

まず【SCM分析4】需要データの階差数列が定常過程になることをADF検定するで、元データは単位根過程であることがわかりました。

単位根過程とは1つ前のデータとの差分を取ったデータ(階差数列)が定常過程になることです。

このことから、p、d、qはすべて1であることが予想できます。

また【SCM分析3】Pythonで出荷データから周期変動と長期トレンドを分離するで、元データは7日間の周期を持つことがわかっています。

このことから、P、D、Qは1、sは7であることが予想できます。

このようにしてパラメータを自分で決めるやり方がSARIMAモデルの手動構築です。

 

学習データとテストデータに分ける

それではこのパラメータでSARIMAモデルを作る前に、半年分の元データを学習用とテスト用に分けておきましょう。

ここでは最初の5ヶ月分のデータで学習し、最後の1ヶ月分のデータでテスト(検証)することとします。

sklearnのtrain_test_splitを使って、次のようにコーディングできます。

from sklearn.model_selection import train_test_split
df_train, df_test = train_test_split(data_df, test_size=31, shuffle=False)

これでdf_trainに学習データが、df_testにテストデータが入りました。

 

SARIMAX関数で予測値を求める

Pythonでは7つのパラメータを指定して学習データを入れれば、SARIMAモデルを構築してくれる関数statsmodels.tsa.statespace.sarimaxがあるので、これをSARIMAXとしてインポートしておきます。

from statsmodels.tsa.statespace.sarimax import SARIMAX

あとは最小限必要な情報を次のように指定することによって、best_modelにSARIMAXモデルが構築されます。

best_model = SARIMAX(df_train,
                    order = (1, 1, 1),
                    seasonal_order = (1, 1, 1, 7))
best_result = best_model.fit()

 

実績値と予測値をグラフで表示する

構築したモデルによる予測値はget_predictionを使って次のようにして得ることができます。

fitting = best_result.get_prediction()# 学習データ期間の予測値をget
fitting_mean = fitting.predicted_mean# 予測値の期待値をget
forecast = best_result.get_prediction('2019-03-01', '2019-03-31')# テストデータ期間の予測値をget
forecast_mean = forecast.predicted_mean# 予測値の期待値をget

テストデータ期間(今まで見ていないデータ)の予測値はget_forecastを使うことにより次のように書くこともできます。

forecast = best_result.get_forecast(len(df_test))# テストデータ期間の予測値をget
forecast_mean = forecast.predicted_mean# 予測値の期待値をget

グラフで表示すると次のようになります。

学習データの期間は予測値と実績値が近くなっていますが、テストデータ期間は少しズレが大きい感じです。

特に減少傾向のトレンドが反映されていません。

 

予測値を評価する

それではテストデータ期間の予測値を二乗平均平方根誤差残差の正規性で評価してみます。

 

二乗平均平方根誤差

これは簡単に次のように計算できます。

np.sqrt(mean_squared_error(df_test.S00239, forecast_mean))
…
357.632117066137

 

残差の正規性

まずは残差をグラフに表示してみます。

df_test['residuals'] = forecast_mean - df_test.S00239
df_test['residuals'].plot()

ゼロ近辺に均等にばらついている(ホワイトノイズ)というには無理があるように見えます。

次にヒストグラムで表示してみます。

xmin, xmax = df_test.residuals.min(), df_test.residuals.max()
mean, std = df_test.residuals.mean(), df_test.residuals.std()
x_axis = np.linspace(xmin, xmax, 100)
pdf = norm.pdf(x_axis, mean, std)

pdf_graph = pd.DataFrame({"X":x_axis, "norm":pdf}).set_index("X")
pdf_graph.plot()
df_test['residuals'].hist(density = True)

200を中心とする正規分布になっているように見えます。

それを確かめるためにKS検定してみます。

para = norm.fit(df_test.residuals)
stats.kstest(df_test.residuals, 'norm', args = para)
…
KstestResult(statistic=0.15481577481145348, pvalue=0.40595494337011895, statistic_location=286.5209028088717, statistic_sign=1)

p値は0.4なので正規分布と見なせます(厳密には正規分布でないということは言えない)が、残差のばらつきの中心はゼロではなく280付近になっています。

つまりホワイトノイズというには無理があります。

 

自動で決めたパラメータで需要予測

次にもっと良いパラメータがないかどうかを調べるために、sを除く6つのパラメータを0、1、2にぞれぞれ変化させて総当りしてみます。

6つのパラメータが3通りの値を取るため、

36=729通り

のSARIMAモデルを構築して最も良いモデルを探すことになります。

なお、7日周期であることは明らかなのでsだけは7で固定します。

 

パラメータのすべての順列を求める

SARIMA関数にp、d、q、P、D、Qの6つのパラメータを入れるために、729通りある6つのパラメータをリストに格納します。

これはitertoolsを使います。

import itertools
p = d = q = P = D = Q = range(0,3)
combinations = list(itertools.product(p, d, q, P, D, Q))
combinations
…
[(0, 0, 0, 0, 0, 0),
 (0, 0, 0, 0, 0, 1),
 (0, 0, 0, 0, 1, 0),
 (0, 0, 0, 0, 1, 1),
 (0, 0, 0, 1, 0, 0),
 (0, 0, 0, 1, 0, 1),
…
 (1, 1, 1, 1, 0, 1),
 (1, 1, 1, 1, 1, 0),
 (1, 1, 1, 1, 1, 1)]

これで729の順列がcombinationsに格納されました。

 

AICが最小になるSARIMAXモデルを求める

良いモデルを決める基準にはAIC(Akaike Information Criterion:赤池情報量規準)を使います。

モデルはパラメータの次数を増やして複雑にするほど学習データに対する適合性は上がりますが、テストデータに対する適合性は下がります(過学習)。

AICは程良く複雑なモデルを判定するための指標といえます。

このAICはSARIMA関数の結果の1つとして得ることができます。

 

まず、先程生成した順列をarima_orders(p、d、q)seasonal_orders(P、D、Q、s)に分けます。

周期を表すsはseasonal_ordersに含めます。

s = 7
arima_orders = [(x[0], x[1], x[2]) for x in combinations]
seasonal_orders = [(x[3], x[4], x[5], s) for x in combinations]

また729通りのパラメータの順列とそれに対応するAICの値を入れるデータフレームを次のように作っておきます。

results = pd.DataFrame(columns=['p', 'd', 'q', 'P', 'D', 'Q', 'AIC'])

そして、for文を使ってすべての順列を順番にSARIMA関数に入れて、結果をこのデータフレームresultsに書き込みます。

for i in range(len(combinations)):
    try:
        model = SARIMAX(df_train,
                    order = arima_orders[i],
                    seasonal_order = seasonal_orders[i])
        result = model.fit()
        results.loc[i, 'p'] = arima_orders[i][0]
        results.loc[i, 'd'] = arima_orders[i][1]
        results.loc[i, 'q'] = arima_orders[i][2]
        results.loc[i, 'P'] = seasonal_orders[i][0]
        results.loc[i, 'D'] = seasonal_orders[i][1]
        results.loc[i, 'Q'] = seasonal_orders[i][2]
        results.loc[i, 'AIC'] = result.aic
        print(i)
    except:
        continue

これでresultsに729通りの結果が書き込まれたので、AICが小さい順に5つの順列を抜き出してみます。

results[results.AIC.rank(ascending=True) < 5]

(p、d、q、P、D、Q)=(1、1、2、1、2、2)の順列のパラメータの時にAICが最小になることがわかりました。

 

実績値と予測値をグラフで表示する

このモデルでできた予測値と実績値をグラフで表示させてみます。

best_model = SARIMAX(df_train,
                    order = (1, 1, 2),
                    seasonal_order = (1, 2, 2, 7))
best_result = best_model.fit()
fitting = best_result.get_prediction()
fitting_mean = fitting.predicted_mean
forecast = best_result.get_forecast(len(df_test))# テストデータ期間の予測値をget
forecast_mean = forecast.predicted_mean# 予測値の期待値をget

plt.plot(fitting_mean, label = 'Predicted_train_data')
plt.plot(df_train, label = 'Actual_train_data')
plt.plot(forecast_mean, label = 'Predicted_test_data')
plt.plot(df_test, label = 'Actual_test_data')
plt.legend()
plt.show()

手動でパラメータを決めたモデルよりもテストデータ期間における予測精度が上がっていることがわかります。

減少トレンドもしっかり反映できているようです。

 

予測値を評価する

二乗平均平方根誤差

np.sqrt(mean_squared_error(df_test.S00239, forecast_mean))
…
278.06369214340816

手動で構築したモデルでは358だったので改善されています。

 

残差の正規性

残差のプロット図は次の通りです。

df_test['residuals'] = forecast_mean - df_test.S00239
df_test['residuals'].plot()

ヒストグラムは次の通りです。

xmin, xmax = df_test.residuals.min(), df_test.residuals.max()
mean, std = df_test.residuals.mean(), df_test.residuals.std()
x_axis = np.linspace(xmin, xmax, 100)
pdf = norm.pdf(x_axis, mean, std)

pdf_graph = pd.DataFrame({"X":x_axis, "norm":pdf}).set_index("X")
pdf_graph.plot()
df_test['residuals'].hist(density = True)

ゼロを中心とする正規分布に近づきました。

KS検定の結果を見ても、そのことが裏付けられます。

para = norm.fit(df_test.residuals)
stats.kstest(df_test.residuals, 'norm', args = para)
…
KstestResult(statistic=0.10494420738598431, pvalue=0.849378740264486, statistic_location=-80.63771642296763, statistic_sign=-1)

これはplot_diagnosticsでモデル診断してみても、良いモデルになっていることがわかります。

best_result.plot_diagnostics(lags=20);

 

 

RMSEでSARIMAXモデルを評価すると。。。

AICを最適化の指標とするSARIMAモデルの自動構築により、手動構築よりも良いモデルを構築することができました。

ではRMSE(二乗平均平方根誤差)を指標として最適モデルを自動構築したらどうなるでしょうか?

試してみたら次のようになりました。

for i in range(len(combinations)):
    try:
        model = SARIMAX(df_train,
                    order = arima_orders[i],
                    seasonal_order = seasonal_orders[i])
        result = model.fit()
        fitting = result.get_prediction()
        fitting_mean = fitting.predicted_mean
        rmse = np.sqrt(mean_squared_error(df_train.S00239, fitting_mean))
        results.loc[i, 'p'] = arima_orders[i][0]
        results.loc[i, 'd'] = arima_orders[i][1]
        results.loc[i, 'q'] = arima_orders[i][2]
        results.loc[i, 'P'] = seasonal_orders[i][0]
        results.loc[i, 'D'] = seasonal_orders[i][1]
        results.loc[i, 'Q'] = seasonal_orders[i][2]
        results.loc[i, 'RMSE'] = rmse
        print(i)
    except:
        continue
results[results.RMSE.rank(ascending=True) < 5]

そして最適パラメータによるモデルで予測値と実績値を表示させたら次のようになりました。

明らかにAICを指標とした最適化モデルより予測精度が劣っています。

残差の二乗誤差が最小になるモデルを構築するよりも、AICを指標とする方が良いモデルになりました。