「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。データ分析・語学に力点を置いています。 →現在、書くことのスケジューリングと書き溜め中です。気長にお待ちください(※2023年11月5日改訂)。

MENU

時系列解析の基礎(24/XX)

 以下の書籍

を中心に時系列解析を勉強していきます。

13. 線形・Gauss型状態空間モデルの逐次解法

13.1. カルマンフィルタ

 線形・\mathrm{Gauss}型状態空間モデルにおいて、推定すべきデータの真の値と点推定値の間の平均二乗誤差を最小にする意味で最適な逐次推定法はカルマンフィルタと呼ばれる。カルマンフィルタは非定常過程でも扱うことができる。カルマンフィルタで行なわれる計算はすべて線形計算で、正規分布の再生性から関心の対象となる分布はすべて正規分布である。

13.1.3 カルマン平滑化

 時点Tまでのカルマンフィルタリングが一旦完了しているとして、カルマン平滑化を考える。

 線形\mathrm{Gauss}形状態空間モデルでは平滑化分布も正規分布になる。そこで平滑化分布を


\begin{aligned}
\mathcal{N}\left(\boldsymbol{s}_t,\boldsymbol{S}_t\right)
\end{aligned}

とおく。
 時点t+1までの平滑化分布を元に時点tでの平滑化分布を求める手続きは以下のとおりである:

    0. 時点t+1での平滑化分布パラメータとして

平均ベクトル\boldsymbol{s}_{t+1},分散共分散行列\boldsymbol{S}_{t+1}を設定する。
    1. 時点tでの更新を以下で行う:

      平滑化利得

  \boldsymbol{A}_t\leftarrow\boldsymbol{C}_t{}^{t}\boldsymbol{G}_{t+1}\boldsymbol{R}_{t+1}^{-1}
      状態の更新

(平均)  \boldsymbol{S}_t\leftarrow\boldsymbol{C}_t+\boldsymbol{A}_t\left(\boldsymbol{S}_{t+1}-\boldsymbol{R}_{t+1}\right){}^{t}\boldsymbol{A}_t
    2. 時点tでの平滑化分布パラメータとして

平均ベクトル\boldsymbol{s}_{t},分散共分散行列\boldsymbol{S}_{t}とする。

このアルゴリズムt=T-1から時間を遡る方向に適用することで、全時点の平滑化分布を求めることができる。

  • 平滑化分布:フィルタリング分布を補正

 上記の平滑化アルゴリズムは、有名な\mathrm{RTS}アルゴリズムである。

##########################
### 自作カルマン平滑化 ###
##########################

library("ggplot2")

### カルマンフィルタリングが完了していることが前提
if(T){
  # データ
  y <- Nile
  t_max <- length(y)
  
  # 1時点分のカルマンフィルタリングを行なう関数
  Kalman_filtering <- function(m_t_minus_1, C_t_minus_1, t){
    # 一期先予測分布
    a_t <- G_t %*% m_t_minus_1
    R_t <- G_t %*% C_t_minus_1 %*% t(G_t) + W_t
    
    # 一期先予測尤度
    f_t <- F_t %*% a_t
    Q_t <- F_t %*% R_t %*% t(F_t) + V_t
    
    # カルマン利得
    K_t <- R_t %*% t(F_t) %*% solve(Q_t)
    
    # 状態の更新
    m_t <- a_t + K_t %*% (y[t]-f_t)
    C_t <- (diag(nrow(R_t)) - K_t %*% F_t) %*% R_t
    
    # フィルタリング分布
    list(m = m_t, C = C_t,
         a = a_t, R = R_t)
  }
  
  # 線形ガウス型状態空間のパラメータを設定
  G_t <- matrix(1, ncol = 1, nrow = 1)
  W_t <- matrix(exp(7.29), ncol = 1, nrow = 1)
  F_t <- matrix(1, ncol = 1, nrow = 1)
  V_t <- matrix(exp(9.62), ncol = 1, nrow = 1)
  m0 <- matrix(0, ncol = 1, nrow = 1)
  C0 <- matrix(1e+7, ncol = 1, nrow = 1)
  
  # フィルタリング分布
  # 状態の領域を確保
  m <- rep(NA_real_, t_max)
  C <- rep(NA_real_, t_max)
  a <- rep(NA_real_, t_max)
  R <- rep(NA_real_, t_max)
  
  # 時点t=1
  KF <- Kalman_filtering(m_t_minus_1 = m0,C_t_minus_1 = C0,t = 1)
  m[1] <- KF$m
  C[1] <- KF$C
  a[1] <- KF$a
  R[1] <- KF$R
  
  #時点t=2-t_max
  for (t in 2:t_max){
    KF <- Kalman_filtering(m_t_minus_1 = m[t-1],C_t_minus_1 = C[t-1],t = t)
    m[t] <- KF$m
    C[t] <- KF$C
    a[t] <- KF$a
    R[t] <- KF$R
  }
  
  # フィルタリングの結果
  filtering_data <- data.frame(time = as.numeric(time(y)),
                               var = "観測値",
                               val = as.numeric(y))
  filtering_data <- rbind(filtering_data,data.frame(time = as.numeric(time(y)),
                                                    var = "平均値",
                                                    val = m))
  filtering_data <- rbind(filtering_data,data.frame(time = as.numeric(time(y)),
                                                    var = "下側95%区間",
                                                    val = mapply(FUN = function(m,v)qnorm(p = 0.95,mean = m,sd = sqrt(v)),m,C)))
  filtering_data <- rbind(filtering_data,data.frame(time = as.numeric(time(y)),
                                                    var = "下側 5%区間",
                                                    val = mapply(FUN = function(m,v)qnorm(p = 0.05,mean = m,sd = sqrt(v)),m,C)))
  
  filtering_data$var <- factor(filtering_data$var,levels = c("観測値","平均値","下側95%区間","下側 5%区間"))
}

