Holt-Winters

中級

04.12.3

Holt-Winters

最終更新 2020-02-26 読了時間 6 分

まとめ

Holt-Winters 法は、Holt 法の水準 \(\ell_t\) と傾き \(b_t\) に、周期 \(m\) 個の季節成分 \(s_t\) を加えた指数平滑法です。季節成分は平滑係数 \(\gamma\) で「同じ月の 1 年前の季節値」を更新し、予測では \(\ell_T + h, b_T + s_{T+h-m}\) として該当する月の季節値を足し戻します。このページでは米国の天然ガス消費量（月次 BCF）を 2020 年 1 月で分け、SES・Holt・ExponentialSmoothing(trend="additive", seasonal="add", seasonal_periods=12) の 3 つの 33 か月予測を 1 枚に重ねます。SES の水平線、Holt の直線に対し、Holt-Winters の予測だけが冬 3300 超・夏 2100 前後の山谷をたどることが図で分かります。

Holtの線形指数平滑法で水準と傾きを \(\alpha\)・\(\beta\) で更新する流れを先に見ておくと、Holt-Winters 法はそこに季節値 12 個を \(\gamma\) で更新する処理を足しただけだと分かります。

直感 #

天然ガス消費量は毎年 1 月に 3000 BCF を超え、7〜8 月に 2000 前後へ落ちる 12 か月周期を持ちます。Holt 法はこの系列を「水準 + 傾き」の直線で表すので、冬と夏の差 1000 BCF がすべて誤差になります。Holt-Winters 法は「1 月はいつも水準より +600、7 月は −300」といった月ごとのズレを 12 個の季節値 \(s_t\) として覚え、予測ではその月の季節値を直線に足し戻します。季節値の更新は \(s_t = \gamma(y_t - \ell_{t-1} - b_{t-1}) + (1 - \gamma) s_{t-m}\) で、12 か月前の同じ月の季節値を \(\gamma\) の重みで最新の観測に寄せる、という SES と同じ形です。加法季節（seasonal="add"）は山谷の高さが年によらず一定とみなし、乗法季節（seasonal="mul"）は水準に比例して山谷も大きくなるとみなします。

水準 ℓ_t を α で更新傾き b_t を β で更新季節値 s_t を γ で 12 か月前の値から更新学習最終月の ℓ_T, b_T と 12 個の季節値予測 ℓ_T + h b_T + s_{T+h−12} 33 か月分の山谷を描く

1. ライブラリ #

1
2
3
4
5
6
7
8
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.holtwinters import Holt
from statsmodels.tsa.holtwinters import SimpleExpSmoothing

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

6 行目ExponentialSmoothing は 3 成分すべてを持つ汎用クラス
ExponentialSmoothing は trend と seasonal の引数で成分を付け外しでき、trend=None, seasonal=None なら SES、trend="add", seasonal=None なら Holt と同じになります。このページでは季節まで付けた Holt-Winters として使います。
7-8 行目SES と Holt も読み込んで 3 本を同じ図に
前の 2 ページと同じ SimpleExpSmoothing と Holt を同じ学習データ・同じ \(\alpha\)、\(\beta\) で動かし、Holt-Winters と重ねます。成分を 1 つずつ足したときに予測線の形がどう変わるかを、3 本の線の違いとして見せる狙いです。
2 行目月初の日付 index が季節周期の前提
seasonal_periods=12 は「12 行で 1 周期」という意味なので、DataFrame が月ごとに 1 行、欠けなく並んでいる必要があります。後で asfreq("MS") を付けるのはこのためです。

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：水準 → 傾き → 季節と要素を足すほど、未来への伸ばし方が豊かになる。

サンプルデータ #

