statsmodelsの季節分解で実装されているアルゴリズム

前回に続いて、statsmodelsの季節分解の話です。
参考(前回の記事): statsmodelsを利用した時系列データの季節分解のやり方

前回の記事は使い方でしたが、今回はソースコードを参照しながらどのような計算方法で季節分解が実装されているのかを見ていきます。
ちなみに、僕の環境のバージョンは以下の通りです。将来のバージョンでは仕様が変わる可能性もあるのでご注意ください。

$ pip freeze | grep statsmodels
> statsmodels==0.13.2

ドキュメントにソースコードが載ってるページがあるので、そこを参照します。
ソース: statsmodels.tsa.seasonal — statsmodels

ソースの先頭に以下のようにコメントが書かれている通り、移動平均を使って実装されています。

"""
Seasonal Decomposition by Moving Averages
"""

利用するデータですが、前回の記事と同じ、このブログのPVです。

# データ件数
print(len(df))
# 140

# 2週間分のデータ
print(df.tail(14))
"""
              pv
date
2022-10-24  2022
2022-10-25  2140
2022-10-26  2150
2022-10-27  1983
2022-10-28  1783
2022-10-29   847
2022-10-30   793
2022-10-31  1991
2022-11-01  2104
2022-11-02  1939
2022-11-03  1022
2022-11-04  1788
2022-11-05   830
2022-11-06   910
"""

それでは順番にやっていきましょう。モデルは加法モデルの方を取り上げます。乗法モデルもトレンド成分を抽出するところまでは一致していますし、その後の処理から異なりますが加法と乗法の違いを考えれば普通に理解できると思います。

そして、ここが重要なのですが、まず周期が奇数の場合を見ていきます。今回は1週間である7日です。なんでこれが重要かというと、偶数の場合は少し挙動が特殊だからです。

季節分解では、元の時系列データからトレンド成分と季節成分を取り出します(残りが残差)が、処理の順番もこの順番になっていて、最初にトレンド成分、次に季節成分を取り出すようになっています。これ順番逆にすると結果変わりますし、それぞれメリットデメリットあるのですがトレンド成分を先にFIXすることを選ばれているようです。

では、トレンド成分の抽出を見ていきましょう。これ、”by Moving Averages”のコメント通り、移動平均です。two_sided=True (デフォルト)の設定場の場合、その対象日を挟むようにして、前後の期間から取った移動平均をその日のトレンド成分の値とします。

7日だったら、3日前,2日前,1日前,当日,1日後,2日後,3日後 の平均を使います。

pandasで計算するなら、rolling()で移動平均取って、shift()でずらすと同じ値が得られます。
モデルで取り出したトレンド成分と、pandasで自分で計算した値見てみましょう。

# モデルで計算
from statsmodels import api as sm


decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=7,
)
print(decompose_result.trend)
"""
date
2022-06-20            NaN
2022-06-21            NaN
2022-06-22            NaN
2022-06-23    1409.857143
2022-06-24    1408.714286
                 ...     
2022-11-02    1495.285714
2022-11-03    1512.000000
2022-11-04            NaN
2022-11-05            NaN
2022-11-06            NaN
Name: trend, Length: 140, dtype: float64
"""

# pandasで計算。7日移動平均をとって、3日シフトする
print(df["pv"].rolling(7).mean().shift(-3))
"""
date
2022-06-20            NaN
2022-06-21            NaN
2022-06-22            NaN
2022-06-23    1409.857143
2022-06-24    1408.714286
                 ...     
2022-11-02    1495.285714
2022-11-03    1512.000000
2022-11-04            NaN
2022-11-05            NaN
2022-11-06            NaN
Name: pv, Length: 140, dtype: float64
"""

データの大部分…て略されてますが、これらは一致します。(極小値の誤差はあり得ますが)

two_sided = False とすると、前後のデータではなくその日含めた過去のデータだけ使われるようになるので、6日前〜当日のデータでの移動平均になります。

