季節ARIMA（SARIMA）｜時系列分析

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：ARMA・ARIMAモデル　|　関連：時系列データと分解・予測の評価指標と時系列CV

要点（BLUF）

SARIMA $(p,d,q)(P,D,Q)_s$ ＝ARIMA に季節周期 $s$ の部品（季節差分・季節AR・季節MA）を掛け合わせたモデル。月次なら $s=12$ 、四半期なら $s=4$ 。
季節差分 $(1-L^s)y_t=y_t-y_{t-s}$ は「前年同期との差」。季節の単位根を除いて定常化し、ラグ $s,2s,\dots$ の高い自己相関を落とします（コードで実証）。
予測は季節の波が乗った点予測＋予測区間。季節素朴予測（前年同期、予測の評価指標と時系列CV）をベースラインに、それを上回って初めて価値があります。

1. SARIMA：ARIMA に季節の部品を掛ける

非季節 ARIMA（ARMA・ARIMAモデル）は隣接ラグ（ $L,L^2,\dots$ ）の依存を扱いました。季節データはそれに加えて周期 $s$ ごとの依存（今年4月と去年4月が似る）を持ちます。SARIMA は、非季節の部品と季節の部品（ラグを $s$ 単位にしたもの）を掛け合わせて両方を同時に表します。

\underbrace{\Phi_P(L^s)\,\phi_p(L)}_{\text{AR（季節 × 非季節）}}\ \underbrace{(1-L)^d(1-L^s)^D}_{\text{差分（季節 × 非季節）}}\ y_t =\underbrace{\Theta_Q(L^s)\,\theta_q(L)}_{\text{MA（季節 × 非季節）}}\ \varepsilon_t

$\phi_p(L)=1-\phi_1L-\cdots-\phi_pL^p$ ：非季節 AR（AR・MAモデル）。 $\Phi_P(L^s)=1-\Phi_1L^s-\cdots-\Phi_PL^{Ps}$ ：ラグ $s$ 単位の季節 AR。
$(1-L)^d$ ：通常差分（トレンド除去）。 $(1-L^s)^D$ ：季節差分（前年同期差・季節の単位根除去）。
$\theta_q(L),\ \Theta_Q(L^s)$ ：非季節 MA と季節 MA。

要するに「非季節 ARIMA の各部品に、ラグを $s$ にした双子をもう一つ用意して掛ける」だけ。 $(P,D,Q)_s$ が季節側の $(p,d,q)$ に対応します。月次データで最頻出はエアライン・モデル SARIMA $(0,1,1)(0,1,1)_{12}$ で、通常差分＋季節差分＋通常MA(1)＋季節MA(1)です。

2. 季節差分が季節ラグの自己相関を落とす

季節性があると、ACF はラグ $s,2s,\dots$ で大きな山を作ります（前年同期と連動）。トレンドだけ取る通常差分ではこの季節の山は残り、季節差分 $(1-L^s)$ が初めてそれを落とします。トレンド＋月別の季節パターン＋AR(1)ノイズという擬似系列で、ラグ12・24の ACF が季節差分の前後でどう変わるかを見ます。

import numpy as np
from statsmodels.tsa.stattools import acf

# 真の季節構造（s=12）: トレンド + 月別固定効果(季節) + AR(1)ノイズ
np.random.seed(3)
s, n = 12, 180
t = np.arange(n)
month_effect = np.array([0, 2, 5, 8, 10, 9, 6, 3, -2, -6, -8, -5.0])  # 12ヶ月の季節パターン(平均≒0)
trend = 0.08 * t
noise = np.zeros(n)
for k in range(1, n):
    noise[k] = 0.5 * noise[k-1] + np.random.normal(0, 1.0)
y = 20 + trend + month_effect[t % s] + noise

d1 = np.diff(y)              # 1階差分（トレンド除去・季節は残る）
seas_diff = y[s:] - y[:-s]   # 季節差分 (1-L^12)
a1 = acf(d1, nlags=24, fft=True)
a2 = acf(seas_diff, nlags=24, fft=True)
print("季節差分の前後で ACF の季節ラグ(12,24)を比較")
print(f"{'lag':>6}{'1階差分のみ':>14}{'季節差分後':>12}")
for lag in (12, 24):
    print(f"{lag:>6}{a1[lag]:>14.3f}{a2[lag]:>12.3f}")

出力：

季節差分の前後で ACF の季節ラグ(12,24)を比較
   lag        1階差分のみ       季節差分後
    12         0.821      -0.447
    24         0.772       0.053

