状態空間モデル

時系列解析の線形モデルの基本として、ARモデルとMAモデルがあるが、それらの組み合わせや、目的変数の次元、説明変数の有無等により、いくつかのバリエーションがある。AR成分とMA成分両方を考慮するならARMAモデル、目的変数が複数のARモデルなら多変量AR(VAR)、ARMAモデルに説明変数を含めるならARMAXモデルなど。

いずれにしても、これらのモデルは、状態空間モデルという一般化された枠組みの中で、カルマンフィルタのアルゴリズムにより解くことができる。すなわち、これらのモデルは状態空間表現を持つ。

以前、状態空間モデルについて調査していた際、ARMAモデルやVAR、ARMAXの状態空間表現は、書籍¹ ²で具体的な表式を見つけられたが、多変量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モデルに対する状態空間表現を求めた。

北川源四郎, 時系列解析入門 ↩
James Durbin, Siem Jan Koopman, Time Series Analysis by State Space Methods ↩