ラグ特徴量と木モデル

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

📎 前提：予測の評価指標と時系列CV（ウォークフォワード・ベースライン）　|　数理：勾配ブースティング・特徴量エンジニアリングと前処理（機械学習テキスト）

要点（BLUF）

時系列予測は、ラグ特徴（ $y_{t-1},y_{t-2},\dots$ ）・移動統計（移動平均・移動標準偏差）・カレンダー特徴（月・曜日）を説明変数にすれば、ふつうの教師あり回帰に変換できます。あとは木モデル（勾配ブースティング／ランダムフォレスト）で当てるだけ。
強みは非線形・交互作用・多数の特徴を自動で扱えること。弱みは外挿できないこと——木は区分定数なので、訓練データのレンジ外（将来のトレンド）を予測できません。差分やトレンド除去で対処します。
検証は必ずウォークフォワード（時間順）。特徴は時点 $t$ までの情報だけで作る（移動平均は shift(1) 必須）——少しでも未来を覗くと情報漏洩で評価が幻になります（予測の評価指標と時系列CV）。

1. 時系列を「教師あり回帰」に変換する

ARIMA や状態空間は「過去の値・ショックの線形結合」という構造を仮定して当てました。機械学習予測は逆で、構造を仮定せず、過去から作った特徴量 → 次の値という写像をデータから学びます。鍵は、時系列を行ごとの $(\text{特徴ベクトル},\ \text{目的変数})$ に作り直すこと：

y_t = f\big(\underbrace{y_{t-1},y_{t-2},\dots,y_{t-p}}_{\text{ラグ}},\ \underbrace{\bar y_{t-1:t-k},\ s_{t-1:t-k}}_{\text{移動統計}},\ \underbrace{\text{month}_t,\text{dow}_t}_{\text{カレンダー}}\big) + \varepsilon_t

ラグ特徴 $y_{t-L}$ ：自己相関を捨う。AR と同じ「過去が効く」を、線形に限らず学習。
移動統計：直近の水準・ばらつきを要約。ただし $t$ 時点の予測に使う移動平均は** $t-1$ までで打ち切る**（shift(1)）。当期の $y_t$ を含めたら未来漏洩です。
カレンダー特徴：月・曜日・祝日。季節性を「いつか」という外生情報で直接渡せる。周期は $\sin,\cos$ にすると連続的に扱えます。

$f$ に木モデル（勾配ブースティング・ランダムフォレスト、勾配ブースティング・バギングとランダムフォレスト）を使うのが定番です。木は特徴のスケールに鈍感で、非線形・交互作用・欠測も扱いやすく、特徴量重要度で解釈もできる——表形式データで強いのと同じ理由です。

flowchart LR
    Y["原系列 y_t"] --> L["ラグ特徴 y_t-1, y_t-2, ..."]
    Y --> R["移動統計（shift(1)で打ち切り）"]
    C["カレンダー（月・曜日・祝日）"] --> T["特徴行列 X（時点tまでの情報のみ）"]
    L --> T
    R --> T
    T --> M["木モデル（勾配ブースティング等）で回帰"]
    M --> P["1期先予測（ウォークフォワードで前進検証）"]

コード①：ラグ・カレンダー特徴＋勾配ブースティングのウォークフォワード予測

線形トレンド＋月次季節＋AR(1)残差を仕込んだ月次系列に、ラグ（1,2,3,12）・移動統計・月の $\sin/\cos$ を特徴として作り、GradientBoostingRegressor でウォークフォワード（毎回過去だけで再学習し1期先）予測します。RMSE を**素朴（前回値）・季節素朴（1周期前）**と比べ、特徴量重要度も出します。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
from sklearn.ensemble import GradientBoostingRegressor

# 真の構造：線形トレンド + 月次季節(周期12) + AR(1)残差
rng = np.random.default_rng(0)
n, period = 180, 12
dates = pd.date_range("2010-01-01", periods=n, freq="MS")
t = np.arange(n)
trend = 0.08 * t
season = 6*np.sin(2*np.pi*t/period) + 3*np.cos(2*np.pi*2*t/period)
ar = np.zeros(n)
for i in range(1, n):
    ar[i] = 0.5*ar[i-1] + rng.normal(0, 1.0)
