この記事は「【リアルタイム電力予測】需要・価格・電源最適化ダッシュボード構築記」シリーズの十日目です
Day1はこちら | 全体構成を見る

今回はLightGBMで価格予測を行っていきます！
需要予測と違って変動要因が複雑なため特徴量設計から丁寧に行っていこうと思います。

特徴量設計

EDAをした結果、電力価格には以下の構造が読み取れました。

1. 長期トレンド

燃料価格の急騰 → 電力価格の基調が上昇する
円安 → 輸入燃料が割高になり、基礎的な電力価格が押し上げられる
この構造は、短期変動では説明できない価格の土台を形成している

2. 短期トレンド

気温 → 空調需要 → 価格上昇（U字構造）
需要 → 電力供給逼迫 → 価格上昇
太陽光 → 供給増 → 価格低下
風力 → 深夜の価格安定に寄与
⇒ 短期的には、「需要（気温）× 再エネ（太陽光・風力）」の組み合わせで価格が決まりそう

3. 時間的な周期性

季節・時間帯・曜日パターンが見られた
→ 年間を通じて繰り返される安定した構造がある

ここで特徴量エンジニアリングとしては以下の方針を考えました。

1. 長期トレンド

エネルギー価格や為替は短期的にはノイズが多いですが、電力価格に影響を与えるのはむしろ中長期的な水準です。
そのため「原系列（raw）」「移動平均（マクロトレンド）」「変化率（直近の傾き）」の3つを同時に入れることで、急騰・急落だけでなく安定期の影響もモデルが捉えやすくすることを狙っています。

石炭（月次）

石炭価格（raw）
90日移動平均
90日前との差

石油・ガス・為替（日次）

それぞれの価格（raw）
30日移動平均
30日前との差

2. 短期トレンド

需要や気温は日々の値だけでなく、直近数日の水準が価格に影響すると考えられます。
例えば「急に寒くなった日」と「1週間ずっと寒い週」では発電計画の余裕度が異なり、価格形成のメカニズムも違ってきます。
このため 7 日移動平均を特徴量として加えています。

需要

需要（raw）
過去7日移動平均

気温

気温（raw）
気温（絶対値差分）
過去7日移動平均

価格

24時間ラグ
48時間ラグ
72時間ラグ
1週間ラグ
※過去の値を使用するのでデータリークにはならない

太陽光・風力

太陽光（raw）
風速（raw）

3. 時間的な周期性

月そのものをカテゴリとして入れるより、日番号（day_of_year）や季節を用いた方が、連続性のある時間構造をモデルが捉えやすい傾向を持っていると考えました。なお、比較のために month も一度特徴量候補として入れています。

一年の中の日番号
月
時間
季節（カテゴリカル）
一日の中の時間帯（カテゴリカル）
祝日フラグ（カテゴリカル）
※ 需要予測の時には年と月も入れていましたが、需要予測ほど年や月による特徴が見えづらいのでカテゴリカルな特徴量を作成しています

補足：ラグと移動平均の違い

移動平均: 過去の一定区間の統計量を計算（トレンドや平滑化を捉える）
ラグ: データを時間的にずらす

特徴量作成

以下の関数で特徴量を一気に作成し、実験ごとに選択しながら使用します。

def make_daypart(hour: int) -> str:
    """時間を6つに分ける"""
    if 5 <= hour <= 7:  return "dawn"
    if 8 <= hour <= 11: return "morning"
    if 12 <= hour <= 13:return "noon"
    if 14 <= hour <= 16:return "afternoon"
    if 17 <= hour <= 19:return "evening"
    return "night"

def make_season(month: int) -> str:
    """月を季節ごとに4つに分ける"""
    if month in (3, 4, 5):
        return "spring"
    if month in (6, 7, 8, 9):
        return "summer"
    if month in (10, 11):
        return "autumn"
    return "winter"

