時系列分解(STL分解)の理論と実装 — トレンド・季節・残差を分離する

毎月の電力消費量のグラフを眺めると、3つの要素が重なっていることに気づきます。まず、年々少しずつ増えている長期的な上昇傾向（トレンド）。次に、夏と冬に消費が増え、春と秋に減るという毎年繰り返される周期的なパターン（季節性）。そして、突然の猛暑や大型イベントによる一時的な変動（残差）。この3つの要素を分離して個別に分析できれば、「長期的な傾向は？」「季節パターンは変化しているか？」「異常な月はどこか？」といった問いに答えることができます。

時系列分解（Time Series Decomposition）は、時系列データをこれらの構成要素に分離する手法であり、あらゆる時系列分析の基礎的なステップです。

時系列分解を理解することで、以下のような場面で活用できます。

1. 異常検知 — 残差成分が異常に大きい時点を検出し、原因を調査する
2. 予測の改善 — 各成分を個別にモデル化して予測し、合成することで精度を向上させる
3. ビジネスインサイトの抽出 — トレンドの変化点を特定し、季節パターンの変化を追跡する

本記事の内容

• なぜ分解するのか — 分析上の利点
• 加法分解と乗法分解
• 古典的分解法（移動平均ベース）
• STL分解の仕組み（LOESSによるロバストな分解）
• STLのパラメータの選び方
• Python（statsmodels）での実装
• 異常検知への応用

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 【時系列解析】自己共分散とは？定義とPythonでの実装で詳しく解説
• SARIMAモデルによる季節時系列予測
• 指数平滑法とHolt-Winters法

なぜ分解するのか

分解の目的

時系列 $x_t$ を分解する最大の理由は、複雑な変動を単純な要素に分けて、各要素を個別に理解・分析できるようにすることです。

例えば、ある小売店の月次売上が前年比10%増えたとします。しかし、この増加は「長期的なトレンドの上昇」なのか、「たまたま季節的に売上が高い月だったから」なのか、「特別なキャンペーンによる一時的な効果」なのか、全体の数字だけでは判断できません。分解によって各成分を取り出すことで、正確な判断が可能になります。

分解の3成分

一般に、時系列は以下の3つの成分に分解されます。

1. トレンド成分 $T_t$ — 長期的な上昇・下降の傾向。数年〜数十年のスケールでの変化
2. 季節成分 $S_t$ — 固定された周期で繰り返されるパターン。月次データなら12か月周期が典型
3. 残差成分 $R_t$ — トレンドでも季節でも説明できない不規則な変動

分解したい対象が定まったところで、次にこれらの成分をどのように組み合わせるか — 加法と乗法の2つのモデルを見ていきましょう。

加法分解と乗法分解

加法モデル

加法モデルでは、観測値は各成分の和として表されます。

$$ \begin{equation} x_t = T_t + S_t + R_t \end{equation} $$

加法モデルが適切なのは、季節変動の振幅が時間とともにほぼ一定の場合です。例えば、気温データは冬と夏の差が毎年ほぼ一定（約30度の差）なので、加法モデルが適しています。

乗法モデル

乗法モデルでは、観測値は各成分の積として表されます。

$$ \begin{equation} x_t = T_t \times S_t \times R_t \end{equation} $$

乗法モデルが適切なのは、季節変動の振幅がトレンドに比例して増減する場合です。例えば、売上が10億円のときの季節変動が1億円で、売上が100億円のときの季節変動が10億円であるような場合です。

対数変換による統一

乗法モデルに対数をとると加法モデルに変換できます。

$$ \log(x_t) = \log(T_t) + \log(S_t) + \log(R_t) $$

そのため、実務では「乗法的な季節性がある場合は対数変換してから加法的に分解する」というアプローチもよく使われます。

判別方法

データが加法的か乗法的かを判断するには、以下の方法が有効です。

1. 目視: 時間とともに季節変動の振幅が変化しているかを確認
2. 対数変換後の確認: 対数変換で振幅が安定すれば乗法的
3. 月別の標準偏差: トレンドの水準別にグループ化し、標準偏差がほぼ一定なら加法的、水準に比例すれば乗法的

加法と乗法の判断方法を理解したところで、次に最も基本的な分解手法である古典的分解法を見ていきましょう。

古典的分解法

アルゴリズム

古典的分解法は、移動平均を使ってトレンドを推定し、そこから季節成分と残差を順に取り出す手法です。加法モデルの場合のアルゴリズムを示します。

