Mímisbrunnr知恵の泉

← 時系列分析 一覧

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

📎 前提：指数平滑法（SES・Holt・Holt-Winters）　|　前方：状態空間モデルの枠組み（状態空間の一般形）・関連：ARMA・ARIMAモデル

要点（BLUF）

1. ETS の分類（Error/Trend/Seasonal）

指数平滑の各部品を、組み合わせ方で名前をつけたのが ETS です。

軸記号意味
Error（誤差）N / A / M残差を水準に加える（加法 A）か掛ける（乗法 M）か
Trend（トレンド）N / A / Ad / M / Mdなし／加法／減衰加法（傾きを ϕ<1\phi<1 で減衰）など
Seasonal（季節）N / A / Mなし／加法／乗法（振幅が水準に比例）

たとえば 指数平滑法（SES・Holt・Holt-Winters） で見たモデルは ETS の言葉でこう呼びます：

誤差を A/M、トレンドを N/A/Ad、季節を N/A/M と組み合わせると 30 モデルになり、ETSModel はこの族を統一的に扱います。

2. 状態空間表現：なぜ予測区間が出せるのか

ETS の本質は、指数平滑の漸化式を観測方程式＋状態方程式に書き直すことです。Holt 線形（ETS(A,A,N)）なら

yt=ℓt−1+bt−1+εt（観測方程式）ℓt=ℓt−1+bt−1+α εt（水準の更新）bt=bt−1+β εt（傾きの更新）\begin{aligned} y_t &= \ell_{t-1} + b_{t-1} + \varepsilon_t &&\text{（観測方程式）}\\ \ell_t &= \ell_{t-1} + b_{t-1} + \alpha\,\varepsilon_t &&\text{（水準の更新）}\\ b_t &= b_{t-1} + \beta\,\varepsilon_t &&\text{（傾きの更新）} \end{aligned}

ポイントは誤差 εt\varepsilon_t が観測にも状態にも同じものが入ること（single source of error, SSOE）。第4章のカルマン（状態空間モデルの枠組み）は観測ノイズと状態ノイズを別々に持つ（複数誤差源）のと対照的です。SSOE では εt∼N(0,σ2)\varepsilon_t\sim N(0,\sigma^2) を仮定すると、hh 期先予測の分布が漸化的に解析計算でき、その分散から予測区間が得られます——これが古典指数平滑（点のみ）との決定的な差です。

flowchart LR
    E["誤差 ε_t（正規ノイズ）"] --> Y["観測 y_t = 前回水準 + 傾き + ε_t"]
    E --> ST["状態更新 ℓ_t, b_t（同じ ε_t）"]
    ST --> FC["h期先の予測分布（平均と分散）"]
    FC --> PI["点予測 + 予測区間"]

コード①：ETSで予測区間つき予測（カバレッジも確認）

指数平滑法（SES・Holt・Holt-Winters） と同じ「確率的にドリフトする水準＋加法季節」系列に、ETSModel(error="add", trend="add", seasonal="add") を当てます。get_prediction(...).summary_frame() で mean・pi_lower・pi_upper を取り出し、24ヶ月ホールドアウトで予測区間が実測を捕らえるか（カバレッジ）と、区間幅がホライズンで広がるかを確かめます。

import warnings
warnings.simplefilter("ignore")
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
from statsmodels.tsa.exponential_smoothing.ets import ETSModel

# 真の構造：確率的にドリフトする水準 + 加法季節(周期12,振幅8) + 観測ノイズ
np.random.seed(3)
n, period = 132, 12
mu = np.zeros(n); mu[0] = 50.0
for k in range(1, n):
    mu[k] = mu[k-1] + 0.4 + np.random.normal(0, 0.8)
t = np.arange(n)
y = mu + 8 * np.sin(2 * np.pi * t / period) + np.random.normal(0, 2.0, n)
y = pd.Series(y, index=pd.date_range("2014-01-01", periods=n, freq="MS"))