Natural Gas Consumption (NATURALGAS) の月次データを使います。数値の単位は BCF（10 億立方フィート）で、2013 年 1 月から 2022 年 9 月までの 117 か月分を辞書で埋め込み、月初の等間隔 index を持つ DataFrame に整えます。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
data = {
    "value": {
        "2013-01-01": 2878.8,
        "2013-02-01": 2567.2,
        "2013-03-01": 2521.1,
        "2013-04-01": 1967.5,
        "2013-05-01": 1752.5,
        "2013-06-01": 1742.9,
        "2013-07-01": 1926.3,
        "2013-08-01": 1927.4,
        "2013-09-01": 1767.0,
        "2013-10-01": 1866.8,
        "2013-11-01": 2316.9,
        "2013-12-01": 2920.8,
        "2014-01-01": 3204.1,
        "2014-02-01": 2741.2,
        "2014-03-01": 2557.9,
        "2014-04-01": 1961.7,
        "2014-05-01": 1810.2,
        "2014-06-01": 1745.4,
        "2014-07-01": 1881.0,
        "2014-08-01": 1933.1,
        "2014-09-01": 1809.3,
        "2014-10-01": 1912.8,
        "2014-11-01": 2357.5,
        "2014-12-01": 2679.2,
        "2015-01-01": 3115.0,
        "2015-02-01": 2925.2,
        "2015-03-01": 2591.3,
        "2015-04-01": 2007.9,
        "2015-05-01": 1858.2,
        "2015-06-01": 1899.9,
        "2015-07-01": 2067.7,
        "2015-08-01": 2052.7,
        "2015-09-01": 1901.3,
        "2015-10-01": 1987.3,
        "2015-11-01": 2249.1,
        "2015-12-01": 2588.2,
        "2016-01-01": 3091.7,
        "2016-02-01": 2652.3,
        "2016-03-01": 2356.3,
        "2016-04-01": 2083.9,
        "2016-05-01": 1965.8,
        "2016-06-01": 2000.7,
        "2016-07-01": 2186.6,
        "2016-08-01": 2208.4,
        "2016-09-01": 1947.8,
        "2016-10-01": 1925.2,
        "2016-11-01": 2159.5,
        "2016-12-01": 2866.3,
        "2017-01-01": 2913.8,
        "2017-02-01": 2340.2,
        "2017-03-01": 2523.3,
        "2017-04-01": 1932.0,
        "2017-05-01": 1892.0,
        "2017-06-01": 1910.4,
        "2017-07-01": 2141.6,
        "2017-08-01": 2093.8,
        "2017-09-01": 1920.5,
        "2017-10-01": 2031.5,
        "2017-11-01": 2357.3,
        "2017-12-01": 3086.0,
        "2018-01-01": 3340.9,
        "2018-02-01": 2710.7,
        "2018-03-01": 2796.7,
        "2018-04-01": 2350.5,
        "2018-05-01": 2055.0,
        "2018-06-01": 2063.1,
        "2018-07-01": 2350.7,
        "2018-08-01": 2313.8,
        "2018-09-01": 2156.1,
        "2018-10-01": 2285.9,
        "2018-11-01": 2715.9,
        "2018-12-01": 2999.5,
        "2019-01-01": 3424.3,
        "2019-02-01": 3019.1,
        "2019-03-01": 2927.8,
        "2019-04-01": 2212.4,
        "2019-05-01": 2134.0,
        "2019-06-01": 2119.3,
        "2019-07-01": 2393.9,
        "2019-08-01": 2433.9,
        "2019-09-01": 2206.3,
        "2019-10-01": 2306.5,
        "2019-11-01": 2783.8,
        "2019-12-01": 3170.7,
        "2020-01-01": 3320.6,
        "2020-02-01": 3058.5,
        "2020-03-01": 2722.0,
        "2020-04-01": 2256.9,
        "2020-05-01": 2072.2,
        "2020-06-01": 2127.9,
        "2020-07-01": 2464.1,
        "2020-08-01": 2399.5,
        "2020-09-01": 2151.2,
        "2020-10-01": 2315.9,
        "2020-11-01": 2442.0,
        "2020-12-01": 3182.8,
        "2021-01-01": 3343.9,
        "2021-02-01": 3099.2,
        "2021-03-01": 2649.4,
        "2021-04-01": 2265.1,
        "2021-05-01": 2117.4,
        "2021-06-01": 2238.4,
        "2021-07-01": 2412.2,
        "2021-08-01": 2433.8,
        "2021-09-01": 2142.3,
        "2021-10-01": 2262.6,
        "2021-11-01": 2693.3,
        "2021-12-01": 3007.3,
        "2022-01-01": 3612.1,
        "2022-02-01": 3064.2,
        "2022-03-01": 2785.4,
        "2022-04-01": 2379.3,
        "2022-05-01": 2247.8,
        "2022-06-01": 2326.9,
        "2022-07-01": 2597.9,
        "2022-08-01": 2566.1,
        "2022-09-01": 2263.3,
    }
}

data = pd.DataFrame(data)
data.rename(columns={"value": "Natural Gas Consumption(BCF)"}, inplace=True)
data.index = pd.to_datetime(data.index)
data = data.asfreq("MS")
data.head()

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

1-2, 123 行目辞書を DataFrame にして CSV なしで再現できるようにする
外側のキー value が列名、内側の日付文字列が index になります。117 行を埋め込んでいるので、FRED からダウンロードしなくてもこのページのコードがそのまま動きます。
124 行目列名に単位 BCF を入れる
sns.lineplot の凡例にこの列名が出るので、図を見た人が 3000 という値の単位を迷いません。モデルの計算には影響しない見た目のための行です。
125-126 行目asfreq("MS") が seasonal_periods=12 の前提を作る
pd.to_datetime で日付型にし、asfreq("MS") で月初の頻度を明示します。Holt-Winters は「12 行前が同じ月」という前提で季節値を更新するので、月が 1 つでも欠けていると季節がずれます。asfreq は欠けた月に NaN を入れるため、欠損があればここで気づけます。頻度を付けないと forecast(33) の index も整数になります。
127 行目先頭 5 行で取り込みを確認
2013 年 1 月の 2878.8 から 5 月の 1752.5 へ下がる 5 行が出れば、日付 index と値が正しく読めています。

	Natural Gas Consumption(BCF)
