STL分解

Timeseries

チェック3: STL で季節性とトレンドを分解する

作成日: 最終更新: 読了時間: 2 分
まとめ
  • STL（Seasonal-Trend decomposition using Loess）を使うと、非線形トレンドや複数の季節性を柔軟に分解できます。
  • サンプルデータを作成し、statsmodels の seasonal_decompose / STL を比較します。
  • トレンド曲線を区分近似（piecewise）することで、トレンド除去後の系列をより滑らかに扱えるようにします。

1. ライブラリの読み込み #

import japanize_matplotlib
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

2. 季節性＋トレンドを持つダミー系列 #

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) を指定してクラシック分解を実行します。ここでは period=14 とします。

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 を用いたロバスト平滑化で、複雑な季節性にも対応できます。

stl = STL(
    df["観測値"],
    period=14,
    seasonal=13,
    robust=True,
)
result = stl.fit()

plt.figure(figsize=(12, 8))
result.plot()
plt.tight_layout()
plt.show()

seasonal / trend / resid を更新した DataFrame も result から取得できます。

5. トレンド成分を区分線形で近似する #

トレンドをそのまま使うと長く尾を引くので、区分線形に落として扱いやすくします。以下は公開されている segments_fit（ruoyu0088 氏）をベースにした実装です。

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)

区分線形の結果を描画します。

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. トレンドを前後に延長しておく #

STL の trend は窓幅の関係で両端が欠けるので、最初と最後の値を複製して補完しておきます。

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. トレンドを引いた残差を確認 #

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()

これで季節性＋ノイズのみが残り、定常化の目処が立ちます。

8. 実務での留意点 #

  • period の決め方
    ACF（自己相関）や FFT のピークから候補を拾い、STL で確認するのが安全。

  • ロバスト設定
    robust=True にすると外れ値の影響を Loess が抑えてくれます。

  • 複数季節性
    本家 STL は単一周期ですが、MSTL（statsmodels.tsa.seasonal.MSTL）を使えば複数周期を同時に扱えます。

  • 再利用
    分解した季節性をダミー特徴として使う／トレンドを差分で除去するなど、後段の前処理で活用できます。

9. まとめ #

ステップ目的備考
STL で分解トレンド・季節性・残差を分けるperiod 選定が肝
トレンドを区分線形に近似なめらかな補間外れ値に強い
トレンドを補完し残差を作成モデル入力を定常化AR 系列やツリーモデルなどに流用

STL 分解は “観察 → 仮説 → 前処理” の橋渡しになる重要ステップです。Notebook でコピペできるようテンプレートを整備しておきましょう。