ステップ1: トレンドの推定

季節周期 $s$ の移動平均でトレンドを推定します。$s$ が偶数（例えば月次データで $s = 12$）の場合、中心化移動平均を使います。

$$ \begin{equation} \hat{T}_t = \frac{1}{s}\left(\frac{1}{2}x_{t-s/2} + x_{t-s/2+1} + \cdots + x_{t+s/2-1} + \frac{1}{2}x_{t+s/2}\right) \end{equation} $$

端の項を $1/2$ にする理由は、$s$ 個の期間にわたる平均を時刻 $t$ に中心化するためです。これにより、移動平均がちょうど季節1周期分をカバーし、季節変動が相殺されて滑らかなトレンドが得られます。

ステップ2: トレンド除去

$$ x_t – \hat{T}_t = \hat{S}_t + \hat{R}_t $$

ステップ3: 季節成分の推定

トレンド除去後のデータを、各季節ごと（例えば、全ての1月のデータ、全ての2月のデータ、…）にグループ化し、各グループの平均を計算します。

$$ \bar{S}_k = \frac{1}{n_k} \sum_{t \in \text{month } k} (x_t – \hat{T}_t), \quad k = 1, 2, \dots, s $$

さらに、季節成分の和が0になるように正規化します。

$$ \hat{S}_k = \bar{S}_k – \frac{1}{s}\sum_{j=1}^{s} \bar{S}_j $$

ステップ4: 残差

$$ \hat{R}_t = x_t – \hat{T}_t – \hat{S}_t $$

Pythonでの古典的分解

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.datasets import get_rdataset
from statsmodels.tsa.seasonal import seasonal_decompose

# AirPassengersデータ
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
passengers = data['passengers'].values
time_index = pd.date_range('1949-01', periods=len(passengers), freq='MS')
ts = pd.Series(passengers, index=time_index)

# 加法モデルと乗法モデルで分解
fig, axes = plt.subplots(4, 2, figsize=(14, 12))

for col, (model_type, title) in enumerate([('additive', 'Additive'),
                                            ('multiplicative', 'Multiplicative')]):
    result = seasonal_decompose(ts, model=model_type, period=12)

    axes[0, col].plot(ts.index, ts, linewidth=1.5)
    axes[0, col].set_title(f'{title} — Observed')

    axes[1, col].plot(result.trend.index, result.trend, linewidth=1.5, color='darkorange')
    axes[1, col].set_title(f'{title} — Trend')

    axes[2, col].plot(result.seasonal.index, result.seasonal, linewidth=1.5, color='green')
    axes[2, col].set_title(f'{title} — Seasonal')

    axes[3, col].plot(result.resid.index, result.resid, linewidth=0.8, color='red')
    axes[3, col].set_title(f'{title} — Residual')

plt.tight_layout()
plt.show()

古典的分解の結果から、加法モデルと乗法モデルの違いが明確にわかります。

1. 加法モデルの残差（左下）: 残差の振幅が時間とともに増大しています。これは、AirPassengersデータの季節変動が乗法的であるにもかかわらず、加法モデルで分解したため、季節成分の一部が残差に漏れていることを示しています。
2. 乗法モデルの残差（右下）: 残差の振幅がほぼ一定で、より「ランダムな」見た目をしています。これは乗法モデルがこのデータに適していることを示唆しています。
3. 季節成分: 加法モデルでは季節成分が時間によらず一定のパターンですが、乗法モデルでは比率として表現されるため、全体の水準が変わっても適切に季節効果を表現できます。

古典的分解法の限界

古典的分解法は概念がわかりやすい反面、いくつかの重要な限界があります。

1. データ端の欠損 — 移動平均により、データの両端 $s/2$ 期分のトレンドが推定できません
2. 季節パターンの固定性 — 季節成分は全期間で同一パターンと仮定されます。実際には季節パターンが緩やかに変化することも多いです
3. 外れ値への脆弱性 — 移動平均は外れ値の影響を受けやすく、トレンドの推定が歪みます

これらの限界を克服するために開発されたのが、次に見るSTL分解です。

STL分解の仕組み

STLとは

STL（Seasonal and Trend decomposition using Loess）は、Cleveland et al.（1990）によって提案された分解手法です。名前の通り、LOESS（LOcally Estimated Scatterplot Smoothing）と呼ばれる局所回帰手法をベースにしています。