出力の意味：通常の1階差分だけでは、ラグ12で $0.821$ 、ラグ24で $0.772$ と季節の自己相関が高いまま残ります（トレンドは消えても季節の波は残る）。季節差分 $(1-L^{12})$ を施すと、ラグ24は $0.053$ までほぼ消え、ラグ12は $-0.447$ に転じます（過剰差分気味の負の自己相関——これは季節MA(1)で吸収するのが定石で、エアライン・モデルの季節MA項の出番）。「季節の山が ACF に残るなら季節差分（ $D=1$ ）を入れる」——これが SARIMA の $D$ を決める基本動作です（定常性と自己相関の ACF 読みの季節版）。

3. コード②：真の季節AR係数を復元する

SARIMA が季節の依存を本当に推定できるかを、真の係数を仕込んだ乗法的季節ARで確かめます。 $(1-\phi L)(1-\Phi L^{12})y_t=\varepsilon_t$ （ $\phi=0.5,\ \Phi=0.6$ ）は展開すると $y_t=\phi y_{t-1}+\Phi y_{t-12}-\phi\Phi y_{t-13}+\varepsilon_t$ という単純な漸化式で生成でき、定常（差分不要）なので係数をきれいに復元できます。

import numpy as np
import warnings
warnings.simplefilter("ignore")
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import acf

# 真の乗法的季節AR: (1 - φL)(1 - Φ L^12) y_t = ε_t
# 展開: y_t = φ y_{t-1} + Φ y_{t-12} - φΦ y_{t-13} + ε_t
phi_true, Phi_true = 0.5, 0.6
s, n = 12, 1200
np.random.seed(8)
eps = np.random.normal(0, 1.0, n)
y = np.zeros(n)
for t in range(13, n):
    y[t] = phi_true*y[t-1] + Phi_true*y[t-12] - phi_true*Phi_true*y[t-13] + eps[t]

res = SARIMAX(y, order=(1, 0, 0), seasonal_order=(1, 0, 0, s), trend="c").fit(disp=False)
phi_hat = res.params[res.param_names.index("ar.L1")]
Phi_hat = res.params[res.param_names.index("ar.S.L12")]
print(f"{'係数':>12}{'真値':>8}{'推定値':>10}")
print(f"{'φ(非季節AR)':>12}{phi_true:>8.3f}{phi_hat:>10.3f}")
print(f"{'Φ(季節AR)':>12}{Phi_true:>8.3f}{Phi_hat:>10.3f}")

a = acf(y, nlags=26, fft=True)
print(f"ACF: lag1={a[1]:.3f}  lag12={a[12]:.3f}（季節スパイク）  lag24={a[24]:.3f}")

出力：

          係数      真値       推定値
    φ(非季節AR)   0.500     0.472
     Φ(季節AR)   0.600     0.605
ACF: lag1=0.472  lag12=0.602（季節スパイク）  lag24=0.376

出力の意味：非季節AR $\hat\phi=0.472$ 、季節AR $\hat\Phi=0.605$ と、仕込んだ $0.5,\ 0.6$ を復元しました（ $\phi$ は標本のばらつきで少しずれますが許容範囲）。seasonal_order=(1,0,0,12) の ar.S.L12 がラグ12の季節AR係数です。ACF はラグ12で $0.602$ とくっきり季節スパイク——隣接ラグだけでなく周期 $s$ ごとに山が立つのが季節データの指紋です。SARIMA はこの季節の山を ar.S.L12 で捉えています。

4. コード③：季節素朴を上回る予測（予測区間つき）

最後に、トレンド＋季節を持つ系列の最後の2年（24ヶ月）をホールドアウトし（時間順分割・予測の評価指標と時系列CV）、SARIMA $(1,1,1)(0,1,1)_{12}$ で予測します。予測区間つきで、季節素朴予測（前年同月）と RMSE を比べます。

import numpy as np
import warnings
warnings.simplefilter("ignore")
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
from statsmodels.tsa.statespace.sarimax import SARIMAX

# 真の季節構造（s=12）: トレンド + 月別固定効果 + AR(1)ノイズ
np.random.seed(3)
s, n = 12, 180
t = np.arange(n)
month_effect = np.array([0, 2, 5, 8, 10, 9, 6, 3, -2, -6, -8, -5.0])
trend = 0.08 * t
noise = np.zeros(n)
for k in range(1, n):
    noise[k] = 0.5 * noise[k-1] + np.random.normal(0, 1.0)
y = 20 + trend + month_effect[t % s] + noise

h = 24                              # 最後の2年をホールドアウト（時間順）
train, test = y[:-h], y[-h:]

res = SARIMAX(train, order=(1, 1, 1), seasonal_order=(0, 1, 1, s)).fit(disp=False)
fc = res.get_forecast(steps=h)
mean = fc.predicted_mean
ci = fc.conf_int(alpha=0.05)