比較できるように、最初の2個の値が見れるようにprintしました。先頭のNaNの数が6個になってますね。表示てませんが、代わりに末尾のNaNは無くなります。

decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=7,
    two_sided=False
)
print(decompose_result.trend.iloc[: 8])
"""
date
2022-06-20            NaN
2022-06-21            NaN
2022-06-22            NaN
2022-06-23            NaN
2022-06-24            NaN
2022-06-25            NaN
2022-06-26    1409.857143
2022-06-27    1408.714286
Name: trend, dtype: float64
"""

さて、トレンド成分の計算方法はわかったので、季節成分の話に戻ります。モデルも two_sided=Trueの最初の例の方に話戻してこちらで進めます。(この記事真似して再現する人は上のtwo_sided=Falseのコードブロックをスキップして実行してください)

季節成分は、トレンド成分を取り除き終わったデータから計算します。(最初と最後の数件のデータはNaNになっているのでこれらは使いません。)
説明が難しいのですが、今回みたいに曜日ごとの周期であれば、月曜の平均、火曜の平均、水曜の平均と、周期分の平均値を取り出し、さらにこうして計算した平均値たちからその平均を引いて標準化します。

言葉で書くとわかりにくいのでコードでやってみます。

import numpy as np


# 先述のトレンド成分の計算
df["pv_trend"] = df["pv"].rolling(7).mean().shift(-3)

# トレンド成分を取り除く
df["pv_detrended"] = df["pv"] - df["pv_trend"]

# 曜日ごとの平均を取る
period_averages = np.array([df["pv_detrended"][i::7].mean() for i in range(7)])
print(period_averages)
"""
[ 243.2556391   397.69172932  341.57894737  257.54285714  139.98496241
 -720.2481203  -675.17293233]
"""

# 曜日ごとの平均から、さらにそれらの平均を引く
period_averages -= period_averages.mean()
print(period_averages)
"""
[ 245.450913    399.88700322  343.77422127  259.73813104  142.18023631
 -718.0528464  -672.97765843]
"""

# 比較用。statsmodelsが算出した季節成分
decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=7,
)
print(decompose_result.seasonal[: 7])
"""
date
2022-06-20    245.450913
2022-06-21    399.887003
2022-06-22    343.774221
2022-06-23    259.738131
2022-06-24    142.180236
2022-06-25   -718.052846
2022-06-26   -672.977658
Name: seasonal, dtype: float64
"""

データの形式が違いますが、モデルの結果と値が一致してるのがわかりますね。

以上で、トレンド成分と季節成分が計算できました。あと、statsmodelsは残差を計算できますがこれは単にトレンド成分と周期成分を元のデータから取り除いてるだけです。 df[“pv_detrended”] はトレンド成分除去済みなのでここから季節成分も引いてみます。

# np.tile で同じ値を繰り返す配列を作って引く
print(df["pv_detrended"] - np.tile(period_averages, 20))
"""
date
2022-06-20           NaN
2022-06-21           NaN
2022-06-22           NaN
2022-06-23     74.404726
2022-06-24     36.105478
                 ...    
2022-11-02     99.940064
2022-11-03   -749.738131
2022-11-04           NaN
2022-11-05           NaN
2022-11-06           NaN
Name: pv_detrended, Length: 140, dtype: float64
"""

# モデルの残差
print(decompose_result.resid)
"""
date
2022-06-20           NaN
2022-06-21           NaN
2022-06-22           NaN
2022-06-23     74.404726
2022-06-24     36.105478
                 ...    
2022-11-02     99.940064
2022-11-03   -749.738131
2022-11-04           NaN
2022-11-05           NaN
2022-11-06           NaN
Name: resid, Length: 140, dtype: float64
"""

以上が、statsmodelsにおける seasonal_decompose の実装の説明になります。

