はじめに

時系列分析を学び始めたときに難しく感じやすいのは、次の 3 点です。

何を最初に確認すればよいのか
モデルごとの前提と使い分けが見えにくいこと
精度改善をチューニングの話だけで捉えてしまいやすいこと

小売の実務では、ここにさらに次の難しさが加わります。

欠品があるので、sales がそのまま需要を表すとは限らない
未来時点で未知の情報を特徴量へ入れると、評価だけ高い情報漏洩モデルになってしまう
モデルより前に、粒度、評価方法、特徴量設計で勝負が決まる

この記事では、小売売上データを例にしながら、次の流れを一通り確認します。

データを読み込み、主キー・期間・欠損週を確認する
可視化、ACF / PACF、STL で系列構造を理解する
ETS、ARIMA / SARIMA、Prophet、XGBoost、LightGBM の考え方を整理する
各モデルの主要パラメータ、設定原則、調整の考え方をまとめる
小売で精度を上げるための特徴量設計、外部データ、モデル融合を整理する
MAE、RMSE、MAPE、WAPE、sMAPE、MASE などの誤差指標を比較する
日本のスーパー・ドラッグストアを例に、予測しやすい商品と難しい商品を整理する
小売の時系列予測を、課題設定から改善まで工程別に俯瞰する

retail-sales-time-series-forecasting

対象データ

今回扱うのは、次の 2 つの CSV です。

data/store_sales.csv
data/promotion_data.csv

主な列は次のとおりです。

store: 店舗 ID
dept: 部門 ID
week: 週の起点日
sales: 売上金額
promotion_sales: 販促施策に関連する売上金額

ここで重要なのは、sales が必ずしも「真の需要」を表すとは限らない点です。たとえば欠品がある週は、需要が弱かったとは限らず、売りたくても売れなかった可能性があります。小売予測では、この違いを無視すると精度評価や施策判断を誤りやすくなります。

準備と読み込み

import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.metrics import mean_absolute_error, mean_squared_error

pd.options.display.max_rows = 400
pd.options.display.max_columns = None
warnings.filterwarnings("ignore")

%matplotlib inline

# 売上データを読み込む
df_sales = pd.read_csv("data/store_sales.csv", parse_dates=["week"])

# 販促関連データを読み込む
df_promotion = pd.read_csv("data/promotion_data.csv", parse_dates=["week"])

df_sales.head(2), df_promotion.head(2)

販促データが存在しない週は、「販促が無かった」ケースと「データが欠落している」ケースを区別して考える必要があります。この教材データでは前者として扱い、promotion_sales の欠損を 0 で埋めます。実務では、意味が不明な欠損を機械的に 0 埋めしない方が安全です。

df_all = pd.merge(
    df_sales,
    df_promotion,
    how="left",
    on=["store", "dept", "week"],
    validate="one_to_one",  # 主キー重複があればここで例外になる
)

df_all["promotion_sales"] = df_all["promotion_sales"].fillna(0)
df_all.head(2)

validate="one_to_one" を入れておくと、結合前に気づけなかった主キー重複を早い段階で検知できます。教材コードでは省かれがちですが、実務では非常に有効な防御です。

最初に確認すべき事項: 粒度・主キー・期間・系列数

時系列データを受け取ったら、モデルを選ぶ前に少なくとも次を確認します。

1 行は何を表しているか
主キーは何か
時間粒度は日次・週次・月次のどれか
データ期間はどのくらいか
系列は 1 本か複数か
売上なのか需要なのか、返品を含むのか、税抜きか税込みか

このデータでは、1 行は「ある週の、ある店舗、ある部門の売上」を表しています。したがって主キー候補は week, store, dept です。

primary_key = ["week", "store", "dept"]

print("主キー候補:", primary_key)
print("売上データ件数:", df_sales.shape[0])
print("期間:", df_sales["week"].min(), "->", df_sales["week"].max())
print("系列数:", df_sales[["store", "dept"]].drop_duplicates().shape[0])

重複チェックもしておきます。主キー重複があると、系列の意味そのものが崩れます。

is_unique = (
    df_sales.shape[0]
    == df_sales[primary_key].drop_duplicates().shape[0]
)
is_unique

系列の本数も確認します。ここでは store × dept の組み合わせ数が系列数です。

series_count = df_sales[["store", "dept"]].drop_duplicates().shape[0]
series_count

このデータでは 3000 系列を超えるため、単系列モデルだけで全体を回すか、複数系列をまとめて扱うモデルを使うか、という設計判断も重要になります。

系列の連続性を確認する

時系列分析では、欠けた時点がないかを最初に調べることが重要です。週次データの一部が抜けていると、次の処理が簡単に壊れます。

ラグ特徴量
移動平均や移動分散
周期解釈
予測結果の時点整合性

各系列について、最小週・最大週・実際の週数を集計し、理論上の週数との差を見ると欠損週を調べられます。

date_summary = (
    df_sales.groupby(["store", "dept"])["week"]
    .agg(["min", "max", "nunique"])
    .reset_index()
)

date_summary["missing_weeks"] = (
    (date_summary["max"] - date_summary["min"]).dt.days // 7 + 1
    - date_summary["nunique"]
)

date_summary[date_summary["missing_weeks"] > 0].head()

(date_summary["missing_weeks"] > 0).sum()

この教材データでは、欠損週を持つ系列がかなり含まれます。したがって「全データをそのままモデルへ渡せばよい」とは言えません。系列ごとの再インデックス、補完方針、あるいは欠損系列の除外が必要です。

以降は、説明を分かりやすくするため store=1, dept=1 の系列を例に使います。まず、この系列が連続しているかを確認します。

df_sample_raw = (
    df_sales.query("store == 1 and dept == 1")
    .sort_values("week")
    .copy()
)

sample_missing_weeks = (
    (df_sample_raw["week"].max() - df_sample_raw["week"].min()).days // 7 + 1
    - df_sample_raw["week"].nunique()
)
sample_missing_weeks

# 欠損週が 0 と確認できたら、週次インデックスを明示する
df_sample = (
    df_sample_raw.set_index("week")
    .asfreq("W-MON")
    .copy()
)

df_sample.head()

もし sample_missing_weeks > 0 なら、asfreq("W-MON") によって欠損週が NaN として挿入されます。その場合は、補完するのか、除外するのか、あるいは欠損そのものを特徴量にするのかを先に決めてから先へ進みます。

予測問題としてどう表すか

時系列予測では、先に記号の意味をそろえておくと、その後のモデル差が理解しやすくなります。

$t$: 時点
$y_t$: 時点 $t$ の観測売上
$d_t$: 欠品が無いと仮定したときの潜在需要
$x_t$: 時点 $t$ の外生変数
$T$: 参照する過去長
$h$: 予測ホライズン

1 ステップ先でも複数ステップ先でも、考え方は同じです。過去の売上と外生変数、さらに未来時点で既に分かっている外生変数を使って未来の売上を推定します。

$$\hat y_{t+1:t+h} = f(y_{t-T+1:t},\ x_{t-T+1:t+h})$$

ここで大事なのは、x_{t+1:t+h} に入れてよいのは「予測時点で既知の情報」だけだということです。

入れてよい例: 祝日、曜日、月末、販促計画、価格計画、チラシ配布予定、気象予報
そのまま入れてはいけない例: 未来の実績販促売上、未来の実績来店客数、未来の実績在庫切れ件数

このデータで特に注意したい情報漏洩

この教材では promotion_sales という列が出てきますが、もしこれが「その週に実際に発生した販促由来売上」なら、未来週の promotion_sales をそのまま説明変数へ入れるのは情報漏洩です。なぜなら、それは予測したい売上と同じ時点で、しかも実績から後付けで計測された値だからです。

したがって、この記事では promotion_sales を使う箇所について、次の前提を明示します。

教材としては「将来の販促強度を表す代替指標」として使う
実務では promotion_flag, discount_rate_plan, coupon_count_plan, display_plan など、予測時点で既知の計画値へ置き換える

この確認を飛ばすと、オフライン評価は高精度でも本番で崩れます。

単一系列を抽出して概観する

まずは単純な折れ線図で、全体の動きを目で確認します。

plt.figure(figsize=(12, 4))
plt.plot(df_sample.index, df_sample["sales"], ".-")
plt.xticks(rotation=90)
plt.title("Store 1 / Dept 1 Weekly Sales")
plt.show()

この段階では、厳密な結論を急ぐ必要はありません。まずは次を大まかに把握します。

長期の上昇傾向や下降傾向があるか
毎年似た波がありそうか
異常に大きい点や小さい点がありそうか
レベルの切り替わりや構造変化がありそうか

この「目で見る」工程は地味ですが重要です。後でモデル誤差が悪化したときに、原因がデータ側なのかモデル側なのかを切り分けやすくなるからです。

季節性を可視化して確認する

週次データでは、年ごとの波形を重ねると季節性が見えやすくなります。そこで、年と年内週番号を作って年別の波形を並べます。

df_plot = df_sample.reset_index().copy()
df_plot["year"] = df_plot["week"].dt.year
df_plot["week_of_year"] = df_plot["week"].dt.isocalendar().week.astype(int)

plt.figure(figsize=(12, 5))
for year in sorted(df_plot["year"].unique()):
    tmp = df_plot[df_plot["year"] == year]
    plt.plot(tmp["week_of_year"], tmp["sales"], ".-", label=str(year))

plt.legend()
plt.title("Seasonal Line Plot")
plt.xlabel("ISO Week")
plt.ylabel("Sales")
plt.show()

この図では、毎年ある程度似た時期に山と谷が現れるかを見ます。小売では次のような論点と結びつきます。

年末年始
ゴールデンウィーク
お盆
気温上昇・低下
新生活需要
花粉、風邪、インフルエンザなどの季節要因

次に、同じ週番号で売上分布がどのくらい散らばるかを箱ひげ図で確認します。

df_boxplot = df_plot.groupby("week_of_year")["sales"].agg(list).reset_index()

plt.figure(figsize=(14, 5))
plt.boxplot(df_boxplot["sales"].values, labels=df_boxplot["week_of_year"].values)
plt.xticks(rotation=90)
plt.title("Seasonal Box Plot")
plt.xlabel("ISO Week")
plt.ylabel("Sales")
plt.show()

箱ひげ図を見ると、次が把握しやすくなります。

毎年高くなりやすい週
毎年低くなりやすい週
特定の週だけばらつきが大きいか
イベントが年ごとに効いたり効かなかったりするか

週次データでは ISO 週 53 が混ざることがあります。実務では、53 週目を 52 週目とまとめるか、年末特殊週として別扱いするかを先に決めると集計が安定します。

年次水準差を確認する

