17th Place Solution (NeuralProphetによる時系列モデリング)

maruyama

17th Place Solution (NeuralProphetによる時系列モデリング)

NeuralProphetによる野菜取引価格の予測

参加者の皆様お疲れ様でした。
また、運営の方々、本コンペを開催してくださりありがとうございました。

17位の解法を共有します。

前に共有したトピック (※1) の考え方をベースに、価格データだけから予測するモデルを作りました。
本トピックでは、NeuralProphet (※2) という時系列モデルを使ってモデリングしています。
NeuralProphetはよくある時系列モデルと同様に、時系列データをトレンド成分＋季節成分＋イベント成分＋自己回帰成分にわけてモデル化します。
Prophet (※3) とのモデル面での違いは、自己回帰成分を扱えるようになったことで予測精度が大幅に向上している点です。

※1 前年同月 or 同年前月の価格を予測値として出力 (https://comp.probspace.com/competitions/vegetable_price/discussions/maruyama-Postc9d85a7c30fdd975aa24)
※2 NeuralProphet (https://neuralprophet.com/contents.html)
※3 Prophet (https://facebook.github.io/prophet/)

セットアップ

import numpy as np
import pandas as pd

from neuralprophet import set_log_level, set_random_seed, NeuralProphet, save, load
set_log_level('ERROR')

train_data = pd.read_csv('data/train_data.csv', index_col=0, parse_dates=[0])
sample_submission = pd.read_csv('data/submission.csv')

前処理

# ds: 年月
# ID: 野菜と地域
# y: 目的変数 (対数変換済み)
df = (
    train_data
    .pipe(lambda df: np.log1p(df))
    .rename_axis(index='ds')
    .reset_index()
    .melt(id_vars=['ds'], var_name='ID', value_name='y')
)

df.head()

	ds	ID	y
0	2016-01-01	えのきだけ_中国	5.752573
1	2016-02-01	えのきだけ_中国	5.726848
2	2016-03-01	えのきだけ_中国	5.505332
3	2016-04-01	えのきだけ_中国	5.429346
4	2016-05-01	えのきだけ_中国	5.451038

# 野菜一覧
veg_list = train_data.columns.str.replace('_.+', '', regex=True).unique().to_list()

学習

# 野菜ごとに分けてモデル化する。
for veg in veg_list:
    # 学習の分散が大きいので、
    # 乱数シードを変えて10回学習と予測を繰り返し、
    # その平均を最終的な予測値とする。
    for seed in range(10):
        print(f'veg: {veg}, seed: {seed}')
        set_random_seed(seed=seed)
        model = NeuralProphet(
            # 月次データのため、年間の季節性のみモデル化する。
            yearly_seasonality=True,
            weekly_seasonality=False,
            daily_seasonality=False,
            # トレンドと季節性を地域別にモデル化する。
            trend_global_local='local',
            season_global_local='local',
            # 自己回帰成分の次数。
            # 昨年度の価格を参照できる12以上の次数を設定すると予測精度が良くなる。
            n_lags=12,
            # 損失関数をMSEにするとコンペの評価指標と一致するため良さそうだが、
            # 試したところ過学習がひどくなったため、デフォルトのHuber lossを採用した。
            # loss_func='MSE'
        )
        df_veg = df.loc[lambda df: df['ID'].str.startswith(veg), :].copy()
        metrics = model.fit(df=df_veg, freq='MS', progress=None, learning_rate=0.01, epochs=20)
        # 野菜の数が多く学習に時間がかかるため、
        # 途中で中断できるよう学習が済んだ野菜のモデルを保存しておく。
        save(model, f'model_20230817_03/model_20230817_03_{veg}_{seed}.np')

予測

forecast = []
for veg in veg_list:
    for seed in range(10):
        print(f'veg: {veg}, seed: {seed}')
        model = load(f'model_20230817_03/model_20230817_03_{veg}_{seed}.np')
        df_veg = df.loc[lambda df: df['ID'].str.startswith(veg), :].copy()
        df_future = model.make_future_dataframe(df_veg, n_historic_predictions=True, periods=1)
        forecast_ = model.predict(df_future)
        forecast_['seed'] = seed
        forecast.append(forecast_)
forecast = pd.concat(forecast, ignore_index=True)

submission = (
    forecast
    # 2019年12月の予測結果だけ抽出する。
    .loc[lambda df: (df['ds'].dt.year == 2019) & (df['ds'].dt.month == 12), ['ID', 'yhat1']]
    # Seed blending
    .groupby('ID', as_index=False)['yhat1'].mean()
    # 提出用にフォーマットを整える。
    .rename(columns={'ID': 'id', 'yhat1': 'y'})
    .merge(sample_submission[['id']], how='right', on='id')
    # 対数変換を元に戻す。
    .assign(y = lambda df: np.expm1(df['y']))
)

submission.to_csv('submission/submission_20230817_03.csv', index=False, header=True)

モデルの解釈

from matplotlib import pyplot as plt
import seaborn as sns
import japanize_matplotlib

veg = 'トマト'
seed = 0

model = load(f'model_20230817_03/model_20230817_03_{veg}_{seed}.np')
model.set_plotting_backend('matplotlib')

df_veg = df.loc[lambda df: df['ID'].str.startswith(veg), :].copy()
df_future = model.make_future_dataframe(df_veg, n_historic_predictions=True, periods=1)
forecast = model.predict(df_future)

Predicting: 0it [00:00, ?it/s]

Predicting: 0it [00:00, ?it/s]

Predicting: 0it [00:00, ?it/s]

Predicting: 0it [00:00, ?it/s]

Predicting: 0it [00:00, ?it/s]

Predicting: 0it [00:00, ?it/s]

Predicting: 0it [00:00, ?it/s]

Predicting: 0it [00:00, ?it/s]

Predicting: 0it [00:00, ?it/s]

関東

まず、関東のトマトの価格の予実をプロットしてみます。細かいずれはありますが、大まかな動きは捉えられていることがわかります。

# 実測値 (黒点) と予測値 (青線)
model.plot(forecast, df_name=f'{veg}_関東', figsize=(12, 4))

次に、モデルの構成要素をプロットしてみます。次のことが読み取れます。

(1段目) 関東のトマトの価格は下降トレンドにある。
(3段目) 春から夏にかけて価格が下がり、秋から冬にかけて価格が上がる、という季節性がある。
(4段目) 同年前月 (AR lag = 1) と前年の同時期 (AR lag = 8-12) が当月の価格と相関している。

一方で、このモデルの課題もいくつか見受けられます。

季節成分が1か月周期で細かく波打っている。
- これはProphetが採用している季節性をフーリエ級数でモデル化するアプローチの限界で、データがない月初以外の期間についても季節性を見出そうとしてしまうためです。
自己回帰成分のうち、偶数月の回帰係数が高くなっている。
- 偶数月だけ高くなる合理的な理由がないため、これは過学習だと思われます。
- もしかしたら前述の季節成分の推定がうまく行っていないことが原因かもしれません。改めて季節成分を見てみると、細かく振動しているせいで奇数月は上昇、偶数月は下降になっています。

model.plot_parameters(df_name=f'{veg}_関東', figsize=(12, 8))

中国

地域を変えて、中国のトマトの価格の予実をプロットしてみます。こちらも関東と同様で、大まかな動きは捉えられています。

# 実測値 (黒点) と予測値 (青線)
model.plot(forecast, df_name=f'{veg}_中国', figsize=(12, 4))

モデルの構成要素をプロットしてみます。季節成分や自己回帰成分の傾向は関東と同じですが、トレンドが関東とは異なり上昇していることがわかります。

model.plot_parameters(df_name=f'{veg}_中国', figsize=(12, 8))