Zenn
🌊

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

2024/09/11に公開
Python
時系列
時系列分析
statsmodels
tech

1. 概要

本記事では、時系列データをトレンド、季節性、残差の3要素に分解する手法であるSTL分解[1](Seasonal Trend decomposition using Loess)について解説します。STL分解は、複雑な時系列データを理解しやすい要素に分解することで、データの特徴を把握するのに役立ちます。

本記事の目的は以下の2点です。

  • STL分解の仕組みとSTLで使われるLoessを理解する
  • STL分解のパラメータの意味を知り、適切に設定できるようになる(statsmodels)

2. Loess

Loess(Locally Estimated Scatterd Smoothing)は局所回帰の一手法で、STL分解において季節性やトレンド等の平滑化に使用されます。以下にLoessの手順を説明します。

2.1 基本的な平滑化手順

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

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

2.2 外れ値への対応

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

  1. 全ての点に対して、平滑化値と元の値との残差を求める
  2. 残差から「ロバスト重み」(2.4参照)を計算する
  3. 近傍重みとロバスト重みの両方を使用して、手順2~4を再実行
  4. 手順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}|}{6median(|r|)}\right) \end{equation}

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

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

3. STL分解のアルゴリズム

3.1 inner loopとouter loop

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

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

outer loopの処理の中にinner loopが含まれる構造となっており、それぞれのループの反復回数はパラメータで指定します。例えばinner loopを2回、outer loopを5回とすると、5回のouter loopの中で2回ずつinner loopが実行されます。このようにinner loopとouter 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)回目の分解の手順を説明し、S_{t}^{(k+1)}、T_{t}^{(k+1)}を計算していきます。

3.2.1 inner loop

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

\begin{equation} Y_{t} - T_{t}^{(k)} \end{equation}

2. 周期成分の平滑化
手順1でトレンドを除去したデータから、季節性の周期の長さごとに点を抽出し、Loess平滑化を適用します。季節性の周期の長さと平滑化のウィンドウサイズはパラメータで指定します。平滑化された周期成分をC_{t}^{(k+1)}とします。詳細は3.3で説明します。

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

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} - S_{t}^{(k+1)} \end{equation}

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

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. 事前に設定したパラメータ「季節性の周期の長さ」ごとに、原系列から点を抽出
  2. 抽出したデータに対して、Loess平滑化

例えば、日次データから曜日成分を分解する場合、曜日ごとに(7ずつ)点を抽出し、Loess平滑化を行います。

4. STL分解のパラメータ

Pythonの統計モデルライブラリStatsmodelsでは、statsmodels.tsa.seaasonal.STL[2]クラスを用いてSTL分解を実行できます。

4.1 statsmodels.tsa.seaasonal.STLクラスのパラメータ

  • period
    季節性の周期の長さ。例えば、月次データから年次成分を抽出したい場合は、12を指定します。

  • seasonal
    季節性をLoess平滑化する際のウィンドウサイズ。

  • trend
    トレンドをLoess平滑化する際のウィンドウサイズ。

  • low_pass
    ローパスフィルタをLoess平滑化する際のウィンドウサイズ。

  • seasonal_deg
    季節性をLoess平滑化する際の次数。

  • trend_deg
    トレンドをLoess平滑化する際の次数。

  • low_pass_deg
    ローパスフィルタで抽出した低周波成分をLoess平滑化する際の次数。

  • robust
    ロバスト重みを使用するかどうかを指定するフラグ。

  • seasonal_jump
    季節性をLoess平滑化する際のスキップ値。seasonal_jump=1とすると全ての点に対して平滑化値を計算し、2とすると1つ飛ばしで平滑化値を計算します。飛ばされた点は、前後の点から線形補間されます。大きな値を設定すると、平滑化値を計算する点が減るため、処理時間が短縮されます。

  • trend_jump
    トレンドをLoess平滑化する際のスキップ値。

  • low_pass_jump
    ローパスフィルタで抽出した低周波成分をLoess平滑化する際のスキップ値。

4.2 fit[3]メソッドのパラメータ

  • inner_iter
    inner loopの反復回数。
  • outer_iter
    outer loopの反復回数。

5. まとめ

本記事では、STL分解の仕組みとパラメータについて解説しました。ライブラリを使うことで簡単に試せるので、ぜひ使ってみてください。

脚注

  1. Robert B. Cleveland, William S. Cleveland, Jeam E. McRae, and Irma Terpenning. 'STL: A Seasonal-Trend Decomposition Procedure Based on LOESS', Journal of Official Statistics, vol6, No.1, pp.3-73, (1990) ↩︎

  2. statsmodels.tsa.seasonal.STL ↩︎

  3. statsmodels.tsa.seasonal.STL.fit ↩︎

Discussion