もし特定の時期に販促やイベントの影響が強そうなら、その区間をいったん除いて年ごとの差を比べるのも有効です。たとえば 19 週から 43 週だけに絞れば、年末年始や大型連休をある程度外した状態で基礎水準を見やすくなります。

df_trend = df_plot.copy()
df_trend = df_trend[
    (df_trend["week_of_year"] >= 19)
    & (df_trend["week_of_year"] <= 43)
]

trend_boxplot = df_trend.groupby("year")["sales"].agg(list).reset_index()

plt.figure(figsize=(8, 4))
plt.boxplot(
    trend_boxplot["sales"].values,
    labels=trend_boxplot["year"].values,
)
plt.title("Trend Box Plot")
plt.xlabel("Year")
plt.ylabel("Sales")
plt.show()

このような図は、季節性よりも年ごとの基礎水準を見たいときに役立ちます。たとえば競合出店や価格改定、店内改装の前後で水準が変わっていないかを疑うきっかけになります。

ACF と PACF で周期と自己相関を確認する

見た目に加えて、自己相関でも周期性を確認します。

ACF（自己相関関数）は、系列とラグ付き系列の相関を見て、どの周期で似た動きが繰り返されているかを調べる方法です。弱定常性を仮定したとき、ラグ $h$ における自己相関は次のように書けます。

$$
\rho_h = \frac{\operatorname{Cov}(y_t, y_{t-h})}{\operatorname{Var}(y_t)}
$$

ここで注目したいのは、一定間隔で山が現れるかどうかです。すべてのラグで相関が高い必要はありません。週次データで 52 ラグ付近に山があれば、1 年周期の可能性が高いと考えられます。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

fig = plot_acf(df_sample["sales"], lags=140)
plt.show()

AR 系の次数ヒントを見るときは PACF も確認します。

fig = plot_pacf(df_sample["sales"], lags=60, method="ywm")
plt.show()

ACF / PACF の読み方は最初から厳密でなくて構いません。実務での第一歩としては次の理解で十分です。

52 ラグ近辺に山があるなら年周期を疑う
差分後でも自己相関が長く残るなら、まだ非定常性が残っている可能性がある
PACF の切れ方は AR 次数の候補を見る手掛かりになる
ACF の切れ方は MA 次数の候補を見る手掛かりになる

ただし、ACF / PACF だけで次数を決め打ちしない方が安全です。最終的には残差、AIC / BIC、時系列 CV、業務上の解釈まで合わせて判断します。

STL で系列を分解する

系列の構造をより明確にしたいときは STL (Seasonal and Trend decomposition using Loess)が便利です。STL は系列を次の 3 つに分けます。

トレンド
季節成分
残差

加法分解として書くと、考え方は次の形です。

$$
y_t = T_t + S_t + R_t
$$

$T_t$: 長期水準の動き
$S_t$: 周期的に繰り返す成分
$R_t$: それ以外の残差

もし季節振幅が水準と一緒に大きくなるなら、加法分解の前に対数変換してから STL をかけることがあります。加法か乗法かは機械的に決めず、「水準が上がると波も大きくなるか」で考えると分かりやすいです。

from statsmodels.tsa.seasonal import STL

# robust=True にすると、外れ値の影響を受けにくくなる
stl = STL(df_sample["sales"], period=52, robust=True)
stl_result = stl.fit()
stl_result.plot();

分解結果を表に持っておくと、後で特徴量として再利用しやすくなります。

df_stl = df_sample.copy()
df_stl["trend"] = stl_result.trend
df_stl["seasonal"] = stl_result.seasonal
df_stl["residual"] = stl_result.resid
df_stl.head()

STL は予測モデルではありませんが、次を確認するのにとても役立ちます。

トレンドが強いか
季節成分が安定しているか
大きな異常が残差に出ていないか
後続モデルへ渡す前にログ変換や外れ値処理が要るか

STL の主要パラメータ

パラメータ	役割	典型設定	どう変えるか
`period`	季節周期長	週次なら `52`	日次なら `7` や `365`、月次なら `12` を候補にする
`robust`	外れ値の影響を抑えるか	`True` を推奨	セールや欠品などで極端値が混じる小売では有効
`seasonal`	季節平滑の窓	既定値で始める	季節波形がギザギザなら広げる
`trend`	トレンド平滑の窓	既定値で始める	トレンドが揺れすぎるなら広げる
`low_pass`	低周波成分の平滑	既定値で始める	周波数分離が不自然なときだけ触る

STL の設定原則

最初は period と robust=True だけ決め、他は既定値でよい
季節周期を間違えると、分解結果全体が不自然になる
外れ値が多い小売系列では robust=True が有効
水準と季節振幅が一緒に大きくなるなら、log1p 変換を検討する

ETS で短期予測の基準を作る

ETS は指数平滑法(Exponential Smoothing)に基づくモデルで、比較的短い系列のベースラインとして使いやすい手法です。直近の観測に比較的大きな重みを置きながら、

水準 level
傾き trend
季節成分 seasonality

を逐次更新していきます。

加法トレンド・加法季節性の予測式は、概念的には次のように書けます。

$$\hat y_{t+h|t} = \ell_t + h b_t + s_{t+h-m(k+1)},\quad k=\left\lfloor\frac{h-1}{m}\right\rfloor$$

$\ell_t$: 時点 $t$ の水準
$b_t$: 時点 $t$ の傾き
$s_t$: 季節成分
$m$: 周期長。週次の年周期なら 52

代表的な加法型更新式は次のとおりです。

$$
\ell_t = \alpha (y_t - s_{t-m}) + (1-\alpha)(\ell_{t-1} + b_{t-1})
$$

$$
b_t = \beta (\ell_t - \ell_{t-1}) + (1-\beta)b_{t-1}
$$

$$
s_t = \gamma (y_t - \ell_t) + (1-\gamma)s_{t-m}
$$

全部を暗記する必要はありません。ETS を「最近の情報を重めに見ながら、水準・トレンド・季節性を滑らかに更新するモデル」と捉えれば十分です。

ここでいう ETS は Error-Trend-Seasonality の略で、モデル構造は ETS(E,T,S) と書きます。

E: 誤差成分。A は加法誤差、M は乗法誤差
T: トレンド成分。N は無し、A は加法トレンド、Ad は減衰付き加法トレンド
S: 季節成分。N は無し、A は加法季節性、M は乗法季節性

たとえば、ETS(A,N,N) は単純指数平滑法、ETS(A,A,N) は Holt の線形トレンドモデル、ETS(A,Ad,A) は減衰トレンド付きの加法季節モデルに対応します。

ARIMA と違って、ETS では ACF / PACF を見て次数を決める発想が中心ではありません。実務では、季節周期 m を先に決めたうえで、ETS(E,T,S) の組み合わせを並べ、各候補で平滑化係数を最尤推定し、AIC / AICc / BIC で比較して絞り込む流れが自然です。

ETS の調整手順

まず季節周期 m を決める。週次の年周期なら 52、月次なら 12、日次なら 7 などが出発点になる
次に ETS(E,T,S) の候補を作る。たとえば誤差は add / mul、トレンドは None / add、必要なら damped_trend=True を足し、季節性は None / add / mul を比べる
各候補について、α、β、γ、必要なら φ を最尤推定する
AIC や、標本が短めのときは AICc を優先して比較し、最後に残差診断と時系列 CV で確かめる

売上に 0 や負値が含まれる場合、乗法誤差や乗法季節性はそのまま使えないことがあります。mul 候補を比較するときは、データの下限を必ず確認します。

from statsmodels.tsa.exponential_smoothing.ets import ETSModel

forecast_origin = pd.to_datetime("2012-07-30")

ets_train = df_sample[df_sample.index <= forecast_origin][["sales"]].copy()

candidate_specs = [
    {"error": "add", "trend": None, "seasonal": None, "damped_trend": False},
    {"error": "add", "trend": "add", "seasonal": None, "damped_trend": False},
    {"error": "add", "trend": "add", "seasonal": "add", "damped_trend": False},
    {"error": "add", "trend": "add", "seasonal": "add", "damped_trend": True},
    {"error": "add", "trend": "add", "seasonal": "mul", "damped_trend": True},
]

ets_compare_rows = []

for spec in candidate_specs:
    try:
        model = ETSModel(
            ets_train["sales"],
            error=spec["error"],
            trend=spec["trend"],
            damped_trend=spec["damped_trend"],
            seasonal=spec["seasonal"],
            seasonal_periods=52 if spec["seasonal"] is not None else None,
        )
        fit = model.fit()
        ets_compare_rows.append(
            {
                **spec,
                "aic": fit.aic,
                "aicc": fit.aicc,
                "bic": fit.bic,
            }
        )
    except Exception:
        continue

pd.DataFrame(ets_compare_rows).sort_values(["aicc", "aic", "bic"])

このように、ETS では「まずモデル構造をいくつか用意し、情報量規準で比べる」という進め方が基本になります。ARIMA のように ACF / PACF を見て p や q を順に決めるというより、候補モデル全体を横並びで比べるイメージです。

from statsmodels.tsa.exponential_smoothing.ets import ETSModel

forecast_origin = pd.to_datetime("2012-07-30")

ets_train = df_sample[df_sample.index <= forecast_origin][["sales"]].copy()

ets_model = ETSModel(
    ets_train["sales"],
    error="add",
    trend="add",
    seasonal="add",
    seasonal_periods=52,
)

ets_fit = ets_model.fit()

まずは学習データへの当てはまりを見ます。

ets_train["fit_sales"] = ets_fit.fittedvalues

ets_plot = ets_train[
    (ets_train.index >= "2011-01-03")
    & (ets_train.index <= "2011-06-27")
]

plt.figure(figsize=(12, 4))
plt.plot(ets_plot.index, ets_plot["sales"], ".-", label="actual")
plt.plot(ets_plot.index, ets_plot["fit_sales"], ".-", label="fit")
plt.legend()
plt.xticks(rotation=90)
plt.title("ETS Fitted Values")
plt.show()

続いて 2 週間先を予測し、実測値と並べます。

h = 2
ets_forecast = ets_fit.forecast(steps=h)

ets_result = pd.DataFrame(
    {
        "week": pd.date_range(
            start=ets_train.index.max() + pd.Timedelta(weeks=1),
            periods=h,
            freq="W-MON",
        ),
        "pred_sales": ets_forecast.values,
    }
)

ets_test = (
    df_sample.reset_index()
    .loc[
        lambda data: (
            (data["week"] > forecast_origin)
            & (data["week"] <= forecast_origin + pd.Timedelta(weeks=2))
        ),
        ["week", "sales"],
    ]
)

ets_result.merge(ets_test, on="week", how="left")

ETS の主要パラメータ

