SARIMAX モデル

最終更新 2020-02-26 読了時間 3 分

まとめ

SARIMAX は ARIMA に季節項（Seasonal）と外生変数（eXogenous）を足したモデルで、order=(p, d, q) の非季節部分に加えて seasonal_order=(P, D, Q, s) で「12 か月前との差分」や「12 か月前の値・誤差との関係」を表します。たとえば seasonal_order=(1, 1, 1, 12) は、季節差分を 1 回取り、12 か月前の値と 12 か月前の誤差を 1 つずつ使う指定です。このページでは、月 1.2 の上昇と振幅 15・周期 12 の季節波を持つ 72 か月の系列を先頭 60 か月で学習し、statsmodels.tsa.statespace.sarimax.SARIMAX で残り 12 か月の予測平均と 95% 信頼区間を描きます。読み終えると、7 つの数の役割と、get_forecast() が返す区間の読み方が分かります。

ARIMA で order=(p, d, q) の 3 つの役割を先に押さえておくと、SARIMAX はそこへ季節版の 4 つの数と外生変数を足したものとして読めます。

SARIMAX は ARIMA に季節項（Seasonal）と外生変数（Exogenous）を拡張したモデルで、周期 12 の月次売上や曜日効果のある日次データのような実務の複雑な時系列に対応できます。

直感 #

月次の売上が年々増えながら毎年 8 月に山を作る系列を ARIMA(1,1,1) だけで追うと、差分でトレンドは消せても「12 か月前と似ている」という情報を使えず、夏の山を毎年外します。SARIMAX はこの系列に対して 2 種類の差分を取ります。d=1 の「今月 − 先月」でトレンドを消し、D=1 の「今月 − 12 か月前」で季節パターンを消すと、残りは平均が一定の系列になります。その残りを、先月の値・先月の誤差（非季節 AR/MA）と、12 か月前の値・12 か月前の誤差（季節 AR/MA）で説明するのが order と seasonal_order の 7 つの数です。外生変数 exog を渡せば、プロモーションや気温のような外部要因を回帰項として同時に推定できます。

原系列 d=1 で先月との差分 D=1 で 12 か月前との季節差分非季節 AR(1)/MA(1) + 季節 AR(1)/MA(1) で説明差分を足し戻して 12 か月先を予測（平均 + 95% 区間）

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：水準 → 傾き → 季節と要素を足すほど、未来への伸ばし方が豊かになる。SARIMAX は季節を係数ではなく季節差分と 12 か月前の値・誤差で表すが、予測線の形は 3 段目と同じ山と谷になる。

式で書くと 7 つの数はどこにあるか #

ラグ演算子 $B$（$B y_t = y_{t-1}$）を使うと、SARIMA(1,1,1)(1,1,1)$_{12}$ は

$$ (1 - \phi B)(1 - \Phi B^{12})(1 - B)(1 - B^{12})\, y_t = (1 + \theta B)(1 + \Theta B^{12})\, \varepsilon_t $$

と書けます。$(1 - B)$ が d=1 の通常差分、$(1 - B^{12})$ が D=1 の季節差分で、この 2 つで 13 か月分の観測が差分のために使われます。$\phi$ と $\theta$ が p=1, q=1 の非季節係数、$\Phi$ と $\Theta$ が P=1, Q=1 の季節係数で、s=12 が $B^{12}$ の 12 です。外生変数を加える場合は左辺の $y_t$ を $y_t - \beta x_t$ に置き換え、$\beta$ も同時に推定します。

Python で 12 か月先の予測と信頼区間を描く #

次のコードは、水準 200・月 1.2 の上昇・振幅 15 の季節波・標準偏差 3 のノイズを持つ 72 か月の系列を先頭 60 か月で学習し、get_forecast(steps=12) で予測平均と 95% 信頼区間を取り出して実測に重ねます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX

rng = np.random.default_rng(32)
dates = pd.date_range("2018-01-01", periods=6 * 12, freq="M")
trend = 1.2 * np.arange(len(dates))
seasonal = 15 * np.sin(2 * np.pi * np.arange(len(dates)) / 12)
noise = rng.normal(0, 3, len(dates))
series = pd.Series(200 + trend + seasonal + noise, index=dates)

train = series.iloc[:-12]
test_index = series.index[-12:]