STLの主な利点は以下の通りです。

1. 季節パターンの変化に対応 — 季節成分が時間とともに緩やかに変化することを許容します
2. ロバスト性 — 外れ値に対してロバストな推定が可能です（反復重み付け）
3. 柔軟なパラメータ制御 — トレンドと季節成分の平滑度を独立に調整できます

LOESSの概要

STLの理解にはLOESS（局所重み付き回帰）の理解が必要です。LOESSは、各点 $x_0$ の近傍にあるデータ点に重みをつけて低次多項式（通常は1次または2次）を当てはめ、$x_0$ における予測値を計算する手法です。

点 $x_0$ におけるLOESSの推定値は、以下の重み付き最小二乗問題の解から得られます。

$$ \begin{equation} \min_{\beta_0, \beta_1} \sum_{i=1}^{N} w_i(x_0) (y_i – \beta_0 – \beta_1 x_i)^2 \end{equation} $$

重み関数 $w_i(x_0)$ は、$x_0$ に近いデータ点ほど大きな重みを持ちます。典型的にはトリキューブ（tricube）カーネルが使われます。

$$ w_i(x_0) = W\left(\frac{|x_i – x_0|}{d_q(x_0)}\right), \quad W(u) = \begin{cases} (1 – u^3)^3 & 0 \leq u < 1 \\ 0 & u \geq 1 \end{cases} $$

ここで $d_q(x_0)$ は $x_0$ から $q$ 番目に近いデータ点までの距離です。$q$（= 帯域幅パラメータ）を大きくするほど多くのデータ点に重みを与え、推定が滑らかになります。

STLアルゴリズムの流れ

STLは外側ループと内側ループの2重ループ構造を持ちます。

外側ループ（ロバスト性のための反復。$n_o$ 回）: – ロバスト重みの更新（外れ値の影響を抑制）

内側ループ（各反復で $n_i$ 回）: 1. トレンドの除去: $x_t – T_t$ で季節+残差を取り出す 2. サブシリーズの平滑化: 同じ季節（例: 全ての1月）のデータをLOESSで平滑化し、季節成分の候補を得る 3. 低周波フィルタ: 季節成分候補から低周波成分を除去（季節成分にトレンドが混入することを防ぐ） 4. 季節成分の確定: ステップ2の結果からステップ3の低周波成分を引く 5. トレンドの更新: $x_t – S_t$ をLOESSで平滑化してトレンドを更新

STLのパラメータ

STLには主に3つの重要なパラメータがあります。

$n_s$ は最も重要なパラメータです。$n_s$ が小さいと季節パターンが急速に変化でき、$n_s$ が大きいと季節パターンが安定します。

STLの理論を理解したところで、実際にPythonで実装してみましょう。

PythonによるSTL分解の実装

statsmodelsのSTL

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.datasets import get_rdataset
from statsmodels.tsa.seasonal import STL

# AirPassengersデータ
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
passengers = data['passengers'].values
time_index = pd.date_range('1949-01', periods=len(passengers), freq='MS')
ts = pd.Series(passengers, index=time_index)

# 対数変換（乗法的季節性→加法的に変換）
log_ts = np.log(ts)

# STL分解
stl = STL(log_ts, period=12, seasonal=13, trend=25, robust=True)
result = stl.fit()

fig, axes = plt.subplots(4, 1, figsize=(12, 10))

axes[0].plot(log_ts.index, log_ts, linewidth=1.5)
axes[0].set_title('Observed (log scale)')
axes[0].set_ylabel('log(passengers)')

axes[1].plot(result.trend.index, result.trend, linewidth=1.5, color='darkorange')
axes[1].set_title('Trend')
axes[1].set_ylabel('log(passengers)')

axes[2].plot(result.seasonal.index, result.seasonal, linewidth=1.5, color='green')
axes[2].set_title('Seasonal')
axes[2].set_ylabel('Seasonal effect')

axes[3].plot(result.resid.index, result.resid, linewidth=0.8, color='red')
axes[3].axhline(y=0, color='gray', linestyle='--', alpha=0.5)
axes[3].set_title('Residual')
axes[3].set_ylabel('Residual')

plt.tight_layout()
plt.show()

# 季節成分の強さ
seasonal_strength = 1 - np.var(result.resid) / np.var(result.seasonal + result.resid)
trend_strength = 1 - np.var(result.resid) / np.var(result.trend + result.resid)
print(f"Seasonal Strength: {seasonal_strength:.4f}")
print(f"Trend Strength: {trend_strength:.4f}")