パラメータ	役割	典型候補	どう変えるか
`error`	誤差の形	`add`, `mul`	レベルが大きいほど誤差振幅も大きいなら `mul` を検討
`trend`	トレンドの有無	`None`, `add`	明確な上昇・下降があれば `add`
`damped_trend`	長期でトレンドを減衰させるか	`True/False`	遠い将来で直線的増加が過大なら `True`
`seasonal`	季節性の形	`None`, `add`, `mul`	水準と一緒に季節幅が増えるなら `mul` を検討
`seasonal_periods`	季節周期	週次なら `52`	月次なら `12`、日次なら `7` など
`initialization_method`	初期値設定	既定値	まず既定値で十分

平滑化係数 `α` `β` `γ` `φ` の見方

error / trend / seasonal の構造を決めたあと、モデル内部では平滑化係数が推定されます。実務では、最初から手で固定するより、まず最尤推定に任せ、その結果を見て構造側を見直す方が進めやすいです。

記号	役割	読み方の目安
`α`	水準の平滑化係数	大きいほど最新値へ強く反応し、小さいほど水準更新が滑らかになる
`β`	トレンドの平滑化係数	大きいほど傾きの変化を素早く追い、小さいほどトレンドを安定的にみる
`γ`	季節成分の平滑化係数	大きいほど季節パターンを早く更新し、小さいほど季節性を固定的にみる
`φ`	減衰係数	`1` に近いほどトレンドを長く維持し、小さいほど将来へ行くほど傾きを弱める

α、β、γ は一般に 0 から 1 の範囲で解釈する
φ は通常 0 より大きく 1 以下で、減衰トレンドを使うときだけ現れる
α や β が極端に大きく、しかも残差が荒いなら、ノイズまで追いかけている可能性がある
γ が大きすぎると、一時的なイベント変動を季節性として抱え込みやすい

減衰トレンドを入れるときは、φ によって将来の傾きが少しずつ弱まります。長期予測で直線的な伸びが不自然に大きくなりやすい系列では、damped_trend=True と φ の組み合わせが効きやすいです。

ETS の設定原則と調整方法

まずは季節周期 m を決める。週次小売なら 52 が第一候補
そのうえで error / trend / seasonal / damped_trend の組み合わせを比べる。ETS は ACF / PACF で決めるより、候補モデル比較で進める方が自然
最初から α β γ φ を細かく手で触るより、まず自動推定に任せる方が安定しやすい
レベルが上がるほど季節幅も大きくなるなら、対数変換か mul を検討する
長期予測でトレンドが過大になりやすいなら damped_trend=True を試す
系列が短めなら AIC より AICc を重視すると比較が安定しやすい
最終判断は情報量規準だけで終えず、残差と時系列 CV でも確かめる

ETS が向く場面

比較的短い単系列
強すぎない外生変数
まず基準モデルを作りたいとき
学習・解釈を速く回したいとき

一方で、価格、天候、販促、競合要因などを直接入れたいときは、後述の Prophet や木系モデルの方が扱いやすいことが多いです。

ARIMA と SARIMA: 自己相関を明示的に使う

ARIMA は系列の自己相関構造を使って予測する古典的な方法です。非定常な系列には差分をとってから、過去値と過去誤差でモデル化します。

まず差分系列を

$$
w_t = (1-B)^d y_t
$$

とおくと、ARIMA$(p,d,q)$ は次のように表せます。

$$
w_t = c + \sum_{i=1}^{p} \phi_i w_{t-i} + \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j}
$$

AR(AutoRegressive): 過去値の影響
I(Integrated): 差分回数
MA(Moving Average): 過去誤差の影響

実務で ARIMA の調整を始めるときは、いきなり多数の次数候補を総当たりするよりも、まず d、次に p と q の当たりを付けてから候補モデルを比較すると進めやすいです。最も基本的な流れは次の 4 段階です。

ADF 検定で系列が定常的かを確認し、差分回数 d の候補を決める
差分後系列の PACF を見て、AR 次数 p の候補を絞る
差分後系列の ACF を見て、MA 次数 q の候補を絞る
複数の ARIMA モデルを当てて AIC / BIC を比較し、より小さいものを優先候補にする

もちろん、これはあくまで出発点です。最終的には、残差がホワイトノイズに近いか、予測性能がローリングオリジン(Rolling Origin、別名: 時系列クロスバリデーション、ロールフォワード) で安定しているかまで確認して初めて「その次数でよい」と言えます。

ここでいうローリングオリジンとは、過去データから未来を予測する検証を、時間軸に沿って予測起点を少しずつ前へずらしながら繰り返す手法です。各回で「その時点までに見えているデータだけ」を使って学習し、その直後の未来区間で検証するため、本番運用にかなり近い形でモデルを評価できます。

ARIMA の次数を決める基本手順

1. ADF 検定で `d` の候補を決める

d は「何回差分をとれば系列をほぼ定常とみなせるか」を表します。ここでよく使うのが ADF (Augmented Dickey-Fuller) 検定です。ざっくり言えば、単位根があるかを調べる検定で、p-value が十分小さければ「少なくともこの系列は非定常とは言い切れない」と判断しやすくなります。

原系列で ADF を行い、定常とみなしにくければ d=1 を候補にする
それでも定常化しなければ 2 階差分も検討する
ただし差分を増やしすぎると系列の構造を壊しやすいので、通常は必要最小限にとどめる

from statsmodels.tsa.stattools import adfuller


def adf_check(series, name="series"):
    result = adfuller(series.dropna())
    print(f"{name} ADF statistic: {result[0]:.4f}")
    print(f"{name} p-value: {result[1]:.4f}")


adf_check(df_sample["sales"], name="original")
adf_check(df_sample["sales"].diff(), name="diff_1")
adf_check(df_sample["sales"].diff().diff(), name="diff_2")

実務では、ADF 検定の結果だけで機械的に d を決めるよりも、時系列プロット、ACF の減衰、差分後の解釈しやすさも合わせて見る方が安全です。たとえば、明らかなトレンドがあり、1 階差分後に平均回帰的な動きへ変わるなら d=1 が自然な候補になります。

2. PACF で `p` の候補を見る

p は AR 項の次数です。差分後系列の PACF を見ると、「何期前までの値が直接効いていそうか」の当たりを付けやすくなります。典型的には、PACF があるラグで目立って、その先で急に小さくなるなら、その付近が p の候補です。

ラグ 1 だけが強ければ p=1 を候補にする
ラグ 1 と 2 が目立つなら p=2 も候補に入れる
きれいに切れないときは 1, 2, 3 程度の小さい値から比較する

ここで大事なのは、PACF は「p を一意に決める魔法の図」ではなく、あくまで候補を絞るための手掛かりだという点です。

3. ACF で `q` の候補を見る

q は MA 項の次数です。差分後系列の ACF を見ると、過去のショックが何期ぶん残っていそうかの目安を得られます。典型的には、ACF が低いラグで大きく、その後すっと弱まるなら、その付近が q の候補です。

ラグ 1 が目立つなら q=1 を候補にする
ラグ 1 と 2 が大きいなら q=2 も候補に入れる
季節ラグでも山があるなら、非季節 q だけでなく季節 MA 項 Q の必要性も疑う

週次小売では、非季節の q を決めるときでも 52 週付近の山を無視しない方がよいです。年周期の自己相関が強いなら、単純な ARIMA より最初から SARIMA の方が素直に当てはまることがあります。

4. AIC / BIC で候補モデルを比較する

PACF と ACF で p と q の候補を絞ったら、複数の ARIMA モデルを実際に当てて AIC / BIC を比較します。どちらも「当てはまり」と「複雑さ」のバランスを見る指標で、基本的には値が小さいほど望ましいと考えます。

AIC は比較的複雑なモデルも許容しやすい
BIC はパラメータ数への罰則がやや強く、より単純なモデルを選びやすい
実務では両方見たうえで、残差診断と時系列 CV も合わせて判断する

from statsmodels.tsa.arima.model import ARIMA

candidate_orders = [
    (0, 1, 1),
    (1, 1, 0),
    (1, 1, 1),
    (2, 1, 0),
    (2, 1, 1),
]

arima_compare_rows = []

for order in candidate_orders:
    try:
        fit = ARIMA(df_sample["sales"], order=order).fit()
        arima_compare_rows.append(
            {
                "order": order,
                "aic": fit.aic,
                "bic": fit.bic,
            }
        )
    except Exception:
        continue

pd.DataFrame(arima_compare_rows).sort_values(["aic", "bic"])

この表で AIC / BIC が小さいモデルを優先候補にし、そのうえで残差の自己相関が残っていないかを確認します。つまり、PACF で p、ACF で q、ADF で d を考え、最後に AIC/BIC で候補比較する、という流れです。

ただし、週次小売では年周期が強いことが多いため、実務では ARIMA より SARIMA / SARIMAX を先に検討することが少なくありません。つまり、非季節成分の $(p,d,q)$ に加え、季節成分 $(P,D,Q,s)$ を入れます。週次年周期なら s=52 です。

まずは原系列と 1 階差分系列の ACF / PACF を見ます。

df_arima = df_sample.reset_index().copy()
df_arima["sales_lag_1"] = df_arima["sales"].shift(1)
df_arima["sales_diff_1"] = df_arima["sales"] - df_arima["sales_lag_1"]

fig = plot_acf(df_sample["sales"], lags=60)
plt.show()

fig = plot_acf(df_arima["sales_diff_1"].dropna(), lags=60)
plt.show()

fig = plot_pacf(df_arima["sales_diff_1"].dropna(), lags=30, method="ywm")
plt.show()

説明用に単純な次数から始めても構いませんが、週次小売なら年周期を無視しない方が自然です。そこでここでは SARIMAX を使って、非季節成分と季節成分を両方入れます。

from statsmodels.tsa.statespace.sarimax import SARIMAX

sarima_train = df_sample[df_sample.index <= forecast_origin][["sales"]].copy()

sarima_model = SARIMAX(
    sarima_train["sales"],
    order=(1, 1, 0),
    seasonal_order=(1, 1, 0, 52),
    enforce_stationarity=False,
    enforce_invertibility=False,
)

sarima_fit = sarima_model.fit(disp=False)
sarima_fit.summary()

予測結果は次のように確認できます。

sarima_forecast = sarima_fit.get_forecast(steps=2)

sarima_result = pd.DataFrame(
    {
        "week": pd.date_range(
            start=sarima_train.index.max() + pd.Timedelta(weeks=1),
            periods=2,
            freq="W-MON",
        ),
        "pred_sales": sarima_forecast.predicted_mean.values,
    }
)

sarima_test = (
    df_sample.reset_index()
    .loc[
        lambda data: (
            (data["week"] > forecast_origin)
            & (data["week"] <= forecast_origin + pd.Timedelta(weeks=2))
        ),
        ["week", "sales"],
    ]
)

