前回：基礎から分かる時系列分析 4.カルマンフィルタ(ローカルトレンドモデル＋周期モデルの実装

逐次データ同化における状態の逐次推定は、次のような分布を推定することだった。

フィルタ分布：

p({\bf x}_t|{\bf y_{1:t}})=\frac{p({\bf y}_t|{\bf x}_t)p({\bf x}_t|{\bf y}_{1:t-1})}{\int p({\bf y}_t|{\bf x}_t)p({\bf x}_t|{\bf y}_{1:t-1})d{\bf x}_t}

予測分布：

p({\bf x}_{t+k}|{\bf y}_{1:t})=\int p({\bf x}_{t+k-1}|{\bf y}_{1:t})p({\bf x}_t|{\bf x}_{t+k-1})d{\bf x}_{t+k-1}

平滑化分布：

p({\bf x}_t|{\bf y}_{1:T})=p({\bf x}_t|{\bf y}_{1:T})\int\frac{p({\bf x}_{t+1}|{\bf x}_t)}{p({\bf x}_{t+1}|{\bf y}_{1:t})}p({\bf x}_{t+1}|{\bf y}_{1:T})d{\bf x}_{t+1}

いずれの分布を計算するにも積分計算を乗り越えなければならない。カルマンフィルタはすべての確率変数が正規分布に従うと仮定することによって積分を行列計算に落とし込み解析的に分布を求める、という手法であった。一方マルコフ連鎖モンテカルロ法(MCMC)は、求めたい分布からのサンプルを直接得て分布を推定する方法である。

MCMCをフィルタリングに用いると各時点でのフィルタ分布からの乱数を毎回発生させる必要があり、現実的ではない。MCMCは逐次的なアプローチよりも平滑化に導入するのが一般的である。

5.1 モデルの定義

英国の自動車運転者死傷者数のデータを用いる。

#---
#　データの準備
#---

y     = UKDriverDeaths
y_log = log(y)
 
oldpar = par(no.readonly = T) 
par(mfrow  = c(2,1)) 
par(oma    = c(0, 0, 0, 0))
par(mar    = c(4, 4, 2, 1))
par(family = "HiraginoSans-W4") 
plot(y)
plot(y_log) # 対数変換必要なさそう
par(oldpar)

ローカルレベル＋周期成分のカルマンフィルタを実行して平滑化まで行う。状態空間モデルは

\begin{cases}
{\bf a}_t=G_t{\bf a}_{t-1}+HH_t{\bf \tau}_t,\ \ \ \ \tau\sim \mathcal{N}({\bf 0},I_{12})\\
y_t=F_t{\bf a}_t+GG_t\epsilon_t,\ \ \ \ \ \ \epsilon\sim \mathcal{N}(0,1)
\end{cases}

状態ベクトルと遷移行列、システムノイズは

\begin{align}
{\bf a}=\begin{pmatrix}x_t\\ s_t\\ s_{t-1}\\ \vdots\\ s_{t-10}\end{pmatrix},\ \  G_t&=\begin{pmatrix}
1 & 0 & \cdots &&&0\\
0 & -1 & -1 &\cdots&& -1\\
0 & 1 & 0 & \cdots && 0\\
0 & 0 & 1 &\ddots&&\vdots&\\
\vdots&\vdots &&\ddots&\\
0 & 0 & \cdots &&1 &0
\end{pmatrix}\\
HH_t&=\begin{pmatrix}
W_x&0&\cdots&0\\
0 & W_s&&\vdots\\
\vdots &&&\\
0 &\cdots&&0
\end{pmatrix}
\end{align}

である。観測モデルの遷移行列とシステムノイズは

F_t=\begin{pmatrix}1 & 1 &0\cdots&0\end{pmatrix},\ \ GG_t=V

である。

#---
#　カルマンフィルタ----
#---

# モデル定義
build_dlm_UKD = function(par){
   return(
     dlmModPoly(order = 1, dW = exp(par[1]), dV = exp(par[2])) + 
     dlmModSeas(frequency = 12, dW = c(exp(par[3]), rep(0, 10)))
   )
}

# パラメータ特定
fit_dlm_UKD = dlmMLE(y = y, parm = rep(0,3), build = build_dlm_UKD)
mod         = build_dlm_UKD(fit_dlm_UKD$par)

# 平滑化
dlmSmoothed_obj = dlmSmooth(y = y, mod = mod)
x               = dropFirst(dlmSmoothed_obj$s[,1])
s               = dropFirst(dlmSmoothed_obj$s[,3])

oldpar     = par(no.readonly = T) 
par(mfrow  = c(3,1))
par(oma    = c(2, 0, 0, 0))
par(mar    = c(2, 4, 1, 1))
par(family = "HiraginoSans-W4") 
plot(y, col = 'lightgray')
lines(x, col = 'red', lty = 'dashed')
legend('topright',
  legend = c('観測値', '平均(MCMC)'),
  lty    = c('solid', 'dashed'),
  col    = c('lightgray', 'red')
)
plot(s, ylab = '周期成分')
mtext(text = 'Time', side = 1, line = 1, outer = TRUE)

5.2 MCMC

MCMCを用いて状態変数の同時事後分布からのサンプリングを行う。状態変数の同時事後分布とは、

p({\bf a}_{0:T}|y_{1:T})\propto p({\bf a}_{0:T})p(y_{1:T}|{\bf a}_{0:T})

のこと。事前分布は、

p({\bf a}_{0:T})=p({\bf a}_0)\Pi_{t=1}^Tp({\bf a}_t|{\bf a}_{t-1})

と書けるのだから、状態の事前分布を

p({\bf a}_0)=p(\begin{pmatrix}x_0\\ s_{0,1}\\ s_{0,2}\\ \vdots\\ s_{0,11}\end{pmatrix})=\mathcal{N}({\bf m}_0,C_0)

と置けば、システムモデルより次のように表現できる。

x_0は平均m_0[1]、分散C_0[1,1]の正規分布に従う。
任意のt(>=1)に対し、x[t]は平均x[t-1],分散W_xの正規分布に従う。
任意のS(>=1)に対し、s_0[S]は平均m_0[S+1],分散C_0[S+1,S+1]の正規分布に従う。
s[1]は平均が「s_0の和の符号を逆にしたもの」、分散W_sの正規分布に従う。
任意の1<t<12に対し、s[t]は平均「s[t-1]からs[1]、s_0[t]からs_0[11]までの和の符号を逆にしたもの」、分散W_sの正規分布に従う。
任意のt(>=12)に対し、s[t]は平均が「s[t-1]からs[t-11]の和の符号を逆にしたもの」、分散W_sの正規分布に従う。

尤度は

p(y_{1:T}|{\bf a}_{0:T})=\Pi_{t=1}^Tp(y_t|{\bf a}_t)

なので、観測モデルより次のように表現できる。

y[t]は平均x[t]+s[t]、分散W_sの正規分布に従う。

以上より、Stanコードは以下の通り。

// model_UKD.stan

data {
  int<lower=1>   t_max; // 時系列長
  vector[t_max]      y; // 観測値

  vector[12]        m0; // 事前分布の平均
  cov_matrix[12]    C0; // 事前分布の共分散行列
}

parameters {
  real              x0;  // 状態 (レベル成分, t = 0)
  vector[11]        s0;  // 状態 (周期成分, t = 0)
  vector[t_max]      x;  // 状態 (レベル成分, t = 1:t_max)
  vector[t_max]      s;  // 状態 (周期成分, t = 1:t_max)
  
  real<lower=0>     Wx;  // システムモデルの分散(レベル成分) 
  real<lower=0>     Ws;  // システムモデルの分散(周期成分)
  cov_matrix[1]      V;  // 観測モデルの分散
}

model {
  // 尤度
  for (t in 1:t_max){
    y[t] ~ normal(x[t]+s[t], sqrt(V[1,1]));
  }
  
  // 事前分布 (レベル成分)
  x0 ~ normal(m0[1], sqrt(C0[1,1]));
  x[1] ~ normal(x0, sqrt(Wx));
  for (t in 2:t_max){
    x[t] ~ normal(x[t-1], sqrt(Wx));
  }
  
  // 事前分布 (周期成分)
  for (S in 1:11){
    s0[S] ~ normal(m0[S+1], sqrt(C0[(S+1),(S+1)]));
  }
  
  s[1] ~ normal(-sum(s0[1:11]), sqrt(Ws));
  
  for (t in 2:11){
    s[t] ~ normal(-sum(s0[t:11])-sum(s[1:(t-1)]), sqrt(Ws));
  }
  for (t in 12:t_max){
    s[t] ~ normal(-sum(s[(t-11):(t-1)]), sqrt(Ws));
  }
}

'model_UKD.stan'を使ってRStanでMCMCを行う。

#---
#　Rstanの準備----
#---

# install.packages("rstan")
library(rstan)
rstan_options(auto_write=TRUE)                     # 呪文
options(mc.cores=parallel::detectCores())          # 呪文
theme_set(theme_get() + theme(aspect.ratio = 3/4)) # 呪文

#---
#　MCMC
#---

set.seed(123)
t_max = length(y)

# モデル
stan_model_out = stan_model(file = 'model_UKD.stan')

# 平滑化
fit_stan = sampling(
  object = stan_model_out,
  data   = list(t_max = t_max, y = y, m0 = mod$m0, C0 = mod$C0),
  pars   = c('x', 's', 'Wx', 'Ws', 'V'),
  seed   = 123
)

traceplot(fit_stan, pars = c('Wx', 'Ws', 'V'), alpha = 0.5)


# サンプリング結果を取り出す(呪文)
stan_mcmc_out_x = rstan::extract(fit_stan, pars = 'x')
stan_mcmc_out_s = rstan::extract(fit_stan, pars = 's')

# 周辺化、平均を求める
s_mcmc_x     = colMeans(stan_mcmc_out_x$x)
s_mcmc_x     = ts(s_mcmc_x, start = c(1969,1), frequency = 12)
s_mcmc_s     = colMeans(stan_mcmc_out_s$s)
s_mcmc_s     = ts(s_mcmc_s, start = c(1969,1), frequency = 12)

oldpar = par(no.readonly = T)
par(mfrow = c(2,1))
par(family = "HiraginoSans-W4") 
plot(y, col = 'lightgray')
lines(s_mcmc_x, col = 'red', lty = 'dashed')
legend('topright',
  legend = c('観測値', '平均(MCMC)'),
  lty    = c('solid', 'dashed'),
  col    = c('lightgray', 'red')
)
plot(s_mcmc_s, ylab = '周期成分')
mtext(text = 'Time', side = 1, line = 1, outer = TRUE)
par(oldpar)

周期成分のシステムノイズがうまく収束できていない。Rhatも値が悪くなっている。