階層的時系列

上級

04.16.1

階層的時系列

最終更新 2020-01-29 読了時間 3 分

まとめ

階層的時系列（HTS）は、全社売上 → 事業部売上 → 製品売上のように「下位の合計が上位になる」構造を持つ複数の系列をまとめて予測する枠組みです。各レベルを別々に ARIMA や ETS で予測すると、製品 A1 と A2 の予測の合計が事業部 A の予測と一致しないので、整合性を保つ調整が必要になります。このページでは、商業動態統計の都道府県 × 商品別販売額を例に階層構造のイメージをつかみ、ボトムアップ・トップダウン・最適調整の 3 つの調整法を比べます。そのうえで製品 4 系列 × 24 か月の合成データを作り、ExponentialSmoothing(trend="add") で末端を 6 か月先まで予測して足し上げるボトムアップ予測を実装し、np.allclose で「全社 = 事業部 A + 事業部 B」が成り立つことを確認します。

ARIMA モデルなど単一系列の予測手法を 1 つ知っていれば十分で、このページではその予測を「複数の系列に当てて合計を合わせる」部分に焦点を当てます。

階層的時系列とは #

時系列データが複数あるとき、そのデータが「下位の合計が上位になる」階層的な構造を持つことがあります。47 都道府県の売上を足すと全国の売上になり、商品別の売上を足すと店舗の売上になる、という関係です。

例：地域・カテゴリごとの販売データ #

経済産業省のサイトでは『商業動態統計』が定期的に公開されています。これは日本の商業を営む企業の販売活動の動向をチェックすることを目的として集められているデータで、たとえば 2022 年 1 月の付表3 百貨店・スーパー販売都道府県別、商品別販売額等には以下のようなデータが含まれています。

百貨店・スーパー販売都道府県別、商品別販売額等のデータ例

このデータには以下のような階層構造が含まれています。全国の販売額は 47 都道府県の合計であり、各都道府県の販売額は衣料品・飲食料品などの商品別の合計でもあるので、地域の軸と商品の軸という 2 つの階層が交差しています。

階層構造の図

階層構造の具体例 #

企業の売上データを例にすると、次のような 3 層の階層が考えられます。

全社売上
├── 事業部A
│   ├── 製品A-1
│   └── 製品A-2
└── 事業部B
    ├── 製品B-1
    └── 製品B-2

このとき、各レベルの時系列には整合性の制約があり、全社売上は事業部 A と事業部 B の売上の和、事業部 A の売上は製品 A-1 と A-2 の売上の和でなければなりません。観測値では定義上この等式が成り立ちますが、7 本の系列を別々に ARIMA や ExponentialSmoothing で予測すると、製品 A-1 と A-2 の予測の合計が事業部 A の予測と一致せず、この整合性が崩れます。

末端 4 系列（製品 A1・A2・B1・B2）を観測足し上げて事業部 A・B と全社を作る各レベルを別々に予測すると合計が合わない調整法（ボトムアップ / トップダウン / 最適調整）で整合させる全レベルで足し算が成り立つ予測

予測の調整手法 #

階層的時系列の予測には、整合性を保つために 3 つの代表的なアプローチがあります。

ボトムアップ（Bottom-Up） #

もっとも細かいレベル（製品単位など）で個別に予測し、上位レベルはその合算で求めます。製品 A-1 の季節性や新製品の立ち上がりといった個別の特性を予測に反映しやすいのが利点で、合計が定義上一致するので調整の計算も要りません。欠点は、末端の系列ほどノイズが大きいことで、製品 1 つの売上が月 10 件程度しかないような系列では個別の予測が不安定になり、集約しても精度が出にくくなります。

トップダウン（Top-Down） #

全社レベルで予測し、過去の構成比率（たとえば直近 12 か月の平均シェア）で下位レベルに配分します。全社売上のような上位の系列はノイズが打ち消し合って安定していることが多く、1 本の予測に集中できるので予測精度が出やすいのが利点です。欠点は、構成比率が変わる場面に弱いことで、新製品の投入や既存製品の販売終了でシェアが月ごとに動くと、過去の比率で配分した下位の予測は実態からずれます。

最適調整（Optimal Reconciliation） #

