[R] 非ガウシアン状態空間対応パッケージ, KFAS の使い方

概要

  1. まだ日本語情報の少ない KFAS を一連の状態空間モデルネタの続きとして紹介する.
  2. KFAS には一番良く使われている dlm パッケージよりも優れた点がいくつもある.
  3. 前回のように, パッケージの理念・構文・具体例を用いた実験を順に紹介していく.

状態空間モデルを扱う Rパッケージの中では dlm が最も有名だが, これは名前の示すように動的線形モデル dynamic linear model, すなわちノイズが正規分布になる, ガウシアン線形状態空間モデルしか扱うことができない.

これに対して KFAS の長所はいくつもあり, 特に正規分布いがいの分布も扱うことができるという点は特筆すべきである. なお, KFAS を日本語で紹介している文献は, 伊東先生の発表スライド,

伊東 (2017, R による状態空間モデリング), そして 野村 (2016, カルマンフィルタ–Rを使った時系列予測と状態空間モデル–) がある.

今回の記事の構成は次のようになる: まず, 次の第2セクションでは KFAS の特徴を俯瞰する. 第3セクションではKFAS の構文を説明する. 第4セクションでは具体的に状態空間モデルを用いて実験する.

KFAS

Helske (2017, KFAS : Exponential Family State Space Models in R) によれば, KFAS は Kalman filtering and smoothing の略であり, カルマンフィルタと平滑化の機能を持つ. KFAS は 2009年ごろからあるパッケージだが, 2011年ごろに大きく改修され, 現在のように使いやすい構文になったという.

最も大きな長所は, ガウシアンだけでなく, 二項分布, ポアソン分布, ガンマ分布などを伴う非ガウシアン線形状態空間モデルにも対応できる. たとえば, 人口とか生物の個体数など, カウントデータの時系列を扱う場合は, 二項分布やポアソン分布を利用したほうがよりフィットしやすいモデルを作ることができるだろう. 実際, 先に挙げた伊東先生のスライドでは「ハヤブサのつがいの数」をポアソン状態空間モデルで推定するモデルの例を挙げている.

技術的に細かいことを書くと, 非ガウシアンの線形状態空間モデルであっても, 未知パラメータの推定は, 多くの場合は一般化線形モデル (GLM) と同じ要領で推定できる. 一方で, フィルタリングするには, 平均値まわりで1次のテーラー展開による近似をおこなう拡張カルマンフィルタ (extended Kalman filter) や, 平均だけでなく高次のモーメント (=標準偏差) でも近似を行うことでより精度を高めた無香カルマンフィルタ (unscented —), あるいはモンテカルロ法を利用した粒子フィルタ (particle —) といったテクニカルな方法が要求される. KFAS においては, ラプラス近似と重点サンプリング (importance sampling) を利用しているらしい.

2つ目の長所は, 先に書いたように, 構文の簡単さである. dlm パッケージと同様, 状態空間モデルに含まれるローカルレベル, ローカルトレンド, ARIMA項などを足していくことで状態空間モデルを記述することができる.

3つ目の長所には, 実行速度の速さが挙げられる. Tusell (2011, Kalman Filtering in R) の実験によると, 計算も, 他のカルマンフィルタを実装したパッケージと比較すると, FKF パッケージを除き, 高速である. 最速の FKF パッケージがサポートするのは, 動的線形モデルの1期先予測の機能のみであり, またモデルの記述もパラメータの行列でしか受け付けないことを考えると, KFAS が使い勝手のよいパッケージであることがわかる.

KFAS の構文

dlm パッケージでは dlmMod オブジェクトに, dlmMod で始まる各種関数を使い状態空間モデルの部品を追加したのと同様に, KFAS では formula オブジェクトの構文を使って同様に状態空間モデルの部品を追加していく. たとえばローカルレベルモデルならば, SSMtrend() 関数を使う.

model <- SSModel(y~SSMtrend(1, Q=NA), H=NA)

のようになる. 状態空間モデルの代入された左辺の modelSSModel というクラスのオブジェクトになる. また, ローカル線形トレンドモデルならば, SSMtrend() の引数を変更し,

