状態空間モデル
時系列解析の線形モデルの基本として、ARモデルとMAモデルがあるが、それらの組み合わせや、目的変数の次元、説明変数の有無等により、いくつかのバリエーションがある。AR成分とMA成分両方を考慮するならARMAモデル、目的変数が複数のARモデルなら多変量AR(VAR)、ARMAモデルに説明変数を含めるならARMAXモデルなど。
いずれにしても、これらのモデルは、状態空間モデルという一般化された枠組みの中で、カルマンフィルタのアルゴリズムにより解くことができる。すなわち、これらのモデルは状態空間表現を持つ。
以前、状態空間モデルについて調査していた際、ARMAモデルやVAR、ARMAXの状態空間表現は、書籍1 2で具体的な表式を見つけられたが、多変量ARMAX(VARMAX)の状態空間表現を見つけられなかった。
この記事では、VARMAXの状態空間表現の導出を行なう。
VARMAXの状態空間表現の導出
$\boldsymbol{y}_n$を$p$次元の目的変数ベクトルとして、$q$次元の説明変数$\boldsymbol{x}_n$の回帰式により
\begin{aligned}
\boldsymbol{y}_n = B\boldsymbol{x}_n + \boldsymbol{\xi}_n \label{regress_arma}
\end{aligned}
で表されるものとする。ここで$B$は$p\times q$次元の回帰係数の行列である。$\boldsymbol{\xi}_n$は残差の$p$次元ベクトルで、次のARMA過程 ARMA$(m, \ell)$ に従うものとする。
\begin{aligned}
\boldsymbol{\xi}_n = \sum_{i=1}^{m} \Phi_i \boldsymbol{\xi}_{n-i} + \boldsymbol{\zeta}_n - \sum_{i=1}^{\ell} \Theta_i \boldsymbol{\zeta}_{n-i},
\ \ \ \ \
\boldsymbol{\zeta}_n \sim N(0,\Sigma_{\zeta}). \label{arma_error}
\end{aligned}
$\Phi_i$は$\phi_{jk}$を($j,k$)成分に持つAR係数行列、$\Theta_i$は$\theta_{jk}$を($j,k$)成分に持つMA係数行列である。
$\boldsymbol{\zeta}_n$は$p$次元白色雑音で、
\begin{aligned}
&E(\boldsymbol{\zeta}_n) = (0,0,\cdots,0)^T, \\
&E(\boldsymbol{\zeta}_n\boldsymbol{\zeta}_n^T) =
\begin{bmatrix}
\sigma_{11} & \cdots & \sigma_{1p} \\
\vdots & \ddots & \vdots \\
\sigma_{p1} & \cdots & \sigma_{pp} \\
\end{bmatrix} =: \Sigma_{\zeta}, \\
&E(\boldsymbol{\zeta}_n\boldsymbol{\zeta}_m^T) = {\rm O}_{p\times p},\ \ \ \ \ n \ne mのとき \\
&E(\boldsymbol{\zeta}_n\boldsymbol{\xi}_m^T) = {\rm O}_{p\times p},\ \ \ \ \ n > mのとき
\end{aligned}
を満たす。ただし、$T$は転置、${\rm O}_{p\times p}$はすべての成分が0の$p\times p$の行列を表し、$\Sigma_{\zeta}$は分散共分散行列である。
ここで、$\tilde{\boldsymbol{\xi}}_{n+i|n-1}$を
\begin{aligned}
\tilde{\boldsymbol{\xi}}_{n+i|n-1} = \sum_{j=i+1}^{m} \Phi_j \boldsymbol{\xi}_{n+i-j} - \sum_{j=i}^{\ell}\Theta_j \boldsymbol{\zeta}_{n+i-j}
\end{aligned}
により定義し、$k = {\rm max}(m, \ell+1)$とすると、次の関係が成り立つ。
\begin{aligned}
\boldsymbol{\xi}_n &= \Phi_1\boldsymbol{\xi}_{n-1} + \tilde{\boldsymbol{\xi}}_{n|n-2} + \boldsymbol{\zeta}_n, \\
\tilde{\boldsymbol{\xi}}_{n+i|n-1} &= \Phi_{i+1}\boldsymbol{\xi}_{n-1} + \tilde{\boldsymbol{\xi}}_{n+i|n-2} - \Theta_i\boldsymbol{\zeta}_n, \\
\tilde{\boldsymbol{\xi}}_{n+k-1|n-1} &= \Phi_{k}\boldsymbol{\xi}_{n-1} - \Theta_{k-1}\boldsymbol{\zeta}_n.
\end{aligned}
ただし、$i>m$のとき$\Phi_i = {\rm O}_{p\times p}$、$i>\ell$のとき$\Theta_i = {\rm O}_{p\times p}$とする。したがって、$\boldsymbol{\alpha}_n$を
\begin{aligned}
\boldsymbol{\alpha_n} &=
(\boldsymbol{\xi}_n^T, \tilde{\boldsymbol{\xi}}_{n+1|n-1}^T, \cdots, \tilde{\boldsymbol{\xi}}_{n+k-1|n-1}^T)^T \\
&=
\begin{bmatrix}
\xi_n(1) \\
\vdots \\
\xi_n(p) \\
\tilde{\xi}_{n+1|n-1}(1) \\
\vdots \\
\tilde{\xi}_{n+1|n-1}(p) \\
\tilde{\xi}_{n+2|n-1}(1) \\
\vdots \\
\vdots \\
\tilde{\xi}_{n+k-1|n-1}(p) \\
\end{bmatrix}
\end{aligned}
の$kp$次元ベクトルとして定義すると、次の状態空間表現が得られる。
\begin{aligned}
\boldsymbol{\alpha}_{n} &= T\boldsymbol{\alpha}_{n-1} + R \boldsymbol{\zeta}_{n}, \ \ (kp次元) \\
\boldsymbol{\xi}_n &= Z\boldsymbol{\alpha}_{n}, \ \ (p次元) \\
T &=
\begin{bmatrix}
\Phi_1 & I_p & & & \\
\Phi_2 & & I_p & & \\
\vdots & & & \ddots & \\
\Phi_{k-1} & & & & I_p \\
\Phi_{k} & {\rm O}_{p\times p} & {\rm O}_{p\times p} & \cdots & {\rm O}_{p\times p} \\
\end{bmatrix}, \ \ (kp\times kp次元) \\
R &=
\begin{bmatrix}
I_p \\
- \Theta_1 \\
\cdots \\
- \Theta_{k-1} \\
\end{bmatrix},\ \ (kp\times p次元) \\
Z &=
\begin{bmatrix}
I_p & {\rm O}_{p\times p} & \cdots & {\rm O}_{p\times p} \\
\end{bmatrix}.\ \ (p\times kp次元)
\end{aligned}
さらに、ここで$\boldsymbol{\beta}_i$を$q$次元ベクトル
\begin{aligned}
\boldsymbol{\beta}_i = (B(i,1), B(i,2), \cdots, B(i, q))^T,\ \ i=1, 2, \cdots, p \label{beta}
\end{aligned}
とし、$\boldsymbol{\alpha}_n^{\ast}$を$(k+q)p$次元ベクトル
\begin{aligned}
\boldsymbol{\alpha^{\ast}_n} &=
\begin{bmatrix}
\boldsymbol{\beta}_1 \\
\vdots \\
\boldsymbol{\beta}_p \\
\boldsymbol{\alpha}_n
\end{bmatrix}, \\
\end{aligned}
と定義すると、次の表式が得られる(行列$Z^{\ast}$の破線は、表現を見やすくするために書いた)。
\begin{aligned}
\boldsymbol{\alpha}^{\ast}_{n} &= T^{\ast}\boldsymbol{\alpha}^{\ast}_{n-1} + R^{\ast} \boldsymbol{\zeta}^{\ast}_{n}, \\
\boldsymbol{y}_n &= Z^{\ast}_n\boldsymbol{\alpha}^{\ast}_{n}, \\
T^{\ast}&=
\begin{bmatrix}
I_{q} & & & & \\
& I_q & & & \\
& & \ddots & & \\
& & & I_q & \\
& & & & T \\
\end{bmatrix},\ \ \ \ \ ((k+q)p\times (k+q)p次元) \\
R^{\ast} &=
\begin{bmatrix}
{\rm O}_{pq\times p} \\
R
\end{bmatrix},\ \ ((k+q)p\times p次元) \\
Z^{\ast} &=
\left[
\begin{array}{c:c}
\begin{matrix}
\boldsymbol{x}^T_n & & & \\
& \boldsymbol{x}^T_n & & \\
& & \ddots & \\
& & & \boldsymbol{x}^T_n \\
\end{matrix}
&
\begin{matrix}
\textit{I}_p & {\rm O}_{p\times p} & \cdots & {\rm O}_{p\times p} \\
\end{matrix}
\end{array}
\right].\ \ (p\times (k+q)p次元)
\end{aligned}
すなわち、VARMAXモデルに対する、状態空間表現が得られた。
p=1の場合 (ARMAX)
この場合のモデルは、回帰係数のベクトルを $\boldsymbol{\beta} = (\beta_1, \beta_2, \cdots, \beta_q)$ として、
\begin{aligned}
y_n &= \boldsymbol{\beta}\cdot\boldsymbol{x}_n + \xi_n \\
\xi_n &= \sum_{i=1}^{m} \phi_i \xi_{n-i} + \zeta_n - \sum_{i=1}^{\ell}\theta_i\zeta_{n-i}
\end{aligned}
となり、状態空間表現は以下により与えられる。
\begin{aligned}
\boldsymbol{\alpha}^{\ast}_{n} &= T^{\ast}\boldsymbol{\alpha}^{\ast}_{n-1} + R^{\ast} \zeta_{n}, \\
y_n &= Z^{\ast}_n\boldsymbol{\alpha}^{\ast}_{n}, \\
\boldsymbol{\alpha}_n^{\ast} &= (\beta_1, \beta_2, \cdots, \beta_q, \xi_n, \tilde{\xi}_{n+1|n-1}, \tilde{\xi}_{n+2|n-1}, \cdots, \tilde{\xi}_{n+k-1|n-1})^T, \\
T^{\ast} &=
\begin{bmatrix}
I_{q} & & & & \\
& \phi_1 & 1 & & & \\
& \phi_2 & & 1 & & \\
& \vdots & & & \ddots & \\
& \phi_{k-1} & & & & 1 \\
& \phi_{k} & & & & \\
\end{bmatrix},\ \ ((k+q)\times (k+q)次元) \\
R^{\ast} &=
\begin{bmatrix}
0 \\
0 \\
\vdots \\
0 \\
1 \\
-\theta_1 \\
-\theta_2 \\
\vdots \\
-\theta_{k-1} \\
\end{bmatrix},\ \ ((k+q)\times 1次元) \\
Z^{\ast}_n &=
\begin{bmatrix}
x_n(1) & x_n(2) & \cdots & x_n(q) & 1 & 0 & \cdots & 0 \\
\end{bmatrix},\ \ (1\times (k+q)次元) \\
\end{aligned}
p=1, q=0の場合 (ARMA)
この場合は、$\boldsymbol{\alpha}_n$を用いた状態空間表現で、$p=1$, $q=0$とすることにより、以下の状態空間表現が得られる。
\begin{aligned}
\boldsymbol{\alpha}_{n} &= T\boldsymbol{\alpha}_{n-1} + R {\zeta}_{n}, \ \ (k次元) \\
{\xi}_n &= Z\boldsymbol{\alpha}_{n}, \ \ (1次元) \\
T &=
\begin{bmatrix}
\phi_1 & 1 & & & \\
\phi_2 & & 1 & & \\
\vdots & & & \ddots & \\
\phi_{k-1} & & & & 1 \\
\phi_{k} & 0 & 0 & \cdots & 0 \\
\end{bmatrix}, \ \ (k\times k次元) \\
R &=
\begin{bmatrix}
1 \\
- \theta_1 \\
\cdots \\
- \theta_{k-1} \\
\end{bmatrix},\ \ (k\times 1次元) \\
Z &=
\begin{bmatrix}
1 & 0 & \cdots & 0 \\
\end{bmatrix}.\ \ (1\times k次元)
\end{aligned}
これは[北川源四郎, 時系列解析入門]で与えられているものと同一の表現である。
q=0, ℓ=0の場合 (VAR)
この場合も同様に、$\boldsymbol{\alpha}_n$を用いた状態空間表現より、$\Theta_i = {\rm O}_{p\times p}$, $k=m$として、以下の状態空間表現が得られる。
\begin{aligned}
\boldsymbol{\alpha}_{n} &= T\boldsymbol{\alpha}_{n-1} + R \boldsymbol{\zeta}_{n}, \ \ (mp次元) \\
\boldsymbol{\xi}_n &= Z\boldsymbol{\alpha}_{n}, \ \ (p次元) \\
T &=
\begin{bmatrix}
\Phi_1 & I_p & & & \\
\Phi_2 & & I_p & & \\
\vdots & & & \ddots & \\
\Phi_{k-1} & & & & I_p \\
\Phi_{k} & {\rm O}_{p\times p} & {\rm O}_{p\times p} & \cdots & {\rm O}_{p\times p} \\
\end{bmatrix}, \ \ (kp\times kp次元) \\
R &=
\begin{bmatrix}
I_p \\
{\rm O}_{p\times p} \\
\cdots \\
{\rm O}_{p\times p} \\
\end{bmatrix},\ \ (mp\times p次元) \\
Z &=
\begin{bmatrix}
I_p & {\rm O}_{p\times p} & \cdots & {\rm O}_{p\times p} \\
\end{bmatrix}.\ \ (p\times mp次元)
\end{aligned}
まとめ
説明変数を考慮した多変量ARMAモデルに対する、状態空間表現の具体的な表式を得た。また、その表現から、ARMAX, ARMA, VARモデルに対する状態空間表現を求めた。