時系列データ解析における状態空間モデルとは

状態空間モデルとは

直接観測できない潜在的な状態の存在を仮定したモデルを状態空間モデルといいます。

状態空間モデルは、次式のように観測方程式と状態方程式によって表されます。

観測方程式

$y_t = h({\bf \alpha}_{t}) + \varepsilon_t$

状態方程式

${\bf \alpha}_t =f({\bf \alpha}_{t-1}) + g({\bf \xi}_t)$

観測方程式から見ていきます。 観測データ $yt$ は、同時点の観測できない状態 $αt$ の関数 $h(αt)$ と観測ノイズ $εt$ の加算で表されます。 状態 $αt$ は任意次元のベクトルで、様々な要素を状態として考慮することができます。

次に状態方程式についてですが、 $αt$ は1時点前の状態 $αt−1$の関数 $f(αt−1)$に状態ノイズの関数 $g(ξt)$が加わる形で表されます。

つまり、状態 $αt$はそれ自体で時系列的に推移していくのに対し、観測 $yt$は状態の関数にノイズが加わることで得られると仮定しています。

関数 $f,g,h$に関しては何も記述していないため、使用するパラメータが推定しうるのであればどのようなものでも構いません。

設計者のアイデアや知見を柔軟に反映することができるため、様々なことができます。 一方で柔軟さの反面、パラメータが本当に推定しうるのか、モデルが妥当であるのか、推定法はどのようなものが良いか、などを把握しコントロールしなくてはならないという難点もあります。

目に見えない「状態」の存在を仮定し観測と分けることのメリットの例として、以下のような点があります

• 適切な状態空間モデルを用いることで、観測できないが本来知りたい状態を推定できる。また、これを用いた現場へのフィードバックが可能。
• 観測ノイズと状態ノイズを分離できる。つまり誤差のうちどの程度が観測によるものか、状態の推移によるものかを分けて考えることができる。
• 欠損値があっても、その時点での状態を推定することで補完ができるため、欠損値の処理に敏感にならなくて良い。

状態空間モデルの目的と推定法

状態空間モデルにおける重要な課題は、データから状態 $αt$を推定することです。

$αt$の推定には、用いるデータ${y1,…,yn}$との関係により以下の3種類があります。

• $n \lt t$：予測（最新のデータよりも先の状態を推定する問題）
• $n=t$：フィルタリング（最新のデータと同時点の状態を推定する問題）
• $n \gt t$：平滑化（最新のデータまでに基づいてそれ以前の状態を推定する問題）

フィルタリングは、予測の「ズレ」を最新のデータを使って修正する作業といえます。次々と送られてくるデータからリアルタイムに位置を推定することは、フィルタリングの具体例といえます。

平滑化は、全てのデータを使った状態の修正であり、過去のデータのみを使うよりも正確な推定ができます。 キャンペーンの効果を振り返るときなど、フィードバックや知見獲得を目的として行われることが多いです。

では、これらの推定作業はどのような方法で行われるのでしょうか。 以下では3種類の状態の推定方法を紹介します。

カルマンフィルタ (Kalman filter)

観測ノイズ $εt$ と状態ノイズ $ξt$ に正規分布を仮定し、関数 $f,g,h$に線形関数を仮定した状態空間モデル（線形ガウス状態空間モデル）は、カルマンフィルタというアルゴリズムで効率的に推定することができます。

カルマンフィルタは逐次的なアルゴリズムで、 $t$ 時点での状態 $αt$の推定に $t−1$ までの値を使う一期先の予測と、 $t$ の値も使って修正を行うフィルタリングを交互に繰り返し、最後に平滑化を行います。

正規性と線形性を仮定するのは、正規分布に従う確率変数同士の線型結合により与えられる変数もまた正規分布に従う、という性質を利用するからです。

平均と共分散さえ推定すれば同時分布が定まるので、計算効率が良くなります。 逆にとらえれば、予測値の分布（予測分布）もフィルタリングにより推定する状態の分布（フィルタ分布）も正規分布になることを仮定しています。

粒子フィルタ (Particle filter)

粒子フィルタは、カルマンフィルタよりも計算量が多くなるものの、正規性と線型性を仮定しなくても良い推定アルゴリズムです。

非正規性と非線型性を許すことで、複雑な予測分布やフィルタ分布を推定することができます。

粒子フィルタもカルマンフィルタと同様に予測とフィルタリングを繰り返す逐次的なアルゴリズムですが、複雑な分布を重み付きの乱数（粒子）で近似します。 このような近似をモンテカルロ近似といいます。

