概要

少し前まではState Space Models(SSMs)というとRの独壇場だったけど、今はPythonでも同様のモデルを構築することができる。ここではガラパゴス語でほとんど情報がないStatsModelsを使ったSSMsのOn-line forecastを試している。

実施期間: 2024年8月
環境：Ubuntu22.04LTS
Python: condaのPython3.10

書籍は何冊か購入したが中でもこれは参考になった。

1. Dataset

いつもの気象庁のサイトから東京における月毎の平均気温、気圧、湿度データを拝借し気温を予測することとする。DL後のcsvファイルから不要な列や列名はエディタで編集しておく。
列名はそれぞれ['temperature', 'pressure', 'humidity']とする。

    df = pd.read_csv('tokyo_temp.csv')
    df['date'] = pd.to_datetime(df['date'])
    df = df.reset_index(drop=True).set_index('date')
    df['pressure_rolling'] = df['pressure'].rolling(5).var().shift(1)
    df['humidity_rolling'] = df['humidity'].rolling(5).var().shift(1)
    df = df.dropna()
    df = df.asfreq('D')

SSMsでは変数をendog(endogenous variable:内生変数)とexog(exogenous variable:外生変数)に区別する。共変変数かどうかの違いはあるが、一般のML的には前者が目的変数で後者が説明変数に相当するイメージ。本例では気温の変化で気圧と湿度は変化しないと仮定し、'temperature'をendog、'pressure'と'humidity'をexogとする。

2. 作戦

その都度、onlineでforecastさせるSSMなのでLSMTのcellように過去の内部状態をstatefullに、つまり過去の時系列データの傾向を引き継ぐように実装する必要がある。もちろんMustではないがこれをしないのであればSSMsを使う意味がない。
簡単のためendodgの'temperature'だけで説明する（時間を追うごとに気温は1℃ずつ上昇する前提）。

表1. 時間を追って入力するendogの変遷

time	fit()	1st forecast	2nd	3rd
1	10
2	11
3	12
4	13
5	14
6	15	15
7		16	16
8		17	17	17
9		?	18	18
10			?	19
11				?

上表の例では、6time step(the longer, the better)の初期時系列データで最初のfit()を実行して状態方程式の変数を取得する。（本例ではこのあと初期化する必要はないのでfit()は行わない。）
online(運用時)においてデータが追加されるたびにその1time step先の?をforecastする流れ。
状態を引き継ぐため必ずデータはどこかでオーバラップさせている。使う時系列データにマルコフ性があれば1点でもオーバラップしていれば十分だと思う。

通常の時系列MLモデルの場合、fit後のモデルを予測のたびに使いまわすが本例のSSMでは各予測ごとに改めてモデルを作成する。この時、初期モデルのパラメータmodel_fit.paramsを使いまわす。また各回のモデルの状態方程式の変数には前回モデルの一部変数を流用する。状態空間を時間的に連続にすることがポイントなので、都度作成するモデルに入れるデータは一部オーバラップさせているというわけ。
秘伝のたれを継ぎ足しているような感じ

3. 状態の継承

statsmodelsのcontributorの方のnotebookを読むと次のような流れとなっている。

初期モデルをfit()したときに返るインスタンス(model_fit)内にはparams(初期パラメータ)と各time stepにおける状態方程式用のpredicted_state(変数), predicted_state_cov(変数の共分散)が入っている。本例では6time step分を入力したので 1(初期変数)+6個の変数が入っている。

6time step分のデータを時間を追いながら各時間の変数と変数の共分散を推論する。何time step分を入れようが各時間のデータ(観測値)で推論されるためそれらは変わらない。つまり各時刻における予測精度に影響しない。MLの場合、Traininig時(=prediction時)のデータ長は精度に大きく影響するが、SSMではモデルに時刻ごとの変数が保持されるためそもそも考え方が異なる。

ちなみにorder=(1, 0, 3)のときの初期パラメータはこんな感じ

rint(model_fit_params)
pressure_rolling   -0.008339
humidity_rolling    0.001281
ar.L1               0.999122
ma.L1              -0.198631
ma.L2              -0.281519
ma.L3              -0.141695
sigma2              4.015282
dtype: float64

下表では簡単のためpredicted_stateとしてスカラ値を書いたが、実際にはAR, I, MAのorderごとの行列が入る。predicted_state_csvも同様の表となるが共分散行列のため次元が１つ多い。実行中にブレイクで止めて中を見てみるとイメージはわかるはず。

