状態空間モデルシリーズでは複数の記事に分けて、時系列データを扱う柔軟なモデルである状態空間モデル(State Space Model) を解説します。
第1回では、時系列データに対するモデリングの一つとして構造的時系列モデル（季節調整モデル）を扱いました。
第2回である本記事では、季節調整モデルを発展させて、状態空間モデルの定式化を行います。具体的には、状態空間モデルの核となる「状態ベクトル」の設計思想を解説し、複雑な時系列の依存関係を確率モデルとして計算可能な形で表現する方法を示します。
Contents
状態ベクトル
ここでは、第1回で扱った構造的時系列モデルを基に、時系列モデルに状態ベクトルを導入する。
なぜ状態ベクトルが必要か？
状態ベクトル表現を導入することで、複数の成分を統一的に扱うことができ、確率モデルに基づいた確率的推論が可能になる。具体的には以下の3つの理由により状態ベクトルを導入する。
- 統一的な表記
- トレンド成分 $\mu_t, \mu_{t-1}$、季節成分 $s_t, s_{t-1}, \cdots, s_{t-6}$ など、複数の変数を一つのベクトル $x_t$ にまとめる
- これにより、複雑な依存関係を簡潔な行列演算で表現できる
- 確率推論の基盤
- 確率モデルとして簡潔に表現できるため、後述する同時分布の分解やマルコフ性を利用することで、フィルタ分布、予測分布、平滑化分布を計算することが可能になる
- そしてこれらが統一的なアルゴリズムで扱えるようになる
- 実装の効率性
- 行列演算としてアルゴリズムが表現できるため、計算が容易（カルマンフィルタ）
ベイズ推論の枠組みに自然に組み込むことができるようになる。
トレンドモデル
まず簡単のために、トレンド成分だけを扱うトレンドモデルの状態ベクトル表現を解説する。
トレンドモデルは以下の式で定義される（参考文献、および、前回記事参照）
$$
\begin{aligned}
y_t &= \mu_t + w_t \\
\mu_t &= 2\mu_{t-1} – \mu_{t-2} + v_t
\end{aligned}
$$
ここで、 $w_t, v_t$ はガウス分布に従うノイズ（正規ホワイトノイズ）とする。
二つ目のトレンド成分を示す式について、状態ベクトル $x_t$ を以下の通り定義する。
$$
x_t =
\begin{bmatrix} \mu_t \\ \mu_{t-1} \end{bmatrix}
$$
$x_t$ を導入すると、トレンド成分は以下のように表現できる。
$$
\begin{aligned}
\begin{bmatrix} \mu_t \\ \mu_{t-1} \end{bmatrix} &=
\begin{bmatrix} 2 & -1 \\ 1 & 0 \end{bmatrix}
\begin{bmatrix} \mu_{t-1} \\ \mu_{t-2} \end{bmatrix} + \begin{bmatrix} 1 \\ 0 \end{bmatrix} v_t \\
x_t &= F x_{t-1} + G v_t
\end{aligned}
$$
同様にして、観測される時系列データ $y_t$ は以下の通り表すことができる。
$$
\begin{aligned}
y_t &= \begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} \mu_t \\ \mu_{t-1} \end{bmatrix} + w_t \\
y_t &= Hx_t + w_t
\end{aligned}
$$
このように、行列 $F, G, H$ と状態ベクトルとの線形結合の形で簡潔に表現することができる。
ここで、状態ベクトル $x_t$ としてなぜ $\mu_{t-1}$ も含めるのか？について考える。
これは、トレンドの更新式が $\mu_t = 2\mu_{t-1} – \mu_{t-2} + v_t$ と2時点前の値を必要とするためである。状態ベクトルに過去の必要な情報をすべて含めることで、1時点前の状態ベクトル $x_{t-1}$ のみから現在の状態 $x_t$ を計算できるようになる。これはとても重要で、この性質（マルコフ性）により、後述する同時分布の分解が可能となり、効率的な計算アルゴリズムが構築できる。
季節調整モデル
次に、季節調整モデルの状態ベクトル表現を扱う。
季節調整モデルは以下の式で定義される（参考文献、および、前回記事参照）
$$
\begin{aligned}
y_t &= \mu_t + s_t + w_t \\
&\mu_t = 2 \mu_{t-1} – \mu_{t-2} + v_{\mu, t} \\
&s_t = -\sum_{l=1}^{T-1} s_{t-l} + v_{s, t}
\end{aligned}
$$
これらの式は上から、観測方程式、トレンド成分、季節成分である。また、$w_t, v_{\mu, t}, v_{s, t}$ はガウス分布に従うノイズ（正規ホワイトノイズ）とする。
季節調整モデルの場合、状態ベクトルとしては、 $\mu_t, \mu_{t-1}, s_t, \cdots, s_{t-5}$ を必要とするため、単純に並べて以下のように定義する。なお、季節成分の周期は $T=7$ (週次の成分)とする。
$$
x_t = [\mu_t, \mu_{t-1}, s_t, s_{t-1}, s_{t-2}, s_{t-3}, s_{t-4}, s_{t-5}]^T
$$
この $x_t$ を使うと、季節調整モデルは以下のように表現できる。なおこの式は行列のインパクトが大きく慣れていないと困惑するかもしれない。しかしよく見ると、季節成分 $s_{t-l}$ の周期分があるだけなので、落ち着いて見ることが重要である。決して難しい話をしているわけではない。
$$
\begin{aligned}
{\left[\begin{array}{c}
\mu_t \\
\mu_{t-1} \\
s_{t-1} \\
s_{t-2} \\
s_{t-3} \\
s_{t-4} \\
s_{t-5}
\end{array}\right]}
& =\left[\begin{array}{cccccccc}
2 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & -1 & -1 & -1 & -1 & -1 & -1 \\
0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0
\end{array}\right]\left[\begin{array}{l}
\mu_{t-1} \\
\mu_{t-2} \\
s_{t-1} \\
s_{t-2} \\
s_{t-3} \\
s_{t-4} \\
s_{t-5} \\
s_{t-6}
\end{array}\right]+\left[\begin{array}{cc}
1 & 0 \\
0 & 0 \\
0 & 1 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0
\end{array}\right]\left[\begin{array}{l}
v_{\mu, t} \\
v_{s, t}
\end{array}\right] \\
x_t &= F x_{t-1} + G v_t
\end{aligned}
$$
最後の式にある通り、行列を $F, G$ と表すことで、先に示したトレンドモデルと同じ形で書けることがわかる。
同様に観測される時系列データ $y_t$ も以下の通り表すことができる。
$$
\begin{aligned}
y_t &= \begin{bmatrix} 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \mu_t \\ \mu_{t-1} \\ s_{t} \\ s_{t-1} \\ s_{t-2} \\ s_{t-3} \\ s_{t-4} \\ s_{t-5} \end{bmatrix} + w_t \\
y_t &= Hx_t + w_t
\end{aligned}
$$
線形ガウス状態空間モデル
上記の通り、トレンドモデル、季節調整モデルはベクトルと行列で下記の通り表すことができる。
$$
\begin{aligned}
x_t &= Fx_{t-1} + Gv_t, \quad v_t \sim N(0, Q) \\
y_t &= Hx_t + w_t, \quad w_t \sim N(0, R)
\end{aligned}
$$
ここで各変数の意味を整理すると以下の通り
- $x_t$: 状態ベクトル（次元: $n \times 1$）
- $y_t$: 観測ベクトル（次元: $m \times 1$）
- $F$: 状態遷移行列（次元: $n \times n$）
- $G$: システムノイズの係数行列（次元: $n \times p$）
- $H$: 観測行列（次元: $m \times n$）
- $v_t$: システムノイズ（次元: $p \times 1$）
- $w_t$: 観測ノイズ（次元: $m \times 1$）
- $Q$: システムノイズの共分散行列（次元: $p \times p$）
- $R$: 観測ノイズの共分散行列（次元: $m \times m$）
これらの式はノイズ成分 $v_t, w_t$ が(多次元)ガウス分布に従うとしているため、全体でガウス分布に従う確率モデルとして書くこともできる。
$$
\begin{aligned}
x_t & \sim N( Fx_{t-1}, GQG^T) \\
y_t & \sim N(Hx_t, R)
\end{aligned}
$$
上記の一つ目の式をシステムモデル（状態方程式）、下段を観測モデル（観測方程式）と呼び、二つを合わせて線形ガウス状態空間モデルと呼ぶ。
(一般)状態空間モデル
次に、上記の線形ガウス状態空間モデルを一般化していくことを考える。
線形ガウス状態空間モデルの名称は以下の2つの特徴に由来する
- 「線形」: 状態遷移と観測の関係が状態ベクトル $x_{t-1}, x_t$ に対して線形である
- $F \cdot x_{t-1}$, $H \cdot x_t$のように定義されている
- 「ガウス」: システムノイズと観測ノイズがガウス分布に従う
- $v_t \sim N(0, Q)$, $w_t \sim N(0, R)$
- この結果、$x_t, y_t$ が全てガウス分布となる
これらの制約を取り除くことで、下記のように一般形式でシステムモデル、観測モデルを表現することができる
$$
\begin{aligned}
x_t & \sim p(x_t | x_{t-1}) \\
y_t & \sim p(y_t | x_t)
\end{aligned}
$$
ここで、
- $p(x_t | x_{t-1})$: 状態遷移の条件付き確率密度関数
- $p(y_t | x_t)$: 観測の条件付き確率密度関数
これらのシステムモデル、観測モデルを合わせて（一般）状態空間モデルと呼ぶ。一般状態空間モデルでは、ガウス分布を仮定せず、一般の何らかの確率モデル $p(\cdot)$ を仮定する。また、線形の関係も仮定せず、$x_t$ と $x_{t-1}$ の関係に対して何らかの関数 $f(x_t | x_{t-1})$、また、 $y_t$と$x_t$ の関係に対しても $g(y_t | x_{t})$ を仮定する。
別の回で詳しく解説するが、線形ガウス状態空間モデルは解析的に解を導出できるのに対し、一般状態空間モデルでは解析的な解を導出することができない。そのため、粒子フィルタ（particle filtering）などの数値的な推論が必要になる。
同時分布の分解
ここでは、今後の推論の基礎になる確率モデルの分解を扱う。
時系列データの確率モデル表現
時系列データは、第1回で述べた通り $y_{1:T} = [y_1, y_2, \cdots, y_{T}]$ のように表される。確率モデルとしては、データ $y_{1:T}$ の同時確率 $p(y_{1:T})$ として表現できる。独立同分布（i.i.d）を仮定できるデータの場合には、各データは独立なので単純に $\prod^T_{i=1} p(y_i)$ と分解できるが、時系列データは順番に意味がありデータ間には依存関係がある。
同時確率の分解は一般に以下のようになる
$$
p(A, B) = p(A | B) p(B)
$$
では、 $p(y_{1:T})$ についても上記の関係を利用して一つずつ分解してみる
$$
p(y_{1:T}) = p(y_T | y_{1:T-1}) p( y_{1:T-1})
$$
これを順々に進めていくと以下のように同時確率の分解公式を得ることができる
$$
\begin{aligned}
p(y_{1:T}) &= p(y_T | y_{1:T-1}) p( y_{1:T-1}) = p(y_T | y_{1:T-1}) p( y_{1:T-1} | y_{1:T-2})p( y_{1:T-2}) = \cdots \\
&= \prod^{T}_{t=1} p(y_t | y_{1:t-1})
\end{aligned}
$$
2変数の同時確率の分解
別の変数 $x_t$ が加わった場合も、上記と同じよう分解することができる。一見、難しく見えるかもしれないが、同時確率 $p(A,B)$ の分解のルールに従って一つずつ進めていけばよい。
$$
\begin{aligned}
p(x_{1:T}, y_{1:T}) &= p(y_T | y_{1:T-1}, x_{1:T}) p( y_{1:T-1}, x_{1:T} ) \\
&= p(y_T | y_{1:T-1}, x_{1:T}) p( x_T | y_{1:T-1}, x_{1:T-1} ) p( y_{1:T-1}, x_{1:T-1} ) \\
&= \cdots \\
&= \prod^T_{t=1} p(y_t | y_{1:t-1}, x_{1:t}) p( x_t | y_{1:t-1}, x_{1:t-1} )
\end{aligned}
$$
1式目は $y_T$ に注目してその他の変数を条件にしただけであり、2式目は残った変数 $y_{1:T-1}, x_{1:T}$ から $x_T$ に注目したものになる。これを繰り返していくことを考えると最後の条件付き確率の総積の形を導くことができる。
状態空間モデルの確率モデル表現
一般状態空間モデルでは、以下の2つの重要な仮定を置く。
マルコフ性（Markov Property）
状態変数 $x_t$ について、現在の状態が与えられれば、未来の状態は過去の状態や観測に依存しないという性質。すなわち
$$
p(x_t∣x_{1:t−1}, y_{1:t−1} ) = p(x_t∣x_{t−1} )
$$
条件付き独立性（Conditional Independence）
観測 $y_t$ について、現在の状態が与えられれば、観測は他の状態や観測に依存しないという性質。すなわち
$$
p(y_t | x_{1:t}, y_{1:t-1}) = p(y_t | x_t)
$$
マルコフ性が成り立つ理由
状態空間モデルでは、状態変数 $x_t$ を「過去の情報を要約したベクトル」として設計する。例えば、
- トレンドモデル: $x_t = [\mu_t, \mu_{t-1}]$ （現在と1期前のトレンド）
- 季節調整モデル: $x_t = [\mu_t, \mu_{t-1}, s_{t-1}, \ldots, s_{t-6}]$ （トレンドと過去6期分の季節成分）
このように、将来を予測するのに必要な過去の情報を状態ベクトルに組み込むことで、マルコフ性を満たすように構築されている。
状態空間モデルでの同時分布
上記の仮定により、一般の二変数の確率モデル（$p(x_{1:T}, y_{1:T})$）の分解式は下記の通り簡素に表現できる。
$$
p(x_{1:T}, y_{1:T}) = \prod^T_{t-1} p(y_t | x_{t}) p( x_t | x_{t-1} )
$$
一般状態空間モデルを振り返ると上記の $p(y_t | y_{1:t-1}, x_{1:t}), p( x_t | y_{1:t-1}, x_{1:t-1} )$ の特殊な形になっていることがわかる。
状態空間モデルの状態変数 $x_t$ は、依存関係があると仮定する過去の成分の範囲をベクトルとして取り込んでいる。この仮定があることで、 $p(y_t | x_{t})$ のように時刻 $t$の観測データは $x_t$ にのみ従うし、$p( x_t | x_{t-1} )$ のように時刻 $t$の状態は $t-1$ 時刻の状態変数 $x_{t-1}$ にのみ依存したのである。このような性質をマルコフ性と呼ぶ。
マルコフ性が成り立つように状態変数を構築することで、同時分布が計算可能になる。
この形により、複雑な時系列の依存関係を扱いやすい形で表現でき、カルマンフィルタや粒子フィルタなどの推定アルゴリズムの基礎となる。
ここまでのまとめ
本記事では、トレンドモデルと季節調整モデルを例に、状態ベクトルの設計から状態空間モデルの確率的定式化までを扱った。
本記事で重要な下記2点である。
- 状態ベクトルの設計思想: 状態変数を「過去の情報を要約したベクトル」として構築することで、複雑な時系列の依存関係をマルコフ性を満たす形で表現できる。
- 確率モデルとしての表現: この設計により、同時分布が以下の扱いやすい形に分解される
$$
p(x_{1:T}, y_{1:T}) = \prod^T_{t=1} p(y_t | x_{t}) p( x_t | x_{t-1} )
$$
次回は、この確率的表現を基にした状態空間モデルの推定・予測アルゴリズムを扱う。