状態空間モデル - チャンドラセカール再帰

[1]:

%matplotlib inline

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

from pandas_datareader.data import DataReader

状態空間モデルに関連するほとんどの操作はカルマンフィルタの再帰を利用していますが、いくつかの特殊な場合では、しばしば「チャンドラセカール再帰」と呼ばれる別の方法を使用することができます。これらは、状態ベクトルの条件付きモーメントを反復的に計算するための代替手段を提供し、場合によってはカルマンフィルタの再帰よりも計算負荷が大幅に軽減されることがあります。詳細については、「Using the 『Chandrasekhar Recursions』 for Likelihood Evaluation of DSGE Models」（Herbst, 2015）を参照してください。ここでは基本的なアイデアのみを概説します。

状態空間モデルとカルマンフィルター

時間不変の状態空間モデルは次のように書けることを思い出してください：

\[\begin{split}\begin{aligned} y_t &= Z \alpha_t + \varepsilon_t, \qquad \varepsilon_t \sim N(0, H) \\ \alpha_{t+1} & = T \alpha_t + R \eta_t, \qquad \eta_t \sim N(0, Q) \\ \alpha_1 & \sim N(a_1, P_1) \end{aligned}\end{split}\]

ここで、\(y_t\) は \(p \times 1\) のベクトルで、\(\alpha_t\) は \(m \times 1\) のベクトルです。

カルマンフィルターの各イテレーション、例えば時刻 \(t\) においては、次の三つの部分に分けることができます：

  1. 初期化：条件付き状態分布 \(\alpha_t \mid y^{t-1} \sim N(a_t, P_t)\) を定義する \(a_t\) と \(P_t\) の指定。

  2. 更新：条件付き状態分布 \(\alpha_t \mid y^{t} \sim N(a_{t|t}, P_{t|t})\) を定義する \(a_{t|t}\) と \(P_{t|t}\) の計算。

  3. 予測：条件付き状態分布 \(\alpha_{t+1} \mid y^{t} \sim N(a_{t+1}, P_{t+1})\) を定義する \(a_{t+1}\) と \(P_{t+1}\) の計算。

もちろん、最初のイテレーションの後、予測部分は次のステップの初期化に必要な値を提供します。

予測ステップに注目すると、カルマンフィルタの再帰式は次のようになります：

\[\begin{split}\begin{aligned} a_{t+1} & = T a_{t|t} \\ P_{t+1} & = T P_{t|t} T' + R Q R' \\ \end{aligned}\end{split}\]

ここで、行列 \(T\) と \(P_{t|t}\) はそれぞれ \(m \times m\) のサイズであり、\(m\) は状態ベクトル \(\alpha\) のサイズを示します。場合によっては、状態ベクトルが非常に大きくなることがあり、その場合、\(P_{t+1}\) を計算するために必要な行列積が計算上の負担になることがあります。

例: 季節的自己回帰モデル

例として、AR(r)モデル（ここでは観測ベクトルの次元としてすでに\(p\)を使ったので、\(r\)を使用します）が状態空間形式に変換できることに注目してください：

\[\begin{split}\begin{aligned} y_t &= \alpha_t \\ \alpha_{t+1} & = T \alpha_t + R \eta_t, \qquad \eta_t \sim N(0, Q) \end{aligned}\end{split}\]

ここで、

\[\begin{split}\begin{aligned} T = \begin{bmatrix} \phi_1 & \phi_2 & \dots & \phi_r \\ 1 & 0 & & 0 \\ \vdots & \ddots & & \vdots \\ 0 & & 1 & 0 \\ \end{bmatrix} \qquad R = \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \qquad Q = \begin{bmatrix} \sigma^2 \end{bmatrix} \end{aligned}\end{split}\]

年次の季節性を示す日次データを使用したARモデルでは、\(r=365\)までのラグを取り入れたモデルを適合させることを検討するかもしれません。この場合、状態ベクトルは少なくとも\(m = 365\)となり、行列\(T\)と\(P_{t|t}\)はそれぞれ\(365^2 = 133225\)の要素を持つことになります。したがって、尤度関数の計算に費やす時間のほとんどは、予測ステップでの行列積に支配されることがあります。

状態空間モデルとチャンドラセカール再帰

チャンドラセカール再帰は、方程式 \(P_{t+1} = T P_{t|t} T' + R Q R'\) を異なる再帰式に置き換えます：