さて、ここからはちょっと補足で、周期が偶数の場合の挙動になります。月次データだと12を使うことが多いので周期が偶数というのはよくあるケースです。
説明をコンパクトにするために、周期(period=)4を例に取り上げます。

周期7の時、3日前から3日後までのデータの平均でトレンド成分を取り出してましたが、こういう期間の取り方ができるのは周期が奇数だったからで、偶数だとこうはいきません。

では、どうするのかというと、例えば周期が4日だったら、前日から翌日までの3日分のデータはそのまま使い、2日前と2日後のデータをそれぞれ0.5倍したものを計算に入れます。

# 周期4の場合のトレンド成分
decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=4,
)
print(decompose_result.trend.iloc[: 5])
"""
date
2022-06-20         NaN
2022-06-21         NaN
2022-06-22    1719.750
2022-06-23    1567.250
2022-06-24    1299.125
Name: trend, dtype: float64
"""

# 以下、 2022-06-22    1719.750 が算出されるまでの計算式
# 元のデータ
df["pv"].iloc[: 5]
"""
date
2022-06-20    1703
2022-06-21    1758
2022-06-22    1732
2022-06-23    1744
2022-06-24    1587
Name: pv, dtype: int64
"""

# 計算 2日前〜2日後の5日分のデータから算出。ただし、最初と最後の日はウェイトが0.5
(1703*0.5 + 1758 + 1732 + 1744 + 1587*0.5)/4
# 1719.75

rolling()とshift()では同じ値を得られなかったので戸惑ったこともありましたが、まぁ、仕掛けがわかってしまえば簡単ですね。その日を挟んだ4日分のデータを考慮する手段としても一定の妥当性があると思います。

で、ここからが僕的には納得がいってないのですが、two_sided=Falseで周期を偶数(今の例では4)にした場合の挙動は少し不思議です。その日を含む過去4日分のデータを使えるので、単純に4日間の平均を取ればいいのに、そういう実装になっておらず、two_sided=Trueの場合の結果を単純にshiftしたものを使っています。要するに4日前の0.5倍と3日前から1日前、そして当日の0.5倍のデータを使ってます。two_sided=Falseにしたコードが以下ですが、上のTrue(省略した場合の値)と値が一致しているのがわかると思います。

decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=4,
    two_sided=False,
)
print(decompose_result.trend.iloc[: 7])
"""
date
2022-06-20         NaN
2022-06-21         NaN
2022-06-22         NaN
2022-06-23         NaN
2022-06-24    1719.750
2022-06-25    1567.250
2022-06-26    1299.125
Name: trend, dtype: float64
"""

ここの挙動だけは将来のバージョンで修正されるのではと思っているのですが、
一旦今はこういう作りになっているということを頭に置いた上で気をつけて使うしかないようです。

statsmodelsを利用した時系列データの季節分解のやり方

要するに、seasonal_decomposeメソッドの使い方の紹介記事です。これもとっくに書いたと思っていたら書いてなかったのでまとめておきます。この記事では季節分解の概要の説明とライブラリの使い方を紹介します。そして、これの次の記事でstatsumodelsがどのような実装で季節分解を行っているのかを解説する予定です。
参考: statsmodels.tsa.seasonal.seasonal_decompose — statsmodels

現実の時系列データは、何かしらの季節性を持っていることが多くあります。季節って単語で言うと、春夏秋冬や、何月、といった粒度のものが想像されやすいですが、1週間の中で見ても曜日の傾向とか、1日の中でも時間帯別の違いなどがあります。

時系列データからこの季節に依存する部分を取り出し、季節成分と、トレンド成分、そして残差へと分解する手法が今回紹介する季節分解です。Wikipediaでは季節調整と書いてあります。また、基本成分など、別の用語を使ってる人もいるようです。(どれが一番メジャーなんだろう。)

