MSTL分解

04.11.5

MSTL分解

最終更新 2020-03-11 読了時間 4 分

まとめ

MSTL は STL 分解（LOESS による季節・トレンド分解）を拡張し、周期 24 と周期 168 のように複数の季節パターンを同時に取り出せる手法です。STL を周期ごとに繰り返し当てる仕組みで、statsmodels 0.14.0 以降の MSTL(series, periods=(24, 168)).fit() の 1 行で使えます。このページでは np.sin と np.cos の 2 本の波に対数トレンドを重ねた 1000 日の人工データを分解し、.seasonal が 2 列になって周期ごとに分かれること、.resid がほぼ 0 に収まることを 4 段のプロットで確かめます。

直感 #

電力消費のデータには「1 日の中で昼に高く夜に低い」24 時間周期と、「平日に高く週末に低い」7 日周期が同時に乗っています。STL は周期を 1 つしか指定できないので、24 時間で分解すると 7 日周期がトレンドか残差に紛れ込み、7 日で分解すると 24 時間周期が残差に残ります。MSTL は「まず短い周期 24 で STL を当てて季節成分を引き、残りに長い周期 168 で STL を当てる」という手順を繰り返して、周期ごとの季節成分を別々の列に取り出します。

観測値周期 24 で STL 季節成分 1 を引く周期 168 で STL 季節成分 2 を引く残りをトレンドと残差に分ける

複数の周期性が重なったデータをトレンド・季節/周期性・残差に分解する #

MSTL は STL 分解（LOESS による季節・トレンド分解）手法を拡張したもので、複数の季節パターンを持つ時系列の分解が可能です。statsmodels の version==0.14.0 以降でのみ使用できます。引数の詳細はドキュメントを確認してください。

K. Bandura, R.J. Hyndman, and C. Bergmeir (2021) MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns. arXiv preprint arXiv:2107.13462.

Anacondaを仮想環境として使用している場合、conda install -c conda-forge statsmodelsでインストールされるものは0.13.Xとなっています(2022/11/1時点)。その場合、作業中の仮想環境の中で以下のコマンドを使用して最近のバージョンをインストールしてください。

pip install git+https://github.com/statsmodels/statsmodels

使うライブラリ #

1
2
3
4
5
6
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

from statsmodels.tsa.seasonal import MSTL

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

6 行目MSTL は statsmodels 0.14 以降にしかない
from statsmodels.tsa.seasonal import MSTL が ImportError になるなら、インストールされている statsmodels が 0.13 以前です。import statsmodels; statsmodels.__version__ で確認し、上の tip の手順で更新します。同じモジュールに単一周期用の STL と古典的な seasonal_decompose も入っています。
1 行目pandas の date_range で日次インデックスを作る
pd.date_range("2018-01-01", periods=1000, freq="D") で 1000 日分の日付を作り、観測値と組にして DataFrame にします。MSTL 自体は Series でも ndarray でも受け取れますが、DatetimeIndex を付けておくと分解結果のプロットの横軸が日付になります。
4 行目seaborn は最初の折れ線にだけ使う
sns.lineplot は元データの形を 1 枚確認するためだけに使います。分解結果の 4 段プロットは matplotlib の subplots で描くので、seaborn を省いても後半のコードは動きます。

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：元データ = トレンド + 季節成分 + 残差。残差に残った大きな山が異常の候補。

サンプルデータを作成 #

周期の異なる三角関数を 2 本足し、さらに i % 14 のノコギリ波と np.log(i**3 + 1) + np.sqrt(i) の上昇トレンドを重ねます。答えを知っているデータなので、分解後の .seasonal 2 列がそれぞれどの波に対応するかを確かめられます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
date_list = pd.date_range("2018-01-01", periods=1000, freq="D")
value_list = [
    10
    + i % 14
    + 2 * np.sin(10 * np.pi * i / 24)
    + 5 * np.cos(2 * np.pi * i / (24 * 7)) * 2
    + np.log(i**3 + 1)
    + np.sqrt(i)
    for i, di in enumerate(date_list)
]