# カルマン平滑化
G_t_plus_1 <- G_t

#
kalman_smoothing <- function(s_t_plus_1, S_t_plus_1, t){
  # 平滑化利得
  A_t <- C[t] %*% t(G_t_plus_1) %*% solve(R[t+1])
  
  # 状態の更新
  s_t <- m[t] + A_t %*% (s_t_plus_1 - a[t+1])
  S_t <- C[t] + A_t %*% (S_t_plus_1 - R[t+1]) %*% t(A_t)
  
  # 平滑化分布の平均と分散を返す
  list(s = s_t, S = S_t)
}

# 平滑化分布の平均と分散を求める

s <- rep(NA_real_, t_max)
S <- rep(NA_real_, t_max)

s[t_max] <- m[t_max]
S[t_max] <- C[t_max]

for (t in (t_max-1):1){
  KS <- kalman_smoothing(s[t+1], S[t+1], t = t)
  s[t] <- KS$s
  S[t] <- KS$S
}


# 平滑化結果を図示
smoothed_data <- data.frame(time = c(as.numeric(time(y))),
                            var = "観測値",
                            val = c(as.numeric(y)))
smoothed_data <- rbind(smoothed_data,data.frame(time = c(as.numeric(time(y))),
                                                var = "平均値",
                                                val = as.numeric(s)))
smoothed_data <- rbind(smoothed_data,data.frame(time = c(as.numeric(time(y))),
                                                var = "下側95%区間",
                                                val = mapply(FUN = function(m,v)qnorm(p = 0.95,mean = m,sd = sqrt(v)),s,S)))
smoothed_data <- rbind(smoothed_data,data.frame(time = c(as.numeric(time(y))),
                                                var = "下側 5%区間",
                                                val = mapply(FUN = function(m,v)qnorm(p = 0.05,mean = m,sd = sqrt(v)),s,S)))

smoothed_data$var <- factor(smoothed_data$var,levels = c("観測値","平均値","下側95%区間","下側 5%区間"))


time_labels <- smoothed_data$time
flg_time_labels <- time_labels%%20==0
flg_time_labels[1] <- T
flg_time_labels[length(flg_time_labels)] <- T
flg_time_labels[unique(time_labels)==1970] <- T

time_labels[!flg_time_labels] <- ""

smoothed_data$time <- factor(smoothed_data$time,levels = unique(smoothed_data$time))


# 
g <- ggplot(data = smoothed_data,mapping = aes(x = time, y = val,group = var, color = var, linetype = var))
g <- g + geom_line(size = 0.6) + theme_bw()
g <- g + geom_line(mapping = aes(x = time, y = val,group = var, color = var, linetype = var),data = smoothed_data,size = 0.6) + theme_bw()
g <- g + labs(title = "ナイル川の流量データに対するカルマン平滑化適用結果",
              x = "年",y = paste0("ナイル川の流量[",expression(10^8),expression(m^3),"]"),group = "",color = "",
              linetype = "")
g <- g + theme(axis.text= element_text(size = 12),
               axis.text.x=element_text(size = 12),
               legend.title = element_blank(),
               legend.text = element_text(size = 12),
               legend.position = "top",
               axis.title = element_text(size = 12),
               axis.title.y = element_text(angle = 90,size = 12),
               plot.title= element_text(size = 16,hjust = 0.5))
g <- g + scale_x_discrete(labels = time_labels)
g <- g + geom_vline(xintercept = (1:length(flg_time_labels))[flg_time_labels], linetype =3)
g <- g + geom_vline(xintercept = (1:length(flg_time_labels))[unique(smoothed_data$time)==1970], linetype = 2)
g <- g + scale_linetype_manual(values = c("solid","solid","dashed","dashed"))
g <- g + scale_color_manual(values = c("black","red","blue","blue"))

plot(g)


参考文献

  • 沖本竜義(2010)「経済・ファイナンスデータの 計量時系列分析」(朝倉書店)
  • 北川源四郎(2020)「Rによる時系列モデリング入門」(岩波書店
  • 柴田里程(2017)「時系列解析」(共立出版)
  • 白石博(2022)「時系列データ解析」(森北出版)
  • 萩原淳一郎,瓜生真也,牧山幸史[著],石田基広[監修](2018)「基礎からわかる時系列分析 Rで実践するカルマンフィルタ・MCMC・粒子フィルタ」(技術評論社)
プライバシーポリシー お問い合わせ