Mímisbrunnr知恵の泉

← 時系列分析 一覧

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

📎 前提：時系列データと分解　|　関連：ARMA・ARIMAモデル（残差に当てる予測モデル）・予測の評価指標と時系列CV

要点（BLUF）

1. STL とは：Loess による分解

時系列データと分解 の古典的分解は、トレンドを中心化移動平均で、季節を月ごとの平均で推定しました。シンプルですが、(1) 1つの外れ値が移動平均を引きずる、(2) 季節は全期間で固定、(3) 端でトレンドが欠ける、という弱点があります。

STL はトレンド・季節を **Loess（locally estimated scatterplot smoothing、局所重み付き回帰）**で推定します。Loess は各点の近傍に重み（距離で減衰）をかけて低次多項式を当てる滑らかなフィットで、STL は

  1. 季節：同じ周期位置の部分系列（例：すべての1月、すべての2月…）を Loess で平滑化 → 季節成分。これにより季節が年々ゆっくり変わってよい。
  2. トレンド：季節を抜いた系列を Loess で平滑化 → トレンド成分。
  3. 上を反復し、ロバスト重み（残差が大きい点を軽くする）を更新 → 外れ値の影響を抑える。

加法分解 yt=Tt+St+Rty_t = T_t + S_t + R_t を前提に、トレンド窓・季節窓・ロバスト反復で滑らかさと頑健さを制御できるのが STL です。

2. seasonal_decompose との違い（外れ値ロバスト性）

コード①：外れ値1点で古典分解とSTLを比べる

トレンド＋加法季節（周期12）に、ある1点だけ大きな外れ値（+18）を仕込みます。seasonal_decompose（移動平均）と STL(robust=True) の両方で分解し、(a) 復元トレンドが真のトレンドにどれだけ近いか（RMSE）、(b) 外れ値がどれだけ残差に落ちたか（トレンド/季節を汚していないか）を比べます。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
from statsmodels.tsa.seasonal import STL, seasonal_decompose

# 真の構造：トレンド + 加法季節(周期12) + ノイズ。さらに1点に大きな外れ値
np.random.seed(8)
n, period = 144, 12
t = np.arange(n)
trend_true = 0.1 * t
seasonal_true = 5 * np.sin(2 * np.pi * t / period)
y = 20 + trend_true + seasonal_true + np.random.normal(0, 0.6, n)
io = 80                      # 外れ値を入れる時点
y[io] += 18.0                # 大きなスパイク（汚染）
y = pd.Series(y, index=pd.date_range("2010-01-01", periods=n, freq="MS"))

# 古典分解（移動平均ベース）と STL（ロバスト）の両方で分解
cls = seasonal_decompose(y, period=period, model="additive")
stl = STL(y, period=period, robust=True).fit()

# トレンド復元の正確さ（外れ値の影響を受けにくいのはどちらか）
mask = ~np.isnan(cls.trend.values)
rmse_trend_cls = np.sqrt(np.mean((cls.trend.values[mask] - (20+trend_true)[mask])**2))
rmse_trend_stl = np.sqrt(np.mean((stl.trend.values[mask] - (20+trend_true)[mask])**2))
# 外れ値はどれだけ残差に落ちたか（大きいほど trend/season を汚していない）
print(f"トレンドRMSE  古典={rmse_trend_cls:.3f}  STLロバスト={rmse_trend_stl:.3f}")
print(f"季節相関(vs真) 古典={np.corrcoef(cls.seasonal.values, seasonal_true)[0,1]:.3f}  "
      f"STL={np.corrcoef(stl.seasonal.values, seasonal_true)[0,1]:.3f}")
print(f"外れ値点の残差 古典={cls.resid.values[io]:.2f}  STL={stl.resid.values[io]:.2f}（注入=+18）")

fig, ax = plt.subplots(3, 1, figsize=(9, 6.5), sharex=True)
ax[0].plot(y.index, y.values, color="gray", lw=1); ax[0].axvline(y.index[io], ls=":", color="r")
ax[0].set_ylabel("原系列"); ax[0].set_title("外れ値1点を含む系列の分解：古典 vs STL(ロバスト)")
ax[1].plot(y.index, cls.trend.values, color="C0", label="古典トレンド")
ax[1].plot(y.index, stl.trend.values, color="C1", label="STLトレンド")
ax[1].plot(y.index, 20+trend_true, color="k", ls="--", lw=1, label="真のトレンド")
ax[1].axvline(y.index[io], ls=":", color="r"); ax[1].set_ylabel("トレンド"); ax[1].legend(loc="upper left")
ax[2].plot(y.index, cls.seasonal.values, color="C0", label="古典季節")
ax[2].plot(y.index, stl.seasonal.values, color="C1", label="STL季節")
ax[2].set_ylabel("季節"); ax[2].legend(loc="upper left"); ax[2].set_xlabel("年月")
plt.tight_layout(); plt.show()

出力：

トレンドRMSE  古典=0.485  STLロバスト=0.140
季節相関(vs真) 古典=0.990  STL=0.993
外れ値点の残差 古典=15.14  STL=19.08（注入=+18）