sarima_result.merge(sarima_test, on="week", how="left")

ARIMA / SARIMA の主要パラメータ

パラメータ	役割	典型候補	どう変えるか
`order=(p,d,q)`	非季節成分	`(1,1,0)`, `(1,1,1)` など	差分後 ACF / PACF、AIC / BIC、残差で調整
`seasonal_order=(P,D,Q,s)`	季節成分	週次なら `s=52`	年周期が見えるならまず入れる
`trend`	定数項・トレンド	`n`, `c`, `t`, `ct`	平均水準や線形トレンドを入れたいとき
`exog`	外生変数	`None` or DataFrame	価格・販促計画など既知の将来変数を入れる
`enforce_stationarity`	定常性制約	`True/False`	収束しづらいとき一時的に `False` を試す
`enforce_invertibility`	可逆性制約	`True/False`	同上

ARIMA / SARIMA の設定原則と調整方法

d と D は「必要最小限」にする。差分しすぎると情報を失う
週次小売で年周期が見えるなら seasonal_order を先に検討する
ACF がゆっくり減衰し PACF が早く切れるなら AR を疑う
PACF がゆっくり減衰し ACF が早く切れるなら MA を疑う
候補モデルは AIC / BIC だけでなく、残差 ACF と時系列 CV で比べる
外生変数を入れるときは SARIMAX を使い、未来 exog の入手可能性を必ず確認する

ARIMA 系が向く場面

単系列で自己相関構造を明示的に使いたい
価格や販促よりも、系列そのものの慣性が強い
データ量がそれほど多くない
解釈しやすい統計モデルを先に作りたい

Prophet でトレンド・季節性・イベントを分けて扱う

Prophet は、トレンド、季節性、休日・イベント効果、追加回帰変数を比較的扱いやすいモデルです。中心の考え方は「時系列を解釈しやすい部品に分ける」ことです。

Prophet の代表的な表現は次のように書けます。

$$
y(t) = g(t) + s(t) + h(t) + x(t)^\top \beta + \varepsilon_t
$$

$g(t)$: トレンド
$s(t)$: 季節性
$h(t)$: 祝日や特定イベントの効果
$x(t)^\top \beta$: 追加回帰変数の効果
$\varepsilon_t$: 残差

線形トレンドなら、たとえば

$$
g(t) = kt + m
$$

と書けます。季節性はフーリエ級数で表すのが Prophet の特徴です。

$$
s(t) = \sum_{n=1}^{N} \left( a_n \cos \frac{2\pi n t}{P} + b_n \sin \frac{2\pi n t}{P} \right)
$$

$P$: 周期長
$N$: フーリエ項数

ここでは promotion_sales を「将来の販促強度を表す代替指標」とみなして説明用に追加します。ただし、先に述べたとおり、もし promotion_sales が実績値なら未来週へそのまま使うのは情報漏洩です。実務では promotion_flag や discount_rate_plan など、予測時点で既知の計画値へ置き換えてください。

from prophet import Prophet

prophet_train = (
    df_all.query("week <= '2012-07-30' and store == 1 and dept == 1")
    .copy()
    .rename(columns={"week": "ds", "sales": "y"})
)

prophet_test = (
    df_all.query("week == '2012-08-06' and store == 1 and dept == 1")
    .copy()
    .rename(columns={"week": "ds", "sales": "y"})
)

prophet_model = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=False,  # 週次集計なので週内季節性は使わない
    daily_seasonality=False,
    seasonality_mode="additive",
    changepoint_prior_scale=0.05,
    seasonality_prior_scale=10.0,
)

# 日本の祝日を追加する
prophet_model.add_country_holidays(country_name="JP")

# 注意:
# `promotion_sales` が実績値なら情報漏洩になる。
# 実務では予測時点で既知の販促計画変数へ置き換える。
prophet_model.add_regressor("promotion_sales", standardize=True)

prophet_model.fit(prophet_train[["ds", "y", "promotion_sales"]])

分解結果と当てはまりは次のように可視化できます。

prophet_fit = prophet_model.predict(prophet_train[["ds", "promotion_sales"]])

prophet_model.plot_components(prophet_fit)
prophet_model.plot(prophet_fit);

1 週先予測は次のように確認できます。

prophet_pred = prophet_model.predict(prophet_test[["ds", "promotion_sales"]])

prophet_test["yhat"] = prophet_pred["yhat"].values
prophet_test[["ds", "y", "promotion_sales", "yhat"]]

複数部門へ横展開する考え方は同じです。ただし、部門ごとに別モデルを立てると運用コストが上がるため、部門数が多いときは後述の LightGBM / XGBoost のような複数系列対応の方が現実的なことも多いです。

Prophet の主要パラメータ

パラメータ	役割	典型候補	どう変えるか
`growth`	トレンド形状	`"linear"`, `"logistic"`	上限・下限があるなら `logistic`
`changepoint_prior_scale`	トレンド変化の柔軟性	`0.01` - `0.5`	小さいと滑らか、大きいと折れやすい
`seasonality_prior_scale`	季節成分の柔軟性	`1` - `20`	季節性を強めたいとき増やす
`holidays_prior_scale`	祝日効果の柔軟性	`1` - `20`	祝日効果を強く出したいとき増やす
`seasonality_mode`	加法 / 乗法	`"additive"`, `"multiplicative"`	水準と一緒に季節幅が増えるなら `multiplicative`
`n_changepoints`	変化点候補数	既定値で可	構造変化が多い系列で増やす
`add_regressor()`	追加回帰変数	価格、販促、気温など	未来時点で既知かを必ず確認する
`add_seasonality()`	独自季節性	月次、四半期、棚替え周期など	標準季節性では足りないとき追加

Prophet の設定原則と調整方法

まずは yearly_seasonality=True と祝日追加から始める
トレンドが波打ちすぎるなら changepoint_prior_scale を下げる
祝日効果が弱すぎるなら holidays_prior_scale を上げる
水準が上がると季節振幅も増えるなら seasonality_mode="multiplicative" を試す
回帰変数は「未来に分かるもの」だけ入れる
追加回帰変数を入れすぎると、線形回帰部分が過学習することがある

たとえば 1 系列に対して changepoint_prior_scale と seasonality_mode をローリングオリジンで比較するなら、次のように書けます。

def wape(y_true, y_pred, eps=1e-8):
    y_true = np.asarray(y_true)
    y_pred = np.asarray(y_pred)
    return 100 * np.sum(np.abs(y_true - y_pred)) / np.maximum(np.sum(np.abs(y_true)), eps)


prophet_df = (
    df_all.query("store == 1 and dept == 1")
    .sort_values("week")
    .rename(columns={"week": "ds", "sales": "y"})
    .copy()
)

# 説明用の代替列。実務では販促計画値など、予測時点で既知の列へ置き換える。
prophet_df["promotion_plan_proxy"] = prophet_df["promotion_sales"]

cutoffs = pd.to_datetime([
    "2011-12-26",
    "2012-03-26",
    "2012-07-30",
])

param_grid = [
    {
        "changepoint_prior_scale": 0.03,
        "seasonality_prior_scale": 5.0,
        "seasonality_mode": "additive",
    },
    {
        "changepoint_prior_scale": 0.10,
        "seasonality_prior_scale": 10.0,
        "seasonality_mode": "additive",
    },
    {
        "changepoint_prior_scale": 0.10,
        "seasonality_prior_scale": 10.0,
        "seasonality_mode": "multiplicative",
    },
]

prophet_cv_rows = []

for params in param_grid:
    fold_scores = []

    for cutoff in cutoffs:
        train = prophet_df.loc[
            prophet_df["ds"] <= cutoff,
            ["ds", "y", "promotion_plan_proxy"],
        ]
        valid = prophet_df.loc[
            (prophet_df["ds"] > cutoff)
            & (prophet_df["ds"] <= cutoff + pd.Timedelta(weeks=4)),
            ["ds", "y", "promotion_plan_proxy"],
        ]

        if valid.empty:
            continue

        model = Prophet(
            yearly_seasonality=True,
            weekly_seasonality=False,
            daily_seasonality=False,
            changepoint_prior_scale=params["changepoint_prior_scale"],
            seasonality_prior_scale=params["seasonality_prior_scale"],
            seasonality_mode=params["seasonality_mode"],
        )
        model.add_country_holidays(country_name="JP")
        model.add_regressor("promotion_plan_proxy", standardize=True)
        model.fit(train)

        pred = model.predict(valid[["ds", "promotion_plan_proxy"]])
        fold_scores.append(wape(valid["y"], pred["yhat"]))

    prophet_cv_rows.append(
        {
            **params,
            "mean_wape": np.mean(fold_scores),
            "std_wape": np.std(fold_scores),
        }
    )

prophet_cv_result = pd.DataFrame(prophet_cv_rows).sort_values("mean_wape")
prophet_cv_result

Prophet が向く場面

トレンド、季節性、休日効果を分けて説明したい
業務側に図で説明しながら進めたい
祝日やイベントが強い
数十系列程度までを個別に見たい

Prophet に外部要因を入れるときのコツ

祝日や大型連休は holidays か add_country_holidays() で入れる
気温や降水量は add_regressor() で入れる
地域イベントは holidays でも regressor でもよい
競合販促や失業率のように未来が未確定なものは、シナリオ値か別モデルの予測値を使う

教材データには気象や人流の実データが含まれていないため、下のコードでは「どのように列を準備してモデルへ渡すか」が分かるよう、週番号から説明用の外部特徴量を作っています。実務ではこの部分を実データ取得テーブルへ置き換えてください。

calendar_helper = prophet_df[["ds"]].drop_duplicates().copy()
calendar_helper["week_of_year"] = calendar_helper["ds"].dt.isocalendar().week.astype(int)
calendar_helper["month"] = calendar_helper["ds"].dt.month

weather_weekly = calendar_helper[["ds"]].copy()
weather_weekly["avg_temp_forecast"] = 15 + 10 * np.sin(2 * np.pi * calendar_helper["week_of_year"] / 52)
weather_weekly["rainfall_forecast"] = 5 + 3 * np.cos(2 * np.pi * calendar_helper["week_of_year"] / 52)

promo_plan_weekly = calendar_helper[["ds"]].copy()
promo_plan_weekly["promo_flag_plan"] = calendar_helper["week_of_year"].isin([1, 13, 26, 39, 48]).astype(int)
promo_plan_weekly["discount_rate_plan"] = promo_plan_weekly["promo_flag_plan"] * 0.10

foot_traffic_weekly = calendar_helper[["ds"]].copy()
foot_traffic_weekly["foot_traffic_data"] = (
    100
    + 12 * np.sin(2 * np.pi * calendar_helper["week_of_year"] / 52)
    + 4 * calendar_helper["month"].isin([12, 1]).astype(int)
)

