周期性分解のアルゴリズムを比較・検証

## はじめに

生成AIが脚光を浴びまくっている最近ですが、生成AIは万能ということはありません。生成AIはまだデータ分析、特に時系列データの扱いが得意ではなく、それは今でも統計モデルの役目です。

時系列データにおいて、**周期性（Seasonality）** は非常に重要な要素です。本記事では、時系列分析で重要な**周期性分解**に着目し、様々な周期性分解の手法を、アルゴリズム・実装方法・計算コストの観点で比較してみます。

## 周期性分解とは

![seasonal-decomp-concept.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3981597/2b79c875-69f2-4606-b5b2-680c4c858f10.png)

時系列 \$y\_t\$ を、周期成分 \$S\_t\$ と残差 \$R\_t\$ に分解する処理です。

```math
\underbrace{y_t}_{\text{オリジナル}} = \underbrace{S_t}_{\text{周期}} + \underbrace{R_t}_{\text{残差}}
```

なお時系列データによってはトレンド成分も重要な要素ですが、本記事では**周期性の抽出**に焦点を当て、トレンドの詳細な議論は割愛します。

## なぜ周期性分解が重要なのか

近年生成AIの性能は飛躍的に向上し、推論なども高いレベルで可能になってきています。しかし、それでもなお**周期性分解のような統計的・ドメイン知識に基づく分析**が不要になるわけではないと思っています。
例えば生成AIに時系列データをそのまま渡して「分析して。」と依頼すれば、「夏に上昇し、冬に下降」*など、それっぽいコメントは返ってくるでしょう。ところが実際の時系列分析では、具体的な数値で**周期成分を切り出し、トレンドや外れ値を定量的に扱う**ことが必要であり、これは依然として統計モデルの領域です。

実際に周期性分解が役立つ場面は次の2つだと思っています。

1. **周期パターンの可視化** — 時系列データから週次・年次などの周期性を取り出してプロットするだけでも、キャンペーン効果や季節要因を把握しやすくなるという大きな価値があります。たとえばECサイトの売上なら「毎週金曜に需要が伸びる」「夏に扇風機がピーク」などの洞察を誰でも一目で把握できます
2. **未来予測モデルの構築** — 時系列をそのまま学習させるよりも、あらかじめ周期成分を除去して残差系列だけをモデルに与えるほうが、モデルが学習すべきパターンが単純化され、精度向上や過学習抑制につながります

## データセット