2013-01-01	2878.8
2013-02-01	2567.2
2013-03-01	2521.1
2013-04-01	1967.5
2013-05-01	1752.5

全期間を描くと、1 月の山（2900 → 3600 BCF）と 7〜8 月の谷（1900 → 2600 BCF）が 12 か月周期で繰り返し、山も谷も 10 年でゆっくり上がっています。山と谷の差は 1000 前後でほぼ一定なので、加法季節が合いそうだと予想できます。

1
2
3
4
plt.figure(figsize=(12, 6))
sns.lineplot(data=data)
plt.grid()
plt.show()

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

1 行目横長にして 10 回の季節周期を数えられるように
figsize=(12, 6) で 117 か月を広げると、1 月の山が 10 回並ぶのが数えられます。山の数 = 周期の回数で、seasonal_periods=12 が妥当かをここで目視確認できます。
2 行目DataFrame を渡すだけで凡例付きの折れ線
sns.lineplot(data=data) は index を横軸、列を線にします。1 列なので 1 本の線になり、凡例に Natural Gas Consumption(BCF) が出ます。
3-4 行目グリッドで山谷の差が一定か確かめる
plt.grid() の横線を基準に、各年の山と谷の差がおよそ 1000 BCF で変わらないことを読みます。差が年々広がるなら seasonal="mul" の乗法季節、一定なら seasonal="add" の加法季節を選ぶ判断材料です。

サンプルデータの図

データの分割 #

データの末尾（2020 年 1 月以降の 33 か月）を予測性能の検証に使い、2013〜2019 年の 84 か月（7 周期分）で学習します。

1
2
3
4
5
6
7
8
data_train = data[data.index < "2020-1-1"]
data_test = data[data.index >= "2020-1-1"]

plt.figure(figsize=(12, 6))
plt.plot(data_train.index, data_train.values, label="train", linewidth=2)
plt.plot(data_test.index, data_test.values, "-.", label="test")
plt.legend()
plt.grid()

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

1-2 行目2020 年 1 月で時間順に切る
data.index < "2020-1-1" で 84 か月を学習、33 か月を検証にします。学習 84 か月は季節周期 12 の 7 周期分で、Holt-Winters が 12 個の季節値を安定して推定するには 2 周期以上が必要なので十分です。ランダム分割は未来を学習に混ぜるので使いません。
5-6 行目学習は太い実線、検証は一点鎖線
linewidth=2 と "-." で線種を分け、次の図で 3 本の予測を重ねたときに「モデルが見た範囲」と「見ていない範囲」を区別できるようにします。
7-8 行目凡例とグリッドは予測図と同じ設定
train・test のラベルと plt.grid() は次の予測図でも同じにしてあり、2 枚を見比べたときに目安がずれません。

データの末尾（2020年以降のデータ）を予測性能の検証に使用するためにデータを分割しますの図

モデルの学習と予測 #

statsmodels.tsa.holtwinters.ExponentialSmoothing のドキュメントを見ながら、trend と seasonal と seasonal_periods を指定します。SES と Holt は前のページと同じ係数で固定し、Holt-Winters は \(\alpha, \beta, \gamma\) を最尤推定に任せます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
ses = SimpleExpSmoothing(data_train)
ses = ses.fit(smoothing_level=0.1)
ses_pred = ses.forecast(33)

holt = Holt(data_train)
holt = holt.fit(smoothing_level=0.1, smoothing_trend=0.2)
holt_pred = holt.forecast(33)

hw = ExponentialSmoothing(
    data_train, trend="additive", seasonal="add", seasonal_periods=12
)
hw = hw.fit()
hw_pred = hw.forecast(33)

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