df = pd.DataFrame(
    {
        "日付": date_list,
        "観測値": value_list,
    }
)

df.head(10)

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

4 行目i % 14 は周期 14 のノコギリ波
0 から 13 まで増えて 0 に戻る成分で、振幅 13 と大きめです。後の periods=(24, 168) には 14 を含めていないので、この成分は季節成分として取り出されず、トレンドか残差のどちらかに混ざります。残差に周期的な山が残っていたらこれが原因です。
5 行目sin の周期は 4.8、振幅 2
np.sin(10 * np.pi * i / 24) は \(2\pi\) で 1 周するので周期は \(24 \times 2 / 10 = 4.8\) 日です。整数周期ではないため、後で指定する periods=24 の中に 5 周期分が入ることになり、MSTL は「24 日の中で繰り返す形」としてこの波を季節成分 1 に取り込みます。
6 行目cos の周期は 168、振幅 10
np.cos(2 * np.pi * i / (24 * 7)) は周期 168（= 24 × 7）で、5 * ... * 2 なので振幅は 10 です。1000 日の中に約 6 周期しか入らない、ゆっくりした大きなうねりで、periods の 2 つ目 24 * 7 がこれに対応します。
7-8 行目log と sqrt で上に凸のトレンドを作る
np.log(i**3 + 1) は序盤で急に増えて後は緩やか、np.sqrt(i) も同じ形で、合わせて 1000 日で約 0 から 52 まで増える曲線です。MSTL はこの成分を .trend に取り出すので、トレンド図が上に凸の曲線になっていれば正しく分かれています。
12-17 行目日本語の列名で DataFrame にする
"日付" と "観測値" は後のプロットの軸ラベルにそのまま出ます。matplotlib に日本語フォントを設定していないと豆腐（□）になるので、japanize_matplotlib を import しておきます。

	日付	観測値
0	2018-01-01	20.000000
1	2018-01-02	24.618006
2	2018-01-03	26.583476
3	2018-01-04	26.587164
4	2018-01-05	28.330645
5	2018-01-06	32.415653
6	2018-01-07	35.578666
7	2018-01-08	35.663289
8	2018-01-09	34.892380
9	2018-01-10	36.617664

折れ線で全体の形を見ます。大きなうねり（周期 168）の上に細かいギザギザ（周期 4.8 と 14）が乗り、全体が右上がりになっています。

1
2
3
4
plt.figure(figsize=(15, 5))
sns.lineplot(x=df["日付"], y=df["観測値"])
plt.grid(axis="x")
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1 行目横幅 15 で 1000 点を潰さない
figsize=(15, 5) にすると 1000 点が 1 ピクセルに 1 点弱で並び、周期 4.8 の細かい波と周期 168 の大きなうねりが同時に見えます。既定サイズだと細かい波が黒い帯になり、2 種類の周期があることが分かりません。
2 行目横軸は datetime なので月ごとの目盛になる
df["日付"] が datetime64 型なので、matplotlib は目盛を月や年単位で自動調整します。周期 168 日のうねりが「およそ半年に 1 回」であることを目盛から読み取れます。
3 行目縦のグリッドで周期を数える
grid(axis="x") の縦線を目印にすると、大きなうねりの山から山までが何か月かを数えられます。実データで periods を決めるときは、この図と ACF から周期の候補を拾います。

サンプルデータを作成 figure

トレンド・季節/周期性・残差に分解する #

periods に取り出したい周期をタプルで渡し、.fit() を呼ぶだけで分解が終わります。周期は短いものから並べます。