表2. 時刻を追って変化する状態方程式用の変数(fit後)

time	temperature	fit()
0	NaN	0.0
1	10	0.1
2	11	0.2
3	12	0.3
4	13	0.4
5	14	0.5
6	15	0.6

この場合、time step6の'15℃'のときの変数が0.6ということ。
表1の1st forecastでendog[15, 16, 17]を入力するが、time step6の'15℃'はオーバラップさせている。これは1st forecastで改めてモデルインスタンスを作成するとき、その初期値として1つ前のモデルの変数を継承したかったためである。

表3. 時刻を追って変化する状態方程式用の変数(1st forecast用モデル用)

time	temperature	fit()	1st forecast	note
0	NaN	0.0
1	10	0.1
2	11	0.2
3	12	0.3
4	13	0.4
5	14	0.5	-> 0.5	新モデル.initialize_known()で初期設定
6	15	0.6	0.6	新モデル.filter()で得られる
7	16		0.7	新モデル.filter()で得られる
8	17		0.8	新モデル.filter()で得られる

これがmodel_fit.predicted_state[:, -2]とmodel_fit.predicted_state_cov[:, ;, -2]の"-2"の理由である。
2nd forecastでendog[16, 17, 18]を入力したときは下表のように1st forecastで得た変数の一部を変数の初期値に継承する。

表4. 時刻を追って変化する状態方程式用の変数(2nd forecast用モデル用)

time	temperature	fit()	1st	2nd	note
0	NaN	0.0
1	10	0.1
2	11	0.2
3	12	0.3
4	13	0.4
5	14	0.5	0.5
6	15	0.6	0.6	-> 0.6	新モデル.initialize_known()で初期設定
7	16		0.7	0.7	新モデル.filter()で得られる
8	17		0.8	0.8	新モデル.filter()で得られる
9	18			0.9	新モデル.filter()で得られる

これがmodel_fit.predicted_state[:, 1]とmodel_fit.predicted_state_cov[:, ;, 1]の"1"の理由である。3rd forecast以降はこの繰り返し。実際にコードで確認するとわかるが、1st[16, 17]と2nd[16, 17]の変数（上例では[0.7, 0.8]）は全く同じになり、初期変数とendogで状態空間を継承することがわかる。

4. Forecast

前述の通り各forecastごとにモデルインスタンスを作成するがfit()はしない。空のモデルに下記を与え、filter()を実行しそのKalman filterで状態方程式を更新する。

・初期パラメータ
・endog
・exog（必要なら）
・前回モデルから流用する初期変数

filter()の結果はFilterResutsオブジェクトとして返される。forecastすることが目的なのでsmoothingは使わない。filteringとsmoothingの違いは前述の書籍に詳しい解説あり、前者がout-of-sampleのforecast用で後者がin-sampleの解析用。

FilterResuts(filtered_results)のattributeのうち、下記は次回のモデルへ継承する変数の参照に使用する。
filtered_results.predicted_state
filtered_results.predicted_state_cov
それと肝心のforecast用に下記の２つが実装されており前者は指定stepsのforecast値、後者は加えて指定信頼区間を出力してくれる。
filtered_results.forecast()
または
filtered_results.get_forecast()
本例では1 time stepずつずらしながらforecastするのでsteps=1にしている。

公式Docではモデルをfit()した結果のSARIMAXResultsオブジェクトでしかforecast系の使用はできないように読めてしまうが、本例のように都度filter()した結果のFilterResultsオブジェクトでも状態を持っているためforecastの実行は可能。
最早、fit() nearly equal filter()な感じ。

5. orderの特定

コレログラムを描いて自分で判断してもよいが、ACF,PACFチャートの読み方に相応の知識が必要なためpmdarima.auto_arimaに最適なorderを探してもらう。しかし複雑な処理が自動化されているのでそれなりに引数は多い。オフィシャルのTips to using auto_arimaが参考になる。
ここでは12ヶ月周期の月次気象情報を使用するため下記のように関数を用意した。

