0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

失業率の予測を「前月の値」と比べる実装:公表タイミングを守るローリング評価と Diebold–Mariano 検定

0
Posted at

完全失業率(季節調整値)の1か月先・3か月先を、ARIMA・SARIMAX・ElasticNet・LightGBM で予測し、「前月の値のまま」と比べた分析の実装をまとめます。結果(1か月先はどのモデルも前月の値に勝てず、3か月先は ElasticNet が誤差を15%減らした)は note に書いたので、この記事は「どう実装したか」に絞ります。

確かめたかったのは次の1点です。

  • いちばん単純な予測「前月の値のまま」に、モデルが本当に勝てるか(そのとき使えた情報だけで)

コードはすべて GitHub の articles/F13_unemployment_forecast/src/ にあるものです。

全体の流れ

ファイル 中身
01_取得.py e-Stat と内閣府からの取得(同梱データを使うときは何もしない)
02_整形.py 労働力調査・有効求人倍率・景気動向指数を月次の1枚の表 data/monthly.csv にする
03_分析.py 特徴量づくり、ローリング予測、評価指標と DM 検定、重要度
03b_頑健性.py 学習期間・特徴量・リーク・原数値を変えた確かめ
04_図.py 図6枚
pip install -r requirements.txt
cd src
python 01_取得.py          # 同梱データを使う(取り直すときは --download)
python 02_整形.py          # data/monthly.csv
python 03_分析.py          # results/forecasts.csv, metrics.csv, importance.csv(約1分)
python 03b_頑健性.py       # results/robustness.csv(約2分)
python 04_図.py            # figures/

乱数の種は 42 に固定しています。

1. データの取り方と整形

使ったのは3つです。

  • 総務省「労働力調査」長期時系列 表1 a-1(完全失業率の季節調整値・原数値)
  • 厚生労働省「一般職業紹介状況」第3表(有効求人倍率、パートを含む一般)
  • 内閣府「景気動向指数」(CI 先行・一致、鉱工業生産指数、労働投入量指数、新規求人数、消費者態度指数)

data/raw/ には 2026-10-26 に取得したファイルを置いています。有効求人倍率の statInfId(e-Stat の表の ID)は毎月変わり、景気動向指数のファイル名の先頭も公表日で変わるので、取り直すときは ID とファイル名を最新に置き換えます。

労働力調査の Excel は「年が書いてある行」と「月の行」が分かれている

長期時系列の Excel は、年が書かれた行と月の行が分かれていて、月の列には「*」の付いた文字が混じります。src/02_整形.py の read_lfs では、年を覚えておきながら行を読み、12月を読んだら年を1つ進めます。値が < > や括弧の付いた文字で入っているセルもあるので、記号を外してから数値にしています。

3つの出所を月の日付(各月1日)でそろえ、outer join で1枚にしています。

2. 特徴量:公表のタイミングをずらし方で表す

この分析でいちばん気をつけたのが、「予測を出す日に、その値はもう公表されていたか」です。予測を出す日を「T月の失業率が公表された日」(T+1月の月末)とすると、

  • 失業率・有効求人倍率・新規求人数・消費者態度指数:T月まで使える(失業率と求人倍率は同じ日に公表)
  • CI・鉱工業生産・労働投入量:T−1月までしか使えない

となります。これを shift の幅(lag_ci=1, lag_job=0)で表しました。あとで「わざとリークさせる」確かめをするときは、この引数を -1 にするだけで済むようにしています。