今回検証に使う時系列データは、Kaggleで提供されている [Store Item Demand Forecasting Challenge](https://www.kaggle.com/competitions/demand-forecasting-kernels-only/data) です。

これは複数の店舗、商品における売り上げデータです。

* 商品: 50 種類 × 店舗: 10 ヶ所 → 500 系列の日次売上
* 期間: 2013‑01‑01 ～ 2017‑12‑31
* 粒度: 1 日
* 周期性：明瞭な 週 × 年 の二重周期

今回はこの内、店舗１・商品１の時系列データを使います。

実際にKaggleのWebページからデータをダウンロードして、`train.csv`を`data/` ディレクトリに配置します。

実際にプロットして眺めてみるだけでも、年の周期性が見て取れますね。

```python
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path

data_path = Path("data/train.csv")

# read csv
raw = pd.read_csv(data_path, parse_dates=["date"])

# store 1, item 1 を選択
sales = (raw.query("store==1 and item==1")
             .set_index("date")["sales"]
             .asfreq("D")
             .fillna(0)
             .rename("sales"))

# plot
ax = sales.plot(figsize=(10,3), title="Store‑1 / Item‑1  Daily Sales")
ax.set_ylabel("Sales")
plt.show()
```

![da16c71b-a755-4661-91b9-39c028caafe8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3981597/da26cac2-e0e5-45c1-adc7-0bc24c778d31.png)


## 今回取り上げる 4 手法

本記事では、私が独断と偏見で選んだ４つの周期性分解手法、Classical、SSA、MSTL、Prophetを比較します。

| No. | 手法                                 | 特徴                                                 | 多重周期対応 |
| --- | --------- | ---------------------------------- | -------------------------------------------------- |
| 1  | Classical Decomposition            | アルゴリズムがシンプルで計算も軽量                                  | 単一     |
| 2  | SSA (Singular Spectrum Analysis) | 固有値分解ベースでトレンドと季節・ノイズを同時抽出                          | 単一     |
| 3   | MSTL                               | STL を階層的に適用し、複数周期性を同時抽出。外れ値にロバスト。季節パターンのゆるやかな歪みに追従 | **複数** |
| 4  | Prophet                            | フーリエ基底で周期性を推定          | **複数** |

### Classical Decomposition

移動平均を用いた周期性分解の手法です。数ある周期性分解の中で、最も古典的で、最もアルゴリズムがシンプルな手法で、計算コストもO(N)と軽量です。

```python
def classical_decompose(series: pd.Series, period: int = 365):
    res = seasonal_decompose(series, model="additive", period=period)
    return res.seasonal, res.resid
```

周期長（例えば１年なら`365`）は明確に指定する必要があり、一度に単一周期しか指定できません。もし分析対象データが複数周期（例えば週次＆年次）をもつ場合でも、どちらかの周期しか指定できません。

また注意点として、一点でも欠損値（NaN）が混ざっているとこの手法は使えなくなります。事前に何らかの方法で欠損値を補完するか、素直に別の手法を検討するしかありません。

### SSA

SSA は時系列を特異値分解で分解し、トレンド・周期・ノイズを抽出する手法です。土木・物理・医用工学など、信号処理的な分野でよく使われるイメージです。

```python
def ssa_decompose(series: pd.Series):
    w = 500
    lower_frequency_bound = 0.005
    lower_frequency_contribution = 0.90

    X = series.values.reshape(1, -1)

    ssa = SingularSpectrumAnalysis(
        window_size=window,
        groups="auto",
        lower_frequency_bound=lower_frequency_bound,
        lower_frequency_contribution=lower_frequency_contribution
    )
    subseries = ssa.fit_transform(X)[0]
    _trend, seasonal, noise = subseries

    seasonal  = pd.Series(seasonal,  index=series.index, name="seasonal")
    resid     = pd.Series(noise,     index=series.index, name="resid")

    return seasonal, resid
```

他の手法と違って、周期長（年次の`365`や週次の`7`など）を指定せずに済む柔軟さが特徴です。

一方で `window_size` や寄与率の閾値（`lower_frequency_bound`, `lower_frequency_contribution`）などハイパーパラメータが多く、これらの選択によって結果が全然変わってくるので、他の手法と比べて正直使いづらいという感覚はあります。実際に現場で使う際はハイパーパラメータのチューニング作業が必要でしょう。

### MSTL

MSTL は STL (Seasonal‑Trend decomposition using Loess) を階層的に適用し、複数周期を同時抽出できる拡張版です。

※ STL は Loess 平滑化でトレンドと季節を分ける軽量な定番アルゴリズムです。

```python
def mstl_decompose(series: pd.Series, periods: tuple = (7, 365)):
    res = MSTL(series, periods=periods).fit()
    seasonal = res.seasonal.sum(axis=1)
    return seasonal, res.resid
```

MSTL は`periods=(7, 365)`のように複数周期をタプルで渡すだけで、それぞれにSTLを階層的に適用し、最終的に各季節成分を合算してくれます。

### Prophet

Prophet は日次粒度のビジネス時系列（売上・アクセス数・在庫推移など）を主なターゲットに設計されたモデルです。

フーリエ基底を用いて周期性の推定を行います。

周期性以外にも、トレンドの変化点検出・休日効果・追加回帰変数も考慮させることができる多機能モデルですが、今回は周期性分解の性能のみに注目していきます。

```python
def prophet_decompose(series: pd.Series,):
    dfp = series.reset_index().rename(columns={series.name: "y", series.index.name or "date": "ds"})
    m = Prophet(seasonality_mode="additive", weekly_seasonality=True, yearly_seasonality=True)
    m.fit(dfp)
    comp = m.predict(dfp)
    trend = pd.Series(comp["trend"].values, index=series.index)
    seas  = pd.Series(comp["additive_terms"].values, index=series.index)
    resid = series - (trend + seas)
    return seas, resid
```

Prophet では `weekly_seasonality=True`, `yearly_seasonality=True` といったフラグを立てるだけで週次・年次の季節性が自動で組み込まれます。周期長を数値で与える必要がなく、`add_seasonality()` を使えばそれ以外の任意の周期も簡単に追加できます。

## ソースコード

全ソースコードはこちらのリポジトリに記載しています：
https://github.com/ShimeiYago/seasonal-decomposition

Python の実行環境は `3.13.2`。Jupyter Notebookで実行してます。

依存ライブラリはリポジトリ内の `requirements.txt` にまとめてあります。

## 結果ハイライト

それでは各手法による周期性分解の結果をいくつかピックアップして紹介します。

### 周期成分プロットの比較

下図は各手法で抽出された365日（年次）周期成分の比較です。

![6c68dee9-b8ef-41c1-9195-384a0605de96.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3981597/457ba291-db8a-40a7-9bde-de2b4f4a60fc.png)

年次パターンを最もわかりやすく捉えているのはClassicalとProphetですね。特にProphetは年次カーブが最も滑らかに見えます。

Prophetの季節成分が太い帯のように見えるのは、年次に加えて週次成分も同時に抽出されているためです。週次の細かいギザギザが密に重なり、全体として幅広いバンド状に表示されます。

MSTLも年次＋週次を一括抽出できるためバンドが太めに見えます。さらにSTL系は季節パターンを Loess で滑らかに追従させるため、年度ごとに少しずつ形が変化し、ClassicalやProphet のように「1 つの固定カーブ」を前提としない点がグラフに表れています。

SSAの周期成分は、周期らしき波形は確認できますが、他と比べてぱっと見「ザラついた」曲線という印象です。先ほども解説した通り、SSAはハイパーパラメータ（`window_size` やグルーピング閾値）が多く、これらの値次第で結果は大きく変わってくるため、これらをチューニングすればもっと滑らかな曲線が得られるかもしれません。とはいえ、その追加コストが結構重いです...。

一方で、SSAの凄いところは、こちらからは周期長（年次は`365`）を指定していなくても周期成分をそれなりの精度で抽出できているところです。周期長がはっきりしないケースではSSAを選択し、はっきりしているケースではClassicalやProphetを選択するのが良さそうですね。

## 実行時間の比較

各手法をそれぞれ複数回実行し、実行時間をざっと比較してみました。

| 手法                                | 平均実行時間 |
| --------------------------------- | ------ |
| Classical Decomposition（周期長`365`） | 0.00 s |
| SSA                               | 1.39 s |
| MSTL                              | 0.33 s |
| Prophet                           | 0.17 s |


Classicalは計算がシンプルなので、実行時間も圧倒的に高速で、丸めると \`0.00s\`になります。その次に早いのがProphetでした。

最も遅いのがSSAでしたが、これは特異値分解の計算量が重いからです。

## 性能評価: 残差 ACF

7 日周期を除去した後の残差に対してACF（自己相関関数）を求めました。ACF（自己相関関数）は、時系列の値とその k ステップ後の値との相関を数値化したものです。値が 0 に近ければ周期性が抜けており、±1 に近いほど周期性が強く残っているということになります。周期性分解の性能を評価する一つの指標として参考になります。

![893159a7-78c1-47f6-bdf7-93a4da72f01f.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3981597/9c171e88-0c62-4c6e-9c1f-ca2ade66f84d.png)


ラグ７の自己相関の絶対値が小さくなるほど「7日の周期性が抜けた」と言えます。

| 手法        | ラグ 7 の自己相関 |
| --------- | ---------- |
| Classical | 0.027      |
| SSA       | 0.200      |
| MSTL      | ‑0.191     |
| Prophet   | 0.031      |

ClassicalとProphetは、残差の自己相関がほぼゼロに近く、良い塩梅で周期を除去できました。一方で、SSAやMSTLは自己相関が大きく残りました。SSAに関してはラグ７以外も多数のラグで自己相関が大きく残る結果となりましたが、SSAはハイパーパラメータが多いのでチューニング次第で結果は大きく変わりそうです。

## 結論まとめ

| 手法        | 可視化の見やすさ | 平均実行時間 | ラグ 7 ACF |
| --------- | -------- | ------ | -------- |
| Classical | ○        | 0.00 s | 0.027    |
| SSA       | △        | 1.39 s | 0.200    |
| MSTL      | △        | 0.33 s | ‑0.191   |
| Prophet   | **◎**    | 0.17 s | 0.031    |

**Prophetを選んでおけば間違いない**というのが今の所の見解です。Prophetは周期性分解としての性能が高く、実行時間も比較的短めで、良いとこ尽くしです。

Classicalも周期性分解としての性能は十分に高く、計算は非常に軽量で簡単に実装できるという点で依然現役だと思います。フロントエンドで軽量に実装したいときや、バルク処理したいときに使用機会がありそうです。ただClassicalは欠損値が一点でもあると使用不可能なので、その場合はやはりProphetを使いましょう。

今回はSSAとMSTLの結果は今ひとつでしたが、それは今回の日次売上データでの話です。データセットが変わるとまた結果は変わってくるかもしれません。
SSA は未知周期や高周波ノイズを含む信号系データに強いですし、MSTL は年度ごとに少しずつ形が変化する季節パターンを Loess で滑らかに追従できるため、複数周期かつゆっくりドリフトする電力需要・問い合わせ件数などのビジネスデータで真価を発揮するはずです。