STLの仕組みを理解してパラメータを適切に設定する【Python】【statsmodels】

1. 概要

1.1 本記事の目的

本記事では、時系列データを「トレンド」「季節性」「残差」の3要素に分解するSTL^[1]（Seasonal-Trend decomposition using Loess）について解説します。
具体的には、以下の2点を理解することで、STL実行時に適切なパラメータを設定できるようになることを目指します。

• STLで使用される平滑化手法Loessのアルゴリズム
• STLのアルゴリズム

1.2 STLとは

STLは、時系列データを以下の3要素に分解する手法です。

• トレンド（Trend）
  時系列の長期的な傾向。
• 季節性（Seasonal）
  一定期間で繰り返す周期的な変動（例：日、週、月、年）。
  季節性の周期の長さは、事前にパラメータで設定する。
• 残差（Remainder）
  トレンドや季節性で説明できない、不規則な変動成分。

時系列データをこれらの要素に分解することで、データの特徴をより明確に理解し、より正確な予測を立てるのに役立ちます。

STLでは、滑らかで解釈しやすいトレンドや季節性を抽出するためにLoessという手法が用いられています。そのため、STLのパラメータにはLoessに関するパラメータも含まれており、Loessについての理解が必要不可欠です。
本記事では、2章でLoess、3章でSTLのアルゴリズム、4章でstatsmodelsでのSTLの実行例とパラメータの決め方について解説します。5章では、複数の季節性を扱えるMSTLを紹介します。

2. Loess

Loess（LOcally Estimated Scatterplot Smoothing）は局所回帰の一手法で、STLにおいて季節性やトレンド等の平滑化に使用されます。以下でLoessのアルゴリズムを説明します。

2.1 基本的な平滑化アルゴリズム

注目する点を1点ずつずらしながら、以下の手順を行います。

1. 注目する点の近傍点に対して、「近傍重み」（2.3参照）を計算する。
   • 使用する近傍点の数は、パラメータ「ウィンドウサイズ」で事前に指定する。
2. 重み付き最小二乗法による回帰を行う。
   • 回帰に用いる多項式の次数は、パラメータで事前に指定する（一般的には1次または2次を使用。STLでは0次または1次）。
3. 注目している点の回帰結果を平滑化値とする。
4. 手順1〜3を全ての点に対して実行する。

2.2 外れ値への対応

基本的な平滑化では、外れ値が回帰線に大きな影響を与える可能性があります。そこで、以下の追加手順を行います。

5. 全ての点に対して、平滑化値と元の値との差（平滑化の残差）を求める。
6. 残差から「ロバスト重み」（2.4参照）を計算する。
7. 近傍重みとロバスト重みの両方を使用して、手順2〜4を再実行する。
8. 手順5〜7を指定回数繰り返す。

ここで2種類の重みは以下の特徴を持っています。

• 近傍重み
  注目する点に近い点ほど大きい重みが与えられる。
• ロバスト重み
  外れ値に小さい重みが与えられる。

以下で、各重みの計算方法を説明します。（数式が苦手な方は読み飛ばして大丈夫です🙆）

2.3 近傍重みの計算方法

$i=1 \sim n$とし、説明変数を$x_{i}$とします。注目する点$x$の近傍点として$q$個の点を選び、そのうち最も遠い点までの距離を$\lambda_{q}(x)$とします。各点の近傍重みは、以下の関数$W(d)$を用いて計算します。

$$\begin{equation} W(d) = \begin{cases} (1-d^{3})^{3} & 0 \leq d < 1 \\ 0 & d \geq 1 \end{cases} \end{equation}$$

ここで、$d$は距離を表しており、距離が短いほど（$d$が0に近いほど）重みが大きくなることがわかります。点$x_i$の重み$\nu_{i}(x)$は以下の式で表されます。

$$\begin{equation} \nu_{i}(x) = W \left(\frac{|x_{i}-x|}{\lambda_{q}(x)}\right) \end{equation}$$

式(2)より、点$x_i$が$x$に近いほど重みが大きくなり、最も遠い近傍点（注目する点から$q$番目に近い点）では重みが0になることがわかります。また、ウィンドウサイズ$q$を大きくすると$\lambda_{q}(x)$が大きくなるため、重み$\nu_{i}(x)$は1に近づき、通常の多項式回帰に近くなります。
つまり、近傍重みは以下の特徴を持つことがわかります。

• 注目する点に近い点ほど大きな重みを与える。
• ウィンドウサイズ$q$が大きいほど、滑らかな回帰線となる。

2.4 ロバスト重みの計算方法

点$x_{i}$における平滑化の残差（平滑化値と元の値との差）を$r_{i}$とします。各点のロバスト重みは、以下の関数$B(r)$を用いて計算します。