定式化しておくと、元の時系列データ$Y_{t}$をトレンド成分$T_t$と季節成分$S_t$、そして残差$e_t$を用いて、
$$
Y_t = T_t + S_t + e_t
$$
と分解することを目指します。上記のは加法モデルと呼ばれる形で、和の代わりに積で分解する乗法モデル、
$$
Y_t = T_t * S_t * e_t
$$
もあります。

季節成分$S_t$は周期性を持っているので、その周期を$p$とすると、$S_{t}=S_{t+p}$を満たします。

具体的に例を見るのが早いと思うので、やっていきましょう。サンプルとして用意したデータはこのブログのpv数です。インデックスを日付けにしていますが、こうしておくとライブラリでplotしたときにx軸に日付が表示されるで便利です。ただ、通常の通し番号のindexでも動きます。

# データ件数
print(len(df))
# 140

# 2週間分のデータ
print(df.tail(14))
"""
              pv
date            
2022-10-24  2022
2022-10-25  2140
2022-10-26  2150
2022-10-27  1983
2022-10-28  1783
2022-10-29   847
2022-10-30   793
2022-10-31  1991
2022-11-01  2104
2022-11-02  1939
2022-11-03  1022
2022-11-04  1788
2022-11-05   830
2022-11-06   910
"""

最後に元のデータのグラフも出るので可視化しませんが、上の例見ても10/29,30や11/5,6など土日にpvが減っていて逆に平日多く、曜日ごと、つまり7日周期がありそうな想像がつきます。詳細省略しますが、自己相関等で評価してもはっきりとその傾向が出ます。

では、やってみましょう。まず、分解自体はライブラリにデータを渡して周期を指定するだけです。

from statsmodels import api as sm


decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=7,  # 周期を指定する
)

結果は以下のプロパティに格納されています。DecomposeResultというデータ型で、ドキュメントはこちらです。
参考: statsmodels.tsa.seasonal.DecomposeResult — statsmodels

順番に表示していきます。

# データの数。
decompose_result.nobs[0]
# 140

# 元のデータ
print(decompose_result.observed[: 10])
"""
date
2022-06-20    1703.0
2022-06-21    1758.0
2022-06-22    1732.0
2022-06-23    1744.0
2022-06-24    1587.0
2022-06-25     654.0
2022-06-26     691.0
2022-06-27    1695.0
2022-06-28    1740.0
2022-06-29    1655.0
Name: pv, dtype: float64
"""
# トレンド成分
print(decompose_result.trend[: 10])
"""
date
2022-06-20            NaN
2022-06-21            NaN
2022-06-22            NaN
2022-06-23    1409.857143
2022-06-24    1408.714286
2022-06-25    1406.142857
2022-06-26    1395.142857
2022-06-27    1370.714286
2022-06-28    1345.285714
2022-06-29    1347.142857
Name: trend, dtype: float64
"""

# 季節成分
print(decompose_result.seasonal[: 14])
"""
date
2022-06-20    245.450913
2022-06-21    399.887003
2022-06-22    343.774221
2022-06-23    259.738131
2022-06-24    142.180236
2022-06-25   -718.052846
2022-06-26   -672.977658
2022-06-27    245.450913
2022-06-28    399.887003
2022-06-29    343.774221
2022-06-30    259.738131
2022-07-01    142.180236
2022-07-02   -718.052846
2022-07-03   -672.977658
Name: seasonal, dtype: float64
"""
# 残差
print(decompose_result.resid[: 10])
"""
date
2022-06-20          NaN
2022-06-21          NaN
2022-06-22          NaN
2022-06-23    74.404726
2022-06-24    36.105478
2022-06-25   -34.090011
2022-06-26   -31.165199
2022-06-27    78.834801
2022-06-28    -5.172718
2022-06-29   -35.917078
Name: resid, dtype: float64
"""

トレンド成分が最初の3項NaNになっているのは、アルゴリズムの都合によるものです。その日を中心とする前後で合計7日(周期分)のデータで移動平均をとっており、要するに、過去の3日、当日、次の3日間、の合計7日分のデータがないと計算できないので最初の3日と、表示していませんが最後の3日間はNaNになっています。この辺の挙動は推定時の引数で調整できます。

