Mímisbrunnr知恵の泉

← 時系列分析 一覧

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

📎 前提：予測の評価指標と時系列CV（時系列CV・ベースライン）・ラグ特徴量と木モデル　|　関連：季節ARIMA（SARIMA）・ETSモデルと状態空間表現

要点（BLUF）

1. バックテストの設計

honest な評価の核心は、「予測時点で実際に使えた情報」だけで予測を作ること。これを破ると（未来の値・全データで作った特徴・テスト込みで決めたハイパラ）、先読みバイアス（look-ahead bias）／情報漏洩で精度が本番より良く見えます。設計の軸は4つ。

flowchart LR
    subgraph WF["ウォークフォワード（時間順・先読みなし）"]
      A["origin t0：[0,t0)で学習 → t0..t0+H を予測"] --> B["origin t0+1：[0,t0+1)で学習 → 予測"]
      B --> C["... 起点を前進 ..."]
    end
    C --> M["各手法の RMSE/MAE/MAPE と 区間スコアを集計"]
    M --> D["ベースライン（素朴・季節素朴）と比較し採否を判断"]

コード①：素朴・季節素朴・SARIMA・ETS・木ML を横並び比較

トレンド＋年次季節＋ノイズの月次系列を、末尾 2424 ヶ月で時間順ホールドアウト。**素朴・季節素朴・SARIMA・ETS・木ML（トレンド除去版）**の5手法で予測し、RMSE/MAE/MAPE を表にします。木MLは ラグ特徴量と木モデル の教訓どおり、線形トレンドを除去してから木で予測し、トレンドを足し戻します。

import warnings; warnings.simplefilter("ignore")
import numpy as np
import pandas as pd
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.exponential_smoothing.ets import ETSModel
from sklearn.ensemble import GradientBoostingRegressor

# 真の構造：線形トレンド + 2調和の年次季節 + ノイズ（月次）
rng = np.random.default_rng(7)
n, period, h = 156, 12, 24
t = np.arange(n)
trend = 0.10 * t
season = 8*np.sin(2*np.pi*t/period) + 4*np.cos(2*np.pi*2*t/period)
y = 50 + trend + season + rng.normal(0, 1.5, n)
train, test = y[:-h], y[-h:]
nt = len(train)

# --- ベースライン：素朴（前回値）・季節素朴（1周期前） ---
fc_naive  = np.repeat(train[-1], h)
fc_snaive = np.array([train[nt - period + (i % period)] for i in range(h)])

# --- SARIMA(1,1,1)(0,1,1,12) ---
sar = SARIMAX(train, order=(1,1,1), seasonal_order=(0,1,1,12), trend="n").fit(disp=False)
fc_sarima = sar.get_forecast(steps=h).predicted_mean

# --- ETS（加法トレンド・加法季節） ---
ets = ETSModel(train, error="add", trend="add", seasonal="add",
               seasonal_periods=period).fit(disp=False)
fc_ets = ets.forecast(h)

# --- 木ML：線形トレンドを除去 → 残差(季節+ノイズ)を木で再帰予測 → トレンドを足し戻す ---
coef = np.polyfit(np.arange(nt), train, 1)         # 訓練でトレンド直線を当てる
trend_fit = np.polyval(coef, np.arange(nt))
detr = train - trend_fit                            # 除トレンド系列（定常的）
lags = [1, 2, 3, 12]
def build(series, lags):
    X = np.column_stack([np.roll(series, L) for L in lags])
    return X[max(lags):], series[max(lags):]
Xd, yd = build(detr, lags)
gbr = GradientBoostingRegressor(n_estimators=300, max_depth=3,
                                learning_rate=0.05, random_state=0).fit(Xd, yd)
hist = list(detr)
det_fc = []
for _ in range(h):
    x = np.array([[hist[-L] for L in lags]])
    p = gbr.predict(x)[0]; det_fc.append(p); hist.append(p)
trend_future = np.polyval(coef, np.arange(nt, nt+h))   # トレンドは外挿（木でなく直線で）
fc_ml = np.array(det_fc) + trend_future

# --- 評価：RMSE / MAE / MAPE を横並び ---
def metrics(p):
    e = p - test
    return (np.sqrt(np.mean(e**2)), np.mean(np.abs(e)), np.mean(np.abs(e/test))*100)

table = {"素朴": fc_naive, "季節素朴": fc_snaive, "SARIMA": fc_sarima,
         "ETS": fc_ets, "木ML(除トレンド)": fc_ml}