model = SARIMAX(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
result = model.fit(disp=False)
forecast = result.get_forecast(steps=12)
pred_mean = forecast.predicted_mean
pred_ci = forecast.conf_int()

fig, ax = plt.subplots(figsize=(7.5, 4))
ax.plot(series.index, series, color="#cbd5f5", label="実測値（全体）")
ax.plot(train.index, train, color="#2563eb", linewidth=1.2, label="学習区間")
ax.plot(test_index, pred_mean, color="#f97316", linewidth=1.6, label="SARIMAX 予測")
ax.fill_between(
    test_index,
    pred_ci.iloc[:, 0],
    pred_ci.iloc[:, 1],
    color="#f97316",
    alpha=0.2,
    label="95% 信頼区間",
)
ax.set_title("SARIMAX モデルの予測")
ax.set_xlabel("年月")
ax.set_ylabel("値")
ax.legend()
ax.grid(alpha=0.3)

fig.tight_layout()
fig.savefig("static/images/timeseries/sarimax.svg")

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

7-11 行目6 周期分の月次データを作る
2018 年 1 月から 72 か月（6 年）の月末日付に、水準 200、月 1.2 の上昇、振幅 15・周期 12 の季節波、標準偏差 3 のノイズを足しています。季節差分 D=1 で 12 か月、通常差分 d=1 で 1 か月分の観測が差分に使われるので、学習に使える 60 か月から 13 か月を除いた 47 か月で係数を推定することになります。周期の 3 倍以上のデータがないと季節係数の推定が不安定になるため、月次なら 36 か月以上を目安にします。
13-14 行目最後の 1 周期を評価用に残す
先頭 60 か月（5 周期）で学習し、末尾 12 か月（1 周期）で山と谷を含めた予測精度を確かめます。評価区間が 1 周期に満たないと、季節パターンのどの部分だけ当たっているのかが判断できません。
16 行目order と seasonal_order の 7 つの数
order=(1, 1, 1) が先月との関係（AR 1、差分 1、MA 1）、seasonal_order=(1, 1, 1, 12) が 12 か月前との関係（季節 AR 1、季節差分 1、季節 MA 1、周期 12）です。D=0 にすると季節差分を取らず、季節パターンを季節 AR の係数だけで追うことになって予測の山が年々鈍ります。s を 52 にすれば週次データの年次季節、7 にすれば日次データの曜日効果になります。
17 行目fit(disp=False) で最適化の途中経過を隠す
fit() はカルマンフィルターで尤度を計算し、$\phi, \theta, \Phi, \Theta$ とノイズ分散を最尤推定します。disp=False は最適化アルゴリズムの反復ログを出さない指定で、省くと十数行の収束メッセージが表示されるだけで結果は変わりません。推定値は result.summary() で、AIC は result.aic で確認できます。
18-20 行目get_forecast は平均と区間を両方返す
forecast() が点予測だけを返すのに対し、get_forecast(steps=12) は予測分布を持つオブジェクトを返し、predicted_mean で平均、conf_int() で既定 95% の区間（alpha=0.05）を取り出せます。区間の幅は先の月ほど広がり、12 か月目では差分を 2 種類取っている分だけ ARIMA 単体より広くなります。conf_int(alpha=0.2) にすれば 80% 区間になります。
26-33 行目fill_between で区間を帯にする
pred_ci の 0 列目（下限）と 1 列目（上限）の間を橙の半透明（alpha=0.2）で塗っています。帯の中に薄い青の実測が 12 か月中 11〜12 か月入っていれば区間は妥当で、半分しか入らないなら分散を過小評価している（ノイズが正規分布から外れている、季節の振幅が年ごとに変わるなど）サインです。

Python コード例の図

実行結果の読み方 #

青の学習区間は 200 付近から 5 年で約 271（200 + 1.2 × 59）まで上がりながら、振幅 15 の季節波を 5 回繰り返します。橙の予測線はその続きとして 6 周期目の山と谷を描き、薄い青の実測とほぼ重なります。周囲の橙の帯が 95% 信頼区間で、1 か月目は標準偏差 3 のノイズに近い幅、12 か月目には差分と係数の不確かさが積み上がって 2〜3 倍に広がります。実測が帯の外に出る月が 12 か月中 2 か月以上あれば、seasonal_order の見直しや残差診断に戻ります。

モデル選択のコツ #

季節周期 $s$ は先に決めます。月次の年次季節なら 12、週次なら 52、日次の曜日効果なら 7 で、これはデータの ACF でラグ 12 や 52 が飛び出すことから確認できます。差分は入れすぎると過差分になってノイズが膨らむので、まず d=1, D=1 から試し、差分後の系列が一方向に傾いていなければそのまま進めます。候補の order・seasonal_order を 3〜4 通り fit() して result.aic を比べ、最後に result.plot_diagnostics() で残差が白色雑音に近いか（ACF が信頼帯に収まり、Q-Q プロットが直線に乗るか）を確認する、という順番が定番です。

活用シーン #

小売の需要予測ではプロモーションの有無や休日フラグを exog に渡して、季節性と施策効果を同時に推定します。エネルギー需要では気温・湿度を外生変数にすると、夏と冬の山を季節項だけで追うより残差が小さくなります。Web トラフィックではイベントカレンダーや広告出稿量を取り込み、週次季節（s=7）と組み合わせるのが典型です。外生変数の具体的な使い方は外生変数付き SARIMAXで扱います。

実装は statsmodels の SARIMAX クラスで足りますが、次数の自動探索が欲しいときは pmdarima.auto_arima を、季節が複数ある大規模データでは Prophet や LightGBM と並べて、同じ 12 か月の評価区間で誤差を比較してください。

ARIMA モデル — 季節性なしの古典的時系列予測
外生変数付き SARIMAX — 外生変数を加えて予測精度を向上
ホルト・ウィンター法 — トレンドと季節性を同時に平滑化して予測