\[P_{t+1} = P_t + W_t M_t W_t'\]

ここで、\(W_t\) は \(m \times p\) の次元を持つ行列で、\(M_t\) は \(p \times p\) の次元を持つ行列です。\(p\) は観測ベクトル \(y_t\) の次元です。これらの行列自体にも再帰的な定式化があります。\(W_t\) と \(M_t\) の計算式については、Herbst (2015) を参照してください。

重要な注意点: カルマンフィルタとは異なり、チャンドラセカール再帰はすべての状態空間モデルに使用することはできません。特に、以下の制限があります（これらはカルマンフィルタに必要なものではありません）：

  • モデルは時間不変である必要があり、時間変動する切片は許可されています。

  • 状態ベクトルの定常初期化を使用しなければなりません（これにより、非定常成分のモデルは除外されます）。

  • 欠損データは許可されていません。

この式がより効率的な計算を意味する理由を理解するために、再び上記のSARIMAXのケースを考えましょう。この場合、\(p = 1\) であり、\(M_t\) はスカラーであるため、チャンドラセカール再帰を次のように書き直すことができます：

\[P_{t+1} = P_t + M_t \times W_t W_t'\]

掛け算される行列 \(W_t\) の次元は \(m \times 1\) であり、\(r=365\) の場合、各行列は \(365\) 要素しか持たず、\(365^2\) 要素を持つのではありません。これにより、予測ステップを完了するために必要な計算量が大幅に減少することを意味します。

収束

パフォーマンスの影響に関する議論を複雑にする要因の一つは、時間不変モデルにおいて、予測される状態共分散行列が定数行列に収束するというよく知られた事実です。これは、ある \(S\) が存在し、すべての \(t > S\) に対して \(P_t = P_{t+1}\) となることを意味します。収束が達成されると、予測ステップから \(P_{t+1}\) の方程式を完全に排除することができます。

AR(r) モデルのような単純な時系列モデルでは、収束は比較的早く達成されるため、Chandrasekhar 再帰を使用することによるパフォーマンス向上の効果が制限されることがあります。Herbst (2015) は、代わりに DSGE（動的確率的一般均衡）モデルに焦点を当てています。これらのモデルは、しばしば大きな状態ベクトルを持ち、収束を達成するために多くの期間が必要となることが多いです。この場合、パフォーマンスの向上は非常に大きくなることがあります。

実践的な例

実践的な例として、明確な季節的要素を持つ月次データを考えます。この場合、消費者物価指数で測定された衣料品のインフレ率を見ていきます。データのグラフは強い季節性を示しています。

[2]:

cpi_apparel = DataReader('CPIAPPNS', 'fred', start='1986')
cpi_apparel.index = pd.DatetimeIndex(cpi_apparel.index, freq='MS')
inf_apparel = np.log(cpi_apparel).diff().iloc[1:] * 1200
inf_apparel.plot(figsize=(15, 5));
../../../_images/examples_notebooks_generated_statespace_chandrasekhar_10_0.png

2つのモデルインスタンスを構築します。最初のものはカルマンフィルターの再帰を使用するように設定し、2 番目のものはチャンドラセカールの再帰を使用するように設定します。この設定は、以下のように ssm.filter_chandrasekhar プロパティで制御されます。

私たちが考えているモデルは季節的な自己回帰モデルであり、過去 15 年間の各年について、最初の6ヶ月を遅延項として、さらに指定された月も遅延項として含めます。これは状態ベクトルが次元 \(m = 186\) を持つことを意味し、この次元が十分大きいため、チャンドラセカールの再帰を使用することで性能向上が期待できるかもしれません。

注: 各モデルで tolerance=0 に設定しています。これにより、フィルターは予測共分散行列が収束したことを認識しなくなります。実際には推奨されません。ここでは、チャンドラセカールの再帰を毎期使用した場合の優れた性能を強調するためにこの設定をしています。後ほど、収束を許可するより現実的な設定での性能を示します。

[3]:

# Model that will apply Kalman filter recursions
mod_kf = sm.tsa.SARIMAX(inf_apparel, order=(6, 0, 0), seasonal_order=(15, 0, 0, 12), tolerance=0)
print(mod_kf.k_states)

# Model that will apply Chandrasekhar recursions
mod_ch = sm.tsa.SARIMAX(inf_apparel, order=(6, 0, 0), seasonal_order=(15, 0, 0, 12), tolerance=0)
mod_ch.ssm.filter_chandrasekhar = True

186

以下のコードを使用して、対数尤度関数の計算時間を計測します：

%timeit mod_kf.loglike(mod_kf.start_params)
%timeit mod_ch.loglike(mod_ch.start_params)

これにより、次のような結果が得られます：

171 ms ± 19.7 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
85 ms ± 4.97 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

この結果の意味は、今回の実験において、チャンドラセカール再帰を使用することでパフォーマンスが約 2 倍向上したことを示しています。

前述のように、前回の実験では予測共分散行列の収束を無効にしたため、そこでの結果は上限値となります。今回は、tolerance=0 引数を削除することで、通常通り収束を許可します。

[4]:

# Model that will apply Kalman filter recursions
mod_kf = sm.tsa.SARIMAX(inf_apparel, order=(6, 0, 0), seasonal_order=(15, 0, 0, 12))
print(mod_kf.k_states)

# Model that will apply Chandrasekhar recursions
mod_ch = sm.tsa.SARIMAX(inf_apparel, order=(6, 0, 0), seasonal_order=(15, 0, 0, 12))
mod_ch.ssm.filter_chandrasekhar = True

186

再度、次のコードを使用して対数尤度関数の計算時間を計測します：

%timeit mod_kf.loglike(mod_kf.start_params)
%timeit mod_ch.loglike(mod_ch.start_params)

これにより、次の結果が得られます：

114 ms ± 7.64 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
70.5 ms ± 2.43 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

チャンドラセカール再帰法は依然としてパフォーマンスを向上させますが、今ではその向上は約 33% にとどまります。これは、収束後には予測共分散行列を計算する必要がなくなるため、収束後の期間では両方のアプローチの計算時間に差がなくなるからです。以下では、収束が達成された期間を確認します：

[ ]:

res_kf = mod_kf.filter(mod_kf.start_params)
print('Convergence at t=%d, of T=%d total observations' %
      (res_kf.filter_results.period_converged, res_kf.nobs))

収束は比較的早く発生したため、すでに半分以上の期間で高価な行列乗算を回避しています。

しかし、前述のように、より大規模な DSGE モデルではサンプルのほとんどまたは全ての期間において収束しない場合があり、そのため最初の例と似たようなパフォーマンス向上を期待できるかもしれません。McAdam と Warne の 2019 年の論文「Euro area real-time density forecasting with financial or labor market frictions」では、彼らのアプローチについて「標準的なカルマンフィルターと比較して、これらの再帰法は、3 つのモデルの対数尤度計算を約 50％ 速くするという経験を持っている」と述べています。これは、私たちが最初の例で見つけた結果とほぼ同じです。

マルチスレッド行列代数ルーチンについての補足

上記のタイミングは、Anaconda 経由でインストールされたNumpyを基にしています。このインストールは Inte lの MKL BLAS および LAPACK ライブラリを使用しており、これらは行列代数の処理を高速化するためにマルチスレッド処理を実装しています。特に、ここで扱っている大きな行列に対する操作に役立ちます。この処理が結果にどのように影響するかを確認するために、ノートブックの最初のセルに以下のコードを追加してマルチスレッド処理を無効にすることができます。

import os
os.environ["MKL_NUM_THREADS"] = "1"

これを行うと、最初の例のタイミングは次のように変わります：

307 ms ± 3.08 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
97.5 ms ± 1.64 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

そして、2番目の例のタイミングは次のように変わります：

178 ms ± 2.78 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
78.9 ms ± 950 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

どちらも遅くなっていますが、典型的なカルマンフィルタはより大きく影響を受けます。

これは予想通りです。単一スレッドとマルチスレッドの線形代数のパフォーマンス差は、典型的なカルマンフィルタのケースでは非常に大きいためです。というのも、チャンドラセカール再帰の目的は行列演算のサイズを削減することだからです。つまり、マルチスレッド線形代数が使用できない場合でも、チャンドラセカール再帰はさらに大きなパフォーマンス向上を提供します。

チャンドラセカール再帰法と単変量フィルタリングアプローチ

チャンドラセカール再帰法を、KoopmanとDurbin（2000）の単変量フィルタリングアプローチと組み合わせることも可能です。これは、Aknouche と Hamdi の 2007 年の論文「周期的チャンドラセカール再帰法」の結果を利用することで実現できます。この組み合わせの初期実装は Statsmodels に含まれています。しかし、実験結果によると、この方法は通常のカルマンフィルタと比較してパフォーマンスが低下する傾向があることが示されています。これは、単変量フィルタリング法における計算の節約と一致しており、節約効果は状態ベクトルが観測ベクトルに比べて小さい場合に最も高いことが示唆されています。

参考文献

Aknouche, Abdelhakim, and Fayçal Hamdi. 「周期的チャンドラセカール帰納法.」 arXivプレプリント arXiv:0711.3857 (2007).

Herbst, Edward. 「DSGEモデルの尤度評価のための「チャンドラセカール帰納法」の使用.」 Computational Economics 45, no. 4 (2015): 693-705.

Koopman, Siem J., and James Durbin. 「多変量状態空間モデルの高速フィルタリングおよび平滑化.」 Journal of Time Series Analysis 21, no. 3 (2000): 281-296.

McAdam, Peter, and Anders Warne. 「金融市場または労働市場の摩擦を考慮したユーロ圏のリアルタイム密度予測.」 International Journal of Forecasting 35, no. 2 (2019): 580-600.

最終更新日: 2025年01月28日