def make_features(d, lag_ci=1, lag_job=0):
    """T月の行に、T月の公表時点で使える特徴量を並べる(lag_job=-1 は『まだ公表されていない T+1 月の求人』を使う=リーク)"""
    f = pd.DataFrame(index=d.index)
    u = d["u_sa"]
    du = u.diff()
    for l in range(6):
        f[f"失業率の変化_{l+1}か月前"] = du.shift(l)
    f["失業率の水準"] = u
    f["失業率_12か月平均との差"] = u - u.rolling(12).mean()
    j = d["jobs_sa"].shift(lag_job)
    no = np.log(d["new_offers"]).shift(lag_job)
    f["有効求人倍率"] = j
    for k in (1, 3, 6):
        f[f"有効求人倍率の変化_{k}か月"] = j - j.shift(k)
    for k in (1, 3, 12):
        f[f"新規求人数の伸び_{k}か月"] = no - no.shift(k)
    cc = d["consumer_conf"].shift(lag_job)
    f["消費者態度指数の変化_3か月"] = cc - cc.shift(3)
    for col, name in (("ci_lead", "CI先行"), ("ci_coin", "CI一致")):
        s = d[col].shift(lag_ci)
        f[f"{name}の変化_1か月"] = s - s.shift(1)
        f[f"{name}の変化_3か月"] = s - s.shift(3)
    for col, name in (("iip", "鉱工業生産"), ("labor_input", "労働投入量")):
        s = np.log(d[col]).shift(lag_ci)
        f[f"{name}の伸び_3か月"] = s - s.shift(3)
    return f

目的変数は水準ではなく変化幅 u(T+h) − u(T) にしました。こうすると「前月の値」はいつも 0 を予測するモデルになり、比べやすくなります。

def targets(u):
    return {h: u.shift(-h) - u for h in H}

3. ローリング予測の実装

予測の対象月が 2005年1月〜2026年8月になるように、予測を出す月 T を1か月ずつ進めます。学習は1985年7月から T までの拡張窓です。ARIMA 系は毎月、ElasticNet と LightGBM は毎年1月に学習し直し、それ以外の月はキャッシュしたモデルで予測します。

落とし穴は学習データの切り方です。h か月先の目的変数 u(T'+h) が分かっているのは T' ≤ T − h の行だけなので、学習に使う行を T - h か月 までに絞っています。

        for kind in [m for m in models if m in ("ElasticNet", "LightGBM")]:
            for h in H:
                key = (kind, h)
                if key not in cache or T.month == refit_month:
                    # 学習に使えるのは、目的変数(T'+h)が T までに分かっている行だけ
                    lo = TRAIN_START if window is None else T - pd.DateOffset(years=window)
                    idx = feats.index[(feats.index >= lo) & (feats.index <= T - pd.DateOffset(months=h))]
                    tr = pd.concat([feats.loc[idx], Y[h].loc[idx].rename("y")], axis=1).dropna()
                    cache[key] = fit_ml(kind, tr.drop(columns="y"), tr["y"])
                x = feats.loc[[T]]
                rec[f"{kind}_{h}"] = float(cache[key].predict(x)[0]) if x.notna().all(axis=None) else np.nan

ElasticNet の正則化の強さは、TimeSeriesSplit を使った交差検証で決めます(ふつうの KFold だと未来のデータで過去を検証してしまうため)。

def fit_ml(kind, X, y):
    if kind == "ElasticNet":
        m = make_pipeline(StandardScaler(), ElasticNetCV(l1_ratio=[.2, .5, .8, 1.0], cv=TimeSeriesSplit(5),
                                                         alphas=30, max_iter=20000, random_state=SEED))
    else:
        m = lgb.LGBMRegressor(n_estimators=300, learning_rate=0.03, num_leaves=15, min_child_samples=20,
                              subsample=0.8, subsample_freq=1, colsample_bytree=0.8, reg_lambda=1.0,
                              random_state=SEED, verbose=-1, n_jobs=2)
    return m.fit(X, y)

SARIMAX は ARIMA(2,1,0) に、1か月前の有効求人倍率と新規求人数の変化を外生変数として足しました。3か月先では外生変数の先の値が分からないので「変化なし(0)」と置いています。

4. 評価:MAE・向き・Diebold–Mariano 検定

評価は平均絶対誤差(MAE)、RMSE、「上がる・横ばい・下がる」の向き(±0.05ポイント未満を横ばい)の一致率です。

「前月の値より良い」が偶然でないかは、Diebold–Mariano 検定(2つの予測の誤差の差が0と言えるかを見る検定)で確かめます。h か月先の予測は誤差が重なり合うので、自己共分散を h−1 次まで足した分散(HAC)を使い、Harvey らの小標本補正をかけています。

def dm_test(e1, e2, h):
    """Diebold–Mariano(絶対誤差、HAC、Harvey らの小標本補正)。e1 が e2 より良ければ統計量は負"""
    dlt = np.abs(e1) - np.abs(e2)
    dlt = dlt[~np.isnan(dlt)]
    n = len(dlt)
    g = [np.mean((dlt[k:] - dlt.mean()) * (dlt[:n - k] - dlt.mean())) for k in range(h)]
    v = (g[0] + 2 * sum(g[1:])) / n
    stat = dlt.mean() / np.sqrt(v)
    stat *= np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)
    p = 2 * stats.t.sf(abs(stat), n - 1)
    return stat, p