def make_features(df: pd.DataFrame, steps_per_day: int = 48) -> pd.DataFrame:
    """
    電力価格予測用の特徴量をまとめて作成する関数
    """
    df = df.copy()

    df["day_of_year"]= df["timestamp"].dt.dayofyear  # 1〜365/366

    # season / daypart
    df["season"]  = df["month"].map(make_season)
    df["daypart"] = df["hour"].map(make_daypart)

    # カテゴリ型にしておく（LightGBM の categorical_feature 用）
    df["season"]    = df["season"].astype("category")
    df["daypart"]   = df["daypart"].astype("category")
    df["is_holiday"]= df["is_holiday"].astype("category")

    df = df.set_index("timestamp")

    # 為替・石油・ガス（日次レベル）
    for col in ["USDJPY", "brent", "henry_hub"]:
        if col not in df.columns:
            continue

        # 30日移動平均
        df[f"{col}_ma30"] = df[col].rolling("30D", min_periods=1).mean()
        # 30日前との差分
        df[f"{col}_chg30"] = df[col] - df[col].shift(freq="30D")

    # 石炭（月次レベル）
    if "coal_aus" in df.columns:

        # 90日移動平均
        df["coal_aus_ma90"] = df["coal_aus"].rolling("90D", min_periods=1).mean()
        # 90日前との差分
        df["coal_aus_chg90"] = df["coal_aus"] - df["coal_aus"].shift(freq="90D")

    # 需要の7日移動平均
    if "demand" in df.columns:
        df["demand_7d_ma"] = df["demand"].rolling("7D", min_periods=1).mean()

    # 気温の特徴量
    if "temperature" in df.columns:
        # 「18.1度からの絶対値差分」や冷暖房負荷指標
        df["temp_abs"] = (df["temperature"] - 18.1).abs()

        # 7日移動平均（気温）
        df["temp_7d_ma"] = df["temperature"].rolling("7D", min_periods=1).mean()

    # 価格ラグ
    price_col = "tokyo_price_jpy_per_kwh"
    if price_col in df.columns:
        # 24, 48, 72時間, 1週間ラグ
        # 30分刻みならそれぞれ 48, 96, 144, 48*7 ステップ
        df["price_lag_24h"] = df[price_col].shift(steps_per_day * 1)
        df["price_lag_48h"] = df[price_col].shift(steps_per_day * 2)
        df["price_lag_72h"] = df[price_col].shift(steps_per_day * 3)
        df["price_lag_1w"]  = df[price_col].shift(steps_per_day * 7)

    df = df.reset_index()

    return df

精度指標

分位点回帰は平均予測ではないのでRMSEや決定係数は通常使用しません。
代わりに、予測分布の評価でよく使われる以下の指標を使用します。

1. PICP (Prediction Interval Coverage Probability)

意味：予測区間 (例: P10–P90) に実際の値が含まれる割合
理想値：設定した信頼区間に近い値（例えば P10–P90なら 80%）
解釈：
- 高すぎる → レンジが広すぎて「安全だが情報量が少ない」
- 低すぎる → レンジが狭すぎて「外れが多い」
- 適切な値 → レンジが現実をよくカバーしている

def coverage(y, ql, qu):  # PICP
    """被覆率"""
    return np.mean((y >= ql) & (y <= qu))

2. MIW (Mean Interval Width)

意味：予測区間の平均幅
解釈：
- 狭い → 予測がシャープだが、カバー率が下がるリスク
- 広い → カバー率は上がるが、情報量が減る
- PICPとのトレードオフを見ながら調整する

def mean_interval_width(ql, qu):
    """平均区間幅"""
    return np.mean(qu - ql)

3. Interval Score

意味：区間の幅とカバー率を同時に評価する指標。小さいほど良い
解釈：
- 幅が狭く、かつカバー率が高いとスコアが良くなる
- 幅が広すぎたり、カバー率が低すぎると悪化する
- PICPとMIWのバランス を数値化したもの

def interval_score(y, ql, qu, alpha=0.2):
    """間隔スコア"""
    # Gneiting & Raftery (2007)
    K = (y < ql).astype(float)
    J = (y > qu).astype(float)
    wid = qu - ql
    penalty = (2/alpha)*((ql - y)*K + (y - qu)*J)
    return np.mean(wid + penalty)

4. Pinball Loss