model <- SSModel(y~SSMtrend(2, Q=list(NA, NA), H=NA)

とする. ARIMA モデルならば,

model <- SSModel(y~SSMarima(ar=c(0, 0), ma=.3, Q=NA), H=0)

というふうに設定する.

ここで, Q, H はそれぞれ, 遷移方程式 (状態方程式) と観測方程式に対応するノイズの分散パラメータである. この分散の大きさが不明で, 推定する必要がある場合, NA を設定する. 上の例では, ローカル線形トレンドモデルは2つの状態方程式で記述されるため, Q にも2つの要素を与える必要がある. また, ARIMA モデルの係数などは NA を設定できないので, ここでは一旦適当な値を与えておく. ARIMAモデルの係数を未知としたい場合の方法は後述する.

SSModel オブジェクト内には, Z とか T とかいう名前の要素が存在する. これは状態空間モデルの行列パラメータであり, KFAS においては

\begin{align}
y_{t}= & \mathbf{Z}\boldsymbol{\alpha}_{t}+\varepsilon_{t},\,\varepsilon_{t}\sim\mathcal{N}(0,H)\\
\boldsymbol{\alpha}_{t+1}= & \mathbf{T}\boldsymbol{\alpha}_{t}+\mathbf{R}\boldsymbol{\eta}_{t},\,\boldsymbol{\eta}_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{Q})\end{align}

に対応しているため, 記述した状態空間モデルが自分の意図したものであるかどうかを確認することができる. さらに, 各 SMM*() 関数には, state_name という引数があり, これで各状態変数にわかりやすい名前をつけることができる.

加えて, SSMCustom() 関数は, これらの行列を直接与えて SSModel オブジェクトを作成することができる.

モデルを記述したら, 未知パラメータを fitSSM() 関数で推定する. 最も簡単な場合では, 作成した SSModel オブジェクトと, 初期値を与えれば計算できる. しかし, 面倒なことに, \mathbf{T}\mathbf{Z} などの分散以外のパラメータを推定したい場合, 推定に用いる関数を新たに定義する必要がある. この方法は後述する.

実験

次の2つのモデルを計算するプログラムを用意した.

  1. 線形トレンド (ドリフト) 付きのARMA(2, 1)

  2. ポアソン線形状態空間モデル

まずは簡単な例として ARMA モデルを計算してみる. しかし単なる ARMA モデルの当てはめは, はわざわざ状態空間モデルを使う意味がないので, もう1つの例として, KFAS の長所を活かせるようにポアソン線形状態空間モデルを作ることにした.

ARMA (2, 1) with drift

Helske (2017, KFAS : Exponential Family State Space Models in R) にあるサンプルプログラムでは, 確定的線形トレンド項 (ドリフト項) を表現するため, 整数列を作成し, 説明変数として与えていたが, これでは予測値を出す際に新たに, 将来のデータを入力した SSModel を作成する必要があり, すこし面倒である. そこで, 代わりに SSMtrend(2, Q=list(0, 0)) を追加する. 遷移方程式の分散がゼロなので, 毎期のトレンドは固定される. そして対応する \mathbf{T} の係数を未知パラメータとすれば, 確定的線形トレンド項の推定になる……はず.

では未知パラメータを推定する際のプログラムを見ていく. 状態空間モデルのパラメータは推定は fitSSM() でおこなう. 今回は分散パラメータ以外にも未知のパラメータがあるため, update.fn() という関数を定義する必要がある. この関数内で SSmodel を定義し, パラメータが未知の部分に引数を与える. 分散パラメータなどは必ず正でなければならないので, 指数変換する. また, 自己回帰モデルの係数は定常性を満たす必要があるが, artransform() 関数を使うことで簡単にこの制約をあたえられる. ちなみに, update.fn() を定義した場合は, モデルのパラメータが上書きされるので, 未知の分散パラメータに NA を与えても与えなくてもどっちでもいいが, どれが未知なのかわかりやすくするために NA と明示したほうがよいだろう.

update.fn <- function(pars, model) {
  model <- SSModel(y ~ SSMtrend(2, Q=list(0, 0)) +
                     SSMarima(ar =artransform(pars[1:2]), ma=pars[3], Q = exp(pars[4])), H = 0)
  model$T["slope", "slope", 1] <- pars[5]
  return(model)
}

次に, fitSSM()最尤推定を行う. この関数に, SSModel と先ほど定義した update.fn(), そして初期値を与える. fitSSM() optim() 関数を用いているので, オプション引数もこれに準じたものが使用できる. 今回のケースでは局所解に陥りやすいので, 最適化アルゴリズム"SANN" を指定した. よって計算にやや時間がかかる.