結果(results/metrics.csv の全期間)は次のとおりです。

1か月先 MAE 前月の値比 3か月先 MAE 前月の値比
前月の値(ナイーブ) 0.083 1.00 0.128 1.00
AR(失業率だけ) 0.085 1.03 0.130 1.02
SARIMAX 0.086 1.04 0.131 1.02
ElasticNet 0.083 1.00 0.109 0.85(p=0.01)
LightGBM 0.088 1.05 0.114 0.89(p=0.09)
平均(3モデル) 0.083 1.00 0.111 0.87(p=0.002)

前月の値と比べた誤差

5. 図の出し方

予測は変化幅で持っているので、図にするときは u(T) に足して水準に戻し、日付を h か月後ろにずらします。

def levels(h, m):
    s = fc["u_T"] + fc[f"{m}_{h}"]
    s.index = s.index + pd.DateOffset(months=h)
    return s

3か月先の予測と実績

重要度は、最後まで学習した LightGBM の SHAP 値(各特徴量が予測をどれだけ動かしたか)の絶対値の平均と、ElasticNet の標準化した係数(e[-1].coef_)を並べました。

重要度

6. つまずいた点と落とし穴

リークを「特徴量を作る関数の引数」で切り替えられるようにした

公表前の値を混ぜると、実際には出せない予測で精度を測ることになります。make_features(df, lag_ci=-1, lag_job=-1) でわざと T+1 月の求人・景気の値を使わせて比べたところ、3か月先の ElasticNet は 0.85 → 0.85 とほとんど変わりませんでした。求人の情報は「その月の値」より「ここ数か月の流れ」として効いているようです。

原数値で比べると、モデルが楽に勝てて見える

季節調整前の原数値で予測すると、前月の値の MAE は 0.213 と大きく、ElasticNet は前月の値比 0.62 と大きく勝てて見えます。ただし ElasticNet の誤差 0.132pt は、季節調整値で前月の値を使った場合(0.128pt)とほぼ同じで、季節の動きを当てているだけでした。原数値では「季節ナイーブ」(前年同月の動きを使う予測)なども比べる相手に加えています。

        nsa.loc[T, f"季節ナイーブ_{h}"] = un.shift(12 - h).loc[T] - un.loc[T]
        nsa.loc[T, f"季節ナイーブ+直近の差_{h}"] = un.shift(12 - h).loc[T] - un.shift(12).loc[T]

1か月先はモデルの問題というより、データの刻み

失業率は0.1ポイント刻みで公表され、2005年以降の月の変化の38%はゼロ、平均の変化幅は0.083ポイントでした。1か月先で前月の値に勝てなかったのは、実装の誤りを疑う前に、この刻みと標本調査の誤差を見ておくべき点でした。

その他

  • 頑健性チェックの 03b_頑健性.py は、03_分析.py を importlib で読み込んで関数を使い回しています。ファイル名が日本語で数字から始まるため、ふつうの import が使えないからです。
  • 3か月先の SARIMAX は外生変数の先を0と置いているので、求人の情報を十分に使えていない可能性があります(ARIMA 系が3か月先でも前月の値に届かなかった一因と考えていますが、推測です)。

言えないこと

  • 予測が当たることと、失業の仕組みが分かったことは別です。重要度は予測に使われたことを示すだけで、因果ではありません。
  • 季節調整値はあとから少し改定されます。2026年10月時点の値で予測しており、当時の公表値で予測した場合とは少し違います。

参考

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?