STL分解の結果から、AirPassengersデータの構造が明確に読み取れます。

1. トレンド成分: 滑らかに上昇しており、航空旅客数の長期的な成長を反映しています。古典的分解と異なり、データの両端でもトレンドが推定されています。
2. 季節成分: 12か月周期の明確なパターンが見られます。STLの利点として、季節パターンが時間とともにわずかに変化していることも捉えられています。初期の季節振幅と後期の季節振幅を比較すると、微妙な違いが観察できます。
3. 残差成分: 概ねランダムな変動で、明確なパターンは見られません。ロバスト推定により外れ値の影響も抑制されています。
4. Seasonal Strength / Trend Strength: 1に近いほどその成分が支配的であることを示します。季節性が非常に強いデータであることが確認できます。

パラメータの影響

STLのパラメータ $n_s$（seasonal）を変えたときの影響を比較しましょう。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.datasets import get_rdataset
from statsmodels.tsa.seasonal import STL

# データ準備
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
passengers = data['passengers'].values
time_index = pd.date_range('1949-01', periods=len(passengers), freq='MS')
log_ts = pd.Series(np.log(passengers), index=time_index)

seasonal_values = [7, 13, 25, 51]

fig, axes = plt.subplots(len(seasonal_values), 2, figsize=(14, 12))

for idx, ns in enumerate(seasonal_values):
    stl = STL(log_ts, period=12, seasonal=ns, robust=True)
    result = stl.fit()

    # 季節成分
    axes[idx, 0].plot(result.seasonal.index, result.seasonal,
                      linewidth=1, color='green')
    axes[idx, 0].set_title(f'Seasonal (n_s={ns})')
    axes[idx, 0].set_ylim(-0.25, 0.25)

    # 残差
    axes[idx, 1].plot(result.resid.index, result.resid,
                      linewidth=0.8, color='red')
    axes[idx, 1].axhline(y=0, color='gray', linestyle='--', alpha=0.5)
    axes[idx, 1].set_title(f'Residual (n_s={ns})')
    axes[idx, 1].set_ylim(-0.15, 0.15)

plt.tight_layout()
plt.show()

$n_s$ の値による分解結果の違いが明確にわかります。

1. $n_s = 7$（小さい値）: 季節パターンが急速に変化できるため、季節成分がデータの不規則な変動の一部も吸収してしまっています。残差は小さくなりますが、過剰適合の懸念があります。
2. $n_s = 13$: バランスの取れた設定です。季節パターンの緩やかな変化は許容しつつ、過度な適合は避けています。
3. $n_s = 25$（中程度）: 季節パターンがかなり安定しています。数年にわたるゆっくりした変化のみを捉えます。
4. $n_s = 51$（大きい値）: 季節パターンがほぼ固定されています。古典的分解に近い結果になります。

一般的には、$n_s$ は季節周期の1〜2倍程度（月次データなら13〜25）を出発点とし、残差にパターンが残っていないかを確認しながら調整します。

パラメータの影響を理解したところで、STL分解の実践的な応用として異常検知への活用を見ていきましょう。

異常検知への応用

STL分解の残差成分を利用して、時系列データ中の異常値を検出する方法を紹介します。

基本的なアイデア

STL分解により $x_t = T_t + S_t + R_t$ と分解された後、残差 $R_t$ がホワイトノイズに従うならば、$|R_t|$ が異常に大きい時点を異常値と判定できます。

具体的には、残差の中央絶対偏差（MAD）を用いてロバストな閾値を設定します。

$$ \text{MAD} = \text{median}(|R_t – \text{median}(R_t)|) $$

正規分布を仮定すると、$\text{MAD} \approx 0.6745 \sigma$ なので、

$$ \hat{\sigma} = \frac{\text{MAD}}{0.6745} $$

$|R_t| > 3\hat{\sigma}$ の時点を異常値とする「3-MAD法」がロバストな異常検知の基本です。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.datasets import get_rdataset
from statsmodels.tsa.seasonal import STL

# AirPassengersデータに人工的な異常値を追加
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
passengers = data['passengers'].values.copy().astype(float)
time_index = pd.date_range('1949-01', periods=len(passengers), freq='MS')

# 異常値の挿入
np.random.seed(42)
anomaly_indices = [30, 72, 100, 130]
for idx in anomaly_indices:
    passengers[idx] *= 1.5  # 50%増