1-3 行目SES は水準だけ、α=0.1 で固定
前の 2 ページと同じ設定で、予測は水平線になります。3 本重ねたときの「トレンドも季節もなし」の基準線です。
5-7 行目Holt は水準 + 傾き、α=0.1・β=0.2 で固定
前のページと同じ設定で、予測は傾いた直線になります。「トレンドあり・季節なし」の基準線として、Holt-Winters との差が季節成分だけから生じるようにしています。
9-11 行目trend と seasonal と seasonal_periods の 3 つで Holt-Winters になる
trend="additive"（"add" と同じ）で傾きを足し、seasonal="add" で山谷の高さが一定の加法季節を選び、seasonal_periods=12 で月次の 1 年周期を指定します。山谷が水準に比例して大きくなるデータなら seasonal="mul" に変えます。seasonal_periods を間違えて 4 や 7 にすると、季節値が実際の月とずれて予測が波打つ位置も外れます。
12 行目fit() に引数を渡さず α・β・γ を最尤推定
SES と Holt では係数を手で固定しましたが、ここでは fit() に何も渡さず、statsmodels が二乗誤差を最小にする \(\alpha, \beta, \gamma\) と初期の水準・傾き・12 個の季節値を推定します。手で決めたいときは fit(smoothing_level=…, smoothing_trend=…, smoothing_seasonal=…) と 3 つ渡します。
13 行目forecast(33) は 12 個の季節値を繰り返しながら伸ばす
\(h\) か月目の予測は \(\ell_T + h, b_T + s_{T+h-12}\) で、33 か月分には季節値が 2 周と 9 か月分繰り返されます。直線の上に山谷が乗った形になり、冬の山と夏の谷が検証期間に現れます。

1
2
3
4
5
6
7
8
plt.figure(figsize=(12, 6))
plt.plot(data_train.index, data_train.values, label="train", linewidth=2)
plt.plot(data_test.index, data_test.values, "-.", label="test")
plt.plot(ses_pred.index, ses_pred.values, "-.", label="prediction(ses)")
plt.plot(holt_pred.index, holt_pred.values, "-.", label="prediction(Holt)")
plt.plot(hw_pred.index, hw_pred.values, "-o", label="prediction(Holt-Winters)")
plt.legend()
plt.grid()

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

2-3 行目学習・検証は前の図と同じ線種
太い実線の train と一点鎖線の test を先に描き、その上に 3 本の予測を重ねます。
4-5 行目SES と Holt は一点鎖線で控えめに
2 つの基準線を "-." で描き、主役の Holt-Winters と見分けやすくしています。SES は水平、Holt は斜めの直線で、どちらも検証データの山谷を横切るだけです。
6 行目Holt-Winters は丸マーカーで強調
"-o" で 33 個の予測点に丸を付けます。丸の列が検証データの一点鎖線の山谷に沿って上下していれば、季節成分が 12 か月周期を正しく捉えています。
7-8 行目凡例で 3 本を見分ける
prediction(ses)・prediction(Holt)・prediction(Holt-Winters) のラベルで線を区別し、plt.grid() の横線で冬の山が 3300 を超えるか、夏の谷が 2100 付近かを読みます。

statsmodels.tsa.holtwinters.ExponentialSmoothingを確認しつつ、パラメータ…の図

実行結果の読み方 #

SES の予測は 2400 BCF 前後の水平線、Holt の予測はそこから緩やかに上がる直線で、どちらも検証データの冬の山（3300〜3600）と夏の谷（2100〜2300）を外します。Holt-Winters の丸マーカーの列だけが、2020 年 1 月・2021 年 1 月・2022 年 1 月に山、各年の夏に谷を作り、検証データの一点鎖線に沿って動きます。山の高さが実測より少し低めに出るなら、加法季節では 10 年で広がった山谷の振幅を表しきれていない兆候で、seasonal="mul" に変えて比べる価値があります。hw.params を表示すると推定された smoothing_level・smoothing_trend・smoothing_seasonal が見られ、hw.fittedvalues で学習期間内の当てはめ線を描くと季節値が 7 周期でどう安定していったかが分かります。

使いどころと落とし穴 #

Holt-Winters 法が向くのは、天然ガス消費量・電力需要・小売売上のように、はっきりした周期（月次なら 12、日次なら 7）とゆるやかなトレンドを持つ系列です。最初の落とし穴は加法と乗法の選び間違いで、山谷の差が水準に比例して広がるデータに seasonal="add" を使うと、遠い未来の山ほど低く予測されます。全期間の図で山と谷の差（約 1000 BCF）が一定か広がるかを先に見て決めます。二つ目は seasonal_periods の誤りで、月次に 4 や日次に 30 を渡しても statsmodels はエラーを出さずに学習するため、予測の山谷が実測とずれる形でしか気づけません。三つ目はデータの短さで、季節値 12 個を推定するには最低 2 周期（24 か月）が必要で、1 年分では fit が失敗するか季節値が観測そのものを丸写しした形になります。季節性が 2 つ重なる日次データ（週次 7 と年次 365）は 1 つの seasonal_periods では表せないので、複数季節性のベンチマークで扱う SARIMAX + フーリエ項などに切り替えます。

Holtの線形指数平滑法 — SESにトレンド成分を加えた拡張
ホルト・ウィンター法 — 状態空間モデルとして Holt-Winters を定式化し予測区間まで出す
ETS モデル — 誤差・トレンド・季節性を分離して予測