意味：分位点予測の誤差を評価する指標。分位点ごとに算出
解釈：
- p10 → 下側リスクの予測精度
- p50 → 中央値予測の精度
- p90 → 上側リスクの予測精度
- 値が小さいほどその分位点の予測が正確
- 分位点ごとの強み・弱みを把握できる

def pinball_loss(y, yq, tau):
    diff = y - yq
    return np.mean(np.maximum(tau*diff, (tau-1)*diff))

5. CRPS (Continuous Ranked Probability Score)

意味：予測分布全体と実際の値の差を評価する指標。小さいほど良い
解釈：
- Pinball Lossを全分位点にわたって統合したような指標
- 分布全体の予測精度を一つの数値で表す
- モデル全体の分布予測の質を比較するのに便利

def crps_tri_quantiles(y, q10, q50, q90, taus=(0.1,0.5,0.9)):
    # 重み付きピンボールの和で近似
    losses = [
        pinball_loss(y, q10, taus[0]),
        pinball_loss(y, q50, taus[1]),
        pinball_loss(y, q90, taus[2]),
    ]
    return np.mean(losses)

上記をまとめた関数

def enforce_monotonic(p10, p50, p90):
    """分位の非交差を直すために各時点で昇順に並べる"""
    stacked = np.vstack([p10, p50, p90]).T
    stacked.sort(axis=1)
    return stacked[:,0], stacked[:,1], stacked[:,2]
    
def evaluate(pred_p10, pred_p50, pred_p90, y_test):
    p10, p50, p90 = enforce_monotonic(pred_p10, pred_p50, pred_p90)

    picp = coverage(y_test, p10, p90)             # 例：P10–P90 の被覆率（理想 0.8 付近）
    miw  = mean_interval_width(p10, p90)          # 平均幅
    is80 = interval_score(y_test, p10, p90, 0.2)  # P80のInterval Score
    ql10 = pinball_loss(y_test, p10, 0.1)         # 分位ごとのPinball
    ql50 = pinball_loss(y_test, p50, 0.5)
    ql90 = pinball_loss(y_test, p90, 0.9)
    crps = crps_tri_quantiles(y_test, p10, p50, p90)

    results = {
    "PICP(P10-P90)": picp,
    "MIW(P10-P90)": miw,
    "IntervalScore_80": is80,
    "Pinball_p10": ql10, "Pinball_p50": ql50, "Pinball_p90": ql90,
    "CRPS_approx": crps
    }
    return results

実験

以前電力需給予測でLightGBMを使用した際は、高レベルAPI（scikit-learnラッパー）を使用していました。scikit-learnの他のモデルと同様にfit(X, y) で学習、predict(X) で予測できたり、GridsearchCVが組み込めることが特徴です。

しかし、今回は低レベルAPI（LightGBM本体）を使っていきます。理由は複数の分位点を同時に学習したいので、lgb.trainを順に呼んで分位点ごとのモデルを比較していきます。

実装コード

def get_categorical_cols(df: pd.DataFrame) -> list[str]:
    """
    DataFrame から dtype='category' の列名だけを拾ってリストで返す
    """
    cat_cols = [
        c for c in df.columns
        if pd.api.types.is_categorical_dtype(df[c].dtype)
    ]
    return cat_cols

def run_quantile_lightGBM(X_train, X_test, y_train, y_test):
    
    categorical_cols = get_categorical_cols(X_train)
    print("categorical_cols:", categorical_cols)
    
    train_data = lgb.Dataset(X_train, label=y_train, categorical_feature=categorical_cols,)
    valid_data = lgb.Dataset(X_test, label=y_test, categorical_feature=categorical_cols,)
    
    params_p50 = {
    'objective': 'quantile',
    'alpha': 0.5,
    'metric': 'quantile',
    'learning_rate': 0.05,
    'num_leaves': 64,
    'min_data_in_leaf': 30,
    'feature_fraction': 0.9,
    'seed': 42,
    'force_row_wise': "True"
}
    model_p50 = lgb.train(
        params_p50,
        train_data,
        valid_sets=[valid_data],
        num_boost_round=1000,
    )

    # 上位90％（高値側）
    params_p90 = params_p50.copy()
    params_p90['alpha'] = 0.9
    model_p90 = lgb.train(params_p90, train_data, valid_sets=[valid_data])
    # 下位10％（安値側）
    params_p10 = params_p50.copy()
    params_p10['alpha'] = 0.1
    model_p10 = lgb.train(params_p10, train_data, valid_sets=[valid_data])
    # 推論
    pred_p10 = model_p10.predict(X_test)
    pred_p50 = model_p50.predict(X_test)
    pred_p90 = model_p90.predict(X_test)
    
    return pred_p10, pred_p50, pred_p90, model_p10, model_p50, model_p90