$$\begin{equation} B(r) = \begin{cases} (1-r^{2})^{2} & 0 \leq r < 1 \\ 0 & r \geq 1 \end{cases} \end{equation}$$

ここで、$r$は残差の大きさを表しており、残差が大きいほど（$r$が大きいほど）重みが小さくなることがわかります。
点$x_i$の重み$\rho_{i}$は、以下の式で表されます。

$$\begin{equation} \rho_{i} = B \left(\frac{|r_{i}|}{6\,\mathrm{median}(|r|)}\right) \end{equation}$$

式(4)より、残差$r_i$が大きいほど重みは小さくなり、残差の中央値の6倍以上大きい場合は重みが0になることがわかります。
つまり、ロバスト重みは以下の特徴を持つことがわかります。

• 外れ値には小さい重みを与える。

3. STLのアルゴリズム

3.1 outer loopとinner loop

STLでは、outer loopとinner loopという2種類のループ処理を組み合わせることで、時系列データの分解を行います。

• outer loop
  inner loopで計算された残差を用いて、ロバスト重みを更新する。
• inner loop
  原系列に対して、トレンド・季節性・残差への分解を繰り返す。

名前の通り、outer loopの処理の中にinner loopが含まれていて、それぞれのループの反復回数はパラメータで指定します。例えば{outer loop: 5回, inner loop: 2回}とすると、5回のouter loopの中で2回ずつinner loopが実行されます。このようにouter loopとinner loopを複数回繰り返すことで、分解の精度を高めていきます。

inner loopとouter loopのフローチャート

3.2 STLの手順

原系列を$Y$、トレンドを$T$、季節性を$S$、残差を$R$とします。この時、時刻$t=1 \sim n$における各系列の点を、$Y_{t}$、$T_{t}$、$S_{t}$、$R_{t}$とします。また、$k$回目の分解で計算されたトレンドを$T_{t}^{(k)}$、季節性を$S_{t}^{(k)}$とします。
以下では$(k+1)$回目の分解の手順を説明し、$T_{t}^{(k+1)}$、$S_{t}^{(k+1)}$を計算していきます。

3.2.1 inner loop

1. トレンドの除去
原系列から前回（$k$回目）の分解で計算されたトレンドを除去します。初回の分解では、$T_{t}^{(0)}=0$とします。

$$\begin{equation} Y_{t}(\text{detrended}) = Y_{t} - T_{t}^{(k)} \end{equation}$$

2. 周期成分の平滑化
手順1でトレンドを除去したデータから、平滑化周期成分$C_{t}^{(k+1)}$を抽出します。この平滑化の方法が少しわかりにくいため、3.3で詳しく説明します。

$$\begin{equation} C_{t}^{(k+1)} = Y_{t}(\text{detrended})\text{の平滑化周期成分} \end{equation}$$

3. 平滑化周期成分から低周波成分を抽出
手順2で計算した平滑化周期成分$C_{t}^{(k+1)}$に対して、移動平均ローパスフィルタを適用することで、低周波成分（平滑化周期成分に混ざった、トレンドのような緩やかな変動）を抽出します。さらに、抽出した低周波成分に対してLoess平滑化を適用します。この結果得られる低周波成分を$L_{t}^{(k+1)}$とします。

$$\begin{equation} L_{t}^{(k+1)} = C_{t}^{(k+1)}\text{の低周波成分} \end{equation}$$

4. 平滑化周期成分から低周波成分を除去
手順2で計算した平滑化周期成分$C_{t}^{(k+1)}$から、手順3で計算した低周波成分$L_{t}^{(k+1)}$を除去します。これにより、平滑化周期成分から低周波成分が除去され、季節性$S_{t}^{(k+1)}$が抽出されます。

$$\begin{equation} S_{t}^{(k+1)} = C_{t}^{(k+1)} - L_{t}^{(k+1)} \end{equation}$$

5. 季節性の除去
原系列から手順4で計算した季節性$S_{t}^{(k+1)}$を除去します。

$$\begin{equation} Y_{t}(\text{deseasonalized}) = Y_{t} - S_{t}^{(k+1)} \end{equation}$$

6. トレンド平滑化
手順5の結果に対してLoess平滑化を適用し、トレンド$T_{t}^{(k+1)}$を抽出します。

$$\begin{equation} T_{t}^{(k+1)} = Y_{t}(\text{deseasonalized})\text{の平滑化結果} \end{equation}$$

7. 手順1〜6を指定回数繰り返し
パラメータで指定された回数だけ、手順1〜6を繰り返します。

3.2.2 outer loop

1. 残差を計算
inner loopによって原系列が季節性とトレンドに分解されたので、残差を計算します。

$$\begin{equation} R_{t} = Y_{t} - T_{t} - S_{t} \end{equation}$$

2. ロバスト重みを計算
手順1で計算した残差$R_{t}$を用いて、各時刻におけるロバスト重みを計算します。

