狀態空間模型 - Chandrasekhar 遞迴¶

[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

雖然大多數與狀態空間模型相關的操作都依賴於卡爾曼濾波遞迴，但在某些特殊情況下，可以使用一種稱為「Chandrasekhar 遞迴」的獨立方法。這些方法提供了一種迭代計算狀態向量條件動量的替代方法，並且在某些情況下，它們的計算強度可能比卡爾曼濾波遞迴低得多。有關完整詳細資訊，請參閱論文「使用「Chandrasekhar 遞迴」進行 DSGE 模型概似評估」(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\)，可以分為三個部分

初始化：指定 \(a_t\) 和 \(P_t\)，它們定義條件狀態分佈，\(\alpha_t \mid y^{t-1} \sim N(a_t, P_t)\)。
更新：計算 \(a_{t|t}\) 和 \(P_{t|t}\)，它們定義條件狀態分佈，\(\alpha_t \mid y^{t} \sim N(a_{t|t}, P_{t|t})\)。
預測：計算 \(a_{t+1}\) 和 \(P_{t+1}\)，它們定義條件狀態分佈，\(\alpha_{t+1} \mid y^{t} \sim N(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) 模型（這裡我們使用 \(r\)，因為我們已經將 \(p\) 用作觀測向量的維度）可以放入狀態空間形式為

\[\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\) 個元素，因此花費在計算概似函數（通過卡爾曼濾波器）的大部分時間可能會被預測步驟中的矩陣乘法所支配。

狀態空間模型和 Chandrasekhar 遞迴¶

Chandrasekhar 遞迴將方程式 \(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)。

重要注意事項：與卡爾曼濾波器不同，Chandrasekhar 遞迴不能用於每個狀態空間模型。特別是，後者有以下限制（前者使用不需要的限制）

除了允許時變截距外，該模型必須是時不變的。
必須使用狀態向量的靜態初始化（這排除了非靜態元件中的所有模型）
不允許遺失資料

要理解為什麼這個公式可以意味著更有效率的計算，請再次考慮上面的 SARIMAX 案例。在這種情況下，\(p = 1\)，因此 \(M_t\) 是一個純量，我們可以將 Chandrasekhar 遞迴重寫為

\[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

我們將建構兩個模型實例。第一個將設定為使用卡爾曼濾波器遞迴，而第二個將設定為使用 Chandrasekhar 遞迴。此設定由 ssm.filter_chandrasekhar 屬性控制，如下所示。

我們心中的模型是一個季節性自迴歸，其中我們包括前 6 個月作為滯後，以及過去 15 年中每個月的給定月份作為滯後。這意味著狀態向量的維度為 \(m = 186\)，這已經足夠大，我們可能會期望通過使用 Chandrasekhar 遞迴看到一些實質性的效能提升。

備註：我們在每個模型中設定 tolerance=0 - 這具有阻止濾波器永遠無法識別預測共變異數矩陣已收斂的效果。這在實務上不建議。我們在這裡這樣做是為了突出 Chandrasekhar 遞迴在每個週期中使用時，而不是典型的卡爾曼濾波器遞迴時的卓越效能。稍後，我們將在更實際的設定中展示效能，我們確實允許收斂。

[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

我們使用以下程式碼來測量對數概似函數的計算時間

%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)

這意味著在此實驗中，Chandrasekhar 遞迴將效能提高了約 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

同樣地，我們計時對數概似函數的計算，使用以下程式碼：

%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)

Chandrasekhar 遞迴仍然可以提升效能，但現在只提升約 33%。原因是在收斂之後，我們不再需要計算預測的共變異數矩陣，因此對於那些收斂後的期間，兩種方法的計算時間將不會有差異。下面我們檢查達到收斂的期間。

[5]:

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))

Convergence at t=186, of T=463 total observations

由於收斂發生得相對較早，我們已經在超過一半的期間避免了昂貴的矩陣乘法。

然而，如上所述，較大的 DSGE 模型可能在樣本中的大多數或所有期間都無法達到收斂，因此我們或許可以期望獲得更類似於第一個範例的效能提升。在他們 2019 年的論文《Euro area real-time density forecasting with financial or labor market frictions》中，McAdam 和 Warne 指出，在他們的應用中，「與標準卡爾曼濾波器相比，我們的經驗是，這些遞迴在三個模型的對數概似計算速度上大約加快了 50%」。這與我們在第一個範例中發現的結果大致相同。

關於多執行緒矩陣代數常式¶

上面的計時是基於透過 Anaconda 安裝的 Numpy，它使用 Intel 的 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)

第二個範例的計時會變為：

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)

兩者都變慢了，但典型的卡爾曼濾波器受到的影響更大。

這並不意外；在典型的卡爾曼濾波器情況下，單執行緒和多執行緒線性代數之間的效能差異要大得多，因為 Chandrasekhar 遞迴的重點是減少矩陣運算的大小。這意味著，如果無法使用多執行緒線性代數，Chandrasekhar 遞迴會提供更大的效能提升。

Chandrasekhar 遞迴與單變量濾波方法¶

也可以將 Chandrasekhar 遞迴與 Koopman 和 Durbin (2000) 的單變量濾波方法結合使用，利用 Aknouche 和 Hamdi 在他們 2007 年的論文《Periodic Chandrasekhar recursions》中的結果。Statsmodels 中包含此組合的初始實現。然而，實驗表明，與即使是通常的卡爾曼濾波器相比，這往往會降低效能。這與單變量濾波方法報告的計算節省相符，這表明當狀態向量相對於觀測向量較小時，節省的幅度最高。

參考文獻¶

Aknouche, Abdelhakim, and Fayçal Hamdi. “Periodic Chandrasekhar recursions.” arXiv preprint arXiv:0711.3857 (2007).

Herbst, Edward. “Using the “Chandrasekhar Recursions” for likelihood evaluation of DSGE models.” Computational Economics 45, no. 4 (2015): 693-705.

Koopman, Siem J., and James Durbin. “Fast filtering and smoothing for multivariate state space models.” Journal of Time Series Analysis 21, no. 3 (2000): 281-296.

McAdam, Peter, and Anders Warne. “Euro area real-time density forecasting with financial or labor market frictions.” International Journal of Forecasting 35, no. 2 (2019): 580-600.

最後更新：2024 年 10 月 03 日