使用する特徴量

feature_cols = [
    # 需要・気象関連
    'demand', 'temperature', 'wind_speed',
    'sunshine_observed', 'solar_radiation',
    'temp_abs', 'demand_7d_ma', 'temp_7d_ma',

    # カレンダー・時間関連
    'is_holiday', 'month', 'hour',
    'day_of_year', 'season', 'daypart',

    # 為替・燃料価格
    'USDJPY', 'brent', 'henry_hub', 'coal_aus',

    # 為替・燃料価格（移動平均・変化率）
    'USDJPY_ma30', 'USDJPY_chg30',
    'brent_ma30', 'brent_chg30',
    'henry_hub_ma30', 'henry_hub_chg30',
    'coal_aus_ma90', 'coal_aus_chg90',

    # 価格ラグ
    'price_lag_24h', 'price_lag_48h',
    'price_lag_72h', 'price_lag_1w'
]

結果

{'PICP(P10-P90)': 0.841828871510931,
 'MIW(P10-P90)': 5.979869014144618,
 'IntervalScore_80': 8.233735669067814,
 'Pinball_p10': 0.3811329615135726,
 'Pinball_p50': 0.8563743893678877,
 'Pinball_p90': 0.4422406053932088,
 'CRPS_approx': 0.5599159854248897}

10%-90%のバンドは84.1％をカバーしています。中央値の予測（Pinball_p50）は 0.856 程度で、まずまずといったところです。（Pinball Loss は値が小さいほど良いです）
Interval Score は「区間の幅」と「カバー率」の両方を同時に評価するため、単独で MIW や PICP を見るよりも総合的な指標になります。
今回は Interval Score が比較的安定しており、分位点間のバランスは悪くないと判断できます。

SHAP では一貫して demand が最も大きく、電力価格の基本構造が「需要ひっ迫度」で大きく左右されていることが分かります
solar_radiation（全天日射）が 3 位に入っており、日中の価格低下を説明する重要な要因となっています
為替の移動平均が燃料価格より上位に来ている点は興味深く、これは為替が燃料コスト全体に広く影響するため、燃料個別よりも説明力が高まった可能性があります
今後方針
- monthとhenry_hub_chg30の貢献度が低いので削除してみます
- is_holidayの貢献度が低いですがリアルタイム予測本番では、現在実測で入っている電力需要も予測値が使用されることから、is_holidayを残しておきます

monthとhenry_hub_chg30を削除

結果

{'PICP(P10-P90)': 0.8525029967463896,
 'MIW(P10-P90)': 6.129476481838045,
 'IntervalScore_80': 8.272437129051259,
 'Pinball_p10': 0.3832302375849058,
 'Pinball_p50': 0.8451305145955873,
 'Pinball_p90': 0.44401347532022,
 'CRPS_approx': 0.5574580758335711}

カバー率が上がりました
10％、90％の予測では Pinball Loss がやや大きくなってしまいましたが、中央値（p50）の損失は小さくなっています

それではこちらのモデルで、テスト期間の実績と予測をplotしてみます。

見えづらいので拡大してみます。
よく当たっていそうなところ（12月）：

こちらを見るとバンドは価格の動きを捉えて予測が収まっています。

外しているところ（5月）：

バンドの範囲外に実測値が存在しています。特にバンドが広くなっているところでも外しているようです。
原因としては、日射量が多いわりに気温が低かった？需要が不安定だった？といった要因が考えられますが、十日目の結果分析にて詳しく見ていこうと思います。

明日

今日はLightGBM単体で予測しましたが、明日はGAMと組み合わせてモデルの精度を上げることを目指します！