prophet_exog = prophet_df[["ds", "y"]].copy()

prophet_exog = (
    prophet_exog
    .merge(weather_weekly, on="ds", how="left", validate="one_to_one")
    .merge(promo_plan_weekly, on="ds", how="left", validate="one_to_one")
    .merge(foot_traffic_weekly, on="ds", how="left", validate="one_to_one")
)

regressor_cols = [
    "avg_temp_forecast",
    "rainfall_forecast",
    "discount_rate_plan",
    "promo_flag_plan",
    "foot_traffic_data",
]

prophet_exog[regressor_cols] = prophet_exog[regressor_cols].fillna(0)

prophet_model = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=False,
    daily_seasonality=False,
    seasonality_mode="multiplicative",
)
prophet_model.add_country_holidays(country_name="JP")

for col in regressor_cols:
    prophet_model.add_regressor(col, standardize=True)

prophet_model.fit(prophet_exog[["ds", "y"] + regressor_cols])

XGBoost / LightGBM でラグ特徴量と外部変数を使う

XGBoost や LightGBM では、系列を表形式データへ落とし込み、ラグ特徴量を作って学習させます。ここが ETS / SARIMA / Prophet との大きな違いです。

1 ステップ予測の考え方は次のように書けます。

$$\hat y_{t+1} = f(x_{t+1})$$

$$x_{t+1} = \left[y_t,\ y_{t-1},\ y_{t-4},\ y_{t-52},\ p_{t+1}^{plan},\ p_t,\ \bar y_{t-4:t-1},\ c_{t+1},\ w_{t+1},\ \ldots\right]$$

ここで $p_{t+1}^{plan}$ は予測時点で既知の販促計画、$c_{t+1}$ はカレンダー要因、$w_{t+1}$ は気象予報などの外生変数を表します。promotion_sales の実績値とは別物です。

木系モデル自体は「時間」を自動で理解しません。したがって、こちら側で時系列らしい特徴量を設計して初めて強みが出ます。特に小売では次の特徴量がよく効きます。

直近ラグ
前年同週ラグ
移動平均、移動分散
価格、販促、カレンダー
天候、イベント、在庫、競合情報
店舗属性、部門属性、地域属性

以下では、まず 1 系列に対して LightGBM 用の特徴量を作ります。

from lightgbm import LGBMRegressor

lgb_sample = (
    df_all.query("store == 1 and dept == 1")
    .sort_values("week")
    .copy()
)

# カレンダー特徴量
lgb_sample["week_of_year"] = lgb_sample["week"].dt.isocalendar().week.astype(int)
lgb_sample["month"] = lgb_sample["week"].dt.month
lgb_sample["week_sin"] = np.sin(2 * np.pi * lgb_sample["week_of_year"] / 52)
lgb_sample["week_cos"] = np.cos(2 * np.pi * lgb_sample["week_of_year"] / 52)

# 売上ラグ特徴量
lgb_sample["sales_lag_1"] = lgb_sample["sales"].shift(1)
lgb_sample["sales_lag_4"] = lgb_sample["sales"].shift(4)
lgb_sample["sales_lag_52"] = lgb_sample["sales"].shift(52)

# 販促の過去情報
lgb_sample["promo_lag_1"] = lgb_sample["promotion_sales"].shift(1)
lgb_sample["promo_lag_52"] = lgb_sample["promotion_sales"].shift(52)

# 直近移動平均
lgb_sample["sales_roll_mean_4"] = lgb_sample["sales"].shift(1).rolling(4).mean()
lgb_sample["sales_roll_mean_13"] = lgb_sample["sales"].shift(1).rolling(13).mean()

# 注意:
# ここでは教材の都合で代替指標として使う。
# 実務では予測時点で既知の販促計画変数に置き換える。
lgb_sample["promotion_plan_proxy"] = lgb_sample["promotion_sales"]

feature_cols = [
    "week_sin",
    "week_cos",
    "month",
    "sales_lag_1",
    "sales_lag_4",
    "sales_lag_52",
    "promo_lag_1",
    "promo_lag_52",
    "sales_roll_mean_4",
    "sales_roll_mean_13",
    "promotion_plan_proxy",
]

lgb_sample = lgb_sample.dropna(subset=feature_cols)

X_train = lgb_sample.loc[lgb_sample["week"] <= "2012-07-30", feature_cols]
y_train = lgb_sample.loc[lgb_sample["week"] <= "2012-07-30", "sales"]

X_test = lgb_sample.loc[lgb_sample["week"] == "2012-08-06", feature_cols]
y_test = lgb_sample.loc[lgb_sample["week"] == "2012-08-06", "sales"]

lgb_model = LGBMRegressor(
    objective="regression",
    n_estimators=500,
    learning_rate=0.05,
    num_leaves=31,
    min_child_samples=20,
    subsample=0.8,
    colsample_bytree=0.8,
    reg_alpha=0.0,
    reg_lambda=1.0,
    random_state=42,
)

lgb_model.fit(X_train, y_train)
y_pred = lgb_model.predict(X_test)

pd.DataFrame({"y_test": y_test.values, "y_pred": y_pred})

同じ特徴量テーブルは、XGBoost にもそのまま渡せます。

from xgboost import XGBRegressor

xgb_model = XGBRegressor(
    objective="reg:squarederror",
    n_estimators=500,
    learning_rate=0.05,
    max_depth=6,
    min_child_weight=1,
    subsample=0.8,
    colsample_bytree=0.8,
    reg_alpha=0.0,
    reg_lambda=1.0,
    random_state=42,
)

xgb_model.fit(X_train, y_train)

複数系列をまとめるときの注意

複数部門や複数店舗をまとめる場合は、系列ごとにラグをずらす必要があります。ここを全体で一括 shift() すると、別系列の値が混ざります。

lgb_all = df_all.sort_values(["store", "dept", "week"]).copy()

lgb_all["sales_lag_1"] = lgb_all.groupby(["store", "dept"])["sales"].shift(1)
lgb_all["promo_lag_1"] = lgb_all.groupby(["store", "dept"])["promotion_sales"].shift(1)

1 店舗だけに絞っているなら groupby(["dept"]) でも足りますが、一般形は groupby(["store", "dept"]) です。

木系モデルの主要パラメータ

パラメータ	LightGBM	XGBoost	役割	調整の勘所
学習率	`learning_rate`	`learning_rate` / `eta`	1 本ごとの更新量	小さくするほど安定するが木本数が必要
木の本数	`n_estimators`	`n_estimators`	学習ステップ数	学習率を下げたら増やす
木の複雑さ	`num_leaves`, `max_depth`	`max_depth`	非線形の強さ	過学習しやすければ先に抑える
子ノード条件	`min_child_samples`	`min_child_weight`	細かい分割をどこまで許すか	小売のノイズに過剰適合するとき上げる
行サンプリング	`subsample`	`subsample`	各木で使う行割合	過学習抑制に効く
列サンプリング	`colsample_bytree`	`colsample_bytree`	各木で使う列割合	特徴量が多いとき有効
正則化	`reg_alpha`, `reg_lambda`	`reg_alpha`, `reg_lambda`	複雑さ抑制	不安定なら強める
目的関数	`objective`	`objective`	最適化する損失	件数や非負値なら `poisson` や `tweedie` も候補

木系モデルの設定原則と調整方法

まず learning_rate を小さめにし、n_estimators を増やす
過学習しやすいときは num_leaves / max_depth を先に抑える
subsample と colsample_bytree は 0.7 - 0.9 から試すと扱いやすい
価格、販促、在庫、気温の相互作用が重要なら木系モデルが有利
カウントに近い非負データなら poisson や tweedie を検討する
早期終了は必ず時系列順の検証期間で行う

実務では、ラグ特徴量を固定したうえで、まず num_leaves、min_child_samples、learning_rate の 3 つから比較すると進めやすくなります。

import lightgbm as lgb


def wape(y_true, y_pred, eps=1e-8):
    y_true = np.asarray(y_true)
    y_pred = np.asarray(y_pred)
    return 100 * np.sum(np.abs(y_true - y_pred)) / np.maximum(np.sum(np.abs(y_true)), eps)


tune_candidates = [
    {"num_leaves": 15, "min_child_samples": 40},
    {"num_leaves": 31, "min_child_samples": 20},
    {"num_leaves": 63, "min_child_samples": 10},
]

valid_start = pd.to_datetime("2012-07-02")
valid_end = pd.to_datetime("2012-07-30")

train_mask = lgb_sample["week"] < valid_start
valid_mask = (lgb_sample["week"] >= valid_start) & (lgb_sample["week"] <= valid_end)

X_train = lgb_sample.loc[train_mask, feature_cols]
y_train = lgb_sample.loc[train_mask, "sales"]
X_valid = lgb_sample.loc[valid_mask, feature_cols]
y_valid = lgb_sample.loc[valid_mask, "sales"]

lgb_cv_rows = []

for candidate in tune_candidates:
    model = LGBMRegressor(
        objective="regression",
        n_estimators=2000,
        learning_rate=0.03,
        num_leaves=candidate["num_leaves"],
        min_child_samples=candidate["min_child_samples"],
        subsample=0.8,
        colsample_bytree=0.8,
        random_state=42,
    )
    model.fit(
        X_train,
        y_train,
        eval_set=[(X_valid, y_valid)],
        eval_metric="l1",
        callbacks=[lgb.early_stopping(100, verbose=False)],
    )

    pred = model.predict(X_valid, num_iteration=model.best_iteration_)
    lgb_cv_rows.append(
        {
            **candidate,
            "best_iteration": model.best_iteration_,
            "valid_wape": wape(y_valid, pred),
        }
    )

pd.DataFrame(lgb_cv_rows).sort_values("valid_wape")

木系モデルが向く場面

系列数が多い
外生変数が多い
非線形な影響を拾いたい
店舗属性や商品属性を一緒に使いたい
Prophet より柔軟に相互作用を扱いたい

一方で、木系モデルは特徴量設計が仕事の中心になります。モデル自体より、どの変数を未来時点で安全に使えるかの設計が精度を左右します。

モデル選択の実務上の目安

ここまで出てきた手法を実務的に整理すると、次のようになります。

手法	向いている場面	強み	主な注意点
STL	まず系列構造を見たい	トレンドと季節性を切り分けやすい	予測モデルではない
ETS	比較的短い単系列	ベースラインとして速い	外生変数を直接入れにくい
SARIMA / SARIMAX	単系列で自己相関と季節性を明示的に使いたい	統計的に解釈しやすい	季節周期や次数選択に慣れが要る
Prophet	トレンド・季節性・祝日・イベントを説明しながら予測したい	可視化しやすく説明しやすい	追加回帰変数は線形、未来既知が前提
XGBoost / LightGBM	多系列、外部変数多め、非線形あり	拡張しやすく精度が出やすい	特徴量設計と漏洩管理が必要