snaive = y[n-h-s:n-s]               # 季節素朴予測（前年同月）
rmse_sarima = np.sqrt(np.mean((mean - test)**2))
rmse_snaive = np.sqrt(np.mean((snaive - test)**2))
print(f"24ヶ月先 RMSE: SARIMA={rmse_sarima:.3f}  季節素朴(前年同月)={rmse_snaive:.3f}")
print(f"改善率 = {(1 - rmse_sarima/rmse_snaive)*100:.1f}%（季節素朴比）")
print(f"予測区間幅: 1ヶ月先={ci[0,1]-ci[0,0]:.2f}  24ヶ月先={ci[h-1,1]-ci[h-1,0]:.2f}")

fidx = np.arange(n-h, n)
plt.figure(figsize=(10, 4.5))
plt.plot(np.arange(n), y, color="gray", lw=1, label="実測")
plt.plot(fidx, mean, color="C1", lw=2, label="SARIMA点予測")
plt.fill_between(fidx, ci[:, 0], ci[:, 1], color="C1", alpha=0.25, label="95%予測区間")
plt.axvline(n-h-1, ls=":", color="k")
plt.xlabel("月インデックス t"); plt.ylabel("y"); plt.legend()
plt.title("SARIMA(1,1,1)(0,1,1)_12 の予測（季節の波が乗り、区間は先ほど広がる）")
plt.tight_layout(); plt.show()

出力：

24ヶ月先 RMSE: SARIMA=1.012  季節素朴(前年同月)=1.481
改善率 = 31.7%（季節素朴比）
予測区間幅: 1ヶ月先=4.07  24ヶ月先=4.88

出力の意味：SARIMA の RMSE は $1.012$ で、季節素朴（前年同月）の $1.481$ より31.7% 良い。季節素朴は「去年と同じ」と仮定するだけなので、トレンドやノイズの変化を取りこぼします。SARIMA は季節差分で季節を、通常差分でトレンドを取り込むので、季節の波に加えて水準の上昇も予測に乗ります（図の橙線がギザギザに上昇）。予測区間は1ヶ月先 $4.07$ から24ヶ月先 $4.88$ へ広がります（和分による不確実性の累積）。ベースライン（季節素朴）を上回って初めて SARIMA に意味がある、という評価姿勢は非季節と同じです。

flowchart LR
    A["原系列 y_t（トレンド + 季節）"] --> B["通常差分 (1-L)^d でトレンド除去"]
    B --> C["季節差分 (1-L^s)^D で季節除去"]
    C --> D["定常系列に 季節ARMA × 非季節ARMA を推定"]
    D --> E["予測 + 予測分散"]
    E --> F["逆差分（通常・季節）で水準へ戻す"]
    F --> G["季節の波が乗った点予測 + 予測区間"]

5. 数式の直観

季節差分＝前年同期との比較： $(1-L^{12})y_t=y_t-y_{t-12}$ は「去年の同じ月からの変化」。季節パターンが毎年ほぼ同じなら、この差で季節が消えて定常になります。通常差分（前月差）と役割が別なので、両方必要なことが多い。
掛け算（乗法）である理由：季節と非季節の依存が相互作用するから。 $\phi(L)\Phi(L^s)$ を展開すると、ラグ $1,\ s,\ s+1,\dots$ に係数が立ち、「隣接の効果と季節の効果が組み合わさったラグ」まで自然に表現できます。
予測の不確実性：季節差分も和分（積分）の一種なので、逆差分で予測を水準に戻すとき分散が累積し、区間が先ほど広がります（ARMA・ARIMAモデルの $\sqrt{h}$ 則の季節版）。

補足：次数の自動選択（要最新確認・任意）

$(p,d,q)(P,D,Q)_s$ の組合せは多く、手で探すのは大変です。pmdarima の auto_arima は AIC を基準に季節次数まで自動探索しますが、API が変わりやすく依存も重いので要最新確認。本ノートのコードは statsmodels だけで完結します。次数選択と妥当性検証の作法は次ノート（モデル選択と残差診断）で扱います。

⚠️ よくある誤解

「季節差分すれば必ず良くなる」ではない：季節が決定論的（毎年完全に同じ）なら季節ダミー回帰の方が素直で、季節差分は過剰差分で人工的な負の季節自己相関（コード①のラグ12が $-0.447$ ）を生みます。 $D$ は 0 か 1 に留め、ACF と季節パターンの安定性で判断します。
「 $s$ は何でもよい」ではない： $s$ はデータの周期そのもの（月次12・四半期4・週次7や52）。間違えると季節を捉えられません。複数周期（週と年など）が混在するとき SARIMA 単独では苦しく、状態空間や Prophet 系の出番（第4章・第7章）。
「季節素朴に勝てば十分」ではない：勝つのは最低条件。残差に季節の山が残っていないか（モデル選択と残差診断の残差ACF・Ljung-Box）まで確認します。
「点予測だけ見ればよい」ではない：季節予測でも予測区間は必須。季節の波は当たっても水準が外れることがあり、区間で不確実性を示します。