戻り値のうち optim.out に含まれているものも optim() のものと同じなので, fitmodel$optim.out$convergence がゼロならば計算が収束したと判定できる. また, -fitmodel$optim.out$value はモデルの対数尤度に等しい.

戻り値のうち, model は, SSModel オブジェクトであり, ここで推定されたパラメータが与えられている. これを用いてフィルタリングなり, 平滑化なりを行える. しかしこのモデルでは, 観測ノイズが存在しないので, 平滑化する意味はないから, 省略する. 20期間の予測を行った結果が以下である. 青い領域は95%信頼区間である. 信頼区間の幅は, predict() 関数で指定できる.

f:id:ill-identified:20170924213531p:plain

ポアソン線形状態空間モデル

非ガウシアンの線形状態空間モデルが KFAS の本領である. 正規分布以外をつかうモデルとして, ポアソン分布の応用例がないか考えてみる. ポアソン分布は件数を表すカウントデータでよく用いられる. たとえば, 先に述べたような動物の個体数, もちろん人間でもいい. あるいは, 何らかのイベントの発生件数. 毎月の事故の発生件数とか, そういったものが考えられる. 前回 bsts パッケージの解説記事

で用いた失業者数の時系列データは, カウントデータと思いきや, 件数を正規化しており, カウントデータにはならない. 原系列に復元するのも大変であるから, 別のデータセットを使うことにした.

使用したのは米国ペンシルバニア州の州都ピッツバーグを含む, アレゲニー郡の交通事故のデータである.

Allegheny County Crash Data - Data.gov

このデータは2004年から2016年までの期間をカバーする. 個別の交通事故が発生時刻まで個別に記録されているが, 今回は月次発生件数データに加工する*1. さらに, 回帰モデルにするための説明変数が欲しい. 事故の予測には, 1日あたりの平均交通量がよく使われるようだ.

しかし, 少し検索してみたところ, 各国の政府自治体が公開している細かいデータはなく, たいていは1年間の平均値などであり, 時系列データとしては不十分である*2. 予測モデルの実用的な観点からは, 事故の後に観測されるような情報は使いたくない. しかし, データを探すのに時間がかかりそうなので, この事故データだけを使うことにする. 具体的には, 死傷者の数といった事後的な情報ではなく, 天候, 道路状況, 飲酒の有無, 運転手の年齢といった, 事故以前に確定しているような情報を説明変数として使うことにする. 使えるかわからない変数が沢山あるので, bsts 使ったほうが良いような気がするが, 今回は KFAS の使用例なので仕方ない.

事故件数を月ごとに集計して件数をプロットしたものが以下になる.

f:id:ill-identified:20170924212938p:plain

毎年, 1月と12月に事故件数が多い. 年末で交通量が増えるからとか, 冬は雪で路面状況が悪いからとか, いろいろ憶測できる.

季節性があるのは明らかなので, モデルは以下の4種類を用意した.

  1. ローカルレベル + 月次トレンド + ガウシアン

  2. ローカルレベル + 月次トレンド + 回帰 + ガウシアン

  3. ローカルレベル + 月次トレンド + ポアソン

  4. ローカルレベル + 月次トレンド + 回帰 + ポアソン

bsts のように予測誤差の累積値をプロットしたものが以下のグラフである.

f:id:ill-identified:20170924214342p:plain

ガウシアンとポアソンとでは大きく差が出ているが, 説明変数で回帰したモデルは, そうでないモデルとあまり差がない. 特にポアソン分布の場合は, 説明変数があるほうが誤差が大きくなっている. やはりいい加減に作成した説明変数はあまり使えないのかもしれない. モデルにもう少し改良の余地があるが, それはまた別の機会にやることにする.

以下が今回のプログラム. パラメータ推定に結構時間がかかることに注意. 読み込むファイル名はダウンロード時のものをそのまま使った.

ARMA モデル


gist54dbf66656d24e963b530c2df8f78cdc


ポアソン モデル


gist9d994625ef9d3f2d843567d8af6a4bae


参考文献


*1:生データには発生場所まで細かく記録されているので, 交通事故防止の施策とか, そういう実用面での観点からも, 単なる時系列でなく空間モデルを適用したほうが良いかもしれない.

*2:もう少しWebスクレイピングの勉強もしようかなと思う.