出力の意味：外れ値があると、古典分解のトレンド RMSE は 0.4850.485 ですが、ロバスト STL は 0.1400.140 と約3.5倍正確です。決め手は残差の行：外れ値点の残差が 古典 15.1415.14／STL 19.0819.08。STL は注入した +18 をほぼ丸ごと残差に落とし（トレンド・季節を汚さない）、古典は外れ値の一部（約3）をトレンド/季節に吸わせてしまっています。図でも、古典トレンド（青）が外れ値の位置（赤点線）で真の直線から盛り上がるのに対し、STL（橙）は真値（破線）に張り付いたまま。外れ値や構造変化を含む実データでは STL のロバスト性が効きます。

Loess は局所回帰（重み付き最小二乗）なので、回帰の重み付き推定の発想を時系列の平滑化に応用したものと見られます。STL の seasonal・trend 窓（奇数）で滑らかさを、robust=True で外れ値耐性を制御します。

3. STLで予測する：季節調整 → 残差にモデル

STL 自体は分解の道具で予測はしません。実務の定石は「STL で季節成分を抜き（季節調整）→ 季節調整後の系列に ARIMA など → 予測に季節を戻す」。statsmodels の STLForecast がこれを一括でやり、予測区間も出します。

flowchart LR
    A["原系列 y_t（トレンド+季節）"] --> B["STLで季節 S_t を推定・除去"]
    B --> C["季節調整後の系列に ARIMA を当てる"]
    C --> D["季節調整系列を予測 + 予測区間"]
    D --> E["季節 S_t を戻す"]
    E --> F["点予測 + 予測区間"]

コード②：STLForecast（STL＋ARIMA）で予測区間つき予測

トレンド＋季節系列を24ヶ月ホールドアウトし、STLForecast(train, ARIMA, model_kwargs=..., period=12, robust=True) で予測します。get_prediction(...).summary_frame() で mean・mean_ci_lower・mean_ci_upper を取り出し、季節素朴予測と RMSE を比べ、区間カバレッジも見ます。

import warnings
warnings.simplefilter("ignore")
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
from statsmodels.tsa.forecasting.stl import STLForecast
from statsmodels.tsa.arima.model import ARIMA

# 真の構造：トレンド + 加法季節(周期12) + ノイズ
np.random.seed(8)
n, period = 144, 12
t = np.arange(n)
y = 20 + 0.1 * t + 5 * np.sin(2 * np.pi * t / period) + np.random.normal(0, 0.6, n)
y = pd.Series(y, index=pd.date_range("2010-01-01", periods=n, freq="MS"))

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

# STL で季節調整 → 季節調整後の系列に ARIMA を当てて予測（区間つき）
stlf = STLForecast(train, ARIMA,
                   model_kwargs=dict(order=(2, 1, 0), trend="t"),
                   period=period, robust=True).fit()
pred = stlf.get_prediction(start=len(train), end=len(train) + h - 1)
fr = pred.summary_frame()   # 列: mean, mean_se, mean_ci_lower, mean_ci_upper

# 季節素朴予測（12ヶ月前の値）と比較
fc_naive = np.resize(train.iloc[-period:].values, h)
rmse_stl = np.sqrt(np.mean((fr["mean"].values - test.values) ** 2))
rmse_naive = np.sqrt(np.mean((fc_naive - test.values) ** 2))
covered = ((test.values >= fr["mean_ci_lower"].values) &
           (test.values <= fr["mean_ci_upper"].values)).sum()
print(f"STL+ARIMA RMSE  = {rmse_stl:.3f}")
print(f"季節素朴予測 RMSE = {rmse_naive:.3f}")
print(f"改善率 = {(1 - rmse_stl/rmse_naive)*100:.1f}%（季節素朴比）")
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(fr.index, fr["mean"], color="C1", lw=2, label="STL+ARIMA点予測")
plt.fill_between(fr.index, fr["mean_ci_lower"], fr["mean_ci_upper"],
                 color="C1", alpha=0.25, label="95%予測区間")
plt.axvline(train.index[-1], ls=":", color="k")
plt.xlabel("年月"); plt.ylabel("値"); plt.legend()
plt.title("STLForecast（STLで季節調整 + ARIMA）による予測区間つき予測")
plt.tight_layout(); plt.show()

出力：

STL+ARIMA RMSE  = 0.886
季節素朴予測 RMSE = 2.084
改善率 = 57.5%（季節素朴比）
95%予測区間カバレッジ = 22/24 = 92%

出力の意味：STL で季節を抜いてから ARIMA を当てた予測の RMSE は 0.8860.886 で、季節素朴予測 2.0842.084 を 57.5% 改善。95% 予測区間は 24 ヶ月中 22 ヶ月（92%）を捕らえ、名目 95% とよく整合（較正されている）。STL が季節を担い、ARIMA が季節調整後のトレンド・自己相関を担う役割分担で、解釈と精度を両立できます。なおここでの区間は「季節調整後の系列に対する ARIMA の予測区間」で、戻す季節成分は確定的に扱われる点に注意（季節の推定誤差ぶんは含みません）。

4. 数式の直観

⚠️ よくある誤解

関連ノート