y = 20 + trend + season + ar
s = pd.Series(y, index=dates)

# 特徴量生成：すべて時点 t-1 までの情報だけを使う（未来を見ない）
def make_features(s):
    df = pd.DataFrame({"y": s.values}, index=s.index)
    for L in [1, 2, 3, 12]:
        df[f"lag{L}"] = df["y"].shift(L)            # 過去の値
    df["roll_mean3"] = df["y"].shift(1).rolling(3).mean()   # shift(1)で漏洩防止
    df["roll_std3"]  = df["y"].shift(1).rolling(3).std()
    df["sin_m"] = np.sin(2*np.pi*df.index.month/12)  # カレンダー（月）
    df["cos_m"] = np.cos(2*np.pi*df.index.month/12)
    return df

feat = make_features(s).dropna()
X = feat.drop(columns="y")
yv = feat["y"]

# ウォークフォワード：過去だけで学習 → 1期先を予測 → 1歩進める
start = 120
pred_ml, actual, pred_naive, pred_snaive = [], [], [], []
for i in range(start, len(feat)):
    Xtr, ytr = X.iloc[:i], yv.iloc[:i]
    model = GradientBoostingRegressor(n_estimators=200, max_depth=3,
                                      learning_rate=0.05, random_state=0)
    model.fit(Xtr, ytr)
    pred_ml.append(model.predict(X.iloc[[i]])[0])
    actual.append(yv.iloc[i])
    pred_naive.append(X.iloc[i]["lag1"])     # 素朴：前回値
    pred_snaive.append(X.iloc[i]["lag12"])   # 季節素朴：1周期前

actual = np.array(actual)
rmse = lambda p: np.sqrt(np.mean((np.array(p)-actual)**2))
print(f"1期先 RMSE  木ML(GBR)={rmse(pred_ml):.3f}  素朴={rmse(pred_naive):.3f}  季節素朴={rmse(pred_snaive):.3f}")
print(f"改善率（季節素朴比）= {(1-rmse(pred_ml)/rmse(pred_snaive))*100:.1f}%")

# 特徴量重要度（最終ステップの学習済みモデル）
imp = pd.Series(model.feature_importances_, index=X.columns).sort_values(ascending=False)
print("特徴量重要度（上位5）:")
for k, v in imp.head(5).items():
    print(f"  {k:>10} = {v:.3f}")

plt.figure(figsize=(9, 4))
times = feat.index[start:]
plt.plot(s.index, s.values, color="0.6", lw=1, label="実測")
plt.plot(times, pred_ml, "o-", ms=3, color="C1", label="木ML 1期先予測")
plt.axvline(times[0], ls=":", color="k")
plt.xlabel("年月"); plt.ylabel("y"); plt.legend()
plt.title("ラグ・季節特徴 + 勾配ブースティングのウォークフォワード予測")
plt.tight_layout(); plt.show()

出力：

1期先 RMSE  木ML(GBR)=1.802  素朴=3.282  季節素朴=2.201
改善率（季節素朴比）= 18.1%
特徴量重要度（上位5）:
       lag12 = 0.878
        lag1 = 0.080
       sin_m = 0.012
   roll_std3 = 0.011
       cos_m = 0.006

出力の意味：木ML の1期先 RMSE は $1.802$ で、素朴（ $3.282$ ）・季節素朴（ $2.201$ ）を上回り、季節素朴比で 18.1% 改善——ベースラインを超えたので「価値がある」と言えます（予測の評価指標と時系列CV）。特徴量重要度を見ると lag12（1周期前）が $0.878$ と圧倒的：この系列は季節が支配的なので、木は自動で「1年前の値」を最重要特徴として選びました。lag1（前回値）が次点。手で「季節成分はこう」と指定せずとも、特徴を並べておけば木が必要なものを選ぶのが機械学習予測の身上です。

2. 木モデルはトレンドを外挿できない（最重要の落とし穴）

