状態空間モデル - R dlm （２）
前回より続く。
線形成長モデルは以下に定義される。
Y[t] = μ[t] + v[t], v[t]～N(0, V)
μ[t] = μ[t-1] + β[t-1] + w1[t], w1[t]～N(0, W1)
β[t] = β[t-1] + w2[t], w2[t]～N(0, W2)
動的線形モデルで表現する場合、各要素は以下になる。
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) 状態
線形成長モデルは、ローカルレベルモデルにトレンドβが加わったもの。
W2=0と場合を考えてみると良い。そうするとβはtに依存しない定数となり、直線的に減少なり増加するトレンドを持つ時系列データを表現できる。初期値β[0]やW2の大きさで、線形の傾向を調整できる。
（実際のモデリングでW2=0とした場合、モデルは可制御でなくなる。可制御であるとはモデルが到達できない領域がなく、可制御でないとはモデルが到達できない領域がある、と解釈する。）
部分部分（ローカル）で線形な傾向を持つデータへ使用されるモデルだと考えることができる。
V = 10、W1 = 6、W2 = 3、μ[0] = 10、β[0] = 1 としてデータを作成してみる。
W2が小さく、β[0]を大きく設定すると、グローバルに直線傾向が強いデータが作成される。
V <- 10
W_mu <- 6
W_trend <- 3
mu.0 <- 10
trend.0 <- 1 # このトレンドの初期値が大きければ、直線的なトレンドの傾向が強くなる
tm <- c(0)
mu <- c(mu.0)
trend <- c(trend.0)
y <- c(NA)
# 40サンプル
for(i in 1:40){
tm <- append(tm, i)
trend_next <- trend[i] + rnorm(1, mean=0, sd=sqrt(W_trend))
trend <- append(trend, trend_next)
mu_next <- mu[i] + trend_next + rnorm(1, mean=0, sd=sqrt(W_mu))
mu <- append(mu, mu_next)
y_next <- mu_next + rnorm(1, mean=0, sd=sqrt(V))
y <- append(y, y_next)
}
df2 <- data.frame(time=tm, trend=trend, mu=mu, y=y)
minY <- min(na.omit(df2$y))-10
maxY <-max(na.omit(df2$y))+10
par(mfrow = c(2, 1))
plot(df2$time, df2$y, ylim=c(minY,maxY),type="o", ylab="y, mu", xlab="time") # 黒 : データ
par(new=TRUE)
plot(df2$time, df2$mu, ylim=c(minY,maxY), type="o", col="blue", ylab="", xlab="") # 青 : mu
plot(df2$time, df2$trend,type="o", ylab="trend", xlab="time") # 黒 : トレンド
abline(h=0, lty=2)
黒 : データy
青 : μ
黒: トレンドβ （グラフは、トレンドは別に描いている。）
では、モデルを関数で定義し、dlmMLE(データ, パラメータ初期値ベクトル, モデル関数)でパラメータを推定してみる。
fn_mod2 <- function(parm) {
V_estim <-parm[1]
W_mu_estim <- parm[2]
W_trend_estim <- parm[3]
return( dlmModPoly( order=2, dV=exp(V_estim), dW=c(exp(W_mu_estim),exp(W_trend_estim)) ) )
}
parm_start <- c(10*runif(1), 10*runif(1), 10*runif(1))
fn_mod2(parm_start)
$FF
[,1] [,2]
[1,] 1 0
$V
[,1]
[1,] 111.6958
$GG
[,1] [,2]
[1,] 1 1
[2,] 0 1
$W
[,1] [,2]
[1,] 1445.989 0.000
[2,] 0.000 1852.898
$m0
[1] 0 0
$C0
[,1] [,2]
[1,] 1e+07 0e+00
[2,] 0e+00 1e+07
状態の初期値m[0]とC[0]を指定しなければ、それらの要素は0と1e+07にデフォルトで指定される。
res_fn_mod2 <- dlmMLE(y=df2$y[-1], parm=parm_start, build=fn_mod2)
dlmMLE()でのパラメータ推定結果は以下。
V_est <- exp(res_fn_mod2$par[1])
9.692269
W_mu_est <- exp(res_fn_mod2$par[2])
3.757845
W_trend_est <- exp(res_fn_mod2$par[3])
7.397736
V= 9.692、W1= 3.757、W2= 7.397となる。真の値V = 10、W1 = 6、W2 = 3と異なる。
関数に渡す初期値を大きくしたり小さくしたりで推定をシミュレーションしてみると、渡す初期値の大きさで結果が大きく変わる。
今回のデータ数は40。データ数を数百程度に増やしてシミュレーションしてみても、データによって結果は近かったり遠かったりとまちまち。
そもそも一系列のデータyから、複数の分散などを推定する訳である。一意に値が定まる訳でもないであろうし、推定自体簡単なことではないことだと想像がつく。推定に関しては今後も調査継続の必要性あり。
観測値1期先予測（f）と、状態のフィルタリング（m）をプロットしてみる。
黒 （○）: データy
青 （○）: μ
黒（×） : 観測値1期先予測
青 （×）: フィルタリングμ
黒（○）: トレンドβ （グラフは、トレンドは別に描いている。）
黒（×）: フィルタリングβ
yやμのフィルタリングはある程度上手く行っているように見えるが、トレンドβの推定に関しては、ちょっと微妙な気がしないでもない。
真の値（V = 10、W1 = 6、W2 = 3）をモデルに指定しての推定も試ししてみたが同様の結果で、上手く行くとは限らない。
継続調査の必要あり。
分散を見てみる。
フィルタリング
dlmSvd2var(fit_mod2_est$U.C, fit_mod2_est$D.C)
これの最後のリスト。
41 [,1] [,2] [1,] 7.495514 4.031254 [2,] 4.031254 13.754984
1期先観測値予測
dlmSvd2var(fit_mod2_est$U.R, fit_mod2_est$D.R)
これの最後のリスト。
40 [,1] [,2] [1,] 33.07085 17.78624 [2,] 17.78624 21.15272
推定したフィルタや予測値の分散共分散行列は対角行列とならない。
モデルの定義では、状態のイノベーション項の分散共分散行列（mod2_est$W）は対角行列とするので、最初、これらの計算結果の出力に違和感があり戸惑ったが、線形成長モデルのμとβ間に相関があることは自明である。
k=10期先予測。
k <- 10 # k=10
mod2_pred <- dlmForecast(fit_mod2_est, nAhead=k)
str(mod2_pred, 1)
List of 4
$ a: num [1:10, 1:2] 94.2 95.3 96.4 97.5 98.6 ...
$ R:List of 10
$ f: num [1:10, 1] 94.2 95.3 96.4 97.5 98.6 ...
$ Q:List of 10
aは状態なので、μとβの2列の行列となる。
Rは状態の分散共分散行列（2×2）。
f、Qは長さ10のベクトルと解釈できる。
a_mu_pred <- unlist(mod2_pred$a)[,1]
a_trend_pred <- unlist(mod2_pred$a)[,2]
f_pred <- unlist(mod2_pred$f)
tm <- c( df2_1$time, length(df2_1$time):(length(df2_1$time)+length(f_pred)-1) )
tr <- c(df2_1$trend, rep(NA,length(f_pred)))
mu <- c(df2_1$mu, rep(NA,length(f_pred)))
y <- c(df2_1$y, rep(NA,length(f_pred)))
m_tr <- c(df2_1$m_trend, rep(NA,length(f_pred)))
m_mu <- c(df2_1$m_mu, rep(NA,length(f_pred)))
a_tr <- c(df2_1$a_trend, a_trend_pred)
a_mu <- c(df2_1$a_mu, a_mu_pred)
f <- c(df2_1$f, f_pred)
df2_2 <- data.frame( time=tm, trend=tr, mu=mu, y=y, m_trend=m_tr, m_mu=m_mu, a_trend=a_tr, a_mu=a_mu, f=f )
time trend mu y m_trend m_mu a_trend a_mu f 1 0 1.00000000 10.0000000 NA 0.00000000 0.000000 NA NA NA 2 1 -0.39189751 11.3792097 15.271752 7.63587107 15.271745 0.00000000 0.00000000 0.00000000 3 2 -0.27076775 11.2985444 7.616363 -7.65529646 7.616393 7.63587107 22.90761608 22.90761608 4 3 -0.86313367 6.4755666 4.722802 -5.03339589 4.091164 -7.65529646 -0.03890386 -0.03890386 5 4 -1.62576596 5.9510299 11.275962 0.30797550 8.824934 -5.03339589 -0.94223215 -0.94223215 6 5 -0.35411010 9.2943873 8.234765 -0.06646885 8.434325 0.30797550 9.13290905 9.13290905 7 6 0.41779153 8.7413314 3.370188 -2.14387961 4.501004 -0.06646885 8.36785609 8.36785609 8 7 -0.98943300 5.4297462 5.423797 -0.86825547 4.728817 -2.14387961 2.35712451 2.35712451 9 8 -0.77556706 7.6265669 8.362854 1.00459543 7.342526 -0.86825547 3.86056160 3.86056160 10 9 0.11205921 10.8966641 13.466506 3.13395497 12.306263 1.00459543 8.34712184 8.34712184 11 10 -0.28963512 8.6545542 9.477838 0.65404066 10.829192 3.13395497 15.44021832 15.44021832 12 11 -1.13876371 9.0680897 16.081112 2.56641211 15.039004 0.65404066 11.48323282 11.48323282 13 12 0.15358498 12.3559981 15.920635 1.86567026 16.302491 2.56641211 17.60541614 17.60541614 14 13 -0.04898526 10.3677927 16.102531 1.00652330 16.570706 1.86567026 18.16816149 18.16816149 15 14 1.44431649 10.4584454 8.183906 -2.90039216 10.312904 1.00652330 17.57722945 17.57722945 16 15 0.37100105 14.2978069 15.276352 0.37037322 13.494011 -2.90039216 7.41251195 7.41251195 17 16 -1.31472734 9.0689054 15.552997 1.07270896 15.170273 0.37037322 13.86438459 13.86438459 18 17 -1.66637109 4.1789927 2.907222 -4.47396319 5.929774 1.07270896 16.24298190 16.24298190 19 18 -5.26607314 -1.1893775 2.527121 -4.02837868 2.284309 -4.47396319 1.45581074 1.45581074 20 19 -4.45189524 -8.2568128 -15.576649 -9.78168978 -12.441492 -4.02837868 -1.74407013 -1.74407013 21 20 -1.25306575 -9.2654350 -10.873195 -5.06094989 -13.445672 -9.78168978 -22.22318215 -22.22318215 22 21 0.11424844 -9.2279975 -15.507593 -3.81357973 -16.187323 -5.06094989 -18.50662166 -18.50662166 23 22 3.82678457 -11.4033086 -13.048461 -0.92188742 -14.624233 -3.81357973 -20.00090297 -20.00090297 24 23 4.10467212 -3.5449498 -3.869699 3.93462463 -6.516162 -0.92188742 -15.54612067 -15.54612067 25 24 3.95499814 0.1072197 4.519471 6.88810939 2.910026 3.93462463 -2.58153764 -2.58153764 26 25 5.71548099 9.2240111 7.336158 5.86411239 7.894166 6.88810939 9.79813530 9.79813530 27 26 6.87021215 20.3360060 14.913178 6.34446366 14.651420 5.86411239 13.75827829 13.75827829 28 27 8.72990299 32.6152238 29.464687 9.86684844 27.545231 6.34446366 20.99588336 20.99588336 29 28 7.33977724 41.4001828 42.137907 11.83243660 41.066797 9.86684844 37.41207992 37.41207992 30 29 6.89292503 48.5144668 52.967099 11.86066321 52.951717 11.83243660 52.89923396 52.89923396 31 30 7.21971539 56.3389234 57.208671 8.69809251 58.932053 11.86066321 64.81238035 64.81238035 32 31 3.81404530 61.1853035 64.224510 7.28160456 64.996398 8.69809251 67.63014524 67.63014524 33 32 5.05958442 60.8899078 60.959617 2.57400845 63.524931 7.28160456 72.27800214 72.27800214 34 33 3.45025926 72.4182725 71.442502 4.79652835 70.231383 2.57400845 66.09893913 66.09893913 35 34 2.57374664 75.5788671 79.736036 6.75475372 78.668939 4.79652835 75.02791110 75.02791110 36 35 1.90856435 80.3943913 80.145944 4.55960767 81.342146 6.75475372 85.42369226 85.42369226 37 36 2.23365459 87.6591308 88.498695 5.63973962 87.910098 4.55960767 85.90175389 85.90175389 38 37 3.05362375 90.0239314 91.289013 4.69940675 91.801429 5.63973962 93.54983746 93.54983746 39 38 1.06568152 90.6666095 89.647751 1.84903939 91.201004 4.69940675 96.50083577 96.50083577 40 39 0.46012447 89.9896574 91.953965 1.39315362 92.202392 1.84903939 93.05004300 93.05004300 41 40 2.34790032 93.2382139 92.920977 1.11258411 93.073868 1.39315362 93.59554516 93.59554516 42 41 NA NA NA NA NA 1.11258411 94.18645220 94.18645220 43 42 NA NA NA NA NA 1.11258411 95.29903632 95.29903632 44 43 NA NA NA NA NA 1.11258411 96.41162043 96.41162043 45 44 NA NA NA NA NA 1.11258411 97.52420454 97.52420454 46 45 NA NA NA NA NA 1.11258411 98.63678866 98.63678866 47 46 NA NA NA NA NA 1.11258411 99.74937277 99.74937277 48 47 NA NA NA NA NA 1.11258411 100.86195688 100.86195688 49 48 NA NA NA NA NA 1.11258411 101.97454100 101.97454100 50 49 NA NA NA NA NA 1.11258411 103.08712511 103.08712511 51 50 NA NA NA NA NA 1.11258411 104.19970922 104.19970922
信頼区間は描かず、予測平均のみを示す。
par(mfrow = c(2, 1))
minY <- min(na.omit(df2_2$y))-20
maxY <-max(na.omit(df2_2$y))+30
plot(df2_2$time, df2_2$y, ylim=c(minY,maxY),type="o", ylab="y, mu, f, m_mu", xlab="time") # 黒 : データ
par(new=TRUE)
plot(df2_2$time, df2_2$mu, ylim=c(minY,maxY), type="o", col="red", ylab="", xlab="") # 赤 : mu
par(new=TRUE)
plot(df2_2$time, df2_2$f, ylim=c(minY,maxY), type="o", pch=4, lty=3, ylab="", xlab="") # 黒× : f
par(new=TRUE)
plot(df2_2$time, df2_2$a_mu, ylim=c(minY,maxY), type="o", pch=4, lty=3, col="red", ylab="", xlab="") # 赤× : a_mu
minTY <- min(na.omit(df2_2$trend))-5
maxTY <- max(na.omit(df2_2$trend))+5
plot(df2_2$time, df2_2$trend, ylim=c(minTY,maxTY), type="o", ylab="trend", xlab="time") # 黒 : トレンド
par(new=TRUE)
plot(df2_2$time, df2_2$a_trend, ylim=c(minTY,maxTY), type="o", pch=4, lty=3, ylab="trend", xlab="time") # ? : a_trend
abline(h=0, lty=2)
黒 （○）: データy
赤 （○）: μ
黒（×） : 観測値1期、k期先予測 f
赤 （×）: 1期、k期先予測 μ
黒（○）: トレンドβ （グラフは、トレンドは別に描いている。）
黒（×）: 1期、k期先予測トレンドβ
トレンドβは最後に1期先予測を行った後、定数となる。μはその時点からベータの累積和となり、直線的に増加または減少する。μと観測値予測fは同じになり、図ではfはμに上書きされている。
