2020/05/23

時系列分析 (5) - 状態空間モデル

AR, MA, ARMA, ARIMA, VAR, VARMA, VARIMA, ... などの時系列モデルは、観測値を直接モデル化するものでした。今回紹介する状態空間モデル (State-space model) の場合は、状態の時系列変化と、その状態から観測される値とに分けてモデル化する手法です。AR モデル、ARMA モデルなども状態空間表現を取ることができますし、状態空間モデルではより多様なシステムが記述できます。例えば、時変的な回帰係数を持つ回帰モデルなどが表現できます。
また、状態空間モデルではカルマンフィルタという強力なアルゴリズムによって、条件付き同時分布が効率よく計算できるため、ARMA のパラメータ推定を行うに際して、ARMA の状態空間表現に対してのカルマンフィルタによる推定が良く行われています。



ガウス型状態空間モデル

ガウス型状態空間モデルの定義

$Y_t$を$l$変量の時系列としたとき、これを以下の2つの方程式によって表現するモデルのことを状態空間モデルと呼びます。$$
\begin{eqnarray*}
X_t &=& F_tX_{t-1} + G_tv_t &\cdots システム方程式 \\
Y_t &=& H_tX_t + w_t &\cdots 観測方程式
\end{eqnarray*}
$$システム方程式は、システムの状態 ($X_t$) が1つ前の時点 ($X_{t-1}$) からどのように更新されるかを示しており、観測方程式は、システムの状態 ($X_t$) から観測値 ($Y_t$) がどのように取り出されるかを示しています。
$X_t$は状態ベクトルと呼ばれています。ここでは$k$変量のベクトルとします ($k$は状態空間モデルのパラメータになります)。$F_t$は回帰係数行列で$k$行$k$列の行列です。システムの過去の状態が現在にどのように影響を与えるかを示しています。$v_t$はシステムノイズと呼ばれ、システムの状態の変化の際に入り込む攪乱を表します。$m$変量の正規ホワイトノイズで、分散共分散行列$Q_t$に従うものとします。$G_t$はシステムノイズがシステムの状態に与える影響を表す$k$行$m$列の行列です。
$H_t$はシステムの状態が観測値にどのように射影されるかを表す行列で、$l$行$k$列の行列です。$w_t$は観測ノイズと呼ばれ、システムの状態を観測する際に入り込む攪乱を表します。$l$変量の正規ホワイトノイズで、分散共分散行列$R_t$に従うものとします。
上記のような条件を置いた状態空間モデルは、ガウス型状態空間モデルとも呼ばれています。
なお、状態空間モデルのシステム方程式、観測方程式は、真のモデルに対して一意ではありません。

AR モデルの状態空間表現

状態空間モデルの例として、$AR(p)$モデルの状態空間表現を見てみましょう。
以下の式で定義される$y_t$は、$AR(p)$モデルです。$$
y_t = \phi_1y_{t-1} + \phi_2y_{t-2} + \ldots + \phi_py_{t-p} + \epsilon_t, \epsilon_t \sim iidN(0, \sigma^2)
$$定数項は$0$とし、攪乱項の分布は正規ホワイトノイズとします。$x_t = (y_t, y_{t-1}, \ldots, y_{t-p+1})'$とおくと、以下の関係式が成り立ちます。$$
x_t = Fx_{t-1} + G\epsilon_t
$$ただし、$$
\begin{eqnarray*}
F &=& \left(\begin{array}{ccccc}
\phi_1 & \phi_2 & \ldots & \phi_p-1 & \phi_p \\
1 & 0 & \ldots & 0 & 0 \\
0 & 1 & \ldots & 0 & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \ldots & 1 & 0
\end{array}\right) (p 行 p 列の行列) \\
G &=& \left(\begin{array}{c}
1 \\
0 \\
0 \\
\vdots \\
0
\end{array}
\right) (p 次元のベクトル)
\end{eqnarray*}
$$です。これが状態空間モデルにおけるシステム方程式にあたります。システムノイズの分散共分散行列$Q$は$Q = \sigma^2$になります。
そして、$H = (1, 0, ..., 0)' (p 次元のベクトル)$と置けば、
$y_t = HX_t$
となり、これが観測方程式にあたります。観測ノイズの分散共分散行列$R$は$R = 0$となります。$F, G, H, Q, R$は、いずれも時不変なので、時不変の状態空間モデルです。

状態の推定

