ill-identified diary

所属組織の見解などとは一切関係なく小難しい話しかしません

[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年ごろに大きく改修され, 現在のように使いやすい構文になったという.

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

技術的に細かいことを書くと, 非ガウシアンの線形状態空間モデルであっても単に未知パラメータを推定するだけなら, 多くの場合は時系列に特有のテクニックは必要ない. 一方でフィルタリングするには, 平均値まわりで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 はモデルの対数尤度に等しい.

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

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

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

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

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

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

https://catalog.data.gov/dataset/allegheny-county-crash-data

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

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

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

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

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

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

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

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

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

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

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

2018-08-11 追記: 結果が不自然なのはパラメータの設定ミスという指摘を受け、結果とその説明を修正

f:id:ill-identified:20180811000928p:plain
各モデルの累積絶対誤差

まず, ガウシアンとポアソンとでは大きく差が出ており, さらにガウシアン/ポアソンの回帰なしモデルどうしと, 回帰ありのモデルどうしが, シフトしたかのように同様な誤差の累積をしている. モデルの分布を変えるだけで, 誤差が大幅に低下することがわかる. さらに回帰変数を追加したモデルは, ないモデルよりも誤差の累積がゆるやかである.

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

ARMA モデル


gist54dbf66656d24e963b530c2df8f78cdc


ポアソン モデル


gist9d994625ef9d3f2d843567d8af6a4bae


参考文献


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

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