def get_orders(endog_df, exdog_df, m=1, with_exog=False):
    if with_exog:
        auto_arima = pm.auto_arima(
                            y=endog_df, 
                            X=exdog_df, 
                            max_order = 5,      # p+q+P+Q
                            m = m,
                            seasonal=True, 
                            stepwise = True,
                            approximation = False,
                            information_criterion='aic',
                            error_action="ignore",
                            suppress_warnings=True)
    else:
        auto_arima = pm.auto_arima(
                            y=endog_df, 
                            max_order = 5,      # p+q+P+Q
                            m = m,
                            seasonal=True, 
                            stepwise = True,
                            approximation = False,
                            information_criterion='aic',
                            error_action="ignore",
                            suppress_warnings=True)
    # 実行
    return auto_arima

解析後のauto_arimaオブジェクトのauto_arima.orderとauto_arima.seasonal_orderを使用する。周期をm=12とすると計算コストがかさんだのでstepwise=Trueでステップワイズ探索アルゴリズムによる計算負荷低減をしてもらっている。もしデータにseasonalityがなければseasonal=Falseとし、SARIMAXへはseasonal_order=(0,0,0,0)を渡せばよい。
pmdarimaはstatsmodelsをラップしているので内生、外生変数の考え方は上述と変わらない。

6. exogenous variable

状態方程式にexogを含めて更新するかどうかは、SARIMAX()モデルの引数であるtime_varying_regressionかmle_regressionで予め定義しておく。外生変数を状態方程式に組み込むオプションなので、両方Falseにするのはもったいない。

引数	解説 (DeepL)
time_varying_regression	説明変数 exog が提供されている場合に使用され、外生回帰変数の係数が時間と共に変化することを許すかどうかを選択する。デフォルトは False
mle_regression	外生変数の回帰係数を最尤推定またはKalman filter(つまり再帰的最小二乗法)で推定するかどうか。 time_varying_regression が True の場合、これは False に設定しなければならない。デフォルトはTrue

内生変数を気温、外生変数を気圧とした時、気圧が気温に与える影響が時間変化すると考えるのであればtime_varying_regression=True、一年を通して影響は一定であると考えるのであればmle_regression=Trueに設定したほうがよい。

コード中のget_forecast()におけるexogは外生変数であるが予測する次time step("t + pred_tw -1: t + pred_tw")の気圧と湿度を入力している点に注意が必要である。簡単のためこのような現実的に不可能なコードにしているが、実際には別モデルで次time stepの気圧と湿度を予測し、それらをfitとforecastに使用しなければならない。この外生変数の精度次第で肝心の内生変数の予測が悪化しかねず、その場合は外生変数を使用しないことがcontributerの方も推奨されている。
次time stepの内生変数をforecastする時、その予測したい時点の外生変数を与えないといけない点がSARIMAXの難点である。

7. code

1time step前の気温をbase lineとしforecast結果をRMSEで比較している。

import numpy as np
import pandas as pd
from tqdm import tqdm
import statsmodels.api as sm
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
import pmdarima as pm
from sklearn.metrics import root_mean_squared_error

def get_orders(endog_df, exdog_df, m=1, with_exog=False):
    if with_exog:
        auto_arima = pm.auto_arima(
                            y=endog_df, 
                            X=exdog_df, 
                            max_order = 5,      # p+q+P+Q
                            m = m,
                            seasonal=True, 
                            stepwise = True,
                            approximation = False,
                            information_criterion='aic',
                            error_action="ignore",
                            suppress_warnings=True)
    else:
        auto_arima = pm.auto_arima(
                            y=endog_df, 
                            max_order = 5,      # p+q+P+Q
                            m = m,
                            seasonal=True, 
                            stepwise = True,
                            approximation = False,
                            information_criterion='aic',
                            error_action="ignore",
                            suppress_warnings=True)
    # 実行
    return auto_arima