rows = {k: metrics(v) for k, v in table.items()}
order = sorted(rows, key=lambda k: rows[k][0])
print(f"{'手法':>16}{'RMSE':>9}{'MAE':>9}{'MAPE(%)':>10}")
for k in order:
    r = rows[k]
    print(f"{k:>16}{r[0]:>9.3f}{r[1]:>9.3f}{r[2]:>10.2f}")
print(f"最良（RMSE）= {order[0]}")

出力：

              手法     RMSE      MAE   MAPE(%)
             ETS    1.437    1.138      1.76
      木ML(除トレンド)    1.458    1.232      1.91
          SARIMA    1.472    1.134      1.74
            季節素朴    3.244    2.661      4.13
              素朴    6.555    5.985      9.44
最良（RMSE）= ETS

出力の意味：「本物の」3手法——ETS 1.4371.437・木ML 1.4581.458・SARIMA 1.4721.472——がほぼ三つ巴で、季節素朴 3.2443.244・素朴 6.5556.555 を大きく引き離しました。勝者が ETS なのは、真の構造が加法トレンド＋加法季節で ETS の想定（ETSモデルと状態空間表現）にぴったり合うから。木ML もトレンドを除去してから当てたので僅差で食らいつきました（生値だと ラグ特徴量と木モデル のとおり外挿できず崩れる）。教訓は2つ：(1)ベースライン（特に季節素朴）を必ず置く——超えて初めて意味がある、(2)勝者はデータの構造で決まるので、複数手法を同じ土俵でバックテストして選ぶ。指標も RMSE だけでなく MAE・MAPE を併記すると（大外れ重視か割合重視か）見え方が変わります。

2. 予測区間の評価：被覆率だけでなく鋭さも

点予測の精度（RMSE 等）に加え、区間予測の質を測らないと不確実性の良し悪しが分かりません。素朴な指標は被覆率（カバレッジ）——95%95\% 区間が実測の何割を捕らえたか。名目 95%95\% に近いほど較正が良い。ただし被覆率は幅を広げれば簡単に上げられるので、それだけでは「無闇に広い区間」を見抜けません。そこで鋭さ（区間幅）も同時に評価する指標を使います。

Lq(y,q^)=max⁡(q (y−q^), (q−1)(y−q^))L_q(y,\hat q) = \max\big(q\,(y-\hat q),\ (q-1)(y-\hat q)\big)

で、これは分位点回帰の損失そのもの。95%95\% 区間なら下限を q=0.025q=0.025、上限を q=0.975q=0.975 で評価し、和を取ります。較正と鋭さを同時に罰する真のスコアリングルール。

Sα=(u−ℓ)+2α(ℓ−y)1{y<ℓ}+2α(y−u)1{y>u}S_\alpha = (u-\ell) + \tfrac{2}{\alpha}(\ell-y)\mathbf 1\{y<\ell\} + \tfrac{2}{\alpha}(y-u)\mathbf 1\{y>u\}

幅 ＋ 区間を外したときの罰。小さいほど良い。狭くて当てれば低スコア、広すぎても外しても高スコアになります。

コード②：被覆率・区間幅・ピンボール・区間スコアで区間の質を測る

同じ系列で SARIMA の 95%95\% 予測区間を点検し、わざと幅を 1.81.8 倍に水増しした区間と比べます。水増しは被覆率を上げますが、鋭さを失うので区間スコアが悪化することを見ます。

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

# トレンド+季節の系列を作り、SARIMAの95%予測区間を点検する
rng = np.random.default_rng(7)
n, period, h = 156, 12, 24
t = np.arange(n)
y = 50 + 0.10*t + 8*np.sin(2*np.pi*t/period) + 4*np.cos(2*np.pi*2*t/period) + rng.normal(0, 1.5, n)
train, test = y[:-h], y[-h:]

sar = SARIMAX(train, order=(1,1,1), seasonal_order=(0,1,1,12), trend="n").fit(disp=False)
fc = sar.get_forecast(steps=h)
ci = fc.conf_int(alpha=0.05)          # 95%区間 [下限, 上限]
lo, hi = ci[:, 0], ci[:, 1]