どれか 1 つだけが常に正解というわけではありません。実務では、

まず系列を理解する
シンプルなベースラインを作る
外生変数や系列数に応じてモデルを広げる
必要ならアンサンブルする

という順番で考えると進めやすくなります。

評価設計: ランダム分割を避け、ローリングオリジンで評価する

時系列では、通常の表形式タスクのようなランダム分割を使いません。ランダム分割では未来情報が訓練に混ざり、過大評価になります。

イメージとしては次の流れです。

ある時点 t_1 までを学習に使い、その先 h 期を検証する
次に起点を t_2 へ進め、再びその時点までで学習し、その先 h 期を検証する
これを複数回繰り返し、各 fold の誤差を平均してモデルを比べる

同じローリングオリジンでも、学習期間を回ごとに広げていく expanding window と、直近一定期間だけを使う sliding window があります。小売の週次予測では、過去をなるべく多く使いたいことが多いため、まずは expanding window から始めると進めやすいです。

原則は次の 3 つです。

訓練期間は検証期間より前である
予測ホライズンに合わせて切り方を決める
本番で使う情報だけで評価する

もっとも基本的なのはローリングオリジン検証です。

df_eval = df_sample.reset_index().copy()

cutoffs = pd.to_datetime([
    "2011-12-26",
    "2012-03-26",
    "2012-07-30",
])

h = 4

for cutoff in cutoffs:
    train = df_eval[df_eval["week"] <= cutoff]
    valid = df_eval[
        (df_eval["week"] > cutoff)
        & (df_eval["week"] <= cutoff + pd.Timedelta(weeks=h))
    ]

    print(cutoff.date(), train.shape, valid.shape)

週次小売では、「1 週先は当たるが 4 週先で崩れる」ということが普通に起こります。したがって、1 つの指標だけでなく、ホライズンごとの誤差を見る方が実務的です。

誤差指標: 何を良くしたいのかを先に決める

モデル改善とは、結局どの損失を下げるかを決めることです。MAPE を下げたいのか、WAPE を下げたいのか、欠品率を下げたいのかで、最適なモデルや特徴量が変わります。

代表的な指標は次のとおりです。

MAE

$$\mathrm{MAE} = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat y_i|$$

絶対誤差の平均です。外れ値に比較的頑健で、金額ベースで誤差を理解しやすいのが利点です。

RMSE

$$\mathrm{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat y_i)^2}$$

大きな誤差をより強く罰します。外れ値や大型売場の誤差を重く見たいときに向きます。

MAPE

$$\mathrm{MAPE} = \frac{100}{n} \sum_{i=1}^{n} \left|\frac{y_i - \hat y_i}{y_i}\right|$$

割合で解釈しやすい一方、真値が 0 に近いと不安定です。低販売 SKU では過大に悪化しやすいので注意が必要です。

WAPE

$$\mathrm{WAPE} = 100 \times \frac{\sum_{i=1}^{n} |y_i - \hat y_i|}{\sum_{i=1}^{n} |y_i|}$$

全体売上に対して誤差がどれだけあるかを見ます。小売では MAPE より安定することが多く、実務で非常によく使われます。WMAPE と呼ばれることも多いです。

予測基準バイアス率

$$\mathrm{PBR} = 100 \times \frac{\sum_{i=1}^{n} (\hat y_i - y_i)}{\sum_{i=1}^{n} |\hat y_i|}$$

このような予測値基準のバイアス率です。標準的な教科書用語として完全に固定されているわけではありませんが、実務では過大予測か過小予測かを把握する補助指標として使えます。値が正なら過大予測、負なら過小予測です。予測値が 0 に近い系列では不安定になりやすいので、WAPE や MAE と併用するのが安全です。

sMAPE

$$\mathrm{sMAPE} = \frac{100}{n} \sum_{i=1}^{n} \frac{2|y_i - \hat y_i|}{|y_i| + |\hat y_i|}$$

MAPE より分母の不安定さをやや抑えますが、解釈の直感性は WAPE よりやや落ちます。

MASE

$$\mathrm{MASE} = \frac{\dfrac{1}{n} \sum_{t=1}^{n} |y_t - \hat y_t|}{\dfrac{1}{n-m} \sum_{t=m+1}^{n} |y_t - y_{t-m}|}$$

単純な季節ナイーブ予測と比べて、どれだけ改善できたかを見る尺度です。$m$ は季節周期で、週次データで年周期を見るなら通常は $m=52$ です。間欠需要や系列比較で役立ちます。

小売でよくある指標の使い分け

指標	向いている場面	長所	弱点
MAE	金額誤差をそのまま見たい	解釈しやすい	規模の違う系列を比較しにくい
RMSE	大外しを強く罰したい	大誤差に敏感	外れ値に引っ張られる
MAPE	業務側へ割合で説明したい	直感的	0 近辺で不安定
WAPE / WMAPE	小売全体の誤差を安定して見たい	実務で使いやすい	小さい系列の重要度が薄まりやすい
予測基準バイアス率	過大予測と過小予測の偏りを見たい	符号で偏りの方向が分かる	予測値が小さい系列では不安定
sMAPE	MAPE を少し安定化したい	比率比較しやすい	やや直感的でない
MASE	ナイーブ予測より良いか見たい	系列間比較に強い	説明には少し慣れが要る

在庫最適化や安全在庫の議論まで行くなら、点予測誤差だけでなく、

予測バイアス
サービスレベル
在庫切れ率
ピンボール損失

まで見ることが多いです。売上予測の最終目的が「補充判断」なら、点予測だけで完結しない点も押さえておくとよいです。

簡単な指標計算関数は次のように書けます。

def regression_metrics(y_true, y_pred, eps=1e-8):
    y_true = np.asarray(y_true)
    y_pred = np.asarray(y_pred)

    mae = mean_absolute_error(y_true, y_pred)
    rmse = np.sqrt(mean_squared_error(y_true, y_pred))
    mape = np.mean(np.abs((y_true - y_pred) / np.maximum(np.abs(y_true), eps))) * 100
    wape = np.sum(np.abs(y_true - y_pred)) / np.maximum(np.sum(np.abs(y_true)), eps) * 100
    pbr = np.sum(y_pred - y_true) / np.maximum(np.sum(np.abs(y_pred)), eps) * 100
    smape = np.mean(
        2 * np.abs(y_true - y_pred) / np.maximum(np.abs(y_true) + np.abs(y_pred), eps)
    ) * 100

    return {
        "MAE": mae,
        "RMSE": rmse,
        "MAPE(%)": mape,
        "WAPE(%)": wape,
        "PBR(%)": pbr,
        "sMAPE(%)": smape,
    }

日本の小売で見るときの誤差レンジの目安

ここで示す数値は、あくまで週次予測のざっくりした目安です。業態、粒度、ホライズン、プロモーション比率、欠品の多さで大きく変わるので、比較用の目安として読んでください。

また、小売では MAPE より WAPE の方が安定することが多いため、以下は主に WAPE の目安で考えます。

予測レベル	良好の目安	実務でよくある範囲	コメント
全社週次売上	`3%` - `8%`	`5%` - `12%`	大きく集約されるので当たりやすい
店舗 × 部門週次	`8%` - `15%`	`10%` - `25%`	店舗差と部門差が効き始める
店舗 × SKU 週次の定番品	`15%` - `25%`	`20%` - `35%`	安定品でも SKU 粒度ではばらつく
店舗 × SKU 週次の生鮮・販促品	`25%` - `40%`	`30%` - `50%` 以上もある	天候、欠品、イベントの影響が強い
低回転 SKU	WAPE 単独では評価しにくい	MASE やサービスレベルを見る	0 売上週が多く、MAPE 系は不安定

「WAPE が 5% 未満なら必ず優秀」「20% 超なら必ず失敗」といった単純な話ではありません。たとえば、生鮮や OTC 医薬品の流行商品は 20% 台でも実用的な場合があり、逆に定番日用品で 20% 台なら改善余地が大きいことがあります。

スーパー・ドラッグストアで予測しやすい商品と難しい商品

ここは実務感覚が非常に大切です。まず結論を言うと、売れ数が少ない商品は一般に予測しやすいどころか、むしろ難しいことが多いです。絶対誤差は小さく見えても、需要発生自体がまばらなため、MAPE 系指標は大きくぶれやすく、間欠需要として別の扱いが必要になります。

日本のスーパーやドラッグストアを例にすると、おおまかには次のように整理できます。

カテゴリ	例	予測しやすさ	主な理由	改善の勘所
定番日用品	ティッシュ、トイレットペーパー、洗剤、歯磨き粉、おむつ	高い	ベース需要が比較的安定している	曜日、月初、特売、価格指数を入れる
定番食品	米、卵、牛乳、食パン、豆腐、納豆	比較的高い	リピート需要が強く週次パターンが出やすい	祝日、気温、価格、欠品補正を入れる
保存食品	カップ麺、レトルト、冷凍食品、飲料ケース	中程度	まとめ買い、特売、天候の影響がある	価格弾力性、前回特売からの日数、在庫を入れる
季節商材	殺虫剤、日焼け止め、カイロ、花粉対策、鍋つゆ	中程度から難しい	季節開始と終了の立ち上がりが急	気温、湿度、降水量、前年同週、季節開始フラグを入れる
OTC 医薬品	風邪薬、咳止め、解熱剤、マスク、整腸剤	難しい	流行、ニュース、気温急変、感染症拡大で跳ねる	流行指数、学校休み、検索量、気象を入れる
青果	葉物、トマト、きゅうり、いちご、みかん	難しい	天候、相場、鮮度、代替関係が強い	気温、降水量、市場価格、入荷量、イベントを入れる
鮮魚・精肉	刺身、切り身、焼肉用、鍋用、惣菜素材	難しい	生鮮特有の週末需要、気温、行事、廃棄制約	気温、曜日、週末、イベント、値引き、在庫を入れる
惣菜・弁当	唐揚げ、弁当、サラダ、揚げ物	かなり難しい	時間帯、天候、来店客数、欠品、廃棄の影響が強い	時間帯別来客、降雨、気温、陳列数、廃棄実績を入れる
化粧品・美容	スキンケア、ヘアケア、限定コフレ	中程度から難しい	キャンペーン、SNS、季節イベントの影響が大きい	キャンペーンフラグ、レビュー、広告投下、連休を入れる
低回転 SKU	高単価サプリ、地域限定商品、輸入食品	低い	間欠需要で 0 売上週が多い	SKU を束ねる、Croston 系、MASE 評価、カテゴリ予測を使う