マルコフ連鎖モンテカルロ法 (Markov Chain Monte Carlo methods: MCMC)

MCMCは、時系列モデルに限らず、複雑な統計モデルのパラメータ推定に用いられます。

粒子フィルタと同様、分布をたくさんの乱数で近似しますが、粒子フィルタが逐次的に予測 → フィルタリング を繰り返すのに対し、MCMCは初めから全時点のデータを使って分布を近似します。そのため平滑化分布のみが計算され、予測分布を得るには追加の計算が必要になります。

状態空間モデルの例

時系列データのモデリングで広く利用されている状態空間モデルの例として、ローカルレベルモデルと構造時系列モデルについて紹介します。

ローカルレベルモデル

最もシンプルな状態空間モデルの例として、ローカルレベルモデルを紹介します。

ローカルレベルモデルは次式で表されます。

$\begin{cases}y_t = \alpha_t + \varepsilon_t,~ &\varepsilon_t\sim N(0, \sigma^2_\varepsilon),~ t=1,\dots,T\\\alpha_t = \alpha_{t-1} + \xi_t,~ &\xi_t\sim N(0, \sigma^2_\xi),~ t=1,\dots,T\end{cases}$

$yt$ は観測値で、 $αt$ は一変数の状態です。 観測方程式では状態 $αt$ にホワイトノイズ $εt$ が加算されることで $yt$ が得られることが表され、状態方程式では状態 $αt$ が時間とともに $ξt$ の変動をともなって推移していくことが表されています。

ここでパラメータが $σ2ε$ と $σ2ξ$ だけであることに着目すると、ローカルレベルモデルは「観測の誤差による変動の大きさと真の水準が推移する変動の大きさに分けるモデル」であると解釈できます。

たとえば、毎日体重を測定して得られたデータがあるとしましょう。 体重の真の水準（状態）は毎日の食事や運動などによって推移します。 一方で、体重計や身につけている衣服・気圧などに依存する測定の誤差も存在します。 ここでローカルレベルモデルを利用すれば真の体重の推移が推定でき、当てはまりが良くなるように真の体重推移の分散 $σ2ε$ と測定誤差の分散 $σ2ξ$を推定できます。

ローカルレベルモデルの他にも、様々なトレンドに関するモデリングが提案されています。 たとえば、トレンドの傾きは直前の傾きと似ているはずだ、というアイデアを盛り込んだローカル線形トレンドモデルなどがあります。

構造時系列モデル

状態空間モデルの枠組みが活用される典型例に構造時系列モデルがあります。 構造時系列モデルは、注目している時系列を、直接解釈可能な複数の時系列の和で表すモデルで、任意の時系列の足し合わせであるため解釈性が高いというメリットがあります。

なかでも、基本的な構造時系列モデルは次の構造で表されます。

時系列データ = トレンド + 周期成分 + ホワイトノイズ

構造時系列モデルを状態空間モデルの枠組みで扱う際は、トレンドや周期の各要素それぞれが状態方程式に従うとし、その状態 $α(1)t,α(2)t,…,α(n)t$を観測方程式において足し合わせます。 トレンドの状態方程式にはローカルレベルモデルやローカル線形トレンドモデルなど、様々なトレンドモデルを用いることができます。

周期成分のモデリングでは、ダミー変数を使うものと三角関数を使うものが一般的です。

まずはダミー変数を使った周期のモデリングについて説明します。

例えば1週間の周期を考えた時、ダミー変数によるモデリングは次式のように行われます。 ここで、各$β$は回帰係数で、$di(t)$は $t=i$であるとき1、それ以外のとき$0$をとるダミー変数です。

$\alpha_{t}^{weekly} = \beta_{月}d_{月}(t) + \beta_{火}d_{火}(t) + \cdots+\beta_{土}d_{土}(t) +\varepsilon_t$

ここでは、日曜日を基準として他の曜日にダミー変数を用意しています。 そのため、同時に推定されるトレンド $αtrendt$は日曜日に対するもので、たとえば火曜日であれば対応する係数$β$火 が加算されます。

続いて、三角関数を使った周期のモデリングについて説明します。

余弦関数 $cos(t)$と正弦関数 $sin(t)$はどちらも周期 $2π$、振幅1です。 異なる周期と振幅を持つ複数の三角関数を足し合わせることで、複雑な周期を近似することができます。

三角関数を足し合わせた関数は滑らかであるため、非整数な周期のモデリングや一年間の周期のモデリングに適しています。 一年の周期でダミー変数を364個用意することは現実的でない他、うるう年を考慮した周期は365.25日になりダミー変数で表現することは困難です。