alpha = 0.05
def evaluate(lo, hi):
    cover = np.mean((test >= lo) & (test <= hi))               # 被覆率
    width = np.mean(hi - lo)                                   # 平均区間幅
    # ピンボール損失（分位点 q）：q=alpha/2（下限）, q=1-alpha/2（上限）
    def pinball(yq, q):
        d = test - yq
        return np.mean(np.maximum(q*d, (q-1)*d))
    pin = pinball(lo, alpha/2) + pinball(hi, 1-alpha/2)
    # 区間スコア（Winkler）：幅 + 区間外への罰（小さいほど良い）
    winkler = np.mean((hi-lo)
                      + (2/alpha)*(lo-test)*(test < lo)
                      + (2/alpha)*(test-hi)*(test > hi))
    return cover, width, pin, winkler

c1, w1, p1, s1 = evaluate(lo, hi)
print(f"SARIMA 95%区間   被覆率={c1*100:.0f}%  平均幅={w1:.2f}  ピンボール={p1:.3f}  区間スコア={s1:.3f}")

# 比較：同じ中心で幅を1.8倍に水増しした区間（被覆は上がるが鋭さを失う）
mid = (lo + hi)/2; half = (hi - lo)/2 * 1.8
lo2, hi2 = mid - half, mid + half
c2, w2, p2, s2 = evaluate(lo2, hi2)
print(f"水増し1.8倍区間   被覆率={c2*100:.0f}%  平均幅={w2:.2f}  ピンボール={p2:.3f}  区間スコア={s2:.3f}")
print(f"被覆は {c2*100:.0f}% > {c1*100:.0f}% と上がるが、区間スコアは {s2:.2f} > {s1:.2f} と悪化（広すぎは減点）")

出力：

SARIMA 95%区間   被覆率=96%  平均幅=5.55  ピンボール=0.162  区間スコア=6.463
水増し1.8倍区間   被覆率=100%  平均幅=9.98  ピンボール=0.250  区間スコア=9.981
被覆は 100% > 96% と上がるが、区間スコアは 9.98 > 6.46 と悪化（広すぎは減点）

出力の意味：SARIMA の 95%95\% 区間は被覆率 96%96\%（名目 95%95\% にほぼ一致＝較正良好）、幅 5.555.55、区間スコア 6.4636.463。これを 1.81.8 倍に水増しすると被覆率は 100%100\% に上がる——けれど幅が 9.989.98 に膨らみ、ピンボール（0.162→0.2500.162\to0.250）も区間スコア（6.46→9.986.46\to9.98）も悪化します。つまり「区間を広げてカバレッジを稼ぐ」のは良い予測ではないと、これらの指標がきちんと見抜く。区間予測は較正（当たる）と鋭さ（狭い）の両立が目標で、被覆率だけ見て満足してはいけない——点予測の RMSE と合わせ、区間スコアまで報告するのが honest な評価です。

3. 数式の直観

4. テキスト全体のまとめ：依存・不確実性・検証の地図

第8章で本テキストは完結します。全8章を貫いてきた問いは一つ——「過去が未来に効く（依存のある）データを、どうモデル化し、不確実性をどう出し、未来をどう検証するか」。各章はこの3つの異なる切り口でした。

flowchart TB
    C1["第1章 基礎：分解・定常性・ACF/PACF・単位根・評価とCV"] --> C2["第2章 ARIMA：差分で定常化し自己相関をAR/MAで"]
    C1 --> C3["第3章 指数平滑/ETS：トレンド・季節へ分解して予測"]
    C2 --> C4["第4章 状態空間/カルマン：隠れ状態を逐次推定（ARIMA/ETSの上位枠）"]
    C3 --> C4
    C2 --> C5["第5章 多変量：VAR・グレンジャー・共和分（系列間の依存）"]
    C2 --> C6["第6章 ボラティリティ：分散の依存をARCH/GARCHで"]
    C4 --> C7["第7章 ベイズ時系列：成分も予測も事後分布で（不確実性を厚く）"]
    C1 --> C8["第8章 機械学習予測：ラグ特徴×木/NN・バックテスト（本章）"]
    C2 --> C8

結びの実務指針：いきなり複雑なモデルに行かず、素朴・季節素朴をベースラインに置き、データの構造（トレンド・季節・自己相関・分散変動・系列間関係）を ACF/分解で見立て、それに合うモデルを複数バックテストで比べる。点予測と**予測区間（被覆＋鋭さ）**を両方報告し、漏洩を徹底排除する。これが、依存のあるデータを予測と不確実性まで含めて扱う作法の全体像です。

⚠️ よくある誤解・落とし穴

関連ノート