次に季節成分は14日分printしましたが、最初の7日間の値が繰り返されて次の7日間でも表示されているのがわかると思います。ずっとこの繰り返しです。

残差は元のデータからトレンド成分と季節成分を引いたものになります。トレンド成分や季節成分に比べて値が小さくなっていて、今回のデータではトレンドと季節である程度分解が綺麗に行えたと考えられます。

さて、データが取れたのでこれを使えばmatplotlib等で可視化できるのですが、大変ありがたいことにこのDecomposeResultが可視化のメソッドを持っています。
少し不便なところは、そのplotメソッドがfigsizeとかaxとかの引数を受け取ってくれないので、微調整とかしにくいのですよね。
個人的な感想ですが、デフォルトでは少しグラフが小さいのでrcParamsを事前に変更してデフォルトのfigsizeを大きくし、それで可視化します。

import matplotlib.pyplot as plt
#  注: DecomposeResult.plot()を実行するだけなら matplotlibのimportは不要。
#       今回画像サイズを変えるためにインポートする

# figure.figsizeの元の設定を退避しておく
figsize_backup = plt.rcParams["figure.figsize"]

# 少し大きめの値を設定
plt.rcParams["figure.figsize"] = [10, 8] 

# 可視化
decompose_result.plot()
plt.show()

# 設定を元に戻す
plt.rcParams["figure.figsize"] = figsize_backup

これで出力される画像が以下です。

一番上が元のpvです。冒頭に書いた通り、生データそのまま見ても周期性が明らかですね。

トレンド成分を見ていくと、お盆の時期にアクセスが減っていますが、その後順調にアクセスが伸びていることがわかりますね。

残差で見ると大きくアクセスが減っているのはそれぞれ祝日に対応しています。
8月に1日だけ異常にアクセス伸びた日がありますがこれは謎です。

説明や例をいろいろ書いてきたので長くなりましたが、基本的には、seasonal_decompose()で分解して、plot()で可視化してそれで完成という超お手軽ライブラリなので、時系列データが手元にあったらとりあえず試す、くらいの温度感で使っていけると思います。

最後に、seasonal_decompose にオプション的な引数が複数あるので使いそうなものを説明してきます。

まず、model= は、 “additive”,”multiplicative” の一方をとり、加法的なモデルか乗法的なモデルを切り替えることができます。デフォルトは、”additive”です。

two_sided= は True/Falseの値をとり、デフォルトはTrueです。これはトレンド成分の抽出方法を指定するもので、Trueの時は、その日を挟むように前後の日付から抽出しますが、Falseの場合は、その日以前の値から算出されます。True or False で並行移動するイメージです。
次回の記事で詳細書こうと思いますが、周期が偶数か奇数かで微妙に異なる挙動をするので注意が必要です。

extrapolate_trend= はトレンド成分の最初と最後の欠損値を補完するための引数です。1以上の値を渡しておくと、その件数のデータを使って最小二乗法を使って線形回帰してトレンドを延長し、NaNをなくしてくれます。使う場合はその回帰が妥当かどうか慎重にみて使う必要がありそうです。

pprintでデータを整形して出力する

前回の記事がtextwrapだったので、文字列の見栄えを整えるつながりで今回はpprintを紹介しようと思います。
参考: pprint — データ出力の整然化 — Python 3.11.0b5 ドキュメント

自分はもっぱらdictやlistの表示に使うのですが、ドキュメントを見ると任意のデータ構造に使えるようなことが書いてありますね。

使い方は簡単で、printすると結果が少しちょっと見にくくなるようなdict等のデータを渡すだけです。値が少し長いデータを使ってprintと見比べてみます。

import pprint