3. inner loopの実行
更新されたロバスト重みを用いて、inner loopを実行します。

3.3 周期成分のLoess平滑化

inner loopの手順2「周期成分の平滑化」では、以下の手順で平滑化を行います。

1. 事前に設定したパラメータ「季節性の周期の長さ」ごとに、inner loopの手順1でトレンドを除去したデータから点を抽出する。
2. 抽出したデータに対して、Loess平滑化を行う。

例えば、日次データから曜日の季節性を抽出する場合、曜日ごとに（6個飛ばし）点を抽出し、Loess平滑化を行います。下の図は、この例を示したものです。

• 上段：トレンドを除去したデータ。破線の丸は、同じ曜日として抽出する点。
• 下段：抽出した点。注目点の近傍点（濃い色の点）に重み付き回帰を行い、赤い回帰線上の値（×）を注目点の平滑化値とする。

4. パラメータの決め方

Pythonの統計モデルライブラリstatsmodelsでは、statsmodels.tsa.seasonal.STL^[2]クラスを用いてSTLを実行できます。
本章では、実行例とパラメータの一覧を示したうえで、結果への影響が大きいパラメータperiod・seasonal・trend・robustの決め方を説明します。本章の内容は、statsmodels 0.14.6で確認しています。

4.1 実行例

以下は、statsmodelsに付属する「大気中のCO2濃度」のデータ（週次）を月次に集計し、年周期（12か月）の季節性で分解する例です。4.3以降の検証も、このデータで行います。

import matplotlib.pyplot as plt
from statsmodels.datasets import co2
from statsmodels.tsa.seasonal import STL

data = co2.load().data["co2"]
data = data.resample("MS").mean().ffill()  # 週次を月次に集計し、欠損を直前の値で補完

res = STL(data, period=12).fit()
res.plot()
plt.show()

図のTrend・Season・Residは、それぞれトレンド・季節性・残差です。分解した各成分は、res.trend・res.seasonal・res.residで取得できます。
まずはパラメータperiodだけを指定し、残りは既定値で試すのがおすすめです。

4.2 パラメータ早見表

パラメータは、STL(...)に渡すものと、fit(...)に渡すものに分かれます。

STLクラスに渡すパラメータは、次のとおりです。

パラメータ	意味	既定値	決め方の目安
period	季節性の周期の長さ。	pandasのデータならインデックスの頻度から推定。 NumPy配列では指定が必須。	月次データの年周期なら12、日次データの週周期なら7。 明示するのがおすすめ（4.3参照）。
seasonal	季節性を平滑化するウィンドウサイズ（何周期分の値を使うか）。	7	3以上の奇数。 大きいほど、季節性のパターンが周期ごとに変わりにくくなる（4.4参照）。
trend	トレンドを平滑化するウィンドウサイズ。	periodとseasonalから自動計算。	periodより大きい奇数。 大きいほどトレンドが滑らかになる（4.5参照）。
low_pass	低周波成分を平滑化するウィンドウサイズ。	periodより大きい最小の奇数。	通常は既定値のまま。
seasonal_deg・trend_deg・low_pass_deg	Loessで使う多項式の次数。	1	0か1。 通常は1。
robust	ロバスト重みを使うか。	False	外れ値を含むデータではTrue（4.6参照）。
seasonal_jump・trend_jump・low_pass_jump	Loessで平滑化値を計算する間隔。 2以上にすると、間の点は前後の点から線形補間される。	1	系列が長く計算が遅いときに大きくする。

fit^[3]メソッドに渡すパラメータは、次のとおりです。

パラメータ	意味	既定値	決め方の目安
inner_iter	inner loopの反復回数。	robust=Trueなら2、Falseなら5。	通常は既定値のまま。
outer_iter	outer loopの反復回数。	robust=Trueなら15、Falseなら0。	通常は既定値のまま。

4.3 period

決め方

• 月次データの年周期なら12、日次データの週周期なら7、時間単位のデータの日周期なら24を指定する。
• 省略すると、pandasのデータではインデックスの頻度から推定される。ただし、週次データでは52になるなど、意図した周期と一致するとは限らないため、明示する。
• 複数の季節性（例：日周期と週周期）を持つデータは、STLでは扱えない。5章で紹介するMSTLを使う。

値を変えたときの影響

• 実際の周期と異なる値にする：周期的な変動を季節性として抽出できず、残差に残る。

仕組み

• 3.3で、季節性を求めるために点を抽出する間隔になる。
  • 例：月次データでperiod=12なら、12か月ごとに同じ月の値を抽出する。

CO2濃度のデータで確認した結果

CO2濃度のデータで、periodを実際の周期（12）からずらして比べました。振幅は、年ごとの季節性の振幅（その年の季節性の最大値と最小値の差）の平均です。