カテゴリ別に見るときの注意

野菜は「日常的に使う」一方で、価格や天候の影響が大きいため、必ずしも予測しやすくありません
殺虫剤やカイロのような季節商材は、季節性がはっきりしていれば当てやすくなる一方、立ち上がり週を外すと大きく崩れます
低回転商品は、売上絶対額が小さくても予測が簡単とは限りません。間欠需要として別設計が必要です

小売で効く外部要因とデータソース

モデルよりも「未来時点で使える説明変数をどこまで作れるか」が精度を左右することがよくあります。特に小売では、売上そのもの以外の情報が効きやすいです。

外部要因は「未来に分かるか」で 4 分類する

完全に未来既知: 曜日、祝日、月末、価格計画、自社販促計画
予報を使える: 気温、降水量、積雪などの気象
予定は分かるが強度は不確実: 地域イベント、コンサート、学校行事
未来が本質的に不明: 競合の直前値引き、突発欠品、急激な景気ショック

Prophet でも木系モデルでも、強い特徴量になるのは主に 1 と 2、条件付きで 3 です。4 は実績値をそのまま使わず、シナリオ、代理変数、別モデル予測で扱います。

小売で使う代表的な外部要因

要因	代表データソース	代表特徴量	Prophet への入れ方	XGBoost / LightGBM への入れ方	注意点
祝日・カレンダー	祝日カレンダー、学校休暇、給料日	祝日フラグ、連休長、月初月末、給料日前後	`holidays`、`add_country_holidays()`	ダミー、前後日数、周期エンコード	未来既知なので最優先
気象	気象庁、民間 API、店舗最寄り観測所	気温、降水量、積雪、平年差	`add_regressor()`	数値特徴量、ラグ、前年差、カテゴリとの相互作用	予報値を使う
自社販促	チラシ計画、クーポン配布、値引き計画	販促フラグ、割引率、掲載面数、エンド陳列	`add_regressor()`	数値、カテゴリ、交互作用	実績販促売上は使わない
価格	価格マスタ、価格改定計画	実価格、割引率、カテゴリ平均比、前年差	`add_regressor()`	価格指数、弾力性特徴、前年差	未来価格計画の確定範囲を確認
競合販促	競合チラシ、立地情報、店舗オープン情報	競合値引き指数、近隣競合数、出店フラグ	既知シナリオなら `add_regressor()`	距離加重特徴、イベントフラグ	未来実績は通常未知
在庫・欠品	ERP、WMS、棚在庫、補充ログ	在庫量、欠品フラグ、在庫消化日数	線形回帰変数としては可能	強い数値特徴量、交互作用	売上と需要の混同に注意
地域イベント	祭り、花火、球場イベント、学校行事	開催フラグ、来場予定数、距離	`holidays` または `add_regressor()`	バイナリ、強度、距離加重	開催有無と集客強度を分ける
人口・人流	国勢調査、モバイル位置情報、観光統計	昼間人口、人流データ、訪日客指数	`add_regressor()`	水準、前年差、移動平均、店舗属性との交互作用	週次粒度へ合わせる必要がある
マクロ経済	GDP、失業率、CPI、ガソリン価格、住宅価格	失業率、可処分所得の代替指標、物価指数	`add_regressor()`	低頻度系列を週次へ結合して使う	公表ラグを必ず考慮する
政策・技術変更	増税、補助金、決済キャンペーン、アプリ刷新	制度変更フラグ、アプリ導入後、配送開始後	`holidays` または `add_regressor()`	構造変化ダミー、前後トレンド、交互作用	構造変化として扱う

Prophet へ外部要因を入れるときの考え方

休日やイベントは holidays が自然
数値の外生変数は add_regressor() が自然
回帰部分は基本的に線形なので、強い非線形相互作用は表現しにくい
未来の値が用意できない変数は入れても本番で使えない

たとえば、天候・販促・価格計画があるなら、次のような考え方になります。

avg_temp_forecast
rainfall_forecast
promo_flag_plan
discount_rate_plan
price_index_plan

を add_regressor() で入れる、という形です。

XGBoost / LightGBM へ外部要因を入れるときの考え方

木系モデルでは、上の変数をそのまま特徴量列へ追加できます。さらに、次のような加工がしやすいのが強みです。

気温 × 季節商材
価格差 × 競合距離
販促フラグ × 曜日
在庫量 × 週末
人流データ × 観光地店舗フラグ

たとえば feature list としては、次のような形が典型です。

avg_temp_forecast
rainfall_forecast
promo_flag_plan
discount_rate_plan
price_index_plan
inventory_on_hand
stockout_flag_lag_1
competitor_campaign_index
foot_traffic_data
unemployment_rate
tax_change_flag
app_launch_after_flag

この柔軟さが、木系モデルが小売の多系列予測で強い理由です。

以下は、外部データを週次表へ結合し、木系モデル用の特徴量へ落とし込む基本例です。ここでも教材データだけで動くよう、説明用の外部特徴量を週番号から生成しています。

weekly_calendar = lgb_sample[["week"]].drop_duplicates().copy()
weekly_calendar["week_of_year"] = weekly_calendar["week"].dt.isocalendar().week.astype(int)
weekly_calendar["month"] = weekly_calendar["week"].dt.month
weekly_calendar["promo_flag_plan"] = weekly_calendar["week_of_year"].isin([1, 13, 26, 39, 48]).astype(int)
weekly_calendar["discount_rate_plan"] = weekly_calendar["promo_flag_plan"] * 0.10
weekly_calendar["avg_temp_forecast"] = 15 + 10 * np.sin(2 * np.pi * weekly_calendar["week_of_year"] / 52)
weekly_calendar["rainfall_forecast"] = 5 + 3 * np.cos(2 * np.pi * weekly_calendar["week_of_year"] / 52)
weekly_calendar["foot_traffic_data"] = (
    100
    + 12 * np.sin(2 * np.pi * weekly_calendar["week_of_year"] / 52)
    + 4 * weekly_calendar["month"].isin([12, 1]).astype(int)
)
weekly_calendar["inventory_on_hand"] = 400 + 40 * np.cos(2 * np.pi * weekly_calendar["week_of_year"] / 52)

external_features = weekly_calendar[
    [
        "week",
        "promo_flag_plan",
        "discount_rate_plan",
        "avg_temp_forecast",
        "rainfall_forecast",
        "foot_traffic_data",
        "inventory_on_hand",
    ]
].copy()

lgb_feature_table = (
    lgb_sample[["week", "sales", "promotion_sales"]].copy()
    .merge(external_features, on="week", how="left", validate="one_to_one")
    .assign(
        week_of_year=lambda d: d["week"].dt.isocalendar().week.astype(int),
        sales_lag_1=lambda d: d["sales"].shift(1),
        sales_lag_52=lambda d: d["sales"].shift(52),
        sales_roll_mean_4=lambda d: d["sales"].shift(1).rolling(4).mean(),
        week_sin=lambda d: np.sin(2 * np.pi * d["week_of_year"] / 52),
        week_cos=lambda d: np.cos(2 * np.pi * d["week_of_year"] / 52),
        discount_x_foot_traffic=lambda d: d["discount_rate_plan"] * d["foot_traffic_data"],
        temp_x_promo=lambda d: d["avg_temp_forecast"] * d["promo_flag_plan"],
    )
)

feature_cols = [
    "sales_lag_1",
    "sales_lag_52",
    "sales_roll_mean_4",
    "promo_flag_plan",
    "discount_rate_plan",
    "avg_temp_forecast",
    "rainfall_forecast",
    "foot_traffic_data",
    "inventory_on_hand",
    "week_sin",
    "week_cos",
    "discount_x_foot_traffic",
    "temp_x_promo",
]

lgb_feature_table = lgb_feature_table.dropna(subset=feature_cols + ["sales"])
lgb_feature_table[feature_cols + ["sales"]].head()

特徴量テーブルを作っただけでは予測モデルにはなりません。次のように学習期間と予測対象週を切り分ければ、追加した特徴量をそのままモデルへ渡せます。

X_train_ext = lgb_feature_table.loc[
    lgb_feature_table["week"] <= "2012-07-30",
    feature_cols,
]
y_train_ext = lgb_feature_table.loc[
    lgb_feature_table["week"] <= "2012-07-30",
    "sales",
]

X_test_ext = lgb_feature_table.loc[
    lgb_feature_table["week"] == "2012-08-06",
    feature_cols,
]
y_test_ext = lgb_feature_table.loc[
    lgb_feature_table["week"] == "2012-08-06",
    "sales",
]

lgb_external_model = LGBMRegressor(
    objective="regression",
    n_estimators=500,
    learning_rate=0.05,
    num_leaves=31,
    min_child_samples=20,
    subsample=0.8,
    colsample_bytree=0.8,
    random_state=42,
)

lgb_external_model.fit(X_train_ext, y_train_ext)
external_pred = lgb_external_model.predict(X_test_ext)

pd.DataFrame({"y_test": y_test_ext.values, "y_pred": external_pred})

精度を高める方法: チューニングより先に見直すこと

「精度を上げたい」という話になると、すぐにハイパーパラメータチューニングへ意識が向きがちです。しかし小売では、その前段で効く改善がたくさんあります。

1. まず売上と需要を区別する

欠品があるとき、観測売上は潜在需要を切り取った値になります。単純化すると、

$$\mathrm{sales}(t) = \min\bigl(\mathrm{demand}(t),\ \mathrm{availableStock}(t)\bigr)$$

と考えられます。

このため、欠品が多いカテゴリでは「売上予測」と「需要予測」を切り分けておかないと、発注判断がぶれます。惣菜や生鮮、人気 OTC 商品では特に重要です。

2. 評価指標を SKU 特性に合わせる

定番品なら WAPE や MAE
大外しを避けたいなら RMSE
低回転商品なら MASE や RMSSE
補充判断まで見るならサービスレベルや欠品率

を使い分けます。

3. 未来既知の変数だけで特徴量を作る

これは精度以前の正しさの条件です。価格実績、販促実績、在庫実績をそのまま未来へ入れると、評価だけ高いモデルになります。

forbidden_future_actuals = {
    "promotion_sales",
    "actual_customer_count",
    "actual_stockout_count",
    "realized_price",
}


def assert_future_known(feature_cols, forbidden_cols):
    leaked = sorted(set(feature_cols) & set(forbidden_cols))
    if leaked:
        raise ValueError(f"情報漏洩の疑いがある列: {leaked}")


candidate_feature_cols = [
    "sales_lag_1",
    "sales_lag_52",
    "promo_flag_plan",
    "discount_rate_plan",
    "avg_temp_forecast",
    "foot_traffic_data",
]

assert_future_known(candidate_feature_cols, forbidden_future_actuals)

