R‎ > ‎

附録 4 状態空間モデル

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

状態空間モデルの定義
Yt を l 変量の時系列としたとき、これを以下の2つの方程式によって表現するモデルのことを状態空間モデルと呼ぶ。
Xt = FtXt-1 + Gtvt        … システム方程式
Yt = HtXt + wt             … 観測方程式
システム方程式は、システムの状態が1つ前の時点からどのように更新されるかを示しており、観測方程式は、システムの状態から観測値がどのように取り出されるかを示している。
Xt は「状態ベクトル」と呼ばれる。ここでは k 変量のベクトルとする (k はモデルのパラメータである)。Ft は回帰係数行列で k 行 k 列の行列である。システムの過去の状態が現在にどのように影響を与えるかを示している。vt はシステムノイズと呼ばれ、システムの状態の変化の際に入り込む攪乱を表す。m 変量の正規ホワイトノイズで、分散共分散行列 Qt に従うものとする。Gt はシステムノイズがシステムの状態に与える影響を表す k 行 m 列の行列である。
Ht はシステムの状態が観測値にどのように射影されるかを表す行列で、l 行 k 列の行列である。wt は観測ノイズと呼ばれ、システムの状態を観測する際に入り込む攪乱を表す。l 変量の正規ホワイトノイズで、分散共分散行列 Rt に従うものとする。
上記のような条件を置いた状態空間モデルは、ガウス型状態空間モデルとも呼ばれる。
状態空間モデルのシステム方程式、観測方程式は、真のモデルに対して一意ではない。

AR モデルの状態空間表現
状態空間モデルの例として、AR(p) モデルの状態空間表現を見てみることにする。
以下の式で定義される yt は、AR(p) モデルである。
yt = φ1yt-1 φ2yt-2 + ... + φpyt-p + εtεt ~ iidN(σ2)
定数項は 0 とし、攪乱項の分布は正規ホワイトノイズとする。xt = (yt, yt-1, ..., yt-p+1)' とおくと、以下の関係式が成り立つ。
xt = Fxt-1 + Gεt
ただし、
F = { (φ1φ2, ...φp-1φp), (1, 0, ..., 0, 0), (0, 1, ..., 0, 0), ..., (0, 0, ..., 1, 0) } (p 行 p 列の行列)
G = (1, 0, ..., 0)' (p 次元のベクトル)
である。これが状態空間モデルにおけるシステム方程式にあたる。システムノイズの分散共分散行列 Q は Q = σ2 になる。
そして、H = (1, 0, ..., 0)' (p 次元のベクトル) と置けば、
yt = HXt
となり、これが観測方程式にあたる。観測ノイズの分散共分散行列 R は R = 0 となる。F, G, H, Q, R のいずれも時不変なので、時不変の状態空間モデルである。

状態の推定
時系列データ Yt = (Y1, Y2, ..., YT) に基づいて状態空間モデルを構築するにあたっては、状態 Xt' = (X1, X2, ..., XT) を推定する必要がある。また状態空間モデルから将来の観測値を予測するには状態 Xt' = (XT+1, XT+2, ...) の推定も必要である。
Xt' (1 ≦ t' < T) の推定を平滑化Xt' (t' = T) の推定をフィルタXt' (t' > T) の推定を予測と呼んでいる。

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

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

非ガウス型状態空間モデル
ガウス型状態空間モデルでは、ノイズが正規分布に従うため、平均値から離れた値を取る可能性が低く、状態や観測値の変位が滑らかになる。従って、観測値の変位が急激な場合、あまり適合が良くない。
そこで、ノイズの分布を正規分布と仮定せず、一般のホワイトノイズであるとの仮定に緩めたものが、非ガウス型状態空間モデルである。ノイズにファットテールな分布を持ってくれば、観測値の急激な変化もモデルに織り込むことができる可能性が出てくる。
定式化すると、以下のようになる。
Xt = FtXt-1 + Gtvt        … システム方程式
Yt = HtXt + wt            … 観測方程式
方程式の形状はガウス型と同じで、各係数行列の条件も同じであるが、vt, wt はそれぞれ、密度関数 q(v), r(w) に従うホワイトノイズであるとする。

推定は、カルマンフィルタの拡張によって行うことができる。ガウス型では条件付き同時分布が正規分布であるという条件を置けたが、非ガウス型ではそうとは限らないため、条件付き同時分布を直接求める必要がある。
1期先予測は、以下のようになる。
p(Xt'|Yt-1) = ∫ p(Xt'|Xt'-1p(Xt'-1|Y1, ..., Yt-1) dXt'-1
フィルタは、以下のようになる。
p(Xt'|Yt) = (p(Yt'|Xt'p(Xt'|Y1, ..., Yt-1)) / p(Yt'|Y1, ..., Yt-1)

一般化状態空間モデル
非ガウス型状態空間モデルでは、状態の推移の線形性を仮定していたが、その仮定を緩め、非線形な推移も許したものが一般化状態空間モデルである。
例えば、確率的ボラティリティモデルのパラメータ推定を行うために、非線形変換を加えた上で、線形のガウス型状態空間モデルで近似して推定を行う方法があるが、一般化状態空間モデルでは確率的ボラティリティモデルを近似無しに表現することができる。
一般化状態空間モデルは、以下のように定式化される。
Xt = F(Xt-1, vt)           … システム方程式
Yt = H(Xt, wt)            … 観測方程式
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) が存在することを要求する。また、初期状態 X0 は密度関数 p0(X) に従うものとする。
モデルの推定は解析的に行うことは一般にはできず、モンテカルロ・フィルタと呼ばれるアルゴリズムによって行う。モンテカルロ・フィルタでは、仮定した分布に従う「粒子」を乱数で生成し、解析的に表現できない同時分布をこの乱数に基づく経験分布関数で近似するという方法をとる。概要としては以下の通り。
  1. 用いる「粒子」の数を m とし、以下、j = 1, ..., m とする
  2. p0(X) に従う m 個の乱数ベクトル f0(j) を生成する (初期状態を乱数で設定)
  3. t = 1, ..., T について、以下の手順を実行する
    1. q(v) に従う m 個の乱数ベクトル vt(j) を生成する (システムノイズを乱数で生成)
    2. 各 j について、pn(j) = F(fn-1(j), vn(j)) を計算する (新しい状態の候補値の算出)
    3. 各 j について、αn(j) = r(G(Yn, pn(j))) |∂G/∂Yn| を計算する (状態 pn(j) が観測値 Yn となる確率の算出)
    4. pn(j) の中から αn(j) に比例する確率で m 回サンプリングを行い、その値を fn(1), ..., fn(m) とする (新しい状態を決定)
サンプリングは、 ランダムサンプリングでも良いが、計算量を減らすため階層化リサンプリング法といったものを用いることもできる。


Comments