状態空間モデル - R dlm （３）
今回は季節要素モデル。
例えば、データが四半期（周期=4）だった場合を考える。
Y[t-1] = α1 + v[t-1]
Y[t] = α2 + v[t]
Y[t+1] = α3 + v[t+1]
Y[t+2] = α4 + v[t+2]
α1+α2+α3+α4 = w[t], w[t] ~ N(0, w_seas)
α1,2,3,4が各四半期成分となり、和が平均0で分散w_seasの正規分布に従うとする。
動的線形モデルで表現する場合、各要素は以下になる。
t=0, θ[0]～N(m[0], C[0])
t>0, Y[t] = Fθ[t] + v[t], v[t]～N(0, V) 観測値
θ[t] = Gθ[t-1] + w[t], w[t]～N(0, W) 状態
ちなみに、
dlmModSeas()によって定義することができる。（dlmModTrig()を使うともっと柔軟に周期成分を捉えることができるようであるが、文献をさっと見ただけではよく理解できなかったので、必要があれば後々調べて行く。）
例えば、V=2、w_seas=4とした場合。
dlmModSeas(frequency=4, dV=2, dW=c(4,0,0))
$FF
[,1] [,2] [,3]
[1,] 1 0 0
$V
[,1]
[1,] 2
$GG
[,1] [,2] [,3]
[1,] -1 -1 -1
[2,] 1 0 0
[3,] 0 1 0
$W
[,1] [,2] [,3]
[1,] 4 0 0
[2,] 0 0 0
[3,] 0 0 0
$m0
[1] 0 0 0
$C0
[,1] [,2] [,3]
[1,] 1e+07 0e+00 0e+00
[2,] 0e+00 1e+07 0e+00
[3,] 0e+00 0e+00 1e+07
今回はこの季節要素モデル単体ではなく、ローカルレベルモデルと組み合わせてみる。ここら辺が状態空間モデルの柔軟なところである。
例えば、ローカルレベルモデルでV=3、W=6とし、季節要素モデル（四半期）でV=2、v_seas=4とした場合、以下のようにdlmModPoly()とdlmModSeas()を組み合わせる。
dlmModPoly( order=1, dV=3, dW=6 ) + dlmModSeas(frequency=4, dV=2, dW=c(4,0,0))
$FF
[,1] [,2] [,3] [,4]
[1,] 1 1 0 0
$V
[,1]
[1,] 5
$GG
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 -1 -1 -1
[3,] 0 1 0 0
[4,] 0 0 1 0
$W
[,1] [,2] [,3] [,4]
[1,] 6 0 0 0
[2,] 0 4 0 0
[3,] 0 0 0 0
[4,] 0 0 0 0
$m0
[1] 0 0 0 0
$C0
[,1] [,2] [,3] [,4]
[1,] 1e+07 0e+00 0e+00 0e+00
[2,] 0e+00 1e+07 0e+00 0e+00
[3,] 0e+00 0e+00 1e+07 0e+00
[4,] 0e+00 0e+00 0e+00 1e+07
この場合、モデルは以下に定義される。
Y[t] = μ[t] + α2 + v[t] v[t] ～N(0, V)
μ[t] = μ[t-1] + w1[t], w1[t]～N(0, w_loc)
α1+α2+α1+α2 = w2[t], w2[t] ~ N(0, w_seas)
動的線形モデルで表現する場合、各要素は以下になる。
t=0, θ[0]～N(m[0], C[0])
t>0, Y[t] = Fθ[t] + v[t], v[t]～N(0, V) 観測値
θ[t] = Gθ[t-1] + w[t], w[t]～N(0, W) 状態
ちなみに、
これまで見てきた各モデルの要素のベクトルや行列を、ブロックとして左や上から詰めたり、行列に対角に挿入したりの形式で表現される。
Vは、ローカルレベルと季節要素モデルそれぞれの観測誤差の和（3+2=5）になっている。
では、各パラメータを指定し、サンプルデータ（40行）を出力してみる。
alpha1 <- 10
alpha2 <- 20
alpha3 <- -5
alpha4 <- -alpha1 -alpha2 -alpha3 # -25
seas <- c(alpha1, alpha2, alpha3, alpha4)
V <- 3+2
w_loc <- 6
w_seas <- 4
tm <- c(0)
mu <- c(100) # μの初期値は100
alpha <- c(NA)
y <- c(NA)
# 40サンプル
for(i in 1:40){
tm <- append(tm, i)
mu_next <- mu[i] + rnorm(1, mean=0, sd=sqrt(w_loc))
mu <- append(mu, mu_next)
alpha_next <- seas[i%%4+1] + rnorm(1, mean=0, sd=sqrt(w_seas))
alpha <- append(alpha, alpha_next)
y_next <- mu_next + alpha_next + rnorm(1, mean=0, sd=sqrt(V))
y <- append(y, y_next)
}
df3 <- data.frame(time=tm, mu=mu, alpha=alpha, y=y)
time mu alpha y 1 0 100.00000 NA NA 2 1 101.03818 17.955605 117.69493 3 2 101.25233 -4.091898 98.77220 4 3 100.06346 -23.062593 72.51658 5 4 96.04669 8.901576 102.28330 6 5 95.93837 20.853973 114.93580 7 6 95.13444 -4.282241 88.50010 8 7 97.25849 -24.034810 71.62968 9 8 91.51833 8.773162 103.92921 10 9 89.47757 23.573238 111.84888 11 10 88.85872 -3.415835 85.99318 12 11 85.27012 -25.549430 60.39846 13 12 83.69630 10.338072 95.27553 14 13 81.79563 19.841001 96.37555 15 14 84.31110 -7.703048 76.44482 16 15 82.14169 -24.087354 62.47521 17 16 83.97272 8.371363 94.25735 18 17 86.51272 22.397388 109.34330 19 18 84.96344 -6.344667 78.87386 20 19 88.88974 -29.737400 56.65816 21 20 86.93938 8.035606 93.87915 22 21 85.26500 17.822354 106.99937 23 22 86.98058 -5.784289 79.34397 24 23 90.86066 -25.089460 64.02854 25 24 88.00473 7.444889 100.68429 26 25 92.70302 16.964471 110.13245 27 26 98.23367 -4.378634 91.80236 28 27 103.50861 -25.343409 75.63039 29 28 101.71510 7.622101 105.03880 30 29 99.40863 20.895265 118.87026 31 30 98.33315 -5.755254 97.06423 32 31 98.25762 -24.370844 72.00850 33 32 96.20491 10.091563 103.16654 34 33 92.37165 23.834555 117.60553 35 34 90.52986 -2.675492 84.79213 36 35 87.80711 -23.419353 62.24392 37 36 87.24032 7.758830 91.15328 38 37 84.09756 20.816570 105.28377 39 38 85.87142 -7.620894 76.79782 40 39 83.18068 -26.800049 57.26111 41 40 79.01976 11.497422 87.86101
プロット。
minY <- min(na.omit(df3$y))-10
maxY <-max(na.omit(df3$y))+10
par(mfrow = c(2, 1))
plot(df3$time, df3$y, ylim=c(minY,maxY),type="o", ylab="y, mu", xlab="time") # 黒 : データ
par(new=TRUE)
plot(df3$time, df3$mu, ylim=c(minY,maxY), type="o", col="blue", ylab="", xlab="") # 青 : mu
abline(h=100, lty=2)
plot(df3$time, df3$alpha, type="o", ylab="seasonal", xlab="time") # 黒 : alpha
abline(h=0, lty=2)
黒 : データy
青 : ローカルレベルモデル部分μ
黒: 季節要素部分α （グラフは、季節要素は別に描いている。）
では、関数定義と推定。
fn_mod3 <- function(parm) {
V <- parm[1]
w_loc <- parm[2]
w_seas <- parm[3]
return( dlmModPoly( order=1, dV=exp(V), dW=exp(w_loc) ) + dlmModSeas(frequency=4, dV=0, dW=c(exp(w_seas),0,0)) )
}
parm_start <- c(10*runif(1), 10*runif(1), 10*runif(1))
res_fn_mod3 <- dlmMLE(y=df3$y[-1], parm=parm_start, build=fn_mod3)
V_est <- exp(res_fn_mod3$par[1])
3.613708
w_loc_est <- exp(res_fn_mod3$par[2])
11.18024
w_seas_est <- exp(res_fn_mod3$par[3])
0.03253725
V= 3.613、w_loc= 11.180、w_seas= 0.032となる。これまでの推定結果で見てきた通り、真の値V = 5、w_loc=6、w_seas=4と異なる結果。
カルマンフィルタの実行。
mod3_est <- fn_mod3( res_fn_mod3$par )
fit_mod3_est <- dlmFilter(df3$y[-1], mod3_est)
str(fit_mod3_est, 1)
fit_mod3_est$m # 41×4の行列。必要なのは最初の2列。
[,1] [,2] [,3] [,4] [1,] 0.00000 0.000000 0.000000 0.0000000 [2,] 29.42375 88.271165 -29.423722 -29.4237217 [3,] 78.80987 19.962313 38.885103 -29.4237080 [4,] 86.52002 -14.003451 12.252181 31.1749761 [5,] 97.81681 4.466474 -25.300197 0.9554599 [6,] 96.28570 18.841810 5.404486 -24.9552884 [7,] 91.08020 -1.571212 20.276052 6.3781730 [8,] 94.54692 -23.586630 -2.546678 19.6680823 [9,] 96.33134 7.245955 -24.060232 -2.8045851 [10,] 93.62298 18.768130 7.954189 -23.7808371 [11,] 90.42758 -3.736562 19.320389 8.2521456 [12,] 86.20605 -24.885758 -3.006170 19.7127197 [13,] 86.81212 8.330082 -24.986893 -3.0549997 [14,] 79.70472 18.238856 9.520411 -24.5146983 [15,] 79.69382 -3.246505 18.240248 9.5210856 [16,] 84.88081 -23.592474 -3.907571 17.9172445 [17,] 84.73424 9.556773 -23.574268 -3.8994040 [18,] 89.55919 18.674913 8.953554 -23.8090888 [19,] 84.56592 -4.521445 19.188505 9.1835815 [20,] 81.61569 -24.265807 -4.218537 19.3252328 [21,] 83.87268 9.476095 -24.493025 -4.3140614 [22,] 86.65627 19.688730 9.193996 -24.5991631 [23,] 84.42293 -4.547578 19.884117 9.2753542 [24,] 87.53267 -24.244067 -4.819311 19.7703504 [25,] 90.37753 9.628814 -24.488764 -4.9157097 [26,] 90.36216 19.773962 9.630145 -24.4882853 [27,] 95.09226 -4.426344 19.408591 9.4886199 [28,] 98.81982 -24.085032 -4.714118 19.2968316 [29,] 96.37653 9.249946 -23.898850 -4.6455587 [30,] 98.76652 19.528847 9.066798 -23.9612860 [31,] 100.96262 -4.429610 19.374948 9.0114099 [32,] 97.21961 -24.305715 -4.167304 19.4693449 [33,] 94.93014 8.790651 -24.147025 -4.1128300 [34,] 97.33823 19.684315 8.622898 -24.2009000 [35,] 90.98728 -4.651370 20.096420 8.7605263 [36,] 87.57253 -24.498573 -4.429692 20.1702737 [37,] 83.67656 8.424338 -24.248019 -4.3496794 [38,] 84.75778 20.263011 8.354492 -24.2690601 [39,] 82.03062 -4.567471 20.429461 8.4060369 [40,] 81.65193 -24.298429 -4.544341 20.4365929 [41,] 79.99277 8.273214 -24.197904 -4.5145095
df3_1 <- data.frame( df3, m_loc=fit_mod3_est$m[,1], m_seas=fit_mod3_est$m[,2],
a_loc=c(NA,fit_mod3_est$a[,1]), a_seas=c(NA,fit_mod3_est$a[,2]), f=c(NA,fit_mod3_est$f) )
time mu alpha y m_loc m_seas a_loc a_seas f 1 0 100.00000 NA NA 0.00000 0.000000 NA NA NA 2 1 101.03818 17.955605 117.69493 29.42375 88.271165 0.00000 0.0000000 0.000000e+00 3 2 101.25233 -4.091898 98.77220 78.80987 19.962313 29.42375 -29.4237219 3.275159e-05 4 3 100.06346 -23.062593 72.51658 86.52002 -14.003451 78.80987 -29.4237081 4.938616e+01 5 4 96.04669 8.901576 102.28330 97.81681 4.466474 86.52002 -29.4237059 5.709632e+01 6 5 95.93837 20.853973 114.93580 96.28570 18.841810 97.81681 19.8782635 1.176951e+02 7 6 95.13444 -4.282241 88.50010 91.08020 -1.571212 96.28570 0.7089925 9.699469e+01 8 7 97.25849 -24.034810 71.62968 94.54692 -23.586630 91.08020 -25.0830134 6.599718e+01 9 8 91.51833 8.773162 103.92921 96.33134 7.245955 94.54692 6.4652266 1.010121e+02 10 9 89.47757 23.573238 111.84888 93.62298 18.768130 96.33134 19.6188627 1.159502e+02 11 10 88.85872 -3.415835 85.99318 90.42758 -3.736562 93.62298 -2.9414819 9.068150e+01 12 11 85.27012 -25.549430 60.39846 86.20605 -24.885758 90.42758 -23.8359724 6.659161e+01 13 12 83.69630 10.338072 95.27553 86.81212 8.330082 86.20605 8.1792081 9.438526e+01 14 13 81.79563 19.841001 96.37555 79.70472 18.238856 86.81212 19.7118112 1.065239e+02 15 14 84.31110 -7.703048 76.44482 79.69382 -3.246505 79.70472 -3.2445687 7.646015e+01 16 15 82.14169 -24.087354 62.47521 84.88081 -23.592474 79.69382 -24.5148283 5.517899e+01 17 16 83.97272 8.371363 94.25735 84.73424 9.556773 84.88081 9.5828005 9.446361e+01 18 17 86.51272 22.397388 109.34330 89.55919 18.674913 84.73424 17.9168992 1.026511e+02 19 18 84.96344 -6.344667 78.87386 84.56592 -4.521445 89.55919 -3.8193778 8.573981e+01 20 19 88.88974 -29.737400 56.65816 81.61569 -24.265807 84.56592 -23.8506414 6.071528e+01 21 20 86.93938 8.035606 93.87915 83.87268 9.476095 81.61569 9.1591109 9.077480e+01 22 21 85.26500 17.822354 106.99937 86.65627 19.688730 83.87268 19.3309907 1.032037e+02 23 22 86.98058 -5.784289 79.34397 84.42293 -4.547578 86.65627 -4.2835635 8.237271e+01 24 23 90.86066 -25.089460 64.02854 87.53267 -24.244067 84.42293 -24.6118931 5.981104e+01 25 24 88.00473 7.444889 100.68429 90.37753 9.628814 87.53267 9.2930274 9.682569e+01 26 25 92.70302 16.964471 110.13245 90.36216 19.773962 90.37753 19.7756597 1.101532e+02 27 26 98.23367 -4.378634 91.80236 95.09226 -4.426344 90.36216 -4.9158219 8.544634e+01 28 27 103.50861 -25.343409 75.63039 98.81982 -24.085032 95.09226 -24.4708669 7.062139e+01 29 28 101.71510 7.622101 105.03880 96.37653 9.249946 98.81982 9.5023191 1.083221e+02 30 29 99.40863 20.895265 118.87026 98.76652 19.528847 96.37653 19.2944622 1.156710e+02 31 30 98.33315 -5.755254 97.06423 100.96262 -4.429610 98.76652 -4.6343585 9.413216e+01 32 31 98.25762 -24.370844 72.00850 97.21961 -24.305715 100.96262 -23.9567479 7.700587e+01 33 32 96.20491 10.091563 103.16654 94.93014 8.790651 97.21961 9.0036746 1.062233e+02 34 33 92.37165 23.834555 117.60553 97.33823 19.684315 94.93014 19.4692037 1.143993e+02 35 34 90.52986 -2.675492 84.79213 90.98728 -4.651370 97.33823 -4.1063123 9.323192e+01 36 35 87.80711 -23.419353 62.24392 87.57253 -24.498573 90.98728 -24.2055755 6.678170e+01 37 36 87.24032 7.758830 91.15328 83.67656 8.424338 87.57253 8.7579920 9.633052e+01 38 37 84.09756 20.816570 105.28377 84.75778 20.263011 83.67656 20.1733599 1.038499e+02 39 38 85.87142 -7.620894 76.79782 82.03062 -4.567471 84.75778 -4.3484429 8.040934e+01 40 39 83.18068 -26.800049 57.26111 81.65193 -24.298429 82.03062 -24.2680267 5.776259e+01 41 40 79.01976 11.497422 87.86101 79.99277 8.273214 81.65193 8.4061765 9.005810e+01
観測値1期先予測（f）と、状態のフィルタリング（m）をプロットしてみる。
minY <- min(na.omit(df3$y))-20
maxY <-max(na.omit(df3$y))+20
par(mfrow = c(2, 1))
plot(df3$time, df3$y, ylim=c(minY,maxY),type="o", ylab="y, mu", xlab="time") # 黒 : データ
par(new=TRUE)
plot(df3$time, df3$mu, ylim=c(minY,maxY), type="o", col="blue", ylab="", xlab="") # 青 : mu
abline(h=100, lty=2)
par(new=TRUE)
plot(df3_1$time, df3_1$f, ylim=c(minY,maxY), type="o", pch=4, lty=3, ylab="", xlab="") # ? : f
par(new=TRUE)
plot(df3_1$time, df3_1$m_loc, ylim=c(minY,maxY), type="o", pch=4, lty=3, col="blue", ylab="", xlab="") # ? : m_mu
minTY <- min(na.omit(df3_1$alpha))-10
maxTY <- max(na.omit(df3_1$alpha))+10
plot(df3$time, df3$alpha, ylim=c(minTY,maxTY), type="o", ylab="seasonal", xlab="time") # 黒 : alpha
abline(h=0, lty=2)
par(new=TRUE)
plot(df3_1$time, df3_1$m_seas, ylim=c(minTY,maxTY), type="o", pch=4, lty=3, ylab="seasonal", xlab="time") # ? : m_trend
abline(h=0, lty=2)
黒 （○）: データy
青 （○）: ローカルレベルモデル部分μ
黒（×） : 観測値1期先予測
青 （×）: フィルタリングμ
黒（○）: 季節要素α （グラフは、季節要素は別に描いている。）
黒（×）: フィルタリングα
前回の線形成長モデルのトレンド成分より、今回の季節要素の方が予測しやすいようである。5順くらいしたところからフィルタリング値が適応していく様子が分かる。
