概要

pythonの状態空間モデルのライブラリといえば、statsmodelsのUnobservedComponents*1だと思います。

www.statsmodels.org

このページを読んでいくと、特に周期変動(seasonal)のモデリングの方法がたくさんありそうだぞということが分かってきます。柔軟に様々な表現をライブラリ内で簡単にカスタマイズできることはありがたいことです。しかしそれぞれが一体どんな挙動になるのか、ドキュメントを読むだけではイメージしにくかったので、実験してみたというのが今回の記事の内容になります。

状態空間モデルってなんだっけ

時系列データの分析や予測に使用される手法です。下図のように、観測値が生成されるプロセスと、実際には観測されない状態が生成されるプロセスを合わせて記述します。

状態を生成する過程を記述する変数として、時間経過によって緩やかに上昇/下降することを表現するトレンドや、曜日変動などを表現する周期変動、外部要因のイレギュラーな影響などを組み込むことができます。外部要因に対する回帰の重みを時間軸方向に共通化しなければ、外部要因の観測値や状態に与える影響が時間経過によって、どのように変化するかを評価できるみたいです。例えば、CMを打った直後から終了までの売り上げに与える影響が、時間が経過すると変化することをモデル上で表現できます。

また、この記事では対象ではないですが、観測値の生成分布なども柔軟に設定することもできます。例えば、観測値が正の値のみの場合、対数正規分布から生成されるように設定することが可能です。もし変数の推定を綺麗に導出できない場合も、MCMCが良い感じに推定してくれます。

このように時系列の生成過程を、人間の直感や知見を柔軟に組み込むことで記述できることが状態空間モデルの利点です。

UnobservedComponentsの周期変動オプションについて

今回対象とする状態空間モデルのライブラリであるUnobservedComponentsには、大きく３つの(私の考える)周期変動を表現するオプションがあります。