h = 24
train, test = y.iloc[:-h], y.iloc[-h:]

# ETS(A,A,A)：加法誤差・加法トレンド・加法季節。状態空間なので予測区間が出せる
fit = ETSModel(train, error="add", trend="add", seasonal="add",
               seasonal_periods=period).fit(disp=False)
pred = fit.get_prediction(start=len(train), end=len(train) + h - 1)
frame = pred.summary_frame(alpha=0.05)   # 列: mean, pi_lower, pi_upper

covered = ((test.values >= frame["pi_lower"].values) &
           (test.values <= frame["pi_upper"].values)).sum()
width1 = frame["pi_upper"].iloc[0] - frame["pi_lower"].iloc[0]
width24 = frame["pi_upper"].iloc[-1] - frame["pi_lower"].iloc[-1]
print(f"alpha={fit.smoothing_level:.3f}  beta={fit.smoothing_trend:.3f}  "
      f"gamma={fit.smoothing_seasonal:.3f}")
print(f"1期先 95%予測区間幅  = {width1:.2f}")
print(f"24期先 95%予測区間幅 = {width24:.2f}")
print(f"95%予測区間カバレッジ = {covered}/{h} = {covered/h*100:.0f}%")

plt.figure(figsize=(9, 4.5))
plt.plot(train.index, train.values, color="gray", lw=1, label="訓練データ")
plt.plot(test.index, test.values, color="k", lw=1.2, label="実測（検証）")
plt.plot(frame.index, frame["mean"], color="C1", lw=2, label="ETS点予測")
plt.fill_between(frame.index, frame["pi_lower"], frame["pi_upper"],
                 color="C1", alpha=0.25, label="95%予測区間")
plt.axvline(train.index[-1], ls=":", color="k")
plt.xlabel("年月"); plt.ylabel("値"); plt.legend()
plt.title("ETS(A,A,A) による予測区間つき予測")
plt.tight_layout(); plt.show()

出力：

alpha=0.388  beta=0.000  gamma=0.000
1期先 95%予測区間幅  = 9.05
24期先 95%予測区間幅 = 19.15
95%予測区間カバレッジ = 24/24 = 100%

出力の意味：推定された α,β,γ\alpha,\beta,\gamma は古典 Holt-Winters（指数平滑法（SES・Holt・Holt-Winters））と同じ 0.388,0,00.388,0,0——同じモデルだからです。違いは出力で、ETS は 95% 予測区間を返します。区間幅は 1 期先 9.059.05 → 24 期先 19.1519.15 とホライズンとともに広がり（先ほど不確実）、24 ヶ月すべて区間内（カバレッジ 24/24）。図では橙の帯が右へラッパ状に開きます。点予測だけでなく区間で語れる——これが状態空間化の恩恵です（ARMA・ARIMAモデル の ARIMA 予測区間と同じ思想）。なおホールドアウトが 24 点と少ないため 100% とやや過カバーですが、区間が広がる挙動と捕捉自体は健全です。

コード②：AICで真の構造を選び、状態成分を復元する

ETS 族の中から、トレンド有無・季節の加法/乗法・減衰の有無を変えた候補を AIC で比べ、真の構造（加法トレンド＋加法季節）が選ばれるかを確かめます。さらに選んだモデルの状態成分 fit.states（level/trend/season）を取り出し、仕込んだ真の水準・季節と相関で照合します。

import warnings
warnings.simplefilter("ignore")
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
from statsmodels.tsa.exponential_smoothing.ets import ETSModel

# 真の構造（前コードと同じ）：加法トレンド + 加法季節 が正解
np.random.seed(3)
n, period = 132, 12
mu = np.zeros(n); mu[0] = 50.0
for k in range(1, n):
    mu[k] = mu[k-1] + 0.4 + np.random.normal(0, 0.8)
t = np.arange(n)
seasonal_true = 8 * np.sin(2 * np.pi * t / period)
y = mu + seasonal_true + np.random.normal(0, 2.0, n)
y = pd.Series(y, index=pd.date_range("2014-01-01", periods=n, freq="MS"))

