このコードは、外れ値を含む時系列データに対してSTL分解を適用し、トレンド・季節性・残差の3成分に分解しています。

import numpy as np, pandas as pd
from statsmodels.tsa.seasonal import STL

np.random.seed(42)
t = np.arange(144.)
data = pd.Series(
    50 + 0.3*t + 10*np.sin(2*np.pi*t/12)
    + np.random.randn(144)*2,
    index=pd.date_range('2010',periods=144,freq='MS'))
data.iloc[30] += 40
data.iloc[80] -= 35

r = STL(data, period=12, robust=True).fit()
print(f"【残差】平均:{r.resid.mean():.2f}"
      f", 標準偏差:{r.resid.std():.2f}")
print(f"【外れ値付近の残差】"
      f"\n  t=30の残差:{r.resid.iloc[30]:+.2f}"
      f"\n  t=80の残差:{r.resid.iloc[80]:+.2f}")

簡単に説明します。

ライブラリのインポート

• numpy：数値計算を行うためのライブラリ。ここでは乱数生成や数値操作に使用しています。
• pandas：時系列インデックス付きのSeriesを作成するために使用しています。
• statsmodels.tsa.seasonal.STL：STL分解を実装したクラスです。

データの生成

• 50 + 0.3*t：ベースレベル50に線形トレンド（傾き0.3）を加えたものです。
• 10*np.sin(2*np.pi*t/12)：振幅10、周期12の季節成分です。
• np.random.randn(144)*2：標準偏差2のノイズです。
• data.iloc[30] += 40、data.iloc[80] -= 35：意図的に2つの外れ値を挿入しています。これにより、外れ値に対する頑健性を検証できます。

STL分解の実行

• STL(data, period=12, robust=True).fit()：STL分解を実行します。
  • period=12：季節の周期を12に設定します。
  • robust=True：頑健な推定を行います。外れ値に低い重みを与えることで、トレンドと季節成分が外れ値に引っ張られるのを防ぎます。
• 戻り値rからr.trend、r.seasonal、r.residでそれぞれの成分を取得できます。

(A) 分解が高速で、大規模データに適している → ×

STL分解はLOESS平滑化を反復的に適用するため、seasonal_decomposeの単純な移動平均よりもむしろ計算コストは高くなります。ただし、実用上問題になるほど遅いわけではありません。

(B) 乗法モデルのみをサポートし、季節成分の精度が高い → ×

STL分解は加法モデルのみを直接サポートしています。乗法モデルが必要な場合は、データの対数を取ってからSTLを適用し、結果を逆変換するのが一般的です。

(C) LOESS平滑化により柔軟な分解が可能で、外れ値に頑健である → ○

正解です。STLの名称自体が「Seasonal and Trend decomposition using LOESS」の略であり、LOESS平滑化が中核技術です。さらにrobust=Trueにより外れ値への頑健性が得られます。

(D) 自動的に最適な周期を検出できる → ×

STLは周期をperiodパラメータで明示的に指定する必要があります。周期の自動検出にはスペクトル解析や自己相関関数などの手法を別途使用します。

seasonal_decompose（classical decomposition）とSTLの主な違いは以下の通りです。

トレンド推定の手法

seasonal_decomposeは単純移動平均でトレンドを推定します。窓幅は周期に固定されるため、トレンドの滑らかさを調整できません。STLはLOESS（局所重み付き回帰）を使用し、パラメータにより平滑化の度合いを柔軟に制御できます。

季節成分の扱い

seasonal_decomposeは全期間で同一の季節パターンを仮定します。つまり、1月のパターンは10年前も今年も同じです。STLは季節成分が時間とともに緩やかに変化することを許容しており、より現実的なモデル化が可能です。

外れ値への頑健性

seasonal_decomposeには外れ値に対する特別な処理がなく、外れ値がトレンドや季節成分を歪めます。STLはrobust=Trueを指定することで、外れ値の影響を抑えた頑健な推定が可能です。

端点の処理

seasonal_decomposeは移動平均の性質上、データの両端にNaN（欠損値）が生じます。STLはLOESS平滑化により端点でも値を推定できるため、欠損が生じません。

robust=Trueの効果を具体的に確認するために、同じデータに対してrobust=Falseと比較してみましょう。

以下、コードです。

for robust in [False, True]:
    r = STL(data, period=12, robust=robust).fit()
    print(f"[robust={robust!s:<5}]"
          f" 残差の標準偏差:{r.resid.std():.2f},"
          f" t=30の残差:{r.resid.iloc[30]:+6.2f},"
          f" t=80の残差:{r.resid.iloc[80]:+6.2f}")

以下、実行結果です。

[robust=False] 残差の標準偏差:3.74, t=30の残差:+27.25, t=80の残差:-22.71
[robust=True ] 残差の標準偏差:4.62, t=30の残差:+38.37, t=80の残差:-35.15

一見するとrobust=Trueの方が残差の標準偏差が大きく、性能が悪いように見えるかもしれません。しかし、注目すべきは外れ値付近の残差の大きさです。

• robust=False：t=30の残差が+27.25（本来+40のはず）。外れ値の影響がトレンドや季節成分に漏れ出しており、トレンドが外れ値に引っ張られています。
• robust=True：t=30の残差が+38.37（本来+40に近い）。外れ値がほぼそのまま残差に封じ込められており、トレンドと季節成分は歪んでいません。

 

ではトレンドと季節成分の推定精度が高く保たれ、外れ値は残差に正しく分離されます。この頑健な推定は、LOESS回帰の重み付けにおいて、残差が大きいデータ点の重みを反復的に低減するBisquare重み関数によって実現されています。

STLにはトレンドと季節成分の平滑化を制御するパラメータがあります。

• seasonal（奇数、デフォルト7）：季節成分のLOESS平滑化の窓幅です。大きくすると季節パターンがより安定的になり、小さくすると季節パターンの時間変化に敏感になります。
• trend（奇数、省略可）：トレンドのLOESS平滑化の窓幅です。指定しない場合は自動計算されます。大きくすると滑らかなトレンドになり、小さくするとトレンドがデータの短期的な変動に追従します。
• robust（デフォルトFalse）：Trueにすると外れ値に頑健な推定を行います。計算時間は増えますが、外れ値がある可能性がある場合は指定することを推奨します。