(数式的なところは、書いても良いのですがただのコピペになるのでやめました。 https://www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.structural.UnobservedComponents.html を見てください。)

1. Seasonal

周期性の見られる期間-1個のパラメータを足し合わせることで、周期性を表現します。 stochastic_seasonal=Falseにすることで、周期変動を表現する式から誤差項がなくなります。

2.Frequency-domain Seasonal

1と大きく異なるのは、三角関数を使って表現する点でしょうか。期間periodの他に、harmonics（フーリエ級数展開の三角関数項の数とほぼ同等と考えています）を設定する必要があります。 Seasonalと違い、月の周期性と週の周期性というような、いくつかの周期性を同時に表現できるみたいです。

3.Cycle

推定する変数が一気に増加した印象があります。学習データが1.5年から12年までの長さであるとうまくいくみたいです。

周期変動を状態空間モデルで表現したい場合、これだけオプションがあるとどれを使ったら良いか分からなくなります(筆者の気持ち)。 Seasonalの方が最も推定コストは低そうですが、Cycleを実際に適用した場合、どのくらいの表現力を持つのかも気になるところです。次の章では、各オプションの挙動について簡単に検証してみます。

実験

使用データ

seabornライブラリから読み込める旅客機乗客数に関するデータを使用しました。データは月毎の乗客数を表現しています。

試しにトレンドや周期性があるかどうかをざっと確認してみると、（オープンデータゆえに綺麗に修正されているからなのか）とても綺麗に時系列を分解できました。トレンドがはっきりと右肩上がりになっているので、モデリングする際はトレンド項を考慮した方が良さそうです。

f:id:green-j:20220227141119p:plain — seasonal_decomposeの出力

上図を見ると1年ごとに周期性が確認できます。念の為周期性成分のみで自己相関を可視化すると、lag=12のところでコレログラムのピークを確認できるため、12ヶ月で１周期とみなして良さそうです。

f:id:green-j:20220227145543p:plain — 周期成分の自己相関

設定

状態空間モデルのモデリングにはstatsmodelsのUnobservedComponentsを使用しました。トレンド成分があるデータであることは確認済みなので、今回はlevel='local linear trend'と設定することにしました。

データは1949年1月から1960年12月まであるので、1959年12月までを学習期間として、1960年の1年を予測期間として検証に使用しました。検証のため、予測期間の誤差をMAE(絶対誤差)で測定しました。

比較手法は次の３つです。

model_1: seasonal=7に設定。seasonalの挙動を知るため。
model_2: freq_seasonal=[{'period':12, 'harmonics':6}]に設定。frequency-domain seasonalの挙動を知るため。
model_3: cycle=Trueに設定。cycleの挙動を知るため。

この３つの手法でそれぞれ学習と予測を実行します。例として、model_1の実行コードを以下に貼ります。

import numpy as np
import pandas as pd

import statsmodels as sm
from statsmodels.tsa.statespace.structural import UnobservedComponents
from statsmodels.datasets import elec_equip
from statsmodels.tsa.seasonal import STL, seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

%matplotlib inline
import matplotlib.pyplot as plt

# pre-processing
df = sns.load_dataset('flights')
df.index = pd.date_range(start='1949-01-01', freq='MS', periods=len(df))

df["mode"] = df["year"].map(lambda x: "train" if x < 1960 else "test")

# 1. seasonal
model_1 = UnobservedComponents(
    df.loc[df["mode"] == "train", "passengers"],
    level="local linear trend",
    seasonal=12
)
res_1 = model_1.fit()
print(res_1.summary())

# prediction
y_pred = res_1.predict(start='1949-01', end='1960-12')
df["pred_1"] = y_pred
mae_1 = np.mean(np.abs(df.loc[df["mode"]=="test", "passengers"] - df.loc[df["mode"]=="test", "pred_1"]))
print("MAE :", mae_1)

この後の結果も含めた、実際に実行したコードはこちらです。

github.com

結果

３つの学習と予測の結果を貼ります。まずはmodel_1(seasonal検証モデル)。

                            Unobserved Components Results                            
=====================================================================================
Dep. Variable:                    passengers   No. Observations:                  132
Model:                    local linear trend   Log Likelihood                -501.715
                   + stochastic seasonal(12)   AIC                           1011.430
Date:                       Sun, 27 Feb 2022   BIC                           1022.546
Time:                               07:39:05   HQIC                          1015.944
Sample:                           01-01-1949                                         
                                - 12-01-1959                                         
Covariance Type:                         opg                                         
====================================================================================
                       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------
sigma2.irregular  4.639e-08     20.916   2.22e-09      1.000     -40.995      40.995
sigma2.level      1.536e-07     51.519   2.98e-09      1.000    -100.976     100.976
sigma2.trend        56.5406     19.161      2.951      0.003      18.985      94.096
sigma2.seasonal     18.7901      4.137      4.542      0.000      10.682      26.898
===================================================================================
Ljung-Box (L1) (Q):                   7.80   Jarque-Bera (JB):                64.85
Prob(Q):                              0.01   Prob(JB):                         0.00
Heteroskedasticity (H):               4.01   Skew:                            -1.04
Prob(H) (two-sided):                  0.00   Kurtosis:                         5.96
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
MAE : 78.7839030141792

続いてmodel_2(frequency domain seasonal検証モデル)。

                                Unobserved Components Results                                
=============================================================================================
Dep. Variable:                            passengers   No. Observations:                  132
Model:                            local linear trend   Log Likelihood                -458.354
                   + stochastic freq_seasonal(12(6))   AIC                            924.707
Date:                               Sun, 27 Feb 2022   BIC                            935.790
Time:                                       07:31:00   HQIC                           929.207
Sample:                                   01-01-1949                                         
                                        - 12-01-1959                                         
Covariance Type:                                 opg                                         
==============================================================================================
                                 coef    std err          z      P>|z|      [0.025      0.975]
----------------------------------------------------------------------------------------------
sigma2.irregular            1.596e-08      9.677   1.65e-09      1.000     -18.967      18.967
sigma2.level                  15.4431      5.667      2.725      0.006       4.335      26.551
sigma2.trend                   0.0154      0.033      0.473      0.637      -0.049       0.080
sigma2.freq_seasonal_12(6)     1.0049      0.221      4.540      0.000       0.571       1.439
===================================================================================
Ljung-Box (L1) (Q):                  19.75   Jarque-Bera (JB):                 1.92
Prob(Q):                              0.00   Prob(JB):                         0.38
Heteroskedasticity (H):               2.28   Skew:                             0.17
Prob(H) (two-sided):                  0.01   Kurtosis:                         2.47
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
MAE : 14.404733712999388

最後にmodel_3(cycle検証モデル)。

                        Unobserved Components Results                         
==============================================================================
Dep. Variable:             passengers   No. Observations:                  132
Model:             local linear trend   Log Likelihood                -629.951
                              + cycle   AIC                           1267.902
Date:                Sun, 27 Feb 2022   BIC                           1279.310
Time:                        07:29:31   HQIC                          1272.537
Sample:                    01-01-1949                                         
                         - 12-01-1959                                         
Covariance Type:                  opg                                         
====================================================================================
                       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------
sigma2.irregular   8.27e-06     96.626   8.56e-08      1.000    -189.384     189.384
sigma2.level      1013.6807    187.769      5.399      0.000     645.660    1381.702
sigma2.trend      3.243e-07      6.161   5.26e-08      1.000     -12.075      12.075
frequency.cycle      0.0436      0.055      0.788      0.431      -0.065       0.152
===================================================================================
Ljung-Box (L1) (Q):                  10.74   Jarque-Bera (JB):                 5.67
Prob(Q):                              0.00   Prob(JB):                         0.06
Heteroskedasticity (H):               7.91   Skew:                            -0.37
Prob(H) (two-sided):                  0.00   Kurtosis:                         3.72
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
MAE : 70.1899506705816

学習期間のAICやBIC、予測期間のMAEのそれぞれにおいてmodel_2が最も良い値になりました。 model_1を単純に複雑化したのがmodel_2ですが、パラメータ数と学習期間における推論能力のトレードオフのような指標であるAIC *2やBICもmodel_2の方が良かったので、モデルのパラメータが増えたことによる過学習の可能性は低そうだと考えられます。 model_3は学習期間のLog likelihoodも最も悪い値になっているので、うまく学習できていないか、今回のタスクに合っていない可能性があります。

数字だけ見ても分からない部分があるので、実際の値とそれぞれの予測値を可視化しました。

f:id:green-j:20220227165912p:plain — 学習＆予測結果

AICやMAEなどが最も良かったmodel_2は、予測期間(点線よりも右側)でもうまく予測できているように見えます。今回のタスクにおいてはfrequency domain seasonalを活用して周期性をモデリングした方が良かったと考えられます。 model_2に比べると、model_1は周期変動を表現しようとはしているものの、予測期間においてトレンドをうまく合わせられていません。またmodel_3は周期変動のようなものをうまく表現できていません。 model_1, model_3が予測に失敗した原因を調査するため、成分ごとに可視化してみました。

f:id:green-j:20220227201128p:plain — model_1の成分ごとの値

f:id:green-j:20220227201208p:plain — model_2の成分ごとの値

f:id:green-j:20220227201234p:plain — model_3の成分ごとの値

model_2とmodel_1を比較すると、model_1の方は1960年の直前でトレンドの値がマイナスになっており、これが予測期間において実際の値よりも小さい値を予測した原因と考えられます。 model_2の方は周期性の表現力が増したおかげなのか、時間の経過とともに周期変動成分(= seasonal component)の振幅が大きくなっています。それゆえトレンドに上下の動きのようなノイズが少なくなるので、緩やかな上昇を表現できるようになります。結果としてmodel_2の予測精度は良くなっています。このことから、周期変動の表現力が不足していると他の成分で表現しようとするあまり、予測に失敗する可能性があることが分かりました。

model_3は予測値でも周期変動をキャッチできていませんでしたが、cycle成分もほぼ横ばいになっていることがわかります。代わりに状態(=Level component)に周期性のようなものが見えてしまっています。ドキュメントには

The cyclical component is intended to capture cyclical effects at time frames much longer than captured by the seasonal component.

と書かれていたのですが、今回はデータが1ヶ月おきに記録されているものだったので、cycleの良さを活かせていない可能性があります。データの記録頻度が日ごとのもので実験するとまた変わるかもしれません*3。

まとめ

状態空間モデルの周期変動を表現するオプションをそれぞれ試しました。超簡単なデータセットでの実験でしたが、周期変動の表現力が不足していると他の成分の表現に悪影響を与える恐れがあり、それゆえ予測に失敗する可能性があるということが分かりました。実際にどのオプションが良いかはデータによると思いますが、周期変動を表現したい場合はseasonalの設定とfrequency domain seasonalの設定の両方を試し、AICやクロスバリデーションでどちらを採用するか決めるのが良いかなと考えます。

オープンデータセットって見つけるの意外と大変ですね...。ダウンロードしてみると非数ばかりのデータや、周期変動の全くないデータなどが多かったです。