木モデルの予測は、葉に割り当てられた訓練目的変数の平均です。つまり予測は必ず $[\min y_{\text{train}},\ \max y_{\text{train}}]$ の内側に収まります——訓練レンジの外（上昇トレンドの将来）を出せない。線形回帰や ARIMA が直線・差分で外挿できるのと対照的で、これは木の構造上の限界です。対処は、トレンドを差分や除去で消してから木に渡すこと。差分系列 $\Delta y_t=y_t-y_{t-1}$ は定常（トレンドが消える）なので木の内挿で扱え、予測後に累積して水準へ戻します（ARMA・ARIMAモデルの I＝差分と同じ発想）。

コード②：木はトレンドを外挿できない → 差分で対処

強い上昇トレンドを持つ系列で、後半（テスト期）は訓練レンジの外になります。(A) 生値でラグ回帰した木を再帰予測（予測を入力に回して多段へ）し、(B) 1階差分で学習した木の予測を累積して水準へ戻します。両者の RMSE と、生値の木の予測が訓練最大値に頭打ちする様子を見ます。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
from sklearn.ensemble import RandomForestRegressor

# 真の構造：強い上昇トレンド + ノイズ（テスト期は訓練レンジの外＝外挿領域）
rng = np.random.default_rng(1)
n, split = 200, 130
t = np.arange(n)
y = 10 + 0.5*t + rng.normal(0, 2.0, n)
lags = [1, 2, 3]

def build(series, lags):
    X = np.column_stack([np.roll(series, L) for L in lags])
    return X[max(lags):], series[max(lags):]   # 先頭のラグ未定義を捨てる

def recursive_forecast(model, history, lags, steps):
    hist = list(history)
    out = []
    for _ in range(steps):
        x = np.array([[hist[-L] for L in lags]])   # 直近のラグから特徴を作る
        p = model.predict(x)[0]
        out.append(p); hist.append(p)              # 予測を履歴に足して前進
    return np.array(out)

# (A) 生値で学習 → 再帰予測（木は訓練の最大値より上を出せない）
Xr, yr = build(y, lags)
mr = RandomForestRegressor(n_estimators=300, random_state=0).fit(Xr[:split-max(lags)], yr[:split-max(lags)])
fc_raw = recursive_forecast(mr, y[:split], lags, n-split)

# (B) 1階差分で学習 → 差分を再帰予測 → 累積して水準へ戻す（トレンドに追随）
dy = np.diff(y)                                  # Δy は定常（トレンドが消える）
Xd, yd = build(dy, lags)
md = RandomForestRegressor(n_estimators=300, random_state=0).fit(Xd[:split-1-max(lags)], yd[:split-1-max(lags)])
dfc = recursive_forecast(md, dy[:split-1], lags, n-split)
fc_diff = y[split-1] + np.cumsum(dfc)            # 逆差分で水準を復元

actual = y[split:]
rmse = lambda p: np.sqrt(np.mean((p-actual)**2))
print(f"訓練データの最大 y       = {y[:split].max():.1f}")
print(f"生値・木の予測の最大値   = {fc_raw.max():.1f}  -> 訓練レンジに頭打ち（外挿できない）")
print(f"再帰予測 RMSE  生値の木 = {rmse(fc_raw):.2f}   差分の木 = {rmse(fc_diff):.2f}")

plt.figure(figsize=(9, 4.5))
plt.plot(t, y, color="0.6", lw=1, label="実測")
plt.plot(t[split:], fc_raw, color="C3", lw=2, label="生値の木（頭打ち＝外挿不可）")
plt.plot(t[split:], fc_diff, color="C0", lw=2, label="差分の木（トレンド追随）")
plt.axvline(split, ls=":", color="k"); plt.xlabel("時点 t"); plt.ylabel("y"); plt.legend()
plt.title("木モデルはトレンドを外挿できない：差分で対処")
plt.tight_layout(); plt.show()

出力：

訓練データの最大 y       = 77.1
生値・木の予測の最大値   = 72.4  -> 訓練レンジに頭打ち（外挿できない）
再帰予測 RMSE  生値の木 = 22.40   差分の木 = 2.42

