以下の書籍
を中心に時系列解析を勉強していきます。
13. 線形・Gauss型状態空間モデルの逐次解法
13.1. カルマンフィルタ
線形・型状態空間モデルにおいて、推定すべきデータの真の値と点推定値の間の平均二乗誤差を最小にする意味で最適な逐次推定法はカルマンフィルタと呼ばれる。カルマンフィルタは非定常過程でも扱うことができる。カルマンフィルタで行なわれる計算はすべて線形計算で、正規分布の再生性から関心の対象となる分布はすべて正規分布である。
13.1.3 カルマン平滑化
時点までのカルマンフィルタリングが一旦完了しているとして、カルマン平滑化を考える。
線形形状態空間モデルでは平滑化分布も正規分布になる。そこで平滑化分布を
とおく。
時点までの平滑化分布を元に時点
での平滑化分布を求める手続きは以下のとおりである：
| 0. | 時点 |
|
| 1. | 時点 |
|
| 平滑化利得 |
||
| 状態の更新(平均) |
||
| 2. | 時点 |
このアルゴリズムをから時間を遡る方向に適用することで、全時点の平滑化分布を求めることができる。
- 平滑化分布：フィルタリング分布を補正
########################## ### 自作カルマン平滑化 ### ########################## 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)