時系列データ$Y_t = (Y_1, Y_2, \ldots, Y_T)$に基づいて状態空間モデルを構築するにあたっては、状態$X_{t'} = (X_1, X_2, \ldots, X_T)$を推定する必要があります。また状態空間モデルから将来の観測値を予測するには、状態$X_{t'} = (X_{T+1}, X_{T+2}, \ldots)$の推定も必要です。
$X_{t'} (1 \le t' < T)$の推定を平滑化、$X_{t'} (t' = T)$の推定をフィルタ、$X_{t'} (t' > T)$の推定を予測と呼んでいます。

状態の推定は、$Y_t$が与えられた条件下での$X_{t'}$の条件付き分布$p(X_{t'}|Y_t)$を求める問題となります。状態空間モデルでは、ノイズが正規分布であるため、$p(X_{t'}|Y_t)$の分布も正規分布となります。つまり、$p(X_{t'}|Y_t)$の平均と分散共分散行列を求めれば、$X_{t'}$が推定できたことになります。

カルマンフィルタ

一般に$Y_t$が与えられた条件下での$X_{t'}$の条件付き同時分布を求める問題は、計算量が膨大になります。しかし、状態空間モデルにおいては逐次的に計算することが可能です。このアルゴリズムをカルマンフィルタと呼んでいます。
$X_{t'}$の条件付き分布を求めるには、前述のようにその平均と分散共分散行列を求めれば良いことが分かっています。条件付き平均、条件付き分散共分散行列をそれぞれ以下のように表記することにします。$$
\begin{eqnarray*}
X_{t'|t} &=& E(X_{t'}|Y_t) \\
V_{t'|t} &=& E((X_{t'} - X_{t'|t})(X_{t'} - X_{t'|t})')
\end{eqnarray*}
$$カルマンフィルタは、1期先予測 ($t = T - 1$の下での$t' = T$における$X_{t'}$の推定) と、フィルタ ($t = T$の下での$t' = T$における$X_{t'}$の推定) の2つの計算から成り立っています。
1期先予測の計算は、以下のようになります。$$
\begin{eqnarray*}
X_{t'|t'-1} &=& F_{t'}X_{t'-1|t'-1} \\
V_{t'|t'-1} &=& F_{t'}V_{t'-1|t'-1}F_{t'}' + G_{t'}Q_{t'}G_{t'}'
\end{eqnarray*}
$$1期先予測の平均は、現時点のフィルタに推移行列を掛けたものです。システムノイズの平均が$0$ですから、直感的にも分かります。分散共分散の第1項は推移行列の影響で、第2項はシステムノイズの影響を表しています。
フィルタの計算は、以下のようになります。$$
\begin{eqnarray*}
K_{t'} &=& V_{t'|t'-1}H_{t'}'(H_{t'}V_{t'|t'-1}H_{t'}' + R_{t'})^{-1} \\
X_{t'|t'} &=& X_{t'|t'-1} + K_{t'}(Y_t - H_{t'}X_{t'|t'-1}) = K_{t'}Y_t + (I - K_{t'}H_{t'})X_{t'|t'-1} \\
V_{t'|t'} &=& V_{t'|t'-1} - K_{t'}H_{t'}V_{t'|t'-1} = (I - K_{t'}H_{t'})V_{t'|t'-1}
\end{eqnarray*}
$$ $K_{t'}$はカルマンゲインと呼ばれています。$Y_t - H_{t'}X_{t'|t'-1}$は$Y_t$の観測誤差、$H_{t'}V_{t'|t'-1}H_{t'}' + R_{t'}$はその分散共分散行列を示しており、カルマンゲインは観測誤差を状態の推定誤差に繰り込む働きをしています。
フィルタの平均の後者の等式 ($X_{t'|t'} = K_{t'}Y_t + (I - K_{t'}H_{t'})X_{t'|t'-1}$) を見ると、フィルタが新しい観測値と1時点前までの情報による1期先予測との加重平均であると見ることもできます。
分散共分散の前者の形 ($V_{t'|t'} = V_{t'|t'-1} - K_{t'}H_{t'}V_{t'|t'-1}$) を見ると、新しい観測値によって状態推定の精度が改善される、と解釈することもできます。

非ガウス型状態空間モデル

ガウス型状態空間モデルでは、ノイズが正規分布に従うため、平均値から離れた値を取る可能性が低く、状態や観測値の変位が滑らかになります。従って、観測値の変位が急激な場合は、あまり適合が良くありません。
そこで、ノイズの分布を正規分布と仮定せず、一般のホワイトノイズであるとの仮定に緩めたものが、非ガウス型状態空間モデルです。ノイズにファットテールな分布を持ってくれば、観測値の急激な変化もモデルに織り込むことができる可能性が出てきます。

非ガウス型状態空間モデルの定義

定式化すると、以下のようになります。$$
\begin{eqnarray*}
X_t &=& F_tX_{t-1} + G_tv_t &\cdots システム方程式 \\
Y_t &=& H_tX_t + w_t &\cdots 観測方程式
\end{eqnarray*}
$$方程式の形状はガウス型と同じで、各係数行列の条件も同じですが、$v_t, w_t$はそれぞれ、密度関数$q(v), r(w)$に従うホワイトノイズであるとします。

状態の推定

推定は、カルマンフィルタの拡張によって行うことができます。ガウス型では条件付き同時分布が正規分布であるという条件を置くことができましたが、非ガウス型ではそうとは限らないため、条件付き同時分布を直接求める必要があります。
1期先予測は、以下のようになります。$$
p(X_{t'}|Y_{t-1}) = \int p(X_{t'}|X_{t'-1}) p(X_{t'-1}|Y_1, \ldots, Y_{t-1}) dX_{t'-1}
$$フィルタは、以下のようになります。$$
p(X_{t'}|Y_t) = \frac{p(Y_{t'}|X_{t'}) p(X_{t'}|Y_1, \ldots, Y_{t-1})}{p(Y_{t'}|Y_1, \ldots, Y_{t-1})}
$$

一般化状態空間モデル

非ガウス型状態空間モデルでは、状態の推移の線形性を仮定していましたが、その仮定を緩め、非線形な推移も許したものが一般化状態空間モデルです。
例えば、確率的ボラティリティモデルのパラメータ推定を行うために、非線形変換を加えた上で、線形のガウス型状態空間モデルで近似して推定を行う方法がありますが、一般化状態空間モデルでは確率的ボラティリティモデルを近似無しに表現することができます。

一般化状態空間モデルの定義

一般化状態空間モデルは、以下のように定式化されます。$$
\begin{eqnarray*}
X_t &=& F(X_{t-1}, v_t) &\cdots システム方程式 \\
Y_t &=& H(X_t, w_t) &\cdots 観測方程式
\end{eqnarray*}
$$ $v, w$は非ガウス型状態空間モデルと同様、密度関数$q(v), r(w)$に従うホワイトノイズであるとします。$F(X, v)$および$H(X, w)$は非線形の関数です。$H(X, w)$については、$Y = H(X, w)$に対して$w = G(Y, X)$を満たす$w$を一意に決定できる、$Y$について微分可能な関数$G(Y, X)$が存在することを要求します。また、初期状態$X_0$は密度関数$p_0(X)$に従うものとします。

モデルの推定

モデルの推定は解析的に行うことは一般にはできません。モンテカルロ・フィルタと呼ばれるアルゴリズムによって行います。モンテカルロ・フィルタでは、仮定した分布に従う「粒子」を乱数で生成し、解析的に表現できない同時分布をこの乱数に基づく経験分布関数で近似するという方法をとります。概要としては以下の通りです。

  1. 用いる「粒子」の数を$m$とし、以下、$j = 1, \ldots, m$とする
  2. $p_0(X)$に従う$m$個の乱数ベクトル$f_0^{(j)}$を生成する (初期状態を乱数で設定)
  3. $t = 1, \ldots, T$について、以下の手順を実行する
    1. $q(v)$に従う$m$個の乱数ベクトル$v_t^{(j)}$を生成する (システムノイズを乱数で生成)
    2. 各$j$について、$p_n^{(j)} = F(f_{n-1}^{(j)}, v_n^{(j)})$を計算する (新しい状態の候補値の算出)
    3. 各$j$について、$\alpha_n^{(j)} = r(G(Y_n, p_n^{(j)})) |\partial G/\partial Y_n|$を計算する (状態$p_n^{(j)}$が観測値$Y_n$となる確率の算出)
    4. $p_n^{(j)}$の中から$\alpha_n^{(j)}$に比例する確率で$m$回サンプリングを行い、その値を$f_n^{(1)}, \ldots, f_n^{(m)}$とする (新しい状態を決定)

サンプリングは、 ランダムサンプリングでも良いのですが、計算量を減らすため階層化リサンプリング法といったものを用いることもできます。