全レベルで独立に予測し、最小分散基準などで整合性を満たすように事後調整します。MinT（Minimum Trace）と呼ばれる方法では、各レベルの予測誤差の共分散行列 \(W\) と集約行列 \(S\) を使って、整合した予測 \(\tilde{y} = S (S^\top W^{-1} S)^{-1} S^\top W^{-1} \hat{y}\) を求めます。各レベルの情報を最大限に使えるので、ボトムアップとトップダウンの中間で精度と整合性を両立できるのが利点です。欠点は、共分散行列 \(W\) の推定が必要で、系列数が 100 を超えると実装も計算も複雑になることです。

Pythonでの実装例 #

scikit-hts や hierarchicalforecast などのライブラリで階層的予測を実装できます。ここではライブラリを使わず、階層の定義・データ生成・ボトムアップ予測の 3 ステップを素の pandas と statsmodels で書いて仕組みを確かめます。

1
2
3
4
5
6
# 階層構造の定義例
hierarchy = {
    "全社": ["事業部A", "事業部B"],
    "事業部A": ["製品A1", "製品A2"],
    "事業部B": ["製品B1", "製品B2"],
}

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

2 行目親をキー、子のリストを値にした辞書で階層を表す
hierarchy は「どの系列がどの系列の合計か」を表す定義で、scikit-hts もこの形式（親 → 子のリスト）で階層を受け取ります。この辞書から集約行列 \(S\) を機械的に作れます。
3 行目最上位の全社は 2 つの事業部の和
"全社": ["事業部A", "事業部B"] は全社 = 事業部 A + 事業部 B という制約です。3 層構造では、この行の子がさらに次の行のキーになっている点が階層のつながりを示します。
4-5 行目末端の製品 4 つが予測の出発点
製品A1〜製品B2 は辞書のキーに現れない（子を持たない）ので末端です。ボトムアップではこの 4 系列だけを予測し、上の 3 系列（事業部 A・B・全社）はこの辞書に沿って足し上げます。

次のコードは、末端の製品 4 系列をランダムウォークで 24 か月分作り、hierarchy の定義どおりに足し上げて事業部と全社の系列を作ります。観測の段階では「合計が一致する」ことが定義上保証されています。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
import pandas as pd

# 人工データの生成
np.random.seed(42)
dates = pd.date_range("2020-01-01", periods=24, freq="MS")

# 末端レベルの系列を生成
a1 = 100 + np.cumsum(np.random.randn(24) * 5)
a2 = 80 + np.cumsum(np.random.randn(24) * 4)
b1 = 120 + np.cumsum(np.random.randn(24) * 6)
b2 = 60 + np.cumsum(np.random.randn(24) * 3)

df = pd.DataFrame({
    "製品A1": a1, "製品A2": a2,
    "製品B1": b1, "製品B2": b2,
}, index=dates)

# 上位レベルを集約で計算
df["事業部A"] = df["製品A1"] + df["製品A2"]
df["事業部B"] = df["製品B1"] + df["製品B2"]
df["全社"] = df["事業部A"] + df["事業部B"]

print(df[["全社", "事業部A", "事業部B"]].head())

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

5-6 行目種を固定して月初の 24 か月を作る
np.random.seed(42) で乱数を固定し、freq="MS" で 2020 年 1 月から 24 か月の月初の日付を作ります。種を固定しないと末端の系列が実行ごとに変わり、後の print の数値が本文と合いません。
9-12 行目末端はランダムウォークで、水準とノイズを製品ごとに変える
np.cumsum(np.random.randn(24) * 5) は標準偏差 5 の変化を積み上げたランダムウォークで、製品 A1 は水準 100・標準偏差 5、B2 は水準 60・標準偏差 3 のように製品ごとに違います。ランダムウォークは単純な平均では当たらないので、次のコードでトレンド付きの指数平滑を使う理由になります。
14-17 行目末端 4 列を 1 つの DataFrame に
日付 index に製品 4 列を並べます。この後、上位レベルも同じ DataFrame の列として足していくので、7 系列が 1 枚の表で揃います。
20-22 行目上位レベルは足し算で作る
hierarchy の定義どおりに 事業部A = 製品A1 + 製品A2、全社 = 事業部A + 事業部B と列を足して作ります。観測データでは全社は必ず事業部の合計に一致し、これが予測でも保ちたい制約です。
24 行目上位 3 系列の先頭 5 か月を確認
print で全社・事業部 A・事業部 B の最初の 5 行を表示し、全社の値が事業部 2 つの和（たとえば 1 月なら 180 + 180 = 360 前後）になっていることを目で確かめます。