1
2
periods = (24, 24 * 7)
mstl = MSTL(df["観測値"], periods=periods).fit()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1 行目periods は短い周期から順に並べる
(24, 168) の 2 つを指定しています。MSTL は内部で短い周期から順に STL を当てるので、順番は昇順にします。ここに 14 を足して (14, 24, 168) にすればノコギリ波も季節成分として取り出せ、残差がさらに 0 に近づきます。周期を 1 つだけ渡すと通常の STL と同じ結果になります。
2 行目fit() の戻り値に 4 つの成分が入る
.fit() が返すオブジェクトは observed・trend・seasonal・resid を属性として持ちます。seasonal は periods の数だけ列を持つ DataFrame になり、列名は seasonal_24 と seasonal_168 です。fit() を忘れると分解前の設定オブジェクトのままで、属性にアクセスすると AttributeError になります。
2 行目windows と iterate で平滑化を調整できる
省略した引数のうち windows は周期ごとの LOESS の窓幅、iterate は STL を繰り返す回数（既定 2）です。季節成分が年々変わるデータでは windows を小さくして追従させ、ノイズが多いデータでは大きくしてなめらかにします。

分解結果のトレンド（.trend）、季節/周期性（.seasonal）、残差（.resid）をそれぞれプロットします。今回は周期の異なる三角関数を 2 本足しているので .seasonal には 2 つの列が含まれています。残差のプロットにところどころ山があるものの、ほとんどの領域で残差が 0 に近く、きれいに分解できていることが確認できます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
_, axes = plt.subplots(figsize=(12, 8), ncols=1, nrows=4, sharex=True)

axes[0].set_title("観測値")
axes[0].plot(mstl.observed)
axes[0].grid()

axes[1].set_title("トレンド")
axes[1].plot(mstl.trend)
axes[1].grid()

axes[2].set_title("季節性")
axes[2].plot(mstl.seasonal.iloc[:, 0], label="周期的な要素１")
axes[2].plot(mstl.seasonal.iloc[:, 1], label="周期的な要素２")
axes[2].legend()
axes[2].grid()

axes[3].set_title("その他の要因・残差")
axes[3].plot(mstl.resid)
axes[3].grid()

plt.tight_layout()
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1 行目4 段を縦に並べ、横軸を共有する
nrows=4, sharex=True で観測値・トレンド・季節性・残差を同じ横軸に揃えます。横軸を共有すると「観測値のこの山は季節成分 2 の山だ」と上下で対応づけて読めます。mstl.plot() でも同じ図が出ますが、季節成分が別々の段に分かれるので、ここでは 2 列を 1 段に重ねるために手で描いています。
8 行目トレンドは log + sqrt の曲線を再現する
.trend はデータ生成で足した np.log(i**3 + 1) + np.sqrt(i) の上に凸の曲線に、i % 14 の平均（約 6.5）が乗った形になります。周期 14 を periods に入れていないため、ノコギリ波の平均部分がトレンドに吸収されています。
12-13 行目seasonal の 2 列を 1 段に重ねる
iloc[:, 0] が周期 24 の成分（中に周期 4.8 の sin が 5 周期入った形、振幅 2 前後）、iloc[:, 1] が周期 168 の成分（cos のうねり、振幅 10 前後）です。2 本を重ねると、短い周期が小さな振幅、長い周期が大きな振幅で分かれているのが一目で分かります。
18 行目残差に残る山は periods に入れなかった周期 14
.resid はほとんど 0 付近ですが、ところどころ周期的な山が見えます。これは i % 14 のノコギリ波のうちトレンドに吸収されなかった揺れで、periods=(14, 24, 168) にすると消えます。実データでは、残差に周期的な山が残っていたら「拾い損ねた周期がある」合図です。

残差のプロットにところどころ山があるものの、ほとんどの領域で残差が0に近い（＝きれいに分解できている）ことが確認できます… figure

使いどころと落とし穴 #

MSTL は電力・来店数・Web アクセスのように「1 日」と「1 週間」の周期が重なるデータで威力を発揮します。周期の候補は ACF や FFT のピークから拾い、periods を昇順に並べて渡します。周期を入れ忘れると、その成分はこのページの i % 14 のようにトレンドや残差に紛れるので、分解後は必ず .resid に周期的な山が残っていないかを確かめてください。1 周期しかないデータなら STL で十分で、statsmodels 0.13 以前の環境では MSTL 自体が import できない点にも注意が必要です。