出力の意味：生値の木の予測は最大 $72.4$ で止まり、訓練最大値 $77.1$ を超えられません——上昇を続ける実測に置いていかれ、RMSE $22.40$ と大破綻。一方、差分で学習した木は $\Delta y$ （定常）を予測して累積するので、トレンドに追随し RMSE $2.42$ ——約9倍の改善です。図では赤線（生値）が水平に張り付き、青線（差分）が実測に沿います。教訓は明快：トレンドのある系列に木を生値で当ててはいけない。必ず差分かトレンド除去で定常化してから渡します（ARIMA の I、回帰のトレンド項と同じ役割）。

3. 多段先予測：再帰（recursive）vs 直接（direct）

$h$ 期先まで一気に予測したいとき、戦略は2つ。

再帰（recursive）：1期先モデルを1つ学習し、予測値を入力に戻して $h$ 回繰り返す。モデルは1個で済むが、予測の誤差が次の入力に伝播して積み上がる。
直接（direct）：ホライズン $h$ ごとに別モデルを学習（ $y_{t+h}$ を特徴から直接当てる）。誤差は伝播しないが、 $h$ 個のモデルが要り、各 $h$ で独立に学ぶので長期の滑らかさを失いやすい。

経験則は「短期は再帰、長期は直接」。1期先（ $h=1$ ）は同じ学習目標なので再帰が有利、 $h$ が伸びると誤差伝播のない直接が追い越します（中間の DirRec・Rectify 等のハイブリッドもあります。要最新確認）。

コード③：再帰 vs 直接をホライズン別 RMSE で比較

定常な AR(2)（トレンドなし＝木の内挿で扱える）で、 $h=1..8$ の多段予測を再帰と直接で行い、複数の予測起点（origin）で平均したホライズン別 RMSEを比べます。

import numpy as np
from sklearn.ensemble import RandomForestRegressor

# 真の構造：定常な AR(2)（トレンド無し＝木の内挿で扱える）。多段予測の戦略だけを比較
rng = np.random.default_rng(2)
n = 400
y = np.zeros(n)
for i in range(2, n):
    y[i] = 0.6*y[i-1] - 0.3*y[i-2] + rng.normal(0, 1.0)

lags = [1, 2, 3, 4, 5, 6]
H = 8                       # 予測ホライズン
ml = max(lags)

def feat_at(series, idx, lags):
    return np.array([series[idx-L] for L in lags])

# 訓練/検証を時間順に分割（origin はすべて訓練の後）
train_end = 320
origins = range(train_end, n-H)

# --- 再帰戦略：1期先モデルを1つ学習し、予測を入力に回して多段へ ---
Xr = np.array([feat_at(y, i, lags) for i in range(ml, train_end)])
yr = np.array([y[i] for i in range(ml, train_end)])
rec = RandomForestRegressor(n_estimators=300, random_state=0).fit(Xr, yr)

# --- 直接戦略：ホライズン h ごとに別モデル（y_{t+h} を直接予測） ---
direct = {}
for h in range(1, H+1):
    Xh = np.array([feat_at(y, i, lags) for i in range(ml, train_end-h)])
    yh = np.array([y[i+h] for i in range(ml, train_end-h)])
    direct[h] = RandomForestRegressor(n_estimators=300, random_state=0).fit(Xh, yh)

err_rec = np.zeros(H); err_dir = np.zeros(H); cnt = 0
for o in origins:
    base = feat_at(y, o, lags)
    # 再帰：1歩ずつ予測を履歴に積む
    hist = list(y[:o+1])
    for h in range(1, H+1):
        x = np.array([[hist[-L] for L in lags]])
        p = rec.predict(x)[0]; hist.append(p)
        err_rec[h-1] += (p - y[o+h])**2
    # 直接：各 h 専用モデルで一発予測（誤差は伝播しない）
    for h in range(1, H+1):
        p = direct[h].predict(base.reshape(1, -1))[0]
        err_dir[h-1] += (p - y[o+h])**2
    cnt += 1