ボトムアップ予測の簡易例 #

末端の 4 系列だけを ExponentialSmoothing(trend="add") で 6 か月先まで予測し、上位は hierarchy に沿って足し上げます。最後に全社の予測が事業部 A と B の予測の和と一致するかを np.allclose で検証します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
from statsmodels.tsa.holtwinters import ExponentialSmoothing

# 末端レベルで個別に予測
forecasts = {}
for col in ["製品A1", "製品A2", "製品B1", "製品B2"]:
    model = ExponentialSmoothing(
        df[col], trend="add", seasonal=None,
    ).fit()
    forecasts[col] = model.forecast(6)

# 上位レベルは合算
forecasts["事業部A"] = forecasts["製品A1"] + forecasts["製品A2"]
forecasts["事業部B"] = forecasts["製品B1"] + forecasts["製品B2"]
forecasts["全社"] = forecasts["事業部A"] + forecasts["事業部B"]

# 整合性の確認
print("全社 == 事業部A + 事業部B:")
print(np.allclose(
    forecasts["全社"],
    forecasts["事業部A"] + forecasts["事業部B"],
))

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

5-9 行目末端 4 系列だけを個別に予測
ループは製品 4 列だけを回し、事業部や全社は予測しません。ボトムアップの定義どおり、予測するのは最下層だけで、上位は合算で決めます。forecast(6) は 6 か月先までの予測を日付 index 付きの Series で返します。
6-8 行目trend="add" の Holt 法を選ぶ理由
末端はランダムウォークで作ったので直近の水準と傾きを追う Holt 法（trend="add"、seasonal=None）が妥当です。seasonal=None にしているのは 24 か月のデータでは 12 個の季節係数を安定して推定できないためで、季節性のある実データで 3 年以上あれば seasonal="add", seasonal_periods=12 に変えます。
12-14 行目上位は予測の足し算で作る
forecasts["事業部A"] = forecasts["製品A1"] + forecasts["製品A2"] のように、観測を作ったときと同じ式を予測に適用します。Series どうしの足し算なので 6 か月分の日付 index が揃っていることが前提で、ここでは全系列が同じ forecast(6) なので問題ありません。
17-21 行目np.allclose で整合性を検証
np.allclose は浮動小数の丸め誤差を許して 2 つの配列が一致するかを返します。ボトムアップでは全社を事業部の和として作っているので必ず True になりますが、全社を別のモデルで独立に予測して同じ検証をすると False になり、それが「調整が必要になる」状況そのものです。

ボトムアップでは定義上、整合性が常に保たれ、np.allclose は True を返します。トップダウンや各レベルの独立予測では、全社の予測と事業部の予測の和がずれるので、構成比での配分や MinT のような事後調整が必要です。

使いどころと落とし穴 #

ボトムアップが向くのは、製品や店舗ごとの個別の動き（新製品の立ち上がりや特定店舗の改装）を予測に反映したい場面で、末端の系列がある程度の量（月数十件以上）を持っているときです。末端が月数件のまばらな系列ばかりなら、ノイズが積み上がって全社の予測が暴れるので、安定した全社の予測を構成比で配分するトップダウンの方が精度が出ます。両者の欠点を避けたいときは MinT を使いますが、共分散行列 \(W\) の推定に全系列の学習期間内残差が必要で、系列数が 100 を超えると hierarchicalforecast のようなライブラリに任せた方が確実です。もう一つの落とし穴は、負の予測値です。末端で trend="add" の予測が下がり続けて 0 を割ると、合算した上位の予測も実態からずれるので、売上のような非負の系列では対数変換や damped_trend=True で予測が発散しないようにします。

ARIMA モデル — 階層各レベルで使える古典的予測モデル
Prophetを使ってみる — Prophetで時系列を分解して予測
STL 分解とフォーキャスト — STL分解後の残差にARIMAを適用して予測