はじめに：この記事でわかること
時系列データを扱う上で、「トレンド」や「季節性」といった変動要因を個別に捉えたい場面は頻繁に訪れます。
そんなときに強力な武器となるのがSTL分解 (Seasonal-Trend decomposition using Loess) である。
STL分解が何をするのか？ それは、以下の画像を見れば一発でわかる。一番上の生の時系列データが、その下に続く3つのシンプルな成分（トレンド、季節性、残差）に見事に分解されていることを理解できるだろう。
この記事では、STL分解がどのような仕組みで動いているのか、その根幹となる考え方を、以下のステップで説明する。
- 時系列分解の基礎知識
- 移動平均を使った簡易的な分解シミュレーション
- STL分解のアルゴリズム解説(論文ベースで)
この記事を読み終える頃には、あなたは以下をわからされている。
- 時系列データを分解するとはどういうことか。
- STL分解の全体像と、そのアルゴリズム。
- なぜSTL分解が多くの場面で強力なツールとなるのか（特に外れ値や欠損値に強い理由）。
STLを直感的に理解するために必要となる基礎知識
STL分解を学ぶ前に、まず「時系列分解」という考え方から。
1-1. 時系列データはさまざまな要素で構成されている
手元にある時系列データは、多くの場合、様々な要因が混ざり合った結果です。例えば、ある商品の月次売上データには、以下のような要素が同時に含まれているかもしれない。
- 会社の成長による長期的な増加（トレンド）
- 景気による数年単位の波（循環）
- 年末商戦による毎年12月の急増（季節性）
- 突発的なキャンペーンやノイズ（残差）
時系列分解とは、このように一体となった動きを、各構成要素に分離するプロセスのことである。
1-2. 時系列を構成する4つの主要要素
時系列データは、主に以下の4つの要素で構成されると考えられている。
- トレンド (Trend, T)
- データの長期的な増加・減少の動き。必ずしも直線とは限らない。
- 循環 (Cycle, C)
- 景気循環のように、固定ではない周期で現れる波。一般的に、季節性よりも周期が長く、振幅も大きい傾向がある。
- 季節性 (Seasonality, S)
- 1年や1週間など、暦に連動した固定周期で繰り返されるパターン。
- 残差 (Remainder, R)
- 上記3つの要素では説明できない、不規則な変動（ノイズ）。
注意点
STL分解の文脈では、循環(C)はトレンド(T)の一部として扱われることが多い。
1-3. 加法モデル
上記の要素の関係は、シンプルな足し算で表現できる。これを加法モデルと呼ぶ。
Y = T + S + R
この記事では、この加法モデルを基本として話を進めます。
※循環(C)はトレンド(T)の一部として扱っている
補足：乗法モデル
また、以下のように掛け算で表現する乗法モデルも考えられる。
Y = T * S * R
これは売上が増えるにつれて季節変動の幅も大きくなるようなケースに適しています。株価などもそうだと考えられる。
ただし、このモデルは以下の式のように両辺の対数を取ることで実質的に加法モデルになるので、基本的には加法モデルを理解しておけば問題ない。
log(Y) = log(T * S * R) = log(T) + log(S) + log(R)
また余談ですが、多くの時系列ライブラリでは、additive（加法）かmultiplicative（乗法）かを引数で指定できます。
STL分解を支える2つのコア技術
以下の２つの技術がコア技術となっている
1, 移動平均による季節性の除去
2, Loess(STL(Seasonal Trend Decompsition Loess)の名前にも使われている。)
2-1. 技術①：移動平均による季節性の除去
移動平均を取ることで季節性を除去、または抽出できる。
移動平均とは、ある時点 t の値を、その前後を含む一定期間（ウィンドウ）の値の平均で置き換える手法である。数式で説明すると以下。
\hat{T}_t = \frac{1}{m} \sum_{j=-k}^{k} y_{t+j}
季節性が偶数の場合
重要なのは、移動平均の計算に用いるwindows幅を季節周期と同じに設定することである。そうすると、季節パターンを効果的に除去したトレンド成分を抽出できる。
ただし、季節周期が偶数（例：4期、12ヶ月）の場合、一つ問題があります。
それは、移動平均を計算したとき、その中心が整数の時点からズレることである。
例えば、以下のような４周期のデータを考える。
これは4周期のデータで以下のような周期となっています。
このデータにウィンドウ幅4の移動平均を適用すると、結果の中心は t=2.5, 3.5, ... のように、時点の中間に来てしまいます。
そこで移動平均の移動平均を取ります。これは中心化移動平均と呼ばれ、2段階の処理で計算します。
まず、windows=4の移動平均を計算する
MA_0 = Nan
MA_1 = Nan
MA_1.5 = (y_0 + y_1 + y_2 + y_3) / 4 = 0 + 3 + 4 + -7 = 0
MA_2.5 = (y_1 + y_2 + y_3 + y_4) / 4 = 3 + 4 + -7 + 0 = 0
MA_3.5 = (y_2 + y_3 + y_4 + y_5) / 4 = 4 + -7 + 0 + 3 = 0
.
.
.
となる。
※今回はトレンド成分は０なので全て移動平均は０になる。
これを抽象化すると移動平均(windows=4)は
MA_t = (y_(t-2) + y_(t-1) + y_(t) + y_(t+1)) / 4
この時、移動平均の移動平均をとる。
windows=2とする。
2-MA_0 = Nan
2-MA_1 = Nan
2-MA_2 = Nan
2-MA_3 = (MA_2 + MA_3) / 2 = (0 + 0) / 2 = 0
2-MA_4 = (MA_3 + MA_4) / 2 = (0 + 0) / 2 = 0
2-MA_5 = (MA_4 + MA_5) / 2 = (0 + 0) / 2 = 0
.
.
.
※⚪︎-MA_tはデータ点tにおける２段階目でかつWindos=⚪︎の移動平均（MovingAverage、略してMA）
これを抽象化すると移動平均の移動平均化(windows=2)は
2-MA_t = (MA_(t-0.5) + MA_(t+0.5)) / 2
これを式展開する
\begin{align*}
\ (\text{2-MA_t}:移動平均の移動平均のこと)
&= \frac{1}{2} ( \text{MA}_{t-0.5} + \text{MA}_{t+0.5} ) \\ \\
&= \frac{1}{2} \left( \frac{1}{4}(y_{t-2} + y_{t-1} + y_{t} + y_{t+1}) + \frac{1}{4}(y_{t-1} + y_{t} + y_{t+1} + y_{t+2}) \right) \\ \\
&= \frac{1}{8} y_{t-2} + \frac{1}{4} y_{t-1} + \frac{1}{4} y_{t} + \frac{1}{4} y_{t+1} + \frac{1}{8} y_{t+2}
\end{align*}
この時、$y_{t-2}$と$y_{t+2}$は季節周期における同じポジションに位置する（上記添付画像のt=0とt=4の関係）。
この式が示すように、2x4中心化移動平均は、実際には5期間の重み付き移動平均と等価である。中心の3点 $(y_{t-1}, y_t, y_{t+1})$には 1/4、両端の $y_{t-2} と y_{t+2}$には 1/8 の重みが与えられます。
このように、中心化移動平均を使うことで、時点のズレを解消し、滑らかなトレンド成分 T を抽出できる。
そして、私たちは時系列に加法モデルY = T + S + Rを想定しているから、元時系列データからこの抽出したトレンドを「引き算」によって季節性の成分を分離することができる。
つまり
Y - \hat{T} = S + R
である。
季節性が奇数の場合
季節周期が奇数（例：7日）の場合は、中心のズレは発生しない。そのため、単純な移動平均を1回計算するだけでよく、移動平均の移動平均は不要。
\begin{align*}
\hat{T}_t &= \frac{1}{m} \sum_{j=-k}^{k} y_{t+j} \\ \\
&= \frac{1}{m} (y_{t-k} + \dots + y_{t-1} + y_t + y_{t+1} + \dots + y_{t+k})
\end{align*}
補足
Y - S のように、ある成分を引き算して取り除くことを「季節調整」と呼びます。STL分解は、この「引き算による成分分離」というアイデアを、より洗練された形で反復的に利用する手法です。
2-2. 技術②：Loess（局所重み付き回帰）
STLの心臓部であり、その優れた特性の源となっているのが Loess である。
Loessは、データ全体に1つの数式を当てはめるのではなく、「ある点の周辺（近傍）のデータだけを使って、その点における値を推定する」というノンパラメトリックな平滑化手法である。具体的には、近傍内の点に「近いほど重い」重みを付け、重み付き最小二乗法で局所的な多項式（直線や二次曲線）を当てはめる。つまり、重み付けしたデータで回帰を行なっている。
Loessは、以下の3ステップで平滑化された値を計算する。
Loessのステップ
ステップ1：近傍の決定
まず、平滑化したい点 x の周辺から、計算に使うデータ点を q 個選ぶ。この q 個のデータ群を「近傍」と呼ぶ。
ステップ2：距離の正規化
次に、近傍内の各点と x との距離を計算します。そして、その距離を「x から最も遠い近傍点までの距離」で割ることで、距離を0から1の範囲に正規化します。この正規化された距離を u とします。
u_i = \frac{|x_i - x|}{\lambda_q(x)} \quad \text{where} \quad \lambda_q(x) = \max_{k \in \text{Neighborhood}} |x_k - x|
ステップ3：近傍点への重み付け
正規化された距離 u をトリキューブ重み関数 (tricube weight function) に入力し、各点の重みを計算します。x に近い点（uが0に近い）ほど重みは1に近づき、遠い点（uが1に近い）ほど重みは0に近づきます。
W(u) =
\begin{cases}
(1 - u^3)^3 & \text{for } 0 \le u < 1 \\
0 & \text{for } u \ge 1
\end{cases}
可視化すると以下のような関数になります。
なぜトリキューブ関数なのか？
主に以下の2つの理由からである。
- 中心付近が平坦: x に最も近い点たちの重みがほぼ均等になるため、推定が安定する。
- 端で滑らかにゼロになる: 近傍の端で重みが滑らかにゼロになるため、最終的に得られる平滑化曲線がカクカクせず、非常に滑らかになる。
ステップ4：重み付き局所回帰
最後に、計算した重みを考慮して、近傍データに多項式（通常は直線 d=1 または二次曲線 d=2）を当てはめます（重み付き最小二乗法）。この局所的にフィットさせた曲線の、点 x における値が、最終的な平滑化された値となります。
Loessのパラメータ
Loessには主に2つのパラメータがあります。
-
q (近傍のサイズ): 大きくするほど、結果は滑らかになる。
- 論文：
-
To use loess, d and q must, of course, be chosen. The choices, in the context of STL,
will be discussed in detail in this section and in Section 3. As q increases, g(x) becomes smoother. As q tends to infinity, the v₁(x) tend to one and g(x) tends to an ordinary least-squares polynomial fit of degree d.
-
d (多項式の次数): データの局所的なカーブが急な場合は d=2、緩やかなら d=1 が適している。
- 論文：
-
Journal of Official Statistics Taking d = 1 is reasonable if the underlying pattern in the data has gentle curvature. But if the pattern has substantial curvature, for example, peaks and valleys, then d = 2 is a better choice.
ここからわかるのはデータが複雑な場合、dは高い次数にして良い。もし、qを大きくするとトレンドが滑らかになる。そしてqを無限大にするとそれは単純な次数dの多項式回帰分析と等価になる。
LoessがSTLにもたらす強み：欠損値への耐性
STLにおいてLoessが果たす役割は、単なる平滑化に留まらない。
Loessは要するに回帰分析なので、存在しないデータに対しても、その欠損値を推定することで補完することができる。これがSTL分解が欠損値に強い理由である。
シミュレーション
データの準備
今回は、以下の加法モデルに従うシミュレーションデータ（Y）を241期間分、生成します。
Y = trend + low_freq_cycle + seasonality + outlier + white_noise
各成分は、以下のように設定した。
trend:Tのこと。右肩上がりの直線的なトレンドlow_freq_cycle:周期96のゆったりとした循環変動seasonality:周期12で、複数の波を合成した複雑な季節パターンoutlier:特定の時点 (t=150) に意図的に挿入した外れ値white_noise:予測不可能なランダムなノイズ（残差）
詳細はコードを参照すること。
以下が生成した時系列データとなる。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# --- Parameters of the simulation ---
# These variables control the shape of the time series.
n_points = 241 # Number of data points (from t=0 to t=240).
noise_amplitude = 5.0 # Controls the magnitude of the random noise.
outlier_amplitude = 40.0 # Controls how extreme the outliers are.
# --- 1. Generate the time axis (t) ---
# Creates a simple numpy array from 0 to 240.
time = np.arange(n_points)
# --- 2. Generate the Trend component ---
# A simple, increasing linear trend using the formula y = 0.5*t + 10.
trend = 0.5 * time + 10
# --- 3. Generate the Low-Frequency Cycle component ---
# A sine wave with a long period of 48, representing a business cycle.
# The formula is Amplitude * sin(2 * pi * t / Period).
low_freq_cycle = 20 * np.sin(2 * np.pi * time / 96)
# --- 4. Generate the Seasonality component ---
# A cosine wave with a shorter, fixed period of 12, representing yearly seasonality.
# --- 1. Add harmonics to the seasonality ---
# Base seasonality with a period of 12
base_seasonality = 15 * np.cos(2 * np.pi * time / 12)
# A half-year harmonic with a period of 6
harmonic_6_months = 5 * np.sin(2 * np.pi * time / 6)
# A quarterly harmonic with a period of 3
harmonic_3_months = 3 * np.cos(2 * np.pi * time / 3)
# Combine them to create a more complex seasonal pattern
seasonality = base_seasonality + harmonic_6_months + harmonic_3_months
# --- 5. Generate White Noise (Remainder) ---
# Generates random numbers from a standard normal distribution (mean=0, std=1).
# np.random.seed(42) ensures that the "random" noise is the same every time the code is run,
# making the simulation reproducible.
np.random.seed(42)
white_noise = np.random.normal(0, 1, n_points) * noise_amplitude
# --- 6. Add Outliers ---
# First, create an array of zeros.
# Then, insert large positive or negative values at specific time points (50, 150, 200)
# to simulate anomalous events.
outliers = np.zeros(n_points)
outlier_indices = [150]
for i in outlier_indices:
sign = np.random.choice([-1, 1]) # Randomly choose if the outlier is positive or negative
outliers[i] = outlier_amplitude * sign
# --- 7. Combine all components ---
# The final time series is created by adding all the individual components together,
# following the additive model: Y = T + Cycle + S + Noise + Outliers.
time_series_final = trend + low_freq_cycle + seasonality + white_noise + outliers
# --- Create a pandas DataFrame for better handling and visualization ---
df_simulation = pd.DataFrame({
'time': time,
'trend': trend,
'low_freq_cycle': low_freq_cycle,
'seasonality': seasonality,
'white_noise': white_noise,
'outliers': outliers,
'time_series_final': time_series_final
})
# --- Display the first 5 rows of the DataFrame ---
print("--- Simulation DataFrame Head ---")
print(df_simulation.head())
# --- Visualize the final, combined time series ---
plt.figure(figsize=(15, 6))
plt.plot(df_simulation['time'], df_simulation['time_series_final'])
plt.title('Simulated Time Series (All Components Combined)', fontsize=16)
plt.xlabel('Time (t)', fontsize=12)
plt.ylabel('Value', fontsize=12)
plt.grid(True)
plt.show()
# --- Visualize each component separately for clarity ---
fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(15, 12), sharex=True)
fig.suptitle('Individual Components of the Simulated Time Series', fontsize=16, y=0.93)
axes[0].plot(df_simulation['time'], df_simulation['trend'], color='blue')
axes[0].set_ylabel('Trend')
axes[1].plot(df_simulation['time'], df_simulation['low_freq_cycle'], color='green')
axes[1].set_ylabel('Cycle (Low Freq)')
axes[2].plot(df_simulation['time'], df_simulation['seasonality'], color='orange')
axes[2].set_ylabel('Seasonality')
axes[3].plot(df_simulation['time'], df_simulation['white_noise'], color='grey', alpha=0.7)
axes[3].set_ylabel('White Noise')
axes[4].stem(df_simulation['time'], df_simulation['outliers'], linefmt='r-', markerfmt='ro', basefmt=' ')
axes[4].set_ylabel('Outliers')
axes[4].set_ylim(-outlier_amplitude-10, outlier_amplitude+10) # Adjust y-axis for better visibility
axes[5].plot(df_simulation['time'], df_simulation['time_series_final'], color='black')
axes[5].set_ylabel('Final Series')
axes[5].set_xlabel('Time (t)')
for ax in axes:
ax.grid(True)
plt.tight_layout(rect=[0, 0, 1, 0.9])
plt.show()
元時系列からトレンドを除去する
加法モデル（Y = T + S + R）を前提としているため、季節性（S）と残差（R）を得るためには、トレンドTを推定し、単純な引き算（Y - T）をすれば良い。
ここでトレンド T を推定するために、中心化移動平均を利用します。今回のデータの季節周期は12なので、2x12中心化移動-平均を計算する。これにより、周期12の季節性がきれいに除去され、滑らかなトレンド成分（と循環成分）が抽出される。下のグラフの赤線がyを中心化移動平均した結果である。見事にトレンドを抽出しているように見える。
各時系列データの説明は以下。
- 灰色の線: 元の時系列データ (Y)
- 赤い線: 中心化移動平均によって抽出されたトレンド成分 (T (+ C))
- 緑の線: 元の系列からトレンド成分を引き算した結果 (Y - T) 。これは、季節成分と残差 (S + R) を表す
補足：トレンドと循環
厳密には、ここで抽出した赤い線には、長期的なトレンドと、景気循環のような低周波の循環の両方が含まれてる。現代の多くの時系列分解（STLを含む）では、この2つを合わせて「トレンド成分」と見なすことが一般的であるらしい。
# To center it, we can shift the result backward by (window-1)/2 = 5.5 periods.
# A simpler way is to use the 'center=True' argument.
# This calculates a 2 x 12 MA automatically when the window is even.
df_simulation['Centered_MA'] = df_simulation['time_series_final'].rolling(window=12, center=True).mean()
# --- ADDITION: Calculate the Detrended Series ---
# This is the new part. Just subtract the Centered_MA from the original series.
df_simulation['Detrended_Series'] = df_simulation['time_series_final'] - df_simulation['Centered_MA']
# --- MODIFICATION: Update the plot to include the new series ---
# Let's plot to compare all three
plt.figure(figsize=(15, 8))
# Plot the original series (lightly, in the background)
plt.plot(df_simulation['time'], df_simulation['time_series_final'], label='Original Series', alpha=0.4, color='grey')
# Plot the Trend-Cycle component
plt.plot(df_simulation['time'], df_simulation['Centered_MA'], label='Centered MA (Trend-Cycle)', color='red', linewidth=2.5)
# Plot the Detrended (Seasonal + Remainder) component
plt.plot(df_simulation['time'], df_simulation['Detrended_Series'], label='Detrended (Seasonal + Remainder)', color='green', linewidth=2)
# Add a horizontal line at y=0 for reference for the detrended series
plt.axhline(0, color='black', linestyle='--', linewidth=1)
plt.title('Decomposition using Centered Moving Average', fontsize=16)
plt.xlabel('Time (t)', fontsize=12)
plt.ylabel('Value', fontsize=12)
plt.legend(fontsize=12)
plt.grid(True)
plt.show()
季節性を推定する
次に、トレンドが除去された系列（S + R）から、季節成分 S を推定する。
最もシンプルな方法は、同じ季節ごとにグループ化し、その平均値を計算することです。これをサイクルサブ系列の平滑化と呼びます。
例えば、12ヶ月周期のデータなら、以下のように処理します。
1, 全ての「1月」のデータを集めてグループを作る
2, 全ての「2月」のデータを集めてグループを作る
3, (...これを12月まで繰り返す)
4, 各グループの平均値を計算する。この平均値が、その月の代表的な季節効果となる
イメージしづらい場合は補足セクションのCycle-Subseries Plotを参照。
この方法で季節成分 S を推定し、さらに (S + R) - S を計算して残差 R を抽出した結果が、以下のグラフです。
- 上のグラフ:
- 緑色の線が、トレンド除去後の系列 (S + R、前のステップで抽出したもの) 。
- 紫色の線が、上記の平均化によって推定された季節成分 (S) 。きれいな繰り返しパターンになっているのが分かる。
- 下のグラフ:
- 青色の線が、上のグラフの緑色の線から紫色の線を引き算した結果、すなわち残差 (R)
これで、元の系列 Y を T, S, R の3つの成分に分解することができた。
# --- This part is the same as before, assuming it has been run ---
df_simulation['month'] = (df_simulation['time'] % 12) + 1
average_seasonal_effect = df_simulation.groupby('month')['Detrended_Series'].mean()
df_simulation['Simplified_Seasonal'] = df_simulation['month'].map(average_seasonal_effect)
# --- ADDITION: Calculate the Remainder component ---
# The remainder is what's left after subtracting the smoothed seasonal
# component from the detrended series. R = (S+R) - S.
df_simulation['Remainder'] = df_simulation['Detrended_Series'] - df_simulation['Simplified_Seasonal']
# --- MODIFICATION: Create a figure with two subplots ---
fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(15, 9), sharex=True)
fig.suptitle('Cycle-Subseries Smoothing and Resulting Remainder', fontsize=18, y=0.96)
# --- Plot 1: Comparison of Seasonal Components ---
axes[0].plot(df_simulation['time'], df_simulation['Detrended_Series'],
label='Detrended (Seasonal + Remainder)', alpha=0.5, color='green')
axes[0].plot(df_simulation['time'], df_simulation['Simplified_Seasonal'],
label='Smoothed Seasonal (Simplified Method)', color='purple', linewidth=2.5)
axes[0].set_ylabel('Value')
axes[0].legend()
axes[0].grid(True)
axes[0].set_title('Step 1: Smoothing the Detrended Series', fontsize=14)
# --- Plot 2: The Remainder Component ---
axes[1].plot(df_simulation['time'], df_simulation['Remainder'],
label='Remainder (Noise + Outliers)', color='blue', alpha=0.8)
axes[1].axhline(0, color='black', linestyle='--', linewidth=1) # Add a zero line
axes[1].set_ylabel('Value')
axes[1].set_xlabel('Time (t)')
axes[1].legend()
axes[1].grid(True)
axes[1].set_title('Step 2: Resulting Remainder after Subtraction', fontsize=14)
plt.tight_layout(rect=[0, 0, 1, 0.93])
plt.show()
# --- Display the DataFrame with the new Remainder column ---
print("\n--- DataFrame with Remainder Component ---")
# Displaying specific columns for clarity
print(df_simulation[['time', 'Detrended_Series', 'Simplified_Seasonal', 'Remainder']].head(15))
分解結果の「答え合わせ」
今回私たちはシミュレーションデータを活用しているので、分解結果と実際の各成分の比較ができる。
比較した結果が以下のグラフになる。
上のグラフでは、青い実線が「真の」成分、オレンジの破線が「私たちが推定した」成分を示している。
グラフの通り、トレンド、季節性、残差のいずれにおいても、オレンジの破線は青い実線に非常によく追従していることがわかる。特に、残差に含まれる3つの大きな外れ値（t＝150付近、-40くらいのところ）も、ほぼ完璧に捉えられている。
この結果から、移動平均を主体としたシンプルな手法でも、時系列の主要な構成要素をかなり高い精度で分解できることがわかる。
そして、以降の章で説明するが、STL分解も根本にある考え方は同じである。
# --- For this comparison, we need the original components.
# Let's assume df_simulation also contains the original component columns for this "ground truth" comparison.
# If not, you would need to regenerate the data with all component columns.
# --- Create comparison plots ---
fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(15, 10), sharex=True)
fig.suptitle('Comparison: Custom Estimate vs. Original Simulation Components', fontsize=18, y=0.95)
# --- 1. Trend Comparison ---
# Note: STL's trend captures both the linear trend and the low-frequency cycle.
original_trend_cycle = df_simulation['trend'] + df_simulation['low_freq_cycle']
axes[0].plot(df_simulation['time'], original_trend_cycle, label='Original Trend + Cycle (Ground Truth)', color='blue', alpha=0.7)
axes[0].plot(df_simulation['time'], df_simulation['Centered_MA'], label='Estimated Trend (Centered MA)', color='orange', linestyle='--')
axes[0].set_ylabel('Trend')
axes[0].legend()
# --- 2. Seasonal Comparison ---
axes[1].plot(df_simulation['time'], df_simulation['seasonality'], label='Original Seasonality (Ground Truth)', color='blue', alpha=0.7)
axes[1].plot(df_simulation['time'], df_simulation['Simplified_Seasonal'], label='Estimated Trend (Centered MA)', color='orange', linestyle='--')
axes[1].set_ylabel('Seasonal')
axes[1].legend()
# --- 3. Remainder/Residual Comparison ---
original_remainder = df_simulation['white_noise'] + df_simulation['outliers']
axes[2].plot(df_simulation['time'], original_remainder, label='Original Remainder (Ground Truth)', color='blue', alpha=0.7)
axes[2].plot(df_simulation['time'], df_simulation['Remainder'], label='STL Estimated Residual', color='orange', linestyle='--')
axes[2].set_ylabel('Residual')
axes[2].set_xlabel('Time (t)')
axes[2].legend()
for ax in axes:
ax.grid(True)
plt.tight_layout(rect=[0, 0, 1, 0.93])
plt.show()
STL分解のアルゴリズム
全体構造：頑健性を実現する「二重ループ」
STLのアルゴリズムは、**「内側ループ」と「外側ループ」**から成る、入れ子構造の二重ループで構成される。
- 内側ループ: トレンド成分と季節成分の分離・更新を担うメインの処理ループ
- 外側ループ: 内側ループ全体を包み込み、外れ値の影響を軽減するための処理を行うループ
STLは、この二重ループを反復計算することで、各成分の推定精度を段階的に向上させていく。
外側ループの役割：外れ値への対処
外側ループは、外れ値問題を解決し、分解の頑健性 (robustness) を確保するために存在する。
（※外れ値に対処する必要性を説明するための説明を補足のセクションに記述している。上記の章で行ったシミュレーションは外れ値に対応していないアルゴリズムであるが、それが外れ値にどのくらい弱いかを示している。「外れ値処理を行わない簡易手法の限界」の部分。）
このプロセスは、各反復において、以下のステップで頑健性重み (robustness weights) \(\rho_v\) を計算・更新することである。
ステップ1：残差 R_v の計算
まず、内側ループによって推定されたトレンド $(T_v)$ と季節性 $(S_v)$ を用い、各時点 v における残差 $(R_v)$ を算出する。
R_v = Y_v - T_v - S_v
この $(R_v)$ の絶対値が大きいほど、そのデータ点は外れ値である可能性が高いと見なされる。
ステップ2：残差スケールhの計算
次に、「どの程度の大きさの残差を異常と見なすか」の基準となるスケール $(h)$ を計算する。$(h)$ は、全残差の絶対値の中央値を用いて、以下のように定義される。この 6 という係数は経験的に定められた定数だと思われます。
h = 6 \times \text{median}(|R_v|)
$(h)$ 自体は中央値を用いているので外れ値の影響を受けにくい指標である。
ステップ3：頑健性重みの計算
各点の残差 $(|R_v|)$ をスケール $(h)$ で正規化し、それをバイスクエア重み関数 (bisquare weight function) $(B(u))$ に入力することで、最終的な頑健性重み $(\rho_v)$ を得る。
\rho_v = B\left(\frac{|R_v|}{h}\right)
ここで、バイスクエア重み関数 $(B(u))$ は以下のような区分関数である。
B(u) =
\begin{cases}
(1 - u^2)^2 & \text{for } 0 \le u < 1 \\
0 & \text{for } u \ge 1
\end{cases}
この関数により、残差が $(h)$ を超える点（外れ値）の重み $(\rho_v)$ は、厳密に 0 となる。
ステップ4：重みを考慮した内側ループの再実行
算出された頑健性重み $(\rho_v)$ は、次の内側ループに渡される。
内側ループにおけるLoess平滑化（ステップ2のサイクルサブ系列平滑化と、ステップ6のトレンド平滑化）では、通常の近傍重み(tricube weight functionで計算したもの) $(v_i(x))$ にこの頑健性重み $(\rho_i)$ が掛け合わされ、新しい重み $(\rho_i \times v_i(x))$ として使用される。
これにより、重みが0の点（明らかな外れ値）は平滑化の計算から完全に無視され、推定に影響を与えなくなる。
STLは、この「内側ループ実行 → 残差チェック → 重み更新」というサイクルを指定された回数だけ繰り返すことで、外れ値に対して非常に頑健な分解結果を生成する。
内側ループの役割：トレンドと季節性の分離
内側ループは、トレンド成分と季節成分を分離・推定するためのメインエンジンであり、以下の6つのステップで構成される。このループが1回実行されるごとに、両成分の推定値が更新される。
ここでは、k回目の反復が完了し、(k+1)回目の更新を行うプロセスを記述する。
Step 1: 脱トレンド (Detrending)
現在のトレンド推定値 $(T_v^{(k)})$ を、元の時系列 $(Y_v)$ から差し引く。
Y_v - T_v^{(k)}
これにより、「季節性 + 残差」(S+R) に相当する、脱トレンド系列が得られる。
なお、最初の反復 (
k=0) では、\(T_v^{(0)} = 0\)として計算を開始する。論文でも、この初期値でうまく機能することが示されていた。
論文：
To carry out Step 1 on the initial pass through the inner loop we need starting values, T, for the trend component. Using T = 0 works quite well. The trend becomes part of the smoothed cycle-subseries, C, but is largely removed in Step 4, the detrending.
Step 2: サイクルサブ系列の平滑化 (Cycle-subseries Smoothing)
Step 1で得られた脱トレンド系列を、周期内の同じ位置にあるデータごとにグループ化する（例：月次データなら1月の系列、2月の系列…）。この各サイクルサブ系列を個別にLoessで平滑化する。
この処理により、ノイズが平滑化された一時的な季節成分 $(C_v^{(k+1)})$ が得られる。
Step 3: ローパスフィルタリング (Low-Pass Filtering)
Step 2で得られた $(C_v^{(k+1)})$ には、季節成分に本来含まれるべきでない低周波の変動（トレンドの漏れ出し）が不純物として混入している可能性がある。
$(C_v^{(k+1)})$ に、複数の移動平均とLoessを組み合わせたローパスフィルタを適用し、この低周波成分 $(L_v^{(k+1)})$ のみを抽出する。
Step 4: 季節成分の純化
一時的な季節成分 $(C_v^{(k+1)})$ から、Step 3で抽出した低周波成分 $(L_v^{(k+1)})$ を差し引く。
S_v^{(k+1)} = C_v^{(k+1)} - L_v^{(k+1)}
これにより、トレンドの影響が取り除かれた、より純粋な季節成分 $(S_v^{(k+1)})$ が得られる。
Step 5: 脱季節化 (Deseasonalizing)
Step 4で更新された季節成分 $(S_v^{(k+1)})$ を、元の時系列 $(Y_v)$ から差し引く。
Y_v - S_v^{(k+1)}
これにより、「トレンド + 残差」(T+R) に相当する、脱季節系列が得られる。
Step 6: トレンドの平滑化 (Trend Smoothing)
Step 5で得られた脱季節系列をLoessで平滑化する。Loessは高周波のノイズ（残差 R）を除去する効果があるため、この処理によって滑らかなトレンド成分 $(T_v^{(k+1)})$ が得られる。
内側ループは、この一連の処理を指定された回数だけ繰り返す。これにより、トレンド成分と季節成分の推定値は、相互に更新されながら徐々に収束していく。
補足
Cycle-Subseries Plot
これを見ると固定された季節性がどのように計算されているのかイメージがつきやすい。
外れ値処理を行わない簡易手法の限界
STL分解のように外れ値に対処するための処理(Outer Loop)を行わないと分解そのものがうまくいかなくなることがある。以下シミュレーションの結果ではあえて外れ値を大きくしている。
我々のシミュレーションのロジックの結果
季節性の推定含めて多くの部分で歪みが生じている。青線が真の値、黄色い破線がシミュレーション結果。
PythonのStatsmodelのSTL分解の結果
外れ値があっても季節性の推定などは問題なさそうである。
参考文献
1,
2,
paper, [A Seasonal-Trend Decomposition Procedure Based on Loess], 1990,
Robert B. Cleveland, William S. Cleveland, Jean E. McRae,2 and Irma Terpenning
3,
コード全体について：
4,
# シミュレーションデータを生成すること
# Assuming 'df_simulation' DataFrame is already loaded and contains the 'time_series_final' column.
# --- 1. Perform STL Decomposition with a Perfectly Fixed Seasonal Pattern ---
# To achieve a perfectly fixed seasonal pattern, we set two key parameters:
#
# seasonal_deg=0:
# This tells the cycle-subseries smoother to fit a constant (a weighted average)
# instead of a line. It prevents the seasonal pattern from having its own trend.
#
# seasonal=1001: (or any large odd number greater than the number of cycles in the data)
# This forces the Loess smoother to use the *entire* time series for its "local" window.
# This turns the smoother into a global one, finding a single, average seasonal
# effect for each period (e.g., one value for all Januaries) that does not change over time.
stl_fixed_seasonality = sm.tsa.STL(df_simulation['time_series_final'],
period=12,
robust=True,
seasonal_deg=0,
seasonal=1001)
result = stl_fixed_seasonality.fit()
# --- 2. Create a Customized Visualization for the Results ---
# This plot uses separate y-axes for each component to ensure details are not lost.
fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(12, 10), sharex=True)
fig.suptitle('STL Decomposition with Perfectly Fixed Seasonality', fontsize=18, y=0.96)
# Plot 1: Original Series (Observed)
axes[0].plot(result.observed.index, result.observed, color='black', label='Observed')
axes[0].set_ylabel('Observed', fontsize=12)
axes[0].legend(loc='upper left')
# Plot 2: Trend Component
axes[1].plot(result.trend.index, result.trend, color='blue', label='Trend')
axes[1].set_ylabel('Trend', fontsize=12)
axes[1].legend(loc='upper left')
# Plot 3: Seasonal Component
# This plot will now show a perfectly repeating pattern.
axes[2].plot(result.seasonal.index, result.seasonal, color='green', label='Seasonal (Fixed)')
axes[2].set_ylabel('Seasonal', fontsize=12)
axes[2].legend(loc='upper left')
# Plot 4: Residual Component
axes[3].plot(result.resid.index, result.resid, color='red', label='Residual', marker='.', linestyle='None')
axes[3].axhline(0, color='black', linestyle='--', linewidth=1)
axes[3].set_ylabel('Residual', fontsize=12)
axes[3].set_xlabel('Time', fontsize=12)
axes[3].legend(loc='upper left')
# Improve the overall layout
for ax in axes:
ax.grid(True)
plt.tight_layout(rect=[0, 0, 1, 0.94])
plt.show()
# --- 3. (Optional) Store the results back into the DataFrame for further analysis ---
df_simulation['stl_trend_fixed'] = result.trend
df_simulation['stl_seasonal_fixed'] = result.seasonal
df_simulation['stl_resid_fixed'] = result.resid
print("\n--- First 15 rows of the decomposition results ---")
print(df_simulation[['time_series_final', 'stl_trend_fixed', 'stl_seasonal_fixed', 'stl_resid_fixed']].head(15))