# 4つの ETS 仕様を AIC で比較（trend, seasonal, damped を変える）
specs = {
    "(A,N,A) トレンド無": dict(trend=None, seasonal="add", damped_trend=False),
    "(A,A,N) 季節無":   dict(trend="add", seasonal=None, damped_trend=False),
    "(A,A,A) 加法季節":  dict(trend="add", seasonal="add", damped_trend=False),
    "(A,A,M) 乗法季節":  dict(trend="add", seasonal="mul", damped_trend=False),
    "(A,Ad,A) 減衰トレンド": dict(trend="add", seasonal="add", damped_trend=True),
}
print(f"{'仕様':<18}{'AIC':>10}")
results = {}
for name, kw in specs.items():
    fit = ETSModel(y, error="add", seasonal_periods=period, **kw).fit(disp=False)
    results[name] = fit
    print(f"{name:<18}{fit.aic:>10.2f}")
best = min(results, key=lambda k: results[k].aic)
print(f"AIC最小 = {best}")

# 状態成分を取り出して真の構造と重ねる（best が加法季節モデル）
fit = results["(A,A,A) 加法季節"]
states = fit.states   # DataFrame: level / trend / seasonal
print(f"水準と真のmuの相関   = {np.corrcoef(states['level'], mu)[0,1]:.3f}")
print(f"季節と真の季節の相関 = {np.corrcoef(states['seasonal'], seasonal_true)[0,1]:.3f}")

fig, ax = plt.subplots(3, 1, figsize=(9, 6.5), sharex=True)
ax[0].plot(y.index, states["level"], color="C0", label="ETS推定 level")
ax[0].plot(y.index, mu, color="k", ls="--", lw=1, label="真の水準 mu")
ax[0].set_ylabel("水準"); ax[0].legend(loc="upper left")
ax[1].plot(y.index, states["trend"], color="C2"); ax[1].set_ylabel("トレンド傾き")
ax[2].plot(y.index, states["seasonal"], color="C1", label="ETS推定 season")
ax[2].plot(y.index, seasonal_true, color="k", ls="--", lw=1, label="真の季節")
ax[2].set_ylabel("季節"); ax[2].legend(loc="upper left")
ax[0].set_title("ETS の状態成分（level / trend / season）と真の構造")
ax[-1].set_xlabel("年月"); plt.tight_layout(); plt.show()

出力：

仕様                       AIC
(A,N,A) トレンド無         645.13
(A,A,N) 季節無           770.73
(A,A,A) 加法季節          637.05
(A,A,M) 乗法季節          661.23
(A,Ad,A) 減衰トレンド       643.44
AIC最小 = (A,A,A) 加法季節
水準と真のmuの相関   = 0.997
季節と真の季節の相関 = 0.996

出力の意味：AIC が最小なのは真の構造 (A,A,A) ＝ 637.05。季節を外した (A,A,N) は 770.73770.73 と大きく悪化し（季節が本質）、季節を乗法にした (A,A,M) は 661.23661.23 と罰せられます（データは加法季節なので乗法は不要に複雑）。減衰トレンド (A,Ad,A) は 643.44643.44 と僅差ですが、非減衰の加法トレンドに負けます。AIC は真の構造に合うモデルを選びました（モデル選択の考え方は モデル選択と残差診断）。さらに状態成分は、水準が真の μ\mu と相関 0.9970.997、季節が真の季節と 0.9960.996 でほぼ完全に復元。ETS は「分解しながら予測する」モデルで、states で中身を覗けるのが解釈上の強みです。

3. ETS と ARIMA の関係・使い分け

一部の ETS は ARIMA と等価です（線形・加法誤差のもの）：

一方、乗法のETS（誤差・季節が M）に対応する ARIMA は存在しません。両者の住み分けは：

4. 数式の直観

⚠️ よくある誤解

関連ノート