if __name__ == '__main__':
    df = pd.read_csv('tokyo_temp_monthly.csv')
    df['date'] = pd.to_datetime(df['date'], format='%Y/%m')
    df['date'] = pd.to_datetime(df['date']).dt.to_period('M').dt.to_timestamp()
    df = df.reset_index(drop=True).set_index('date')
    df = df.drop(columns='sunshine_hours')
    # df['pressure_rolling'] = df['pressure'].rolling(5).var()    #.shift(1)
    # df['humidity_rolling'] = df['humidity'].rolling(5).var()    #.shift(1)
    # df = df.dropna()
    df = df.asfreq('MS')

    with_exog = True
    initial_tw = 12 * 30
    m = 12       # the number of periods in each season

    # infer orders
    endog_data = df.iloc[:initial_tw]['temperature']
    # exog_data = df.iloc[:initial_tw][['pressure_rolling', 'humidity_rolling']]
    exog_data = df.iloc[:initial_tw][['pressure', 'humidity']]
    auto_arima = get_orders(endog_data, exog_data, m, with_exog)
    print(auto_arima.summary())

    # kwargs = {'order': auto_arima.order, 
    ##           'trend':'c', 
    #           'seasonal_order': (0, 0, 0, 0), 
    #           'initialization': 'approximate_diffuse', 
    #           'enforce_stationarity': False, 
    #           'enforce_invertibility': False}
    kwargs = {'order': auto_arima.order, 
    #          'trend':'c', 
              'seasonal_order': auto_arima.seasonal_order, 
              'initialization': 'approximate_diffuse', 
              'enforce_stationarity': False, 
              'enforce_invertibility': False}

    if with_exog:
        model = SARIMAX(endog_data, exog=exog_data, **kwargs)
    else:
        model = SARIMAX(endog_data, **kwargs)
    model_fit = model.fit(disp=False)
    model_fit_params = model_fit.params

    # 初期状態方程式行列 最後の1ts分が逐次更新の初回の最初の1tsとダブっているから
    predicted_state_estimates = model_fit.predicted_state[:, -2]
    predicted_state_covariances = model_fit.predicted_state_cov[:, :, -2]

    pred_tw = 6
    is_1st2 = True
    for t in tqdm(range(initial_tw, len(df)-pred_tw-1)):
        endog_data = df.iloc[t - 1: t + pred_tw - 1]['temperature']
        # exog_data = df.iloc[t - 1: t + pred_tw - 1][['pressure_rolling', 'humidity_rolling']]
        exog_data = df.iloc[t - 1: t + pred_tw - 1][['pressure', 'humidity']]
        if with_exog:
            model = SARIMAX(endog_data, exog=exog_data, **kwargs)
        else:
            model = SARIMAX(endog_data, **kwargs)
        model.initialize_known(predicted_state_estimates, predicted_state_covariances)

        filtered_results = model.filter(model_fit_params)

        # forecast結果取得
        # prediction = filtered_results.get_forecast(steps=1, exog=exog_data.iloc[-1])
        if with_exog:
            # prediction = filtered_results.get_forecast(steps=1, exog=df.iloc[t + pred_tw -1: t + pred_tw][['pressure_rolling', 'humidity_rolling']])
            prediction = filtered_results.get_forecast(steps=1, exog=df.iloc[t + pred_tw -1: t + pred_tw][['pressure', 'humidity']])
        else:
            prediction = filtered_results.get_forecast(steps=1)
        if is_1st2:
            result_df = pd.concat([df.iloc[t+pred_tw-1: t+pred_tw]['temperature'], 
                                    prediction.predicted_mean, 
                                    prediction.conf_int(alpha=0.10)], axis=1)
            result_df = result_df.rename(columns={0: 'pred_mean'})
            is_1st2 = False
        else:
            wk_df = pd.concat([df.iloc[t+pred_tw-1: t+pred_tw]['temperature'], 
                                prediction.predicted_mean, 
                                prediction.conf_int(alpha=0.10)], axis=1)
            wk_df = wk_df.rename(columns={0: 'pred_mean'})
            result_df = pd.concat([result_df, wk_df])

        # 次mmodel用に変数を確保
        predicted_state_estimates = filtered_results.predicted_state[:, 1]
        predicted_state_covariances = filtered_results.predicted_state_cov[:, :, 1]

    result_df['base_line'] = result_df['temperature'].shift(1)
    result_df = result_df.dropna()
    result_df.to_csv('result_df1.csv')

    wk_result_lst = [[root_mean_squared_error(result_df['temperature'], result_df['base_line']), 
                        root_mean_squared_error(result_df['temperature'], result_df['pred_mean'])]]
    
    print(f'base_acc: {round(float(wk_result_lst[0][0]),5)}, acc: {round(float(wk_result_lst[0][1]),5)}')

気温にトレンドはないのでtrend='c'はコメントアウト。温暖化の影響がトレンドとして認められるのならtrend='ct'などとすべきか。
外生変数は湿度と気圧をそのまま使用したが、各移動分散で試したときのコードをコメントアウトしている。

以上