まとめSTL(Seasonal-Trend decomposition using Loess)は、時系列を「トレンド + 季節成分 + 残差」の 3 本に分ける手法で、LOESS 平滑化を使うため直線でないトレンドや形の変わる季節性にも追従します。このページでは周期 14 の季節性と、100 日目で上がり 200 日目で下がる折れ線トレンドを持つ 365 日の人工データを作り、古典的な seasonal_decompose(period=14) と STL(period=14, seasonal=13, robust=True) の結果を見比べます。さらに STL のトレンドを segments_fit で 4 区間の折れ線に近似し、両端を 7 日ずつ延長してから観測値から引くことで、季節性とノイズだけが残った定常に近い系列を作ります。
トレンド成分を扱いやすくする で np.polyfit によるトレンド推定と「トレンドを引いてから学習する」戦略を先に見ておくと、このページで折れ線トレンドを多項式でなく STL で取り出す理由がはっきりします。
直感
#
2 週間ごとに繰り返す売上の波があり、その上に 100 日目から始まったキャンペーンで売上が伸び、200 日目で終わって落ちていく店を想像してください。多項式 1 本ではこの「途中で折れる」トレンドを追えず、周期 14 の波を平均で捉えることもできません。STL は窓の中で局所的に重みつき回帰(LOESS)を繰り返し、トレンドは「近所の点だけ見てなめらかに」、季節成分は「14 日おきの同じ位置の点だけ見て」推定するので、折れたトレンドも形の変わる季節性も別々に取り出せます。
観測値→
周期 14 ごとの同位置を LOESS で平滑化→
季節成分→
観測値 − 季節成分を LOESS で平滑化→
トレンド→
残り→
残差
1. ライブラリの読み込み
#
1
2
3
4
5
6
| import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.seasonal import STL
|
何が起きているか(ステップで見る)
#
コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。
見るポイント:元データ = トレンド + 季節成分 + 残差。残差に残った大きな山が異常の候補。
2. 季節性+トレンドを持つダミー系列
#
周期 14・28・7 の成分と、100 日目から上がり 200 日目から下がる折れ線トレンド、一様ノイズを足した 365 日のデータを作ります。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
| date_list = pd.date_range("2021-01-01", periods=365, freq="D")
value_list = [
10
+ i % 14
+ np.log(1 + i % 28 + np.random.rand())
+ np.sqrt(1 + i % 7 + np.random.rand()) * 2
+ (((i - 100) / 10)) * (i > 100)
- ((i - 200) / 7) * (i > 200)
+ np.random.rand()
for i, di in enumerate(date_list)
]
df = pd.DataFrame(
{
"日付": date_list,
"観測値": value_list,
}
)
df.head(10)
|
3. 古典的な seasonal_decompose
#
周期 period=14 を指定して古典分解を実行します。トレンドは 14 日の中心化移動平均、季節成分は 14 日ごとの同位置の平均です。
1
2
3
4
5
6
7
8
| dr = seasonal_decompose(
df["観測値"], period=14, extrapolate_trend="freq"
)
plt.figure(figsize=(12, 8))
dr.plot()
plt.tight_layout()
plt.show()
|
trend が長期成分、seasonal が 14 日周期の季節性、resid がそれ以外の残差です。ただし seasonal_decompose は季節成分を全期間で同一の形と仮定するため、複数周期や形の変わる季節性には弱い点に注意してください。
4. STL による分解
#
STL クラスは LOESS(局所重みつき回帰)を使ったロバスト平滑化で、形が変わる季節性や折れたトレンドにも追従します。
1
2
3
4
5
6
7
8
9
10
11
12
| stl = STL(
df["観測値"],
period=14,
seasonal=13,
robust=True,
)
result = stl.fit()
plt.figure(figsize=(12, 8))
result.plot()
plt.tight_layout()
plt.show()
|
result.seasonal・result.trend・result.resid を列として足した DataFrame を作っておくと、後段の特徴量作成や異常検知にそのまま使えます。
5. トレンド成分を区分線形で近似する
#
STL のトレンドは 365 点の曲線で、そのまま使うと細かい揺れまで引き継いでしまいます。折れ点の数を 4 区間と決めて折れ線に落とすと、「100 日目から上がり 200 日目から下がる」という構造が数字で読めるようになります。以下は公開されている segments_fit(ruoyu0088 氏)をベースにした実装です。
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
| from scipy import optimize
from scipy.special import huber
def segments_fit(X, Y, count, is_use_huber=True, delta=0.5):
"""
original: https://gist.github.com/ruoyu0088/70effade57483355bbd18b31dc370f2a
"""
xmin = X.min()
xmax = X.max()
seg = np.full(count - 1, (xmax - xmin) / count)
px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax]
py_init = np.array(
[Y[np.abs(X - x) < (xmax - xmin) * 0.01].mean() for x in px_init]
)
def func(p):
seg = p[: count - 1]
py = p[count - 1 :]
px = np.r_[np.r_[xmin, seg].cumsum(), xmax]
return px, py
def err(p):
px, py = func(p)
Y2 = np.interp(X, px, py)
if is_use_huber:
return np.mean(huber(delta, Y - Y2))
else:
return np.mean((Y - Y2) ** 2)
r = optimize.minimize(err, x0=np.r_[seg, py_init], method="Nelder-Mead")
return func(r.x)
|
区分線形の結果を STL トレンドに重ねて描きます。
1
2
3
4
5
6
7
8
9
10
| trend_data = result.trend.dropna()
x = np.arange(0, len(trend_data))
px, py = segments_fit(x, trend_data, 4)
plt.figure(figsize=(12, 4))
plt.plot(x, trend_data, "-x", label="STLトレンド")
plt.plot(px, py, "-o", label="区分線形近似")
plt.legend()
plt.ylim(trend_data.min() - 2, trend_data.max() + 2)
plt.show()
|
6. トレンドを前後に延長しておく
#
ここで使う trend_data は観測値と長さが同じですが、seasonal_decompose のように両端が窓幅の半分(7 日)欠ける分解結果にも対応できるよう、最初と最後の値を 7 個ずつ複製して前後に延ばしておきます。
1
2
3
4
5
6
7
8
| period = 14
trend_data_first, trend_data_last = trend_data.iloc[0], trend_data.iloc[-1]
for _ in range(int(period / 2)):
trend_data[trend_data.index.min() - 1] = trend_data_first
trend_data[trend_data.index.max() + 1] = trend_data_last
trend_data = trend_data.sort_index()
|
7. トレンドを引いた残差を確認
#
観測値から延長したトレンドを引き、季節性とノイズだけの系列を作ります。
1
2
3
4
5
6
7
| plt.figure(figsize=(12, 4))
plt.title("観測値とトレンド除去後の値")
plt.plot(result.observed, label="観測値")
plt.plot(result.observed - trend_data, label="トレンド差し引き")
plt.grid()
plt.legend()
plt.show()
|
これで季節性とノイズのみが残り、ADF 検定にかけて定常化の目処を立てられます。
8. 実務での留意点
#
period は ACF(自己相関)や FFT のピークから候補を拾い、STL で分解して残差に周期的な山が残らないかで確認するのが安全です。robust=True にすると外れ値の影響を LOESS が抑えてくれるので、異常検知に使うなら必ず付けます。本家の STL は単一周期ですが、日次と週次が重なるようなデータでは MSTL(statsmodels.tsa.seasonal.MSTL)で複数周期を同時に扱えます。分解した result.seasonal を曜日ダミーの代わりの特徴量にしたり、result.trend を差分で除去したりと、後段の前処理での再利用先も多いので、DataFrame の列として保存しておくと便利です。
9. まとめ
#
| ステップ | 目的 | 備考 |
|---|
| STL で分解 | トレンド・季節性・残差を分ける | period 選定が肝 |
| トレンドを区分線形に近似 | なめらかな補間 | 外れ値に強い |
| トレンドを補完し残差を作成 | モデル入力を定常化 | AR 系列やツリーモデルなどに流用 |
STL(period=14, seasonal=13, robust=True) は、LOESS で局所的に平滑化するため、100 日目・200 日目で折れるトレンドも、周期ごとに形の変わる季節性も、seasonal_decompose より鋭く分けられます。取り出したトレンドは segments_fit で 4 区間の折れ線に近似すれば折れ目の位置が数字で読め、観測値から引けば季節性とノイズだけの系列が残ります。period の選定とロバスト設定の 2 点を押さえた上で、この一連の手順を Notebook のテンプレートにしておくと、新しい系列を受け取ったときに 5 分で分解まで進めます。