# サンブルデータ作成
sample_data = {
    1: "1つ目のキーの値",
    2: "2つ目のキーの値",
    3: "3つ目のキーの値",
    4: "4つ目のキーの値",
    5: "5つ目のキーの値",
    6: "6つ目のキーの値",
    7: "7つ目のキーの値",
}
# printした結果
print(sample_data)
"""
{1: '1つ目のキーの値', 2: '2つ目のキーの値', 3: '3つ目のキーの値', 4: '4つ目のキーの値', 5: '5つ目のキーの値', 6: '6つ目のキーの値', 7: '7つ目のキーの値'}
"""

# pprintした結果
pprint.pprint(sample_data)
"""
{1: '1つ目のキーの値',
 2: '2つ目のキーの値',
 3: '3つ目のキーの値',
 4: '4つ目のキーの値',
 5: '5つ目のキーの値',
 6: '6つ目のキーの値',
 7: '7つ目のキーの値'}
"""

pformat というメソッドもあって、こちらを使うと整形したものをprintするのではなく文字列として返してくれます。一応試しますが、文字列で戻ってきてるのをみないといけないので一旦変数に格納して通常のprintで出力します。

p_str = pprint.pformat(sample_data)
# 結果確認
print(p_str)
"""
{1: '1つ目のキーの値',
 2: '2つ目のキーの値',
 3: '3つ目のキーの値',
 4: '4つ目のキーの値',
 5: '5つ目のキーの値',
 6: '6つ目のキーの値',
 7: '7つ目のキーの値'}
"""

さて、このpprintですが、基本的にはそのまま使えば十分なのですが細かい調整ができるようにいろんな引数を取れます。

例えば、 indent= (デフォルト1)でインデントの文字数を指定できますし、width= (デフォルト80)で、横幅の文字数の最大値を指定できます。ただしwidthはベストエフォートでの指定なので、データによっては収めることできずにはみ出します。ちょっとwidthの指定によって結果が変わる例も見ておきましょう。さっきのdictはwidthが大きくても改行されたので、もう少しコンパクトなのを使います。

sample_data_mini = {
    1: '1つ目のキーの値',
    2: '2つ目のキーの値',
    3: '3つ目のキーの値',
}

# 80文字に収まるので、width未指定だと1行で出力
pprint.pprint(sample_data_mini, indent=4)
"""
{1: '1つ目のキーの値', 2: '2つ目のキーの値', 3: '3つ目のキーの値'}
"""

# width   が小さいと収まるように改行される。
pprint.pprint(sample_data_mini, indent=4, width=30)
"""
{   1: '1つ目のキーの値',
    2: '2つ目のキーの値',
    3: '3つ目のキーの値'}
"""

また、データの構造によっては、辞書やリスト、タプルの入れ子になっていることもあると思います。そのようなとき、depthという引数を指定することにより何階層目まで出力するか指定することもできます。オーバーした分は省略記号… になります。ドキュメントのサンプルでちょっとやってみます。

tup = ('spam', ('eggs', ('lumberjack', ('knights',
       ('ni', ('dead', ('parrot', ('fresh fruit',))))))))

# depth未指定
pprint.pprint(tup, width=20)
"""
('spam',
 ('eggs',
  ('lumberjack',
   ('knights',
    ('ni',
     ('dead',
      ('parrot',
       ('fresh '
        'fruit',))))))))
"""

# depth=3を指定
pprint.pprint(tup, width=20, depth=3)
"""
('spam',
 ('eggs',
  ('lumberjack',
   (...))))
"""

何かAPIとか叩いて巨大なJSONが帰ってきたとき、中身を確認するのに先立って上の階層のkeyだけちょっと見たい、って場面で非常に便利です。

このほかにも、辞書の出力をするときにkeyでソートしてくれるsort_key= (デフォルトでTrue)や、widthの範囲に収まるならばできるだけ1行にまとめてくれるcompact= (デフォルトでTrue)などのオプションもあります。正直のこの二つはわざわざFalseを指定することはないかなと追うので結果は省略します。