4. 系列を束ねるか分けるかを見直す

SKU 単位では疎すぎるならカテゴリ単位へ集約する
全商品共通モデルより、定番品 / 季節品 / 低回転品で分けた方がよいことがある
店舗特性が強いなら、商圏タイプ別モデルに分ける

5. ラグだけでなく変化率と相互作用を作る

小売では「値そのもの」より、次が効くことがよくあります。

前年同週差
前回販促からの経過週数
価格前年差
競合との差
在庫量の急減

6. 異常値を放置しない

システム障害
臨時休業
台風や大雪
棚替え
データ欠損

などは、ノイズとして片づけず、「意味のある異常」として扱うべきケースです。削除するのか、フラグ化するのか、別の構造変化として扱うのかを決めます。

df_anomaly = df_sample.copy()
df_anomaly["rolling_median_13"] = df_anomaly["sales"].shift(1).rolling(13, min_periods=4).median()
df_anomaly["rolling_iqr_13"] = (
    df_anomaly["sales"].shift(1).rolling(13, min_periods=4).quantile(0.75)
    - df_anomaly["sales"].shift(1).rolling(13, min_periods=4).quantile(0.25)
)
df_anomaly["anomaly_flag"] = (
    (df_anomaly["sales"] > df_anomaly["rolling_median_13"] + 3 * df_anomaly["rolling_iqr_13"])
    | (df_anomaly["sales"] < df_anomaly["rolling_median_13"] - 3 * df_anomaly["rolling_iqr_13"])
).fillna(False)

df_anomaly.loc[df_anomaly["anomaly_flag"], ["sales"]].head()

7. 間欠需要には別の考え方を使う

低回転 SKU は、通常の回帰モデルだけでは扱いにくいことがあります。次のような方針が有効です。

Croston 系や SBA、TSB を候補にする
需要発生有無と需要量を分けてモデル化する
SKU 単体でなく、上位カテゴリ予測を下ろす
MASE や在庫 KPI を軸に評価する

df_intermittent = df_sample.copy()
df_intermittent["is_positive_sale"] = (df_intermittent["sales"] > 0).astype(int)

zero_ratio = 1 - df_intermittent["is_positive_sale"].mean()
avg_positive_sale = df_intermittent.loc[
    df_intermittent["is_positive_sale"] == 1,
    "sales",
].mean()

pd.Series(
    {
        "zero_week_ratio": zero_ratio,
        "average_positive_sale": avg_positive_sale,
    }
)

8. 集計レベルと整合性を見直す

チェーン全体、地域、店舗、部門、SKU の予測がバラバラだと、経営数値と現場数値が噛み合いません。小売では階層予測と整合化が重要です。

bottom-up
top-down
middle-out
予測整合化

のどれで運用するかを決める必要があります。

9. モデルを分業させる

ベースラインは ETS / SARIMA
構造要因は Prophet
非線形補正は LightGBM / XGBoost

という分業が有効なことがあります。1 モデルで全てを背負わせる必要はありません。

10. 予測値に業務制約を入れる

負値を 0 に切る
最大陳列量や供給上限を反映する
発注単位に丸める
SKU 廃番や新商品導入を別扱いにする

数理モデルの出力を、そのまま業務値として使わない方が安全です。

アンサンブルとモデル融合

単一モデルより、複数モデルの組み合わせの方が安定することがよくあります。考え方は単純で、異なる誤り方をするモデルを組み合わせるということです。

もっとも基本的な重み付き平均は次の形です。

$$\hat y_t^{ens} = \sum_{k=1}^{K} w_k \hat y_t^{(k)},\quad w_k \ge 0,; \sum_{k=1}^{K} w_k = 1$$

小売でよく使われる融合パターンは次のとおりです。

手法	何をするか	向いている場面	注意点
単純平均	複数モデルの平均をとる	まず安定させたい	重み最適化なしでも効くことがある
重み付き平均	CV 成績に応じて重みをつける	モデルごとの強弱がある	同じ検証期間で比較する
セグメント別最良モデル採用	商品群ごとに最良モデルを使う	定番品と季節品で性質が違う	運用が複雑になる
残差学習	先に Prophet などで予測し、残差を GBDT で学習する	構造 + 非線形補正を両取りしたい	残差学習でも漏洩を避ける
スタッキング	複数モデル予測をさらに上位モデルで学習する	データ量が十分ある	out-of-fold 予測で学習しないと漏洩する

残差モデリングの考え方は、たとえば次のように書けます。

$$y_t = \hat y_t^{\mathrm{Prophet}} + r_t$$

$$\hat y_t = \hat y_t^{\mathrm{Prophet}} + \hat r_t^{\mathrm{GBDT}}$$

この形は、小売のように「カレンダー・季節性ははっきりあるが、価格や在庫の非線形も強い」ケースで相性がよいです。

重み付き平均の基本形は次のとおりです。

$$
\hat y_t^{ens} = \sum_{k=1}^{K} w_k \hat y_t^{(k)}, \qquad w_k \ge 0, \quad \sum_{k=1}^{K} w_k = 1
$$

まずは、実装が簡単で効果も確認しやすい重み付き平均から始めるのが安全です。下のコードでは 2012-08-06 の 1 週先予測を例に流れを示します。実務では、重みを同じ週の実績誤差で合わせず、必ず別の検証期間で決めてください。

def wape(y_true, y_pred, eps=1e-8):
    y_true = np.asarray(y_true)
    y_pred = np.asarray(y_pred)
    return 100 * np.sum(np.abs(y_true - y_pred)) / np.maximum(np.sum(np.abs(y_true)), eps)


pred_table = pd.DataFrame(
    {
        "week": [pd.Timestamp("2012-08-06")],
        "y_true": prophet_test["y"].values,
        "ets_pred": ets_result.loc[
            ets_result["week"] == pd.Timestamp("2012-08-06"),
            "pred_sales",
        ].values,
        "sarima_pred": sarima_result.loc[
            sarima_result["week"] == pd.Timestamp("2012-08-06"),
            "pred_sales",
        ].values,
        "prophet_pred": prophet_test["yhat"].values,
        "lgb_pred": y_pred,
    }
)

model_cols = ["ets_pred", "sarima_pred", "prophet_pred", "lgb_pred"]

score_df = pd.DataFrame(
    {
        "model": model_cols,
        "wape": [wape(pred_table["y_true"], pred_table[col]) for col in model_cols],
    }
)
score_df["raw_weight"] = 1 / score_df["wape"]
score_df["weight"] = score_df["raw_weight"] / score_df["raw_weight"].sum()

weight_map = dict(zip(score_df["model"], score_df["weight"]))

pred_table["ensemble_pred"] = sum(
    pred_table[col] * weight_map[col]
    for col in model_cols
)

score_df, wape(pred_table["y_true"], pred_table["ensemble_pred"])

小売の時系列予測で特に難しいところを工程順に整理する

小売の時系列予測は、モデル選択だけの問題ではありません。むしろ難所の多くは、その前段にあります。難しい順に並べると、おおむね次のようになります。

難しさ順位	工程	なぜ難しいか	典型的な失敗	改善の方向
1	課題設定と目的定義	売上最大化なのか、欠品削減なのか、廃棄削減なのかで正解が変わる	指標だけよくても業務効果が出ない	予測の利用先を先に決める
2	データ収集と統合	POS、価格、販促、在庫、気象、競合、会員、人流が別系統	店舗キーや週定義がずれて結合が壊れる	粒度とキーを先に固定する
3	需要と売上の切り分け	欠品や陳列不足で売上が上限制約で切り詰められる	売れなかったと需要が無いを混同する	欠品補正、在庫特徴量、対象定義を見直す
4	特徴量設計	予測時点で既知かどうか、非線形相互作用、季節性が複雑	実績値を未来特徴量へ入れて漏洩する	予測時点で既知の変数だけで設計する
5	評価設計	どのホライズンで何を良しとするかが曖昧	ランダム分割や MAPE 一本化で誤る	rolling CV と複数指標で見る
6	商品・店舗のセグメンテーション	定番品と季節品、都市店と郊外店では性質が違う	全系列に 1 モデルを押し付ける	セグメント別モデルや階層化を使う
7	モデル選択	統計モデル、Prophet、GBDT のどれが効くかは一律でない	モデル名だけで選ぶ	ベースラインから順に比べる
8	学習とハイパーパラメータ調整	調整空間が広い	先に大規模探索して時間を失う	指標とベースラインを固定してから触る
9	アンサンブルと整合化	系列ごとに強いモデルが違う	予測値が階層で不整合になる	重み付け、最良モデル採用、予測整合化を使う
10	運用と監視	店舗改装、競合出店、制度変更で前提が変わる	一度作って終わりにする	drift 監視と定期再学習を組む

見て分かるとおり、モデル学習そのものは全体の一部です。課題設定、データ設計、特徴量設計、評価設計の方が難しいことが多いです。

小売時系列予測の実務フロー

最後に、課題設定から改善までの流れを、実務で迷いにくい順番に整理します。

何のための予測かを決める
対象粒度を決める
売上を予測するのか需要を予測するのか決める
主キー、期間、欠損週、欠品の有無を確認する
カレンダー、価格、販促、在庫、天候などの予測時点で既知の変数を整理する
可視化、ACF / PACF、STL で系列構造を理解する
ETS や SARIMA で基準モデルを作る
Prophet で祝日・イベント・追加回帰変数を確認する
多系列や非線形が強ければ XGBoost / LightGBM を試す
rolling CV と複数誤差指標で評価する
商品群や店舗群で分割して改善点を探す
必要ならアンサンブル、整合化、業務制約を入れる
本番では drift 監視と定期再学習を行う

まとめ

この記事では、小売売上データを例に、基礎的な確認から代表的な予測モデル、誤差指標、特徴量設計、精度改善、アンサンブル、そして小売実務の難所までを一通り見てきました。

特に大事なのは、いきなりモデル選定から始めないことです。まずは次の順で進めると、学習でも実務でもぶれにくくなります。

粒度、主キー、欠損週、欠品の有無を確認する
未来時点で使ってよい情報だけを明確にする
可視化と STL で系列構造を把握する
ETS や SARIMA を基準にする
祝日・イベント・説明性が大事なら Prophet を使う
多系列や非線形が大きいなら XGBoost / LightGBM を使う
指標、特徴量、アンサンブル、階層整合まで含めて改善する

時系列分析は、モデル名を覚えるだけでは身につきにくい分野です。実際に 1 本の系列を見て、季節性、欠損、外部要因、評価指標の意味を手で確かめながら進めると、理解しやすくなります。特に小売では、モデルそのものよりも「何を予測したいのか」「未来に何が分かるのか」を定義する力が精度を大きく左右します。