• period=12：残差の標準偏差は0.21、振幅は約6.1になる。
• period=11・period=13：残差の標準偏差が約0.5に増え、振幅は約4.8に小さくなる。
• period=6：残差の標準偏差が約1.1に増え、振幅は約1.6まで小さくなる。

4.4 seasonal

決め方

• 3以上の奇数で指定する。statsmodelsのドキュメントでは、通常は7以上を推奨している。
• 季節性のパターンが年々変わると考えられる場合は小さく、ほぼ一定と考えられる場合は大きくする。

値を変えたときの影響

• 小さくする：季節性のパターンが周期ごとに変化しやすくなる。
• 大きくする：季節性のパターンが周期ごとに変化しにくくなる。季節性に含まれなくなった変動は残差に移る。

仕組み

• 3.3の周期成分の平滑化で使うウィンドウサイズ。「何周期分の値を使って平滑化するか」を表す。
  • 例：月次データでseasonal=7なら、前後3年を含む7年分の同じ月の値から季節性を求める。

CO2濃度のデータで確認した結果

CO2濃度のデータで、trend=23に固定してseasonal=7とseasonal=35を比べました。trendの既定値はseasonalによって変わるため、seasonalだけの影響を見るように固定しています。図は、年ごとの季節性の振幅（その年の季節性の最大値と最小値の差）です。

• seasonal=7：振幅が約5.1〜6.9の間で年ごとに上下する。
• seasonal=35：振幅が約5.4から約6.5へなだらかに増え続ける。
• 残差の標準偏差は0.21から0.25に増える。

4.5 trend

決め方

• periodより大きい奇数で指定する。
• トレンドが細かく波打って見える場合は大きく、急な変化を捉えたい場合は小さくする。

値を変えたときの影響

• 大きくする：トレンドが滑らかになる。トレンドに含まれなくなった変動は残差に移る。
• 小さくする：急な変化を捉えやすくなる。ただし、小さくしすぎると、トレンドが短期的な変動まで拾って波打つ。

仕組み

• 3.2.1の手順6で、トレンドをLoess平滑化する際のウィンドウサイズ。
• パラメータseasonalとは異なり、観測点の数を表す。
• 省略すると、$1.5 \times \text{period} / (1 - 1.5 / \text{seasonal})$以上の最小の奇数が使われる。4.1の例（period=12、seasonal=7）では23になる。

CO2濃度のデータで確認した結果

CO2濃度のデータで、trend=23（既定値）とtrend=101を比べました。図は、見やすさのため1980〜1990年を切り出しています。

• trend=101のトレンドは、trend=23よりも直線に近い。
• trend=101の残差には、数年単位の緩やかな波が現れる。
• データ全体では、trend=101にすると残差の標準偏差が0.21から0.40に増える。

4.6 robust

決め方

• 外れ値を含む可能性があるデータではTrueにする。

値を変えたときの影響

• Trueにする：外れ値の影響が季節性やトレンドに漏れにくくなり、外れ値は残差に現れる。ただし、計算時間が増える。

仕組み

• outer loop（3.1参照）でロバスト重み（2.4参照）を計算し、inner loopの平滑化で使う。
• 既定では、outer loopを15回繰り返す（Falseでは0回）。

CO2濃度のデータで確認した結果

CO2濃度のデータのうち、1966年7月・1979年1月・1991年7月の3点に+15の外れ値を加えて比較しました。図の赤い破線は、外れ値を加えた時点です。

• robust=False：外れ値の影響が季節性に漏れ出し、外れ値を加える前の季節性との差が最大3.95になる。
  • 季節性は同じ月の前後の年の値から求めるため、外れ値を加えた年だけでなく、前後の年の季節性にも乱れが広がる。
• robust=True：差が最大0.87に収まる。外れ値はほぼそのまま（約15）残差に現れるため、残差から外れ値を見つけやすい。

5. STLの拡張

本記事では、時系列データを「トレンド」「季節性」「残差」に分解する手法であるSTLについて解説しました。しかし、STLは単一の季節性しか扱えないため、複数の季節性を持つ時系列データに対しては十分な分解性能を発揮できません。
この課題を解決するために、STLを拡張したMSTLという手法が考案されています。MSTLを用いることで、時系列データを「トレンド」「複数の季節性」「残差」に分解することができます。
MSTLの詳細については、以下の記事で解説しているので、ぜひご参照ください。

脚注

1. Robert B. Cleveland, William S. Cleveland, Jean E. McRae, and Irma Terpenning. 'STL: A Seasonal-Trend Decomposition Procedure Based on Loess', Journal of Official Statistics, Vol. 6, No. 1, pp. 3-73, (1990)

2. statsmodels.tsa.seasonal.STL

3. statsmodels.tsa.seasonal.STL.fit