log_ts = pd.Series(np.log(passengers), index=time_index)

# STL分解
stl = STL(log_ts, period=12, seasonal=13, robust=True)
result = stl.fit()

# MADベースの異常検知
residuals = result.resid.dropna()
mad = np.median(np.abs(residuals - np.median(residuals)))
sigma_hat = mad / 0.6745
threshold = 3 * sigma_hat

anomalies = np.abs(residuals) > threshold

fig, axes = plt.subplots(3, 1, figsize=(12, 9))

# 元データと異常値
axes[0].plot(log_ts.index, log_ts, linewidth=1.5)
for idx in anomaly_indices:
    axes[0].axvline(x=time_index[idx], color='red', alpha=0.3, linestyle='--')
axes[0].set_title('Observed (with injected anomalies)')

# 残差と閾値
axes[1].plot(residuals.index, residuals, linewidth=0.8)
axes[1].axhline(y=threshold, color='red', linestyle='--', label=f'+3σ = {threshold:.3f}')
axes[1].axhline(y=-threshold, color='red', linestyle='--', label=f'-3σ = {-threshold:.3f}')
axes[1].scatter(residuals[anomalies].index, residuals[anomalies],
                color='red', s=50, zorder=5, label='Detected anomalies')
axes[1].set_title('Residuals with Anomaly Detection')
axes[1].legend()

# 検出結果の確認
detected_months = residuals[anomalies].index
true_anomaly_months = [time_index[i] for i in anomaly_indices]
axes[2].plot(log_ts.index, log_ts, linewidth=1, alpha=0.5)
for m in detected_months:
    axes[2].axvline(x=m, color='red', alpha=0.5, linewidth=2)
for m in true_anomaly_months:
    axes[2].axvline(x=m, color='blue', alpha=0.3, linewidth=2, linestyle='--')
axes[2].set_title('Detected (red) vs True (blue dashed) Anomalies')

plt.tight_layout()
plt.show()

print(f"MAD = {mad:.4f}")
print(f"Estimated σ = {sigma_hat:.4f}")
print(f"Threshold (3σ) = {threshold:.4f}")
print(f"\nDetected anomalies: {len(detected_months)} points")
print(f"True anomalies: {len(anomaly_indices)} points")
for m in detected_months:
    true = "✓" if m in true_anomaly_months else "×"
    print(f"  {m.strftime('%Y-%m')} (True anomaly: {true})")

異常検知の結果から、STL分解ベースの手法の有効性が確認できます。

1. ロバスト性: STL分解は robust=True により外れ値に対してロバストな分解を行うため、外れ値がトレンドや季節成分に影響を与えにくく、残差に適切に分離されます。
2. 検出精度: 人工的に挿入した4つの異常値のうち、多くが残差の3-MAD閾値を超えて検出されます。50%の増加は、対数スケールでは約0.4の変化に相当し、通常の残差変動より明らかに大きいです。
3. 偽陽性: 正常なデータ点が誤って異常と検出される（偽陽性）は少なく、閾値の設定が適切であることが示唆されます。

古典的分解法とSTLの比較

ここまでの内容を踏まえて、両手法の特徴を比較しましょう。

実務では、データの特性に応じて使い分けます。探索的なデータ分析では古典的分解で概要を掴み、より精密な分析ではSTLを使うというアプローチが有効です。

まとめ

本記事では、時系列分解の理論と実装について詳しく解説しました。

• 加法分解と乗法分解: 季節変動の振幅がトレンドに依存するかどうかでモデルを選択する
• 古典的分解法: 移動平均でトレンドを推定し、季節成分を平均で求める。シンプルだが季節パターンの変化や外れ値に弱い
• STL分解: LOESSベースのロバストな分解手法。季節パターンの緩やかな変化に対応でき、外れ値にも頑健
• パラメータ $n_s$: 季節成分の平滑度を制御する最重要パラメータ。小さいと柔軟、大きいと安定
• 異常検知への応用: STL残差のMADベースの閾値判定により、季節性を考慮した異常検知が可能

時系列分解は、予測モデルの前処理としても非常に重要です。分解によってデータの構造を理解した上で、各成分に適したモデルを選択することで、予測精度の向上が期待できます。

次のステップとして、以下の記事も参考にしてください。

• 指数平滑法とHolt-Winters法
• Facebook Prophetによる時系列予測
• SARIMAモデルによる季節時系列予測