rmse_rec = np.sqrt(err_rec/cnt); rmse_dir = np.sqrt(err_dir/cnt)
print(f"検証 origin 数 = {cnt}")
print(f"{'h':>3}{'再帰RMSE':>12}{'直接RMSE':>12}")
for h in range(1, H+1):
    print(f"{h:>3}{rmse_rec[h-1]:>12.3f}{rmse_dir[h-1]:>12.3f}")
print(f"平均RMSE  再帰={rmse_rec.mean():.3f}  直接={rmse_dir.mean():.3f}")

出力：

検証 origin 数 = 72
  h      再帰RMSE      直接RMSE
  1       1.042       1.150
  2       1.256       1.145
  3       1.268       1.135
  4       1.238       1.176
  5       1.258       1.182
  6       1.180       1.167
  7       1.235       1.166
  8       1.259       1.167
平均RMSE  再帰=1.217  直接=1.161

出力の意味： $h=1$ では再帰が勝つ（ $1.042$ vs $1.150$ ）——1期先は再帰モデルそのものの学習目標だから。ところが $h\ge2$ では直接が一貫して有利で、平均 RMSE も直接 $1.161$ < 再帰 $1.217$ 。再帰は予測誤差が入力に積み上がる（ $h$ が伸びるほど効く）のに対し、直接は各ホライズンを独立に当てるので伝播がありません。実務では「短期は再帰で手早く、長期や重要なホライズンは直接（やハイブリッド）で」と使い分けます。なお、データの性質（強い自己相関・滑らかなトレンド）次第で勝敗は変わるので、両方をバックテストで比べるのが正解です（バックテストと予測の評価）。

4. 数式の直観

教師あり回帰への変換＝自己回帰の一般化：AR( $p$ ) は $y_t=\sum\phi_i y_{t-i}+\varepsilon_t$ という線形写像。ラグ特徴＋木は同じラグを入力にしつつ、写像を非線形・交互作用込みに拡張したもの。線形で十分なら ARIMA が効率的、非線形・多特徴・閾値構造があるなら木が活きます。
木の予測は訓練ラベルの凸包の内側：葉の値は訓練目的変数の平均なので、出力は $[\min,\max]$ に閉じる。だからトレンドは外挿不能——差分（ARMA・ARIMAモデル）やトレンド除去（STL分解）で「定常な部分」だけを木に任せ、トレンドは別途（直線・ARIMA）で外挿します。
再帰の誤差伝播：1期先誤差 $e_1$ が次の入力に入り、 $h$ 期先には $e_1,\dots,e_{h-1}$ が非線形に効いて分散が膨らむ。直接はこれを断つ代わりに、ホライズンごとに学習データが減り（ $y_{t+h}$ の組が少ない）バリアンスが増える——バイアス・バリアンスのトレードオフ（汎化と過学習・バイアスバリアンス分解、機械学習テキスト）です。

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

「木は何でも当たる」ではない：木は外挿できません。トレンドのある系列に生値で当てると将来に頭打ち（コード②で RMSE $22.40$ ）。必ず差分・トレンド除去で定常化を。
「移動平均特徴は当期を含めてよい」ではない： $t$ 時点の予測に使う移動平均が $y_t$ を含むと未来漏洩。shift(1) で $t-1$ までに打ち切る。ターゲットエンコーディングや将来カレンダー以外の「未来由来の特徴」も厳禁です（予測の評価指標と時系列CV）。
「k-fold でチューニングしてよい」ではない：時系列の交差検証はシャッフル禁止。TimeSeriesSplit／ウォークフォワードで時間順に（訓練・検証・テストと交差検証、機械学習テキスト）。普通の k-fold は未来で過去を当てる漏洩になります。
「特徴量重要度＝因果」ではない：重要度は「予測にどれだけ使われたか」で、相関・交絡を含む。介入効果は別問題（グレンジャー因果や因果推論テキストへ）。
「多段予測は再帰一択」ではない：長期は直接やハイブリッドが有利なことが多い（コード③）。ホライズンと目的に応じて選び、必ずバックテストで確認します。