【Pythonで需要予測②】SARIMAモデルの自動構築で出荷数量を予測してみた
前回は重回帰分析を使って、物流センターからの出荷数量の需要予測を行いました。
今回は時系列データの予測手法として有名なSARIMAモデルを使って、前回と同じ６か月分の出荷データから需要予測を行ってみました。
最初の５か月分のデータを使ってSARIMAモデルを構築し、残り１か月分のデータで予測値の精度を検証しました。
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モデルのパラメータp、d、qと、季節性を考慮するためのパラメータP、D、Q、sの合計７つがあります。
これらのパラメータは元データである出荷データを分析することによって予め当たりを付けることができます。
まず【SCM分析4】需要データの階差数列が定常過程になることをADF検定するで、元データは単位根過程であることがわかりました。
単位根過程とは１つ前のデータとの差分を取ったデータ（階差数列）が定常過程になることです。
このことから、p、d、qはすべて１であることが予想できます。
また【SCM分析3】Pythonで出荷データから周期変動と長期トレンドを分離するで、元データは７日間の周期を持つことがわかっています。
このことから、P、D、Qは１、ｓは７であることが予想できます。
このようにしてパラメータを自分で決めるやり方がSARIMAモデルの手動構築です。
学習データとテストデータに分ける
それではこのパラメータでSARIMAモデルを作る前に、半年分の元データを学習用とテスト用に分けておきましょう。
ここでは最初の５ヶ月分のデータで学習し、最後の１ヶ月分のデータでテスト（検証）することとします。
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では７つのパラメータを指定して学習データを入れれば、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)
２００を中心とする正規分布になっているように見えます。
それを確かめるために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)ｐ値は0.4なので正規分布と見なせます（厳密には正規分布でないということは言えない）が、残差のばらつきの中心はゼロではなく２８０付近になっています。
つまりホワイトノイズというには無理があります。
自動で決めたパラメータで需要予測
次にもっと良いパラメータがないかどうかを調べるために、ｓを除く６つのパラメータを０、１、２にぞれぞれ変化させて総当りしてみます。
６つのパラメータが３通りの値を取るため、
36=729通り
のSARIMAモデルを構築して最も良いモデルを探すことになります。
なお、７日周期であることは明らかなのでｓだけは７で固定します。
パラメータのすべての順列を求める
SARIMA関数にp、d、q、P、D、Qの６つのパラメータを入れるために、７２９通りある６つのパラメータをリストに格納します。
これは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)]これで７２９の順列がcombinationsに格納されました。
AICが最小になるSARIMAXモデルを求める
良いモデルを決める基準にはAIC（Akaike Information Criterion：赤池情報量規準）を使います。
モデルはパラメータの次数を増やして複雑にするほど学習データに対する適合性は上がりますが、テストデータに対する適合性は下がります（過学習）。
AICは程良く複雑なモデルを判定するための指標といえます。
このAICはSARIMA関数の結果の１つとして得ることができます。
まず、先程生成した順列を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]また７２９通りのパラメータの順列とそれに対応する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に７２９通りの結果が書き込まれたので、AICが小さい順に５つの順列を抜き出してみます。
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手動で構築したモデルでは３５８だったので改善されています。
残差の正規性
残差のプロット図は次の通りです。
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を指標とする方が良いモデルになりました。