R で状態空間モデル: {dlm} の対数尤度計算について

状態空間モデルを扱う R パッケージである {dlm} では、モデルのパラメータ推定を最尤法で行なう。その際に使われる対数尤度の計算式はパッケージによって違いがあるようで、推定結果が自分の参照しているテキストとずれる。その違いを確かめたい。

当該のテキストはこちら。今回はこれを持っている方向けの内容。

状態空間時系列分析入門

作者: J.J.F.コマンダー,S.J.クープマン,Jacques J.F. Commandeur,Sime Jan Koopman,和合肇
出版社/メーカー: シーエーピー出版
発売日: 2008/09
メディア: 単行本
購入: 2人クリック: 4回
この商品を含むブログを見る

参考サイト「状態空間時系列分析入門」をRで再現するはこのテキストを {dlm}, {KFAS} で再現しており、対数尤度の差異についても記載がある。

テキストでは単変量状態空間モデルの時刻 ${ t }$ における観測値 ${ y_t }$ と状態 ${ a_t }$ を以下の式であらわしている。

${ y_t = z'_t a_t + \epsilon_{t} \ \ \ \ \ \ \epsilon_t \sim NID(0, \sigma^{2}_{ \epsilon } ) }$

${ a_{ t+1 } = T_t a_t + R_t \eta_t \ \ \ \ \ \ \eta_t \sim NID(0, Q_t) }$

1期先の状態を更新するカルマンフィルタは、

${ a_{ t+1 } = T_t a_t + K_t (y_t - z'_t a_t ) }$

${ K_t = \frac{P_t}{F_t} }$

${ K_t }$ : カルマンゲイン
${ P_t }$ : 状態推定誤差分散
${ F_t }$ : カルマンフィルタでの1期先予測誤差分散

これらを用いて、対数尤度関数は式 (8.7) で以下のように定義されている。

$\log L_d = - \frac{n}{2} log(2 \pi) - \frac{1}{2} \sum_{t=d+1}^n ( \log F_t + \frac{ v_{t}^{2} }{ F_t })$

${ n }$ : 観測時系列の長さ
${ d }$ : 推定する状態数
${ v_t }$ : カルマンフィルタの1期先予測誤差

この式で計算される対数尤度は {dlm} のものと異なる。差異がある理由として記載する 2 点がわかったが、これらを解消してもまだ結果はあわない。

1. テキストの記載による差異

まず、テキストの各モデルに対数尤度として記載されている数値は上式の計算結果そのものでない。P162 に添付されているプログラムに対数尤度を系列の長さ ${ n }$ で割る処理がしれっと入っており、処理後の値が対数尤度として記載されている。そのため、対数尤度に戻すためには各モデルに記載されている数値に系列の長さ ${ n }$ をかける必要がある。

例えば 2.1 ローカルレベルモデル (確定的レベルモデル) に対数尤度として記載されている数値 0.3297597 に、系列の長さ 192 をかけた 63.31386 が式 (8.7) の計算結果にあたる。

2. `{dlm}`の対数尤度計算式との差異

テキストで利用しているデータ利用して、{dlm} との差異を確かめたい。まずデータを読み込む。データは著者のサポートサイト、An Introduction to State Space Time Series Analysis から入手できる。

ukdrivers <- read.table('UKdriversKSI.txt', skip = 1)
ukdrivers <- ts(ukdrivers[[1]], start = c(1969, 1), frequency = 12)
ukdrivers <- log(ukdrivers)

y <- ukdrivers

パラメータを受け取ってローカルレベルモデルを返す関数 buildFun を定義し、dlm::dlmMLE で最尤推定する。結果はリストになり、推定されたパラメータ par や対数尤度 (を代替するもの/後述) value がプロパティとして含まれている。

library(dlm)
parm <- c(var(y))
buildFun <- function(x) {
  dlmModPoly(order = 1, dV = parm, dW = 0)
}
dlmFit <- dlmMLE(y, parm = parm, build = buildFun)
class(dlmFit)
# [1] "list"
names(dlmFit)
# [1] "par"         "value"       "counts"      "convergence" "message"    
dlmFit$par
# [1] 0.02935256

続けて、最尤推定したパラメータ dlmFit$par を使ってモデルを作成する。プロパティの意味は vignette の記載がわかりやすい。テキストと式/添字をあわせると

FF: ${ z'_t }$ に対応。
V: ${ \sigma^{2} }$ に対応。
GG: ${ T_t }$ に対応。
W: ${ R_t Q_t }$ に対応。
m0: ${ a_1 }$ の事前分布の平均。
C0: ${ a_1 }$ の事前分布の分散。
J...: 各パラメータが時間変化するとき利用。

dlmMod <- buildFun(dlmFit$par)
class(dlmMod)
# [1] "dlm"
names(dlmMod)
#  [1] "m0"  "C0"  "FF"  "V"   "GG"  "W"   "JFF" "JV"  "JGG" "JW" 
V(dlmMod)
#            [,1]
# [1,] 0.02935256

ここで、{dlm} の内部で最尤推定がどのように行なわれているのか確かめたい。 dlm::dlmMLE の定義をみると、最尤推定は dlm::dlmLL を optim で最適化することで行っている。

dlmMLE
# function (y, parm, build, method = "L-BFGS-B", ..., debug = FALSE) 
# {
#     logLik <- function(parm, ...) {
#         mod <- build(parm, ...)
#         return(dlmLL(y = y, mod = mod, debug = debug))
#     }
#     out <- optim(parm, logLik, method = method, ...)
#     return(out)
# }
# <environment: namespace:dlm>

dlm::dlmLL はその名前から対数尤度のようにも見えるが、対数尤度とは符号が逆 ( optim で最適化するため、尤もらしいほど値が小さくなる )。

dlmLL(y = y, mod = buildFun(dlmFit$par))
# [1] -230.7721

補足パラメータ推定時の dlm::dlmLL の値は dlm::dlmMLE の結果にも含まれている。

dlmFit$value
# [1] -230.7721

dlm::dlmLL の処理は少し条件分岐が多いため省略するが、上のモデルについて必要な処理は以下のようなものであることがわかる。カルマンフィルタ dlm::dlmFilter を適用した結果から予測誤差、予測誤差標準偏差を求めて対数尤度の代替となる値を算出している。

# 予測誤差、予測誤差標準偏差を計算
residuals_dlmcompat <- function (object) {
  mod <- object$mod
  # yは観測値 (入力), a は状態の予測値
  res <- drop(object$y - object$a %*% t(mod$FF))
  f <-  function(i) {
    # 予測誤差標準偏差
    diag(crossprod(object$D.R[i, ] * t(mod$FF %*% object$U.R[[i]])) + mod$V)
  }
  SD <- drop(t(sqrt(sapply(seq(along = object$U.R), f))))
  return(list(res = res, sd = SD))
}

# dlmLLと同じ値を計算
LL_dlmcompat <- function(y, mod) {
  filtered <- dlmFilter(y, mod)
  resid <- residuals_dlmcompat(filtered) 
  # 予測誤差分散
  F <- as.numeric(resid$sd)^2
  v <- as.numeric(resid$res)
  return(0.5 * sum(log(F) + v^2/F))
}

LL_dlmcompat(y = y, mod = buildFun(dlmFit$par))
# [1] -230.7721

これを上式 (8.7) と比べると以下の差異がある。

符号が逆 ( optim での処理のため )
第1項目がない ( パラメータ非依存のため optim には影響ない )
第2項目の総和の範囲が異なる

これをテキストに準じて書き換えると、

LL_ckbookcompat <- function(y, mod) {
  n <- length(y)
  filtered <- dlmFilter(y, mod)
  resid <- residuals_dlmcompat(filtered) 
  F <- (as.numeric(resid$sd)^2)[2:n]
  v <- as.numeric(resid$res)[2:n]
  return(-0.5 * (n * log(2 * pi) + sum(log(F) + v^2/F)))
}
LL_ckbookcompat(y = y, mod = buildFun(dlmFit$par))
# [1] 62.39492

となり、テキストの対数尤度 63.31386 に結構近づいた (が、あわない、、)。対数尤度の計算結果はパッケージによって異なる、ということは覚えておいたほうがよさそうだ。

補足: このモデルにおいては、 {KFAS} の結果も上記修正後の結果と一致している。

library(KFAS)
kfasMod <- SSModel(y ~ SSMtrend(degree=1, Q=list(matrix(0))), H = matrix(NA))
kfasFit <- fitSSM(kfasMod, inits = parm, method = "BFGS")

kfasFit$optim.out$value
# [1] -62.39492
KFS(kfasFit$model)$logLik
# [1] 62.39492