状態空間モデル - R dlm （４）
ARIMA(p,d,q)モデルを動的線形モデルで表現してみる。
ARIMAモデルの状態空間表現は一意に決定されないらしい、ここでは文献に従った形式を見て行く。
I(d)は差分を表し、通常のモデリングアプローチでは、データがトレンドを持ち非定常な場合は、その差分を取って定常に変換してからARMA(p,q)をあてはめるといった処理を行うが、動的線形モデルではこれまでに説明してきたような線形成長成分をモデルに含めることにより同様の処理ができる、よって必要であれば、線形成長+ARMA(p,q)モデルをあてはめでは使用するとのこと。
ではまず、AR(p=2)を例で見てみる。
モデルは以下に定義される。
Y[t] = φ1 Y[t-1] + φ2 Y[t-2] + ε[t], ε[t] ~ N(0, σ^2)
動的線形モデルで表現する場合、各要素は以下になる。
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) 状態
この場合、状態式は、
θ1[t] = φ1 θ1[t-1] + θ2[t-1] + ε[t]
θ2[t] = φ2 θ1[t-1]
となり、あわせて、
θ1[t] = φ1 θ1[t-1] +φ2 θ1[t-2] + ε[t]
と書ける。
観測式より、
Y[t] =θ1[t]
となり、AR(p=2)が動的線形モデルで表現されているのが分かる。
dlmModARMA()関数により、φ1=0.5、φ2=-0.3、σ^2=1を定義する場合、以下になる。
dlmModARMA(ar=c(0.5,-0.3), sigma2=1)
$FF
[,1] [,2]
[1,] 1 0
$V
[,1]
[1,] 0
$GG
[,1] [,2]
[1,] 0.5 1
[2,] -0.3 0
$W
[,1] [,2]
[1,] 1 0
[2,] 0 0
$m0
[1] 0 0
ARMA(p=2,q=1)の場合は。
Y[t] = φ1 Y[t-1] + φ2 Y[t-2] + ε[t] + ψ1 ε[t-1], ε[t] ~ N(0, σ^2)
ARMAの場合も観測式は変わらず、状態式が以下に変わって定義される。
θ[t] = Gθ[t-1] + w[t], w[t]～N(0, W) 状態
最後のWであるが、
と計算できる。
このように、MA成分は、誤差が過去の誤差と相関を持つと考えられる場合に含める。
dlmModARMA()関数により、φ1=0.5、φ2=-0.3、ψ1=0.4, σ^2=1を定義する場合、以下になる。
dlmModARMA(ar=c(0.5,-0.3), ma=c(0.4), sigma2=1)
$FF
[,1] [,2]
[1,] 1 0
$V
[,1]
[1,] 0
$GG
[,1] [,2]
[1,] 0.5 1
[2,] -0.3 0
$W
[,1] [,2]
[1,] 1.0 0.40
[2,] 0.4 0.16
$m0
[1] 0 0
$C0
[,1] [,2]
[1,] 1e+07 0e+00
[2,] 0e+00 1e+07
では、dlmModARMA(ar=c(0.5,-0.3), sigma2=5)とdlmModARMA(ar=c(0.5,-0.3), ma=c(0.4), sigma2=5)の、サンプルデータを出力してみる。
（sigma2=5（σ^2=5）と上の例から変更していることに注意）
これらのモデルよりの出力の違いはMA(q=1)部分（ψ1 ε[t-1]）のみ。
ar1 <- 0.5
ar2 <- -0.3
ma1 <- 0.4
inv <- 5
# AR(p=2)要素があるので、time=-1,0のパラメータ初期値を与えておく
tm <- c(-1, 0)
ar_2 <- c(100, 100)
eta_ar <- c(NA, NA) # MA(0)部分
eta_arma <- c(NA, 0) # MA(1)部分
y_ar <- c(NA, NA)
y_arma <- c(NA, NA)
# 60サンプル
for(i in (1+2):(60+2)){
tm <- append(tm, (i-2))
eta_ar_next <- rnorm(1, mean=0, sd=sqrt(inv))
eta_arma_next <- ma1*eta_arma[i-1] + eta_ar_next
ar_2_next <- ar1*ar_2[i-1] + ar2*ar_2[i-2]
eta_ar <- append(eta_ar, eta_ar_next)
eta_arma <- append(eta_arma, eta_arma_next)
ar_2 <- append(ar_2, ar_2_next)
y_ar <- append(y_ar, ar_2_next+eta_ar_next)
y_arma <- append(y_arma, ar_2_next+eta_arma_next)
}
df4_ar_arma <- data.frame(time=tm, ar_2=ar_2, eta_ar=eta_ar, eta_arma=eta_arma, y_ar=y_ar, y_arma=y_arma)
# time<=10を削る
df4_ar_arma_2 <- df4_ar_arma[13:nrow(df4_ar_arma),]
60サンプル作成し、最初の10データを切り落とした。データdf4_ar_arma_2。
これは、AR成分からのデータが定常と見えるようになるのに、いくらかステップを必要とするためである。（指定したパラメータ（φ1=0.5、φ2=-0.3）は理論上定常条件を満たす。だたし、シミュレーションでデータを出力する際は、最初の数ステップそうではないデータが出力される。）
プロット。
minY <- -10
maxY <- 10
plot(df4_ar_arma_2$time, df4_ar_arma_2$y_ar, ylim=c(minY,maxY), type="o", ylab="AR(p=2), ARMA(p=2, q=1)", xlab="time")
par(new=TRUE)
plot(df4_ar_arma_2$time, df4_ar_arma_2$y_arma, ylim=c(minY,maxY), type="o", col="blue", ylab="", xlab="")
abline(h=0, lty=2)
黒 : AR(p=2)
青 : ARMA(p=2, q=1)
当然だがψ1を大きく設定するほどと、ARとARMA系列の違いが大きくなる。
ここで、AR系列（y_ar）、ARMA系列（y_arma）に対して自己相関、偏自己相関を求めてみる。
AR(p=2)。
# 自己相関の確認(ar)
par(mfrow = c(2, 1))
acf(df4_ar_arma_2$y_ar, plot=TRUE)
pacf(df4_ar_arma_2$y_ar, plot=TRUE)
理論と異なりPartial ACFのLag=1,2が見られない。
ARMA(p=2, q=1)。
# 自己相関の確認(arma)
par(mfrow = c(2, 1))
acf(df4_ar_arma_2$y_arma, plot=TRUE)
pacf(df4_ar_arma_2$y_arma, plot=TRUE)
理論と異なりPartial ACFのLag=1,2が見られない。ただし、ACFのLag=1は出ている。
では、もうちょっと複雑なデータを作成し、複数のモデルをあてはめて比べてみる。
ローカルレベルモデル＋季節＋季節要素モデル＋ARMA(p=2, q=1)データを作成してみる。
dlmで定義する場合は以下。
dlmModPoly( order=1, dV=3, dW=6 ) + dlmModSeas(frequency=4, dV=2, dW=c(4,0,0)) + dlmModARMA(ar=c(0.5,-0.3), ma=c(0.2), sigma2=5)
$FF
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 1 0 0 1 0
$V
[,1]
[1,] 5
$GG
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0.0 0
[2,] 0 -1 -1 -1 0.0 0
[3,] 0 1 0 0 0.0 0
[4,] 0 0 1 0 0.0 0
[5,] 0 0 0 0 0.5 1
[6,] 0 0 0 0 -0.3 0
$W
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 6 0 0 0 0 0.0
[2,] 0 4 0 0 0 0.0
[3,] 0 0 0 0 0 0.0
[4,] 0 0 0 0 0 0.0
[5,] 0 0 0 0 5 1.0
[6,] 0 0 0 0 1 0.2
$m0
[1] 0 0 0 0 0 0
$C0
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00
[2,] 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00
[3,] 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00
[4,] 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00
[5,] 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00
[6,] 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07
データ作成。
## dlmModPoly( order=1, dV=3, dW=6 ) ##
V_loc <- 3
w_loc <- 6
## dlmModSeas(frequency=4, dV=2, dW=c(4,0,0)) ##
alpha1 <- 5
alpha2 <- 10
alpha3 <- -2
alpha4 <- -alpha1 -alpha2 -alpha3 # -13
seas <- c(alpha1, alpha2, alpha3, alpha4)
V_seas <- 2
w_seas <- 4
## ARMA(p=2, q=1), dlmModARMA(ar=c(0.5,-0.3), ma=c(0.2), sigma2=5) のデータ作成 ##
ar1 <- 0.5
ar2 <- -0.3
ma1 <- 0.4
sig2 <- 5
## Dataベクトル
tm <- c(-1, 0)
mu <- c(NA, 100) # μ(線形成長部分)の初期値は100
alpha <- c(NA, NA) # Seasonal
arma <- c(0, 0)
ips <- c(NA, rnorm(n=1, mean=0, sd=sqrt(sig2)))
y <- c(NA, NA)
# 70サンプル
for(i in (1+2):(70+2)){
tm <- append(tm, (i-2))
# Trend
mu_next <- mu[i-1] + rnorm(1, mean=0, sd=sqrt(w_loc))
mu <- append(mu, mu_next)
# Seasonal
alpha_next <- seas[(i-2)%%4+1] + rnorm(1, mean=0, sd=sqrt(w_seas))
alpha <- append(alpha, alpha_next)
# ARMA(p=2, q=1)
ips_next <- rnorm(1, mean=0, sd=sqrt(sig2)) + ma1*ips[i-1]
arma_next <- ar1*arma[i-1] + ar2*arma[i-2] + ips_next
arma <- append(arma, arma_next)
ips <- append(ips, ips_next)
# y
y_next <- mu_next + alpha_next + arma_next
y <- append(y, y_next)
}
df_0 <- data.frame(time=tm, mu=mu, alpha=alpha, arma_part=arma, ips=ips, y=y)
# time<=10を削る
df_1 <- df_0[13:nrow(df_0),]
ARMA過程はデータ全体（y）に対するものでなく、独立に作成されるARMA要素であることに注意。
先ほどと同様、最初の10データを切り落とした。
モデルを3つほど作り、推定とあてはめを行ってみる。
また、本当のパラメータをあてたモデルも作成してみる。
モデルA、データをシミュレーション作成したすべての項を含む。
dlmModPoly( order=1, dV=V, dW=w_loc ) + dlmModSeas(frequency=4, dW=c(w_seas,0,0)) + dlmModARMA(ar=c(ar1,ar2), ma=c(ma1), sigma2=sig2)
パラメータは7つ。
状態は以下。
ちなみに、
モデルB、ARMA項はAR(p=1)のみとする。
dlmModPoly( order=1, dV=V, dW=w_loc ) + dlmModSeas(frequency=4, dW=c(w_seas,0,0)) + dlmModARMA(ar=c(ar1), sigma2=sig2)
パラメータは5つ。
モデルC、ARMA項はなし。
dlmModPoly( order=1, dV=V, dW=w_loc ) + dlmModSeas(frequency=4, dW=c(w_seas,0,0))
パラメータは3つ。
また、モデルAに本当のパラメータ値を設定したものを、モデルTrueとする。
まず、モデルA、モデルB、モデルCに対し関数を定義し、dlmMLE()によるパラメータ推定を行う。
モデルA。
# モデル関数定義 - すべて
# dlmModPoly( order=1, dV=3, dW=6 )
# + dlmModSeas(frequency=4, dV=2, dW=c(4,0,0))
# + dlmModARMA(ar=c(0.5,-0.3), ma=c(0.2), sigma2=5)
fn_mod5A <- function(parm) {
V <- parm[1]
w_loc <- parm[2]
w_seas <- parm[3]
para_ar <- parm[4:5]
para_ma <- parm[6]
sig2 <- parm[7]
return(
dlmModPoly( order=1, dV=exp(V), dW=exp(w_loc) ) + dlmModSeas( frequency=4, dV=0, dW=c(exp(w_seas),0,0) ) + dlmModARMA( ar=para_ar, ma=para_ma, sigma2=exp(sig2) )
)
}
# parm = c( log(V), log(w_loc), log(w_seas), para_ar1, para_ar2, para_ma1, log(sigma2) )
## 推定
( parm_start <- c( 3*runif(1), 3*runif(1), 3*runif(1), 10*runif(1)-5, 10*runif(1)-5, 10*runif(1)-5, 3*runif(1)) )
( res_fn_mod5A <- dlmMLE(y=df_1$y, parm=parm_start, build=fn_mod5A) ) # たまにエラーが出る
res_fn_mod5A$value # 何回か試し、最大のものを採用
mod5A_est <- fn_mod5A( res_fn_mod5A$par )
$FF
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 1 0 0 1 0
$V
[,1]
[1,] 7.5373
$GG
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0.0000000 0
[2,] 0 -1 -1 -1 0.0000000 0
[3,] 0 1 0 0 0.0000000 0
[4,] 0 0 1 0 0.0000000 0
[5,] 0 0 0 0 1.5101608 1
[6,] 0 0 0 0 -0.9808136 0
$W
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3.255962 0.000000e+00 0 0 0.000000e+00 0.000000e+00
[2,] 0.000000 3.221217e-08 0 0 0.000000e+00 0.000000e+00
[3,] 0.000000 0.000000e+00 0 0 0.000000e+00 0.000000e+00
[4,] 0.000000 0.000000e+00 0 0 0.000000e+00 0.000000e+00
[5,] 0.000000 0.000000e+00 0 0 9.692737e-10 9.920913e-09
[6,] 0.000000 0.000000e+00 0 0 9.920913e-09 1.015446e-07
$m0
[1] 0 0 0 0 0 0
$C0
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1e+07 0e+00 0e+00 0e+00 0e+00 0e+00
[2,] 0e+00 1e+07 0e+00 0e+00 0e+00 0e+00
[3,] 0e+00 0e+00 1e+07 0e+00 0e+00 0e+00
[4,] 0e+00 0e+00 0e+00 1e+07 0e+00 0e+00
[5,] 0e+00 0e+00 0e+00 0e+00 1e+07 0e+00
[6,] 0e+00 0e+00 0e+00 0e+00 0e+00 1e+07
モデルB。
# モデル関数定義 - AR(1)まで、AR(2)とMA(1)はなし
# dlmModPoly( order=1, dV=3, dW=6 )
# + dlmModSeas(frequency=4, dV=2, dW=c(4,0,0))
# + dlmModARMA(ar=c(0.5), sigma2=5)
fn_mod5B <- function(parm) {
V <- parm[1]
w_loc <- parm[2]
w_seas <- parm[3]
para_ar <- parm[4]
sig2 <- parm[5]
return(
dlmModPoly( order=1, dV=exp(V), dW=exp(w_loc) )
+ dlmModSeas( frequency=4, dV=0, dW=c(exp(w_seas),0,0) )
+ dlmModARMA( ar=para_ar, sigma2=exp(sig2) )
)
}
# parm = c( log(V), log(w_loc), log(w_seas), para_ar1, log(sigma2) )
## 推定
( parm_start <- c( 3*runif(1), 3*runif(1), 3*runif(1), 10*runif(1)-5, 3*runif(1)) )
( res_fn_mod5B <- dlmMLE(y=df_1$y, parm=parm_start, build=fn_mod5B) )
res_fn_mod5B$value # 何回か試し、最大のものを採用
mod5B_est <- fn_mod5B( res_fn_mod5B$par )
$FF
[,1] [,2] [,3] [,4] [,5]
[1,] 1 1 0 0 1
$V
[,1]
[1,] 0.9031128
$GG
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 0.000000e+00
[2,] 0 -1 -1 -1 0.000000e+00
[3,] 0 1 0 0 0.000000e+00
[4,] 0 0 1 0 0.000000e+00
[5,] 0 0 0 0 2.072311e-06
$W
[,1] [,2] [,3] [,4] [,5]
[1,] 17.36985 0.00000000 0 0 0.0000000
[2,] 0.00000 0.01695153 0 0 0.0000000
[3,] 0.00000 0.00000000 0 0 0.0000000
[4,] 0.00000 0.00000000 0 0 0.0000000
[5,] 0.00000 0.00000000 0 0 0.6254222
$m0
[1] 0 0 0 0 0
$C0
[,1] [,2] [,3] [,4] [,5]
[1,] 1e+07 0e+00 0e+00 0e+00 0e+00
[2,] 0e+00 1e+07 0e+00 0e+00 0e+00
[3,] 0e+00 0e+00 1e+07 0e+00 0e+00
[4,] 0e+00 0e+00 0e+00 1e+07 0e+00
[5,] 0e+00 0e+00 0e+00 0e+00 1e+07
モデルC。
# モデル関数定義 - ARMAはなし
# dlmModPoly( order=1, dV=3, dW=6 )
# + dlmModSeas(frequency=4, dV=2, dW=c(4,0,0))
fn_mod5C <- 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 = c( log(V), log(w_loc), log(w_seas) )
## 推定
( parm_start <- c( 3*runif(1), 3*runif(1), 3*runif(1) ) )
( res_fn_mod5C <- dlmMLE(y=df_1$y, parm=parm_start, build=fn_mod5C) )
res_fn_mod5C$value # 何回か試し、最大のものを採用
mod5C_est <- fn_mod5C( res_fn_mod5C$par )
$FF
[,1] [,2] [,3] [,4]
[1,] 1 1 0 0
$V
[,1]
[1,] 1.986954
$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,] 16.69477 0.000000000 0 0
[2,] 0.00000 0.006144571 0 0
[3,] 0.00000 0.000000000 0 0
[4,] 0.00000 0.000000000 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[1:T]（Y[1], Y[2], … , Y[T]）が与えられた元でのフィルタリング分布、θ[T]|Y[1:T] ～N(m[T]=s[T], C[t]=S[T])、から始まり、後ろ向きに状態を推定して行く。
Y[1:T] が与えられた元でのt+1時点の状態の分布をθ[t+1]|Y[1:T]～N(s[t+1], S[t+1])、とすると、t時点の状態分布、θ[t]|Y[1:T]、は以下になる。（t < T）
θ[t]|Y[1:T]～N(s[t], S[t])
s[t] = m[t] + C[t] G[t+1]’ R[t+1]^-1 (s[t+1] – a[t+1])
S[t] = C[t] – C[t] G[t+1]’ R[t+1]^-1 (R[t+1] – S[t+1]) R[t+1]^-1 G[t+1] C[t]
t+1, t, t-1, t-2, t-3, … と後ろ向きに推定して行く。
m[t]、C[t]と言うのは、Y[1:t]が与えられた元でのフィルタリング分布θ[t]|Y[1:t] ～ N(m[t], C[t])から来ている。
式を見ると、フィルタリング分布、θ[t]|Y[1:t]、に、t+1時点のスムージング分布、θ[t+1]|Y[1:T]、とt+1状態予測分布、θ[t]|Y[1:t-1]、の差で調整が行われている形式で式が与えられていることが分かる。
dlmでは、dlmSmooth()関数を使用。この関数は平滑化推定した状態を返す。
モデルA、モデルB、モデルC、モデルTrueに対して平滑化を行う。
モデルA。
fit_mod5A_est <- dlmSmooth(df_1$y, mod5A_est)
#str( fit_mod5A_est, 1 )
df_temp <- dropFirst(fit_mod5A_est$s) # 60*6
# 1: Trend
# 2,3,4： Seasonal
# 5, 6: Theta1, 2
A_Trend <- df_temp[,1]
A_Seasonal <- df_temp[,2]
A_Theta <- df_temp[,5]
A_y <- A_Trend + A_Seasonal + A_Theta
必要なのは、1、2、5行目のみ。（状態の定義を確認しておく）
モデルB。
fit_mod5B_est <- dlmSmooth(df_1$y, mod5B_est)
#str( fit_mod5B_est, 1 )
df_temp <- dropFirst(fit_mod5B_est$s) # 60*5
# 1: Trend
# 2,3,4： Seasonal
# 5: Theta1
B_Trend <- df_temp[,1]
B_Seasonal <- df_temp[,2]
B_Theta <- df_temp[,5]
B_y <- B_Trend + B_Seasonal + B_Theta
モデルC。
fit_mod5C_est <- dlmSmooth(df_1$y, mod5C_est)
#str( fit_mod5C_est, 1 )
df_temp <- dropFirst(fit_mod5C_est$s) # 60*4
# 1: Trend
# 2,3,4： Seasonal
C_Trend <- df_temp[,1]
C_Seasonal <- df_temp[,2]
C_y <- C_Trend + C_Seasonal
本当のパラメータ値を与えたTrue modelを作成してみる。
# パラメータ c( log(V=5), log(w_loc=6), log(w_seas=4), para_ar1=0.5, para_ar2=-0.3, para_ma1=0.4, log(sigma2=5) )
mod5True_est <- fn_mod5A( c( log(5), log(6), log(4), 0.5, -0.3, 0.4, log(5) ) )
fit_mod5True_est <- dlmSmooth(df_1$y, mod5True_est)
#str( fit_mod5True_est, 1 )
df_temp <- dropFirst(fit_mod5True_est$s) # 60*6
# 1: Trend
# 2,3,4： Seasonal
# 5, 6: Theta1, 2
True_Trend <- df_temp[,1]
True_Seasonal <- df_temp[,2]
True_Theta <- df_temp[,5]
True_y <- True_Trend + True_Seasonal + True_Theta
全体データのdata.frameを作成。
df_2 <- data.frame( df_1, A_Trend, A_Seasonal, A_Theta, A_y, B_Trend, B_Seasonal, B_Theta, B_y, C_Trend, C_Seasonal, C_y, True_Trend, True_Seasonal, True_Theta, True_y)
横長になり見にくいので、観測、Trend、Seasonal、ARMA別の4つにdata.frameを作成しておく。
# y
df_2_y <- data.frame( time=df_2$time, y=df_1$y, A_y, B_y, C_y, True_y)
# Trend
df_2_Trend <- data.frame( time=df_2$time, mu=df_1$mu, A_Trend, B_Trend, C_Trend, True_Trend )
# Seasonal
df_2_Seasonal <- data.frame( time=df_2$time, alpha=df_1$alpha, A_Seasonal, B_Seasonal, C_Seasonal, True_Seasonal )
# ARMA
df_2_ARMA <- data.frame( time=df_2$time, alpha=df_1$arma_part, A_Theta, B_Theta, True_Theta )
プロットしてみる。
par(mfrow = c(4, 1))
# - Trend
yRange <- c( min(na.omit(df_2$mu))-10, max(na.omit(df_2$mu))+10 )
plot(df_2$time, df_2$mu, ylim=yRange, type="o", ylab="Trend", xlab="time")
par(new=TRUE)
plot(df_2$time, df_2$A_Trend, ylim=yRange, type="o", col="blue", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$B_Trend, ylim=yRange, type="o", col="green", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$C_Trend, ylim=yRange, type="o", col="red", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$True_Trend, ylim=yRange, type="o", col="orange", ylab="", xlab="")
abline(h=100, lty=2)
# - Seasonal
yRange <- c( min(na.omit(df_2$alpha))-5, max(na.omit(df_2$alpha))+5 )
plot(df_2$time, df_2$alpha, ylim=yRange, type="o", ylab="Seasonal", xlab="time")
par(new=TRUE)
plot(df_2$time, df_2$A_Seasonal, ylim=yRange, type="o", col="blue", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$B_Seasonal, ylim=yRange, type="o", col="green", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$C_Seasonal, ylim=yRange, type="o", col="red", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$True_Seasonal, ylim=yRange, type="o", col="orange", ylab="", xlab="")
abline(h=0, lty=2)
# - ARMA
yRange <- c( min(na.omit(df_2$arma_part))-5, max(na.omit(df_2$arma_part))+5 )
plot(df_2$time, df_2$arma_part, ylim=yRange, type="o", ylab="ARMA", xlab="time")
par(new=TRUE)
plot(df_2$time, df_2$A_Theta, ylim=yRange, type="o", col="blue", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$B_Theta, ylim=yRange, type="o", col="green", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$True_Theta, ylim=yRange, type="o", col="orange", ylab="", xlab="")
abline(h=0, lty=2)
# - y
yRange <- c( min(na.omit(df_2$y))-10, max(na.omit(df_2$y))+10 )
plot(df_2$time, df_2$y, ylim=yRange, type="o", ylab="y", xlab="time")
par(new=TRUE)
plot(df_2$time, df_2$A_y, ylim=yRange, type="o", col="blue", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$B_y, ylim=yRange, type="o", col="green", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$C_y, ylim=yRange, type="o", col="red", ylab="", xlab="")
par(new=TRUE)
plot(df_2$time, df_2$True_y, ylim=yRange, type="o", col="orange", ylab="", xlab="")
abline(h=100, lty=2)
青：モデルA
緑：モデルB
赤：モデルC
オレンジ：モデルTrue
Trend（ローカルレベルモデル部分）、ARMA部分が多少ずれている感じだが、Seasonalとyは精度良く重なっているように見える。
以下、それぞれのデータを出力しておく。
df_2_y
time y A_y B_y C_y True_y 1 11 73.28358 75.99096 73.40635 73.57516 73.43901 2 12 95.63303 94.79082 95.59404 95.53243 95.31520 3 13 98.53241 98.32929 98.75898 99.03478 99.15421 4 14 96.08600 92.27635 95.73925 95.31300 95.92870 5 15 82.63808 83.18041 82.69760 82.78869 82.25243 6 16 102.75459 104.70546 102.83121 102.88708 103.37537 7 17 106.84606 107.37548 106.67923 106.50316 105.87153 8 18 96.05072 96.16221 95.93277 95.77066 95.88189 9 19 79.82388 81.34860 80.08475 80.40366 80.10964 10 20 101.96427 98.70060 101.63623 101.24588 101.79034 11 21 97.91004 98.84474 98.15980 98.46897 98.88447 12 22 88.57271 90.07115 88.72245 88.90063 89.62626 13 23 79.23368 81.70584 79.50323 79.84233 79.61283 14 24 107.79603 106.36573 107.52452 107.19261 107.44635 15 25 112.00800 111.94567 111.97638 111.92955 111.54963 16 26 105.10427 103.71238 104.82302 104.48072 103.99548 17 27 89.20956 90.50919 89.38792 89.59617 89.07500 18 28 109.48524 107.96741 109.35587 109.19515 109.21430 19 29 108.95314 107.55349 108.91044 108.87185 109.12747 20 30 95.86211 96.45140 96.05580 96.28315 96.49760 21 31 84.82737 85.60520 84.94196 85.10225 85.33067 22 32 106.88476 108.16270 107.12502 107.38884 107.62509 23 33 117.51555 114.25135 117.04068 116.49414 116.31739 24 34 105.26862 106.23093 105.32864 105.38024 104.81727 25 35 93.14478 94.48602 93.14190 93.13511 93.30988 26 36 111.94797 112.69907 111.92647 111.89902 111.64500 27 37 113.09740 112.70658 112.99325 112.88288 112.92492 28 38 99.60180 101.18561 99.97215 100.39820 100.60706 29 39 93.27292 89.71142 92.95177 92.58032 93.22450 30 40 109.94548 110.05342 109.94667 109.97050 109.97343 31 41 108.87800 114.48105 109.51061 110.24554 110.67722 32 42 110.63067 109.94920 110.48961 110.33122 109.79586 33 43 107.36584 101.93530 106.69605 105.89477 105.73129 34 44 119.54284 121.09874 119.79842 120.08210 119.50202 35 45 121.50668 121.68850 121.51563 121.52466 120.98641 36 46 111.50938 109.77431 111.46001 111.40298 110.99179 37 47 98.99019 95.32858 98.72339 98.40222 98.84376 38 48 110.73148 112.12661 111.04847 111.42779 111.56197 39 49 113.03435 113.88980 113.03452 113.05115 113.28914 40 50 102.47221 105.65022 102.74417 103.07611 103.49664 41 51 96.31936 96.80014 96.27437 96.20235 96.40401 42 52 119.98334 118.51142 119.66790 119.31267 118.86787 43 53 119.35296 120.41398 119.45288 119.53649 119.64156 44 54 110.02348 109.48797 109.79517 109.55648 109.21425 45 55 92.94468 95.41333 93.29195 93.67943 93.88187 46 56 114.69503 114.20701 114.80487 114.94458 115.06694 47 57 122.16771 116.72011 121.65915 121.06935 120.75076 48 58 105.33282 106.85882 105.73141 106.19403 106.57378 49 59 96.56283 96.62765 96.58396 96.60242 96.20744 50 60 118.56346 117.87485 118.43300 118.30015 117.87052 51 61 120.36601 120.37190 120.47268 120.55232 120.49276 52 62 113.94532 109.59657 113.31656 112.60442 111.92860 53 63 90.80152 93.33706 91.10714 91.42967 91.72766 54 64 106.11965 108.92360 106.27266 106.48679 106.77951 55 65 107.32811 108.73927 107.36043 107.40850 107.83465 56 66 96.40241 98.79657 96.77677 97.22078 97.88382 57 67 92.79161 89.23014 92.41752 91.97532 92.05985 58 68 110.61652 110.47621 110.69833 110.79809 110.62977 59 69 114.01988 114.20505 113.99161 113.94294 114.01415 60 70 104.70104 105.35985 104.67579 104.65095 104.32458
df_2_Trend
time mu A_Trend B_Trend C_Trend True_Trend 1 11 87.12561 92.73568 86.79976 86.76771 94.73273 2 12 90.91686 93.90524 89.16125 89.21775 94.91931 3 13 90.63124 94.71098 90.77275 90.82252 94.72449 4 14 94.44074 95.42898 96.74195 96.64829 95.27583 5 15 95.19361 94.50129 96.04186 95.97913 95.63841 6 16 93.72475 93.80788 96.48650 96.57548 95.53820 7 17 96.50994 93.95720 98.40483 98.28510 96.18294 8 18 94.93648 94.33522 97.11462 97.11365 95.65824 9 19 96.44080 94.76141 93.55581 93.58908 94.93094 10 20 96.72234 95.84623 95.01442 94.93594 94.54656 11 21 93.59476 95.52122 90.16362 90.24676 93.95346 12 22 95.19088 95.59998 90.11656 90.25378 94.52967 13 23 95.73103 96.32604 92.94949 93.01604 96.37015 14 24 98.23800 98.12002 100.96690 100.89230 98.66561 15 25 98.64091 99.29615 103.76236 103.69853 100.54146 16 26 100.87863 100.44534 105.94969 105.84563 101.86726 17 27 102.75371 100.99328 102.72781 102.75356 101.86252 18 28 102.53568 102.10263 102.93636 102.90984 101.69631 19 29 102.33942 102.55630 100.65659 100.62872 101.20496 20 30 106.04527 102.40536 97.55569 97.66461 100.92281 21 31 103.44166 102.50898 98.18008 98.23811 101.40324 22 32 102.20754 102.94860 101.00859 101.12126 102.48765 23 33 104.01406 103.94027 108.45816 108.23979 104.46046 24 34 105.83723 103.52187 106.77435 106.77626 104.99547 25 35 106.13540 103.51918 106.24489 106.25060 104.98886 26 36 105.25531 104.09588 105.66006 105.64375 105.18037 27 37 101.51205 104.99703 104.66158 104.62555 105.00832 28 38 103.74558 105.72936 101.65996 101.80489 104.62929 29 39 107.27372 107.14586 105.78141 105.67569 105.45657 30 40 108.11166 107.02386 103.72609 103.72715 106.22575 31 41 105.73444 106.94850 101.69371 101.98885 107.02846 32 42 110.01180 109.29354 111.82855 111.74086 109.99025 33 43 113.51907 111.34420 119.25034 118.97685 111.95026 34 44 115.63737 111.04898 113.78991 113.85261 111.94881 35 45 114.20380 111.42587 113.24503 113.25932 111.89837 36 46 112.49679 111.88131 112.87235 112.81709 111.22361 37 47 110.83887 111.58723 111.55012 111.48084 109.92774 38 48 108.50605 109.71142 105.09641 105.20432 108.45616 39 49 109.04362 108.43827 104.73955 104.77838 107.98117 40 50 108.92396 107.53467 104.38594 104.49360 107.81192 41 51 111.42104 108.00390 109.26313 109.28297 108.87200 42 52 114.74495 108.68081 113.27488 113.08915 110.03366 43 53 111.52418 108.72189 111.21976 111.26020 109.85675 44 54 110.26100 109.22130 111.08634 110.97329 110.02617 45 55 107.69165 109.48939 106.56166 106.76258 109.22451 46 56 110.33002 110.82388 108.71624 108.72537 109.54748 47 57 110.78981 111.94756 112.98347 112.78490 110.31674 48 58 110.74215 110.71799 107.46955 107.61589 109.38567 49 59 111.37021 110.14762 109.62185 109.68288 109.94375 50 60 112.63455 109.60525 112.18069 112.08258 110.07536 51 61 112.39075 108.76541 112.23050 112.26994 109.37545 52 62 108.06545 107.92812 114.33195 114.02274 108.82763 53 63 101.66834 105.21225 104.34030 104.50904 105.85975 54 64 101.23011 103.59169 100.22668 100.27318 104.00323 55 65 100.92076 103.18238 99.05595 99.12206 102.93854 56 66 103.00017 103.38266 98.50688 98.64636 102.48171 57 67 104.43158 104.61717 105.15783 105.04670 103.80255 58 68 106.73214 104.31319 104.61367 104.58836 104.24529 59 69 106.02372 103.94861 105.64305 105.65562 104.70392 60 70 108.38857 103.66402 106.12875 106.07646 105.15568
df_2_Seasonal
time alpha A_Seasonal B_Seasonal C_Seasonal True_Seasonal 1 11 -13.36256699 -13.218408 -13.308385 -13.192554 -15.22011446 2 12 5.17355551 6.279094 6.405786 6.314685 8.38132976 3 13 8.44608722 8.304729 8.143135 8.212256 7.12305756 4 14 0.04760801 -1.365415 -1.242839 -1.335288 -0.40862301 5 15 -14.09916173 -13.218408 -13.303044 -13.190440 -14.71714949 6 16 8.66657308 6.279094 6.397764 6.311608 7.25100681 7 17 8.98996455 8.304729 8.158881 8.218064 8.49805433 8 18 0.59489880 -1.365416 -1.263531 -1.342990 -0.97358871 9 19 -14.02022332 -13.218408 -13.290399 -13.185413 -15.20200264 10 20 6.39373344 6.279094 6.394634 6.309944 8.20207879 11 21 6.77721872 8.304729 8.169140 8.222209 7.95224371 12 22 -4.30859705 -1.365416 -1.290416 -1.353157 -1.25767154 13 23 -15.19242573 -13.218408 -13.259588 -13.173713 -14.95542934 14 24 8.07927546 6.279094 6.369604 6.300317 7.86671328 15 25 8.67145494 8.304729 8.192122 8.231017 8.26182296 16 26 0.13101631 -1.365416 -1.321428 -1.364906 -0.93894737 17 27 -13.49776330 -13.218408 -13.216373 -13.157396 -14.66529427 18 28 6.22057484 6.279094 6.329915 6.285317 7.03522535 19 29 8.64530608 8.304729 8.224293 8.243135 9.00478322 20 30 -2.65547851 -1.365416 -1.365751 -1.381465 -1.91993373 21 31 -12.85039838 -13.218408 -13.158755 -13.135863 -13.48667036 22 32 6.84012772 6.279094 6.282823 6.267579 5.73840056 23 33 11.62320506 8.304729 8.253653 8.254345 9.73505022 24 34 -1.31846980 -1.365416 -1.404153 -1.396019 -2.42625885 25 35 -14.90706177 -13.218408 -13.104981 -13.115490 -12.60340554 26 36 2.75753311 6.279094 6.251513 6.255267 6.18198685 27 37 11.27941576 8.304729 8.259551 8.257334 8.31707979 28 38 -0.47805002 -1.365416 -1.431335 -1.406693 -2.82830642 29 39 -8.63046568 -13.218407 -13.052037 -13.095371 -10.85250821 30 40 6.60639917 6.279094 6.221405 6.243346 6.14670762 31 41 6.42913659 8.304729 8.254991 8.256685 6.06131829 32 42 -0.56436973 -1.365416 -1.436631 -1.409637 -1.44521418 33 43 -11.59419067 -13.218407 -13.018123 -13.082076 -10.00566072 34 44 1.54597072 6.279094 6.185495 6.229492 4.75551093 35 45 8.41748316 8.304729 8.276804 8.265342 7.32980631 36 46 -0.16818281 -1.365416 -1.446524 -1.414112 -1.52956476 37 47 -12.44755222 -13.218407 -13.011503 -13.078621 -11.07358138 38 48 4.79046549 6.279094 6.171584 6.223467 5.02285080 39 49 7.77073833 8.304729 8.295083 8.272772 8.21259479 40 50 -5.52502745 -1.365416 -1.453431 -1.417483 -1.90870516 41 51 -12.91564519 -13.218407 -13.019921 -13.080627 -12.62610519 42 52 7.16855449 6.279094 6.174577 6.223516 6.53229046 43 53 10.50594694 8.304729 8.302314 8.276285 8.01910346 44 54 -0.84858308 -1.365416 -1.449286 -1.416809 -0.92030022 45 55 -14.08696249 -13.218407 -13.029213 -13.083150 -13.97036895 46 56 5.27121128 6.279094 6.164698 6.219210 5.95838757 47 57 8.96362883 8.304729 8.323501 8.284451 9.82713360 48 58 -7.07166428 -1.365416 -1.462107 -1.421863 -2.20730405 49 59 -13.95720918 -13.218407 -13.023245 -13.080456 -13.46526984 50 60 5.42283601 6.279094 6.161969 6.217572 6.36334987 51 61 4.29264261 8.304729 8.316055 8.282389 8.07775193 52 62 -0.44433216 -1.365416 -1.450816 -1.418316 -0.09090592 53 63 -13.59436263 -13.218407 -13.021515 -13.079366 -13.96720730 54 64 7.18261358 6.279094 6.151945 6.213607 5.84252114 55 65 8.55606165 8.304729 8.326863 8.286440 8.69889453 56 66 -4.51813711 -1.365416 -1.470867 -1.425582 -2.04356950 57 67 -10.18812438 -13.218407 -12.999383 -13.071380 -11.90183563 58 68 5.14654833 6.279094 6.141321 6.209722 5.23132574 59 69 9.18859305 8.304729 8.328986 8.287323 8.41749164 60 70 0.79856678 -1.365416 -1.470449 -1.425510 -1.44580978
df_2_ARMA
time alpha A_Theta B_Theta True_Theta 1 11 -0.4794681 -3.5263139 -0.0850285006 -6.07360211 2 12 -0.4573817 -5.3935164 0.0270037193 -7.98544421 3 13 -0.5449209 -4.6864201 -0.1569037302 -2.69334086 4 14 1.5976522 -1.7872135 0.2401354160 1.06148534 5 15 1.5436321 1.8975249 -0.0412172489 1.33117340 6 16 0.3632648 4.6184910 -0.0530613392 0.58616134 7 17 1.3461566 5.1135459 0.1155271784 1.19053631 8 18 0.5193396 3.1923975 0.0816833493 1.19723785 9 19 -2.5966944 -0.1944019 -0.1806573296 0.38070583 10 20 -1.1518024 -3.4247252 0.2271767560 -0.95829166 11 21 -2.4619386 -4.9812138 -0.1729639327 -3.02123490 12 22 -2.3095767 -4.1634166 -0.1036981970 -3.64574397 13 23 -1.3049305 -1.4017861 -0.1866730342 -1.80188361 14 24 1.4787524 1.9666134 0.1880219929 0.91402104 15 25 4.6956335 4.3447933 0.0218973117 2.74634464 16 26 4.0946239 4.6324553 0.1947647583 3.06716718 17 27 -0.0463816 2.7343198 -0.1235160722 1.87777232 18 28 0.7289911 -0.4143127 0.0895948354 0.48276135 19 29 -2.0315891 -3.3075370 0.0295654899 -1.08226932 20 30 -7.5276731 -4.5885491 -0.1341336622 -2.50527696 21 31 -5.7638933 -3.6853695 -0.0793622637 -2.58589985 22 32 -2.1629076 -1.0649890 -0.1663868208 -0.60096201 23 33 1.8782878 2.0063561 0.3288581981 2.12188387 24 34 0.7498567 4.0744759 -0.0415628899 2.24805487 25 35 1.9164429 4.1852524 0.0019933021 0.92443052 26 36 3.9351280 2.3241025 0.0148943374 0.28263753 27 37 0.3059342 -0.5951841 0.0721247332 -0.40047196 28 38 -3.6657260 -3.1783352 -0.2564738914 -1.19391850 29 39 -5.3703340 -4.2160325 0.2224013231 -1.37956587 30 40 -4.7725878 -3.2495325 -0.0008265259 -2.39902633 31 41 -3.2855754 -0.7721744 -0.4380949412 -2.41255762 32 42 1.1832368 2.0210784 0.0976868770 1.25082104 33 43 5.4409647 3.8095125 0.4638392587 3.78669184 34 44 2.3594988 3.7706751 -0.1769891233 2.79769632 35 45 -1.1146063 1.9579038 -0.0062005831 1.75822990 36 46 -0.8192210 -0.7415799 0.0341897975 1.29774590 37 47 0.5988702 -3.0402436 0.1847650153 -0.01039476 38 48 -2.5650444 -3.8639049 -0.2195242129 -1.91703685 39 49 -3.7800111 -2.8532052 -0.0001177280 -2.90462139 40 50 -0.9267270 -0.5190279 -0.1883410003 -2.40657855 41 51 -2.1860299 2.0146470 0.0311609106 0.15811413 42 52 -1.9301608 3.5515104 0.2184451626 2.30192816 43 53 -2.6771621 3.3873583 -0.0691926554 1.76570705 44 54 0.6110683 1.6320857 0.1581121346 0.10837659 45 55 -0.6600073 -0.8576555 -0.2404946915 -1.37227381 46 56 -0.9062060 -2.8959695 -0.0760681287 -0.43892484 47 57 2.4142625 -3.5321792 0.3521812147 0.60688167 48 58 1.6623399 -2.4937521 -0.2760306149 -0.60458708 49 59 -0.8501737 -0.3015568 -0.0146386640 -0.27104378 50 60 0.5060746 1.9905068 0.0903408265 1.43180166 51 61 3.6826219 3.3017563 -0.0738709121 3.03955618 52 62 6.3242019 3.0338666 0.4354262236 3.19187368 53 63 2.7275425 1.3432185 -0.2116450932 -0.16488412 54 64 -2.2930748 -0.9471818 -0.1059628809 -3.06623936 55 65 -2.1487083 -2.7478438 -0.0223839425 -3.80278458 56 66 -2.0796168 -3.2206771 -0.2592452779 -2.55431885 57 67 -1.4518376 -2.1686175 0.2590683629 0.15913381 58 68 -1.2621713 -0.1160771 -0.0566570487 1.15315296 59 69 -1.1924385 1.9517145 0.0195759577 0.89273554 60 70 -4.4860967 3.0612526 0.0174882250 0.61470634
状態空間モデル - 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順くらいしたところからフィルタリング値が適応していく様子が分かる。
状態空間モデル - 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はμに上書きされている。
状態空間モデル - R dlm （１）
状態空間モデル（正確には、時系列モデルの状態空間表現）の学習メモ。Rのdlmパッケージを使う。
参考書は ”和合 2013”。
単純なモデル（ローカルレベルモデル）を状態空間表現し、カルマンフィルタを実行してみる。
動的線形モデル（Dynamic Linear Model） / 線形ガウス状態空間モデル
t=0, θ[0]～N(m[0], C[0])
t>0, Y[t] = F[t] θ[t] + v[t], v[t]～N(0, V[t]) 観測値
θ[t] = G[t]θ[t-1] + w[t], w[t]～N(0, W[t]) 状態
Yはm次元のベクトル、θはp次元のベクトルとするが、多変量時系列でなければm=1とシンプルになる。
m=1とすると、それぞれの要素の次元は、Fが1×p、v、Vはスカラー、Gはp×p、wはp×1、Wはp×pとなる。
式のパラメータ（F[t]、V[t]、G[t]、W[t]）は[t]を付けた時間依存を表す形式で表記したが、多くの場合は時間依存とすることはない。
状態（θ）とは観測不能で、あるシステムを構成する要素、観測値（Y）は状態より発生する現象から実際に観測できるものを解釈する。
カルマンフィルタにより与えられたデータからθとYの分布を推定することができる。
イメージ図。
θ[0]から始まり、θ[1]を予測し、Y[1]を予測する。そして新たなデータy[1]が観測され、その情報を用いてθ[1]が修正される（フィルタリング）。
そのθ[1]からθ[2]を予測し、Y[2]を予測する。そして新たなデータy[2]が観測され、その情報を用いてθ[2]が修正される
これを繰り返しにより、θ、Yの分布が求められる。
各予測値とフィルタリング値の分布は以下に定義される。動的線形モデルは正規分布を仮定しているので、平均と分散により分布が特定される。
t=0時点の状態の分布、θ[0]～N(m[0], C[0])、が与えられているとする。
(1) Y[1:t-1]（Y[1],Y[2], … , Y[t-1]）が与えられた元での1期先状態予測分布
θ[t]|Y[1:t-1] ～ N(a[t], R[t])
a[t] = G[t] m[t-1]
R[t] = G[t] C[t-1] G[t]’ + W[t]
(2) Y[1:t-1]が与えられた元での1期先観測予測分布
Y[t]|Y[1:t-1] ～ N(f[t], Q[t])
f[t] = F[t] a[t]
Q[t] = F[t] R[t] F[t]’ + V[t]
(3) Y[1:t]が与えられた元でのフィルタリング分布
θ[t]|Y[1:t] ～ N(m[t], C[t])
m[t] = a[t] + R[t] F[t]’ Q[t]^-1 e[t] e[t] = Y[t] – f[t]
C[t] = R[t] - R[t] F[t]’ Q[t]^-1 F[t] R[t]
カルマンフィルタは基本的にY[1:t]が与えられたときのθ[t]のフィルタリング分布の推定を目的としている。Y[1:t]が与えられたときのθ[t+k]やY[t+k]のk期先予測の分布に関しては後述（k>1）。
この分野でもっともシンプルなモデルとして用いられるローカルレベルモデルで、状態空間モデルの定義やカルマンフィルタを試してみる。
ローカルレベルモデル （ランダムウォーク＋ノイズモデル）では、F=1、G＝1と定義される。
t=0, θ[0]～N(m[0], C[0])
t>0, Y[t] = θ[t] + v[t], v[t]～N(0, V) 観測値
θ[t] =θ[t-1] + w[t], w[t]～N(0, W) 状態
θ[t]が時刻tでの真の値、Y[t]がθ[t]の観測値、v[t]が観測誤差。θ[t]は一期先の値θ[t-1]にランダムな値w[t]が加わって変化するシステムだと考える単純なモデル。
θをランダムウォークとすると、Yはそれに誤差vを加えたもの。
# m[0]=10、c[0]=5
θ[0]=10、V=3、W=6と設定して、シミュレーションで長さ20データを作成してみる。
library(dlm)
V <- 3
W <- 6
theta.0 <- 10
time <- c(0)
theta <- c(theta.0)
y <- c(NA)
for(i in 1:20){
time <- append(time, i)
theta_next <- theta[i] + rnorm(1, mean=0, sd=sqrt(W))
theta <- append(theta, theta_next)
y_next <- theta_next + rnorm(1, mean=0, sd=sqrt(V))
y <- append(y, y_next)
}
df1 <- data.frame(time=time, theta=theta, y=y)
minY <- min(df1$theta)-10
maxY <-max(df1$theta)+10
plot(df1$time, df1$theta, ylim=c(minY,maxY),type="o", col="blue", ylab="val(theta, y)")
par(new=TRUE)
plot(df1$time, df1$y, ylim=c(minY,maxY), type="o", ylab="")
青：θ、黒：Y
では、ここから逆に、このデータ（Yのみ）が与えられたとして、カルマンフィルタにより状態と観測値の推定を行う。
当てはめの際は、まずモデルを定義する。
モデルは上でデータを作成したとおり、ローカルレベルモデルとし、m[0]=10、C[0]=50と仮定する。また、V=3、W=6は既知とする。
（状態の事前分布N(m[0], C[0])に関する情報が不確かな場合、C[0]を大きく仮定するらしい。）
dlm()関数を使用した場合は以下。
m.0 <- 10
C.0 <- 50
mod1 <- dlm( m0=m.0, C0=C.0, FF=1, V=3, GG=1, W=6 )
unlist(mod1)
m0 C0 FF V GG W 10 50 1 3 1 6
class(mod1) # "dlm"
typeof(mod1) # "list"
dlmModPoly()を使ってもローカルレベルモデルを定義できる。
mod2 <- dlmModPoly( order=1, m0=m.0, C0=C.0, dV=3, dW=6 )
unlist(mod2)
m0 C0 FF V GG W 10 50 1 3 1 6
dlmFilter(データ, 当てはめるモデル)により、カルマンフィルタを実行し、1期先状態予測、1期先観測値予測、フィルタリングの値を計算する。
fit_mod1 <- dlmFilter(df1$y[-1], mod1)
str(fit_mod1, 1)
List of 9
$ y : num [1:20] 11.48 14.89 16.27 15.19 7.64 ...
$ mod:List of 6
..- attr(*, "class")= chr "dlm"
$ m : num [1:21] 10 11.4 14 15.7 15.3 ...
$ U.C:List of 21
$ D.C: num [1:21, 1] 7.07 1.69 1.5 1.48 1.48 ...
$ a : num [1:20] 10 11.4 14 15.7 15.3 ...
$ U.R:List of 20
$ D.R: num [1:20, 1] 7.48 2.97 2.87 2.86 2.86 ...
$ f : num [1:20] 10 11.4 14 15.7 15.3 ...
- attr(*, "class")= chr "dlmFiltered"
fit_mod1$a # 1期先状態予測
fit_mod1$f # 1期先観測値予測
fit_mod1$m # フィルタリング
データと、それら計算結果をまとめたデータフレームを作成。
df_res1 <- data.frame(time=df1$time, y=df1$y, theta=df1$theta, a=c(NA,fit_mod1$a), f=c(NA,fit_mod1$f), m=fit_mod1$m)
time y theta a f m 1 0 NA 10.00000 NA NA 10.000000 2 1 11.480221 11.12716 10.000000 10.000000 11.404956 3 2 14.887411 14.39119 11.404956 11.404956 14.005587 4 3 16.268663 14.81166 14.005587 14.005587 15.664657 5 4 15.192051 14.39024 15.664657 15.664657 15.318650 6 5 7.640275 11.47224 15.318650 15.318650 9.697648 7 6 11.918582 12.24209 9.697648 9.697648 11.323486 8 7 11.739846 13.41268 11.323486 11.323486 11.628283 9 8 19.019994 18.91138 11.628283 11.628283 17.039391 10 9 21.572069 19.41500 17.039391 17.039391 20.357542 11 10 20.391132 17.97181 20.357542 20.357542 20.382131 12 11 15.116908 18.52915 20.382131 20.382131 16.527721 13 12 19.366015 21.61766 16.527721 16.527721 18.605496 14 13 21.751131 19.29930 18.605496 18.605496 20.908260 15 14 16.585866 17.66393 20.908260 20.908260 17.744048 16 15 17.432607 16.43332 17.744048 17.744048 17.516058 17 16 22.007343 19.82772 17.516058 17.516058 20.803907 18 17 18.873734 17.81004 20.803907 20.803907 19.390922 19 18 19.547199 19.00535 19.390922 19.390922 19.505325 20 19 17.828754 19.51896 19.505325 19.505325 18.277990 21 20 23.217935 22.09485 18.277990 18.277990 21.894281
plot(df_res1$time, df_res1$theta, ylim=c(minY,maxY),type="o", col="blue", ylab="val") # 状態
par(new=TRUE)
plot(df_res1$time, df_res1$y, ylim=c(minY,maxY), type="o", ylab="") # 観測値
par(new=TRUE)
plot(df_res1$time, df_res1$a, ylim=c(minY,maxY), type="l", col="green", ylab="") # 1期先状態予測
par(new=TRUE)
plot(df_res1$time, df_res1$f, ylim=c(minY,maxY), type="l", col="red", ylab="") # 1期先観測値予測
par(new=TRUE)
plot(df_res1$time, df_res1$m, ylim=c(minY,maxY), type="l", col="purple", ylab="") # フィルタリング
青（点と線）：θ、黒（点と線）：Y
緑：1期先状態予測、赤：1期先観測値予測、紫：フィルタリング
ローカルレベルモデルの場合は、1期先状態予測と1期先観測値予測は同じ（a[t] = f[t]）となる。よって、図では緑が赤で上書きされている。
また、1期先状態予測は1期前のフィルタリングと同じ（a[t] = m[t-1]）となる。
分散を見てみる。
時普遍のモデルの場合、C[t]はtが増えるに従って極値に近づくらしい。
以下の指定で、フィルタリング結果から分散（RとCのみ）を取り出せる。
R <- unlist(dlmSvd2var(fit_mod1$U.R, fit_mod1$D.R))
C <- unlist(dlmSvd2var(fit_mod1$U.C, fit_mod1$D.C))
df_res1 <- data.frame(df_res1, R=c(NA,R), C=C)
time y theta a f m R C 1 0 NA 10.00000 NA NA 10.000000 NA 50.000000 2 1 11.480221 11.12716 10.000000 10.000000 11.404956 56.000000 2.847458 3 2 14.887411 14.39119 11.404956 11.404956 14.005587 8.847458 2.240343 4 3 16.268663 14.81166 14.005587 14.005587 15.664657 8.240343 2.199313 5 4 15.192051 14.39024 15.664657 15.664657 15.318650 8.199313 2.196379 6 5 7.640275 11.47224 15.318650 15.318650 9.697648 8.196379 2.196169 7 6 11.918582 12.24209 9.697648 9.697648 11.323486 8.196169 2.196154 8 7 11.739846 13.41268 11.323486 11.323486 11.628283 8.196154 2.196153 9 8 19.019994 18.91138 11.628283 11.628283 17.039391 8.196153 2.196152 10 9 21.572069 19.41500 17.039391 17.039391 20.357542 8.196152 2.196152 11 10 20.391132 17.97181 20.357542 20.357542 20.382131 8.196152 2.196152 12 11 15.116908 18.52915 20.382131 20.382131 16.527721 8.196152 2.196152 13 12 19.366015 21.61766 16.527721 16.527721 18.605496 8.196152 2.196152 14 13 21.751131 19.29930 18.605496 18.605496 20.908260 8.196152 2.196152 15 14 16.585866 17.66393 20.908260 20.908260 17.744048 8.196152 2.196152 16 15 17.432607 16.43332 17.744048 17.744048 17.516058 8.196152 2.196152 17 16 22.007343 19.82772 17.516058 17.516058 20.803907 8.196152 2.196152 18 17 18.873734 17.81004 20.803907 20.803907 19.390922 8.196152 2.196152 19 18 19.547199 19.00535 19.390922 19.390922 19.505325 8.196152 2.196152 20 19 17.828754 19.51896 19.505325 19.505325 18.277990 8.196152 2.196152 21 20 23.217935 22.09485 18.277990 18.277990 21.894281 8.196152 2.196152
t>=3 くらいからC[t]はほぼ極値となったようである。
状態の1期先予測の分散であるR[t]はY[t]が与えられた後に計算されるC[t]より大きくなる。
k期先予測（k>1）を考える。
Y[1:t]が与えられたときのθ[t+k]やY[t+k]の分布の予測に関しては後述（k>1）。
t時点の状態の分布、θ[t]|Y[t]～N(m[t], C[t])、が与えられているとする。
m[t] = a[t]、C[t] = R[t]と置き換え、カルマンフィルタで見たように、状態と観測値の推定が繰り返される。
(1) Y[1:t]が与えられた元でのk期先状態予測分布
θ[t+k]|Y[1:t] ～ N(a[t+k], R[t+k])
a[t+k] = G[t+k] a[t+k-1]
R[t+k] = G[t+k] R[t+k-1] G[t+k]’ + W[t+k]
k期先の状態の予測分布においては、1期先で見た式と比べると、m、Cがa、Rに置き換わっただけである。要するに、推定されないフィルタリング分布を自身の状態予測分布で置き換え、逐次推定を繰り返している。
(2) Y[1:t]が与えられた元でのk期先観測予測分布
Y[t+k]|Y[1:t] ～ N(f[t+k], Q[t+k])
f[t+k] = F[t+k] a[t+k]
Q[t+k] = F[t+k] R[t+k] F[t+k]’ + V[t+k]
以上が定義される。
dlmForecast(あてはめたモデル, 予測期間)により、k期先状態予測、k期先観測値予測、フィルタリングの値を計算する。
k=10とする。
k <- 10 # k=10
pred_mod1 <- dlmForecast(fit_mod1, nAhead=k)
str(pred_mod1, 1)
List of 4 $ a: num [1:10, 1] 21.9 21.9 21.9 21.9 21.9 ... $ R:List of 10 $ f: num [1:10, 1] 21.9 21.9 21.9 21.9 21.9 ... $ Q:List of 10
df_pred1 <- data.frame(time=21:30, a=pred_mod1$a, R=unlist(pred_mod1$R), f=pred_mod1$f, Q=unlist(pred_mod1$Q))
# df_res1とdf_pred1の必要項目を結合
tm <- c(df_res1$time, df_pred1$time)
tt <- c(df_res1$theta, rep(NA,nrow(df_pred1)))
y <- c(df_res1$y, rep(NA,nrow(df_pred1)))
m <- c(df_res1$m, rep(NA,nrow(df_pred1)))
a <- c(df_res1$a, df_pred1$a)
R <- c(df_res1$R, df_pred1$R)
f <- c(df_res1$f, df_pred1$f)
Q <- c(rep(NA,nrow(df_res1)), df_pred1$Q)
df_all1 <- data.frame(time=tm, theta=tt, y=y, m=m, a=a, R=R, f=f, Q=Q)
time theta y m a R f Q 1 0 10.00000 NA 10.000000 NA NA NA NA 2 1 11.12716 11.480221 11.404956 10.000000 56.000000 10.000000 NA 3 2 14.39119 14.887411 14.005587 11.404956 8.847458 11.404956 NA 4 3 14.81166 16.268663 15.664657 14.005587 8.240343 14.005587 NA 5 4 14.39024 15.192051 15.318650 15.664657 8.199313 15.664657 NA 6 5 11.47224 7.640275 9.697648 15.318650 8.196379 15.318650 NA 7 6 12.24209 11.918582 11.323486 9.697648 8.196169 9.697648 NA 8 7 13.41268 11.739846 11.628283 11.323486 8.196154 11.323486 NA 9 8 18.91138 19.019994 17.039391 11.628283 8.196153 11.628283 NA 10 9 19.41500 21.572069 20.357542 17.039391 8.196152 17.039391 NA 11 10 17.97181 20.391132 20.382131 20.357542 8.196152 20.357542 NA 12 11 18.52915 15.116908 16.527721 20.382131 8.196152 20.382131 NA 13 12 21.61766 19.366015 18.605496 16.527721 8.196152 16.527721 NA 14 13 19.29930 21.751131 20.908260 18.605496 8.196152 18.605496 NA 15 14 17.66393 16.585866 17.744048 20.908260 8.196152 20.908260 NA 16 15 16.43332 17.432607 17.516058 17.744048 8.196152 17.744048 NA 17 16 19.82772 22.007343 20.803907 17.516058 8.196152 17.516058 NA 18 17 17.81004 18.873734 19.390922 20.803907 8.196152 20.803907 NA 19 18 19.00535 19.547199 19.505325 19.390922 8.196152 19.390922 NA 20 19 19.51896 17.828754 18.277990 19.505325 8.196152 19.505325 NA 21 20 22.09485 23.217935 21.894281 18.277990 8.196152 18.277990 NA 22 21 NA NA NA 21.894281 8.196152 21.894281 11.19615 23 22 NA NA NA 21.894281 14.196152 21.894281 17.19615 24 23 NA NA NA 21.894281 20.196152 21.894281 23.19615 25 24 NA NA NA 21.894281 26.196152 21.894281 29.19615 26 25 NA NA NA 21.894281 32.196152 21.894281 35.19615 27 26 NA NA NA 21.894281 38.196152 21.894281 41.19615 28 27 NA NA NA 21.894281 44.196152 21.894281 47.19615 29 28 NA NA NA 21.894281 50.196152 21.894281 53.19615 30 29 NA NA NA 21.894281 56.196152 21.894281 59.19615 31 30 NA NA NA 21.894281 62.196152 21.894281 65.19615
ローカルレベルモデルなので、a[t+k] = f[t+k] = m[t]となる。ランダム項が与えられるのみで構成されるので、一度推定されたm[t]と、k期先状態と観測値予測は同じとなる。
分散は、期間が先になるほど広がって行くのが分かる。
状態と観測値予測をプロットしてみる。
# aとfの95%CI
ci_a <- cbind(df_all1$a+df_all1$R*qnorm(0.025), df_all1$a+df_all1$R*qnorm(0.975))
ci_f <- cbind(df_all1$f+df_all1$Q*qnorm(0.025), df_all1$f+df_all1$Q*qnorm(0.975))
minY <- min(na.omit(df_all1$theta))-30
maxY <-max(na.omit(df_all1$theta))+30
plot(df_all1$time, df_all1$theta, xlim=c(0,30), ylim=c(minY,maxY),type="o", col="blue", xlab="time", ylab="val") # 状態
par(new=TRUE)
plot(df_all1$time, df_all1$y, xlim=c(0,30), ylim=c(minY,maxY), type="o", xlab="", ylab="") # 観測値
par(new=TRUE)
plot(df_all1$time, df_all1$a, xlim=c(0,30), ylim=c(minY,maxY), type="l", col="green", xlab="", ylab="") # 状態予測
par(new=TRUE)
plot(df_all1$time, ci_a[,1], xlim=c(0,30), ylim=c(minY,maxY), type="l", lty=2, col="green", xlab="", ylab="")
par(new=TRUE)
plot(df_all1$time, ci_a[,2], xlim=c(0,30), ylim=c(minY,maxY), type="l", lty=2, col="green", xlab="", ylab="")
par(new=TRUE)
plot(df_all1$time, df_all1$f, xlim=c(0,30), ylim=c(minY,maxY), type="l", col="red", xlab="", ylab="") # 観測値予測
par(new=TRUE)
plot(df_all1$time, ci_f[,1], xlim=c(0,30), ylim=c(minY,maxY), type="l", lty=2, col="red", xlab="", ylab="")
par(new=TRUE)
plot(df_all1$time, ci_f[,2], xlim=c(0,30), ylim=c(minY,maxY), type="l", lty=2, col="red", xlab="", ylab="")
青（点と線）：θ、黒（点と線）：Y
緑：状態予測、赤：観測値予測
緑（破線）：状態予測の95%信頼区間、赤（破線）：観測値予測の95%信頼区間
ローカルレベルモデルの場合は、状態予測値と観測値予測値は同じとなるので、図では緑が赤で上書きされている。
もちろん、信頼区間は異なる。ただし単に観測値誤差（V=3）を足し引きしただけ。
観測値予測の信頼区間は予測区間（time>=21）からのみ。
では、最後にパラメータの推定を考える。
これまでは、モデルのパラメータは与えられたものとして見てきた。
通常のモデリングプロセスでは、パラメータを推定し、カルマンフィルタを実行し状態や予測値の推定をする流れとなるだろう。
ローカルレベルモデルの場合は、V、Wが推定するパラメータとなる。
推定するパラメータを引数とする関数を作成する。
その関数に対し、dlmMLE(データ, パラメータ初期値, 関数)を適用する。dlmMLE()は後ろでoptim()が動いている。
使い方は、関数を定義し、optim()を適用という流れと同じ。
前回記事参照。
関数ではdV=exp(V_estim)のように、expを取った。こうすると、最適値計算において負の値を探索しないらしい。
返り値はV_estim=log(aV)となるので、推定したパラメータの対数が返ってくる。
これからカルマンフィルタを実行するモデルにこれら推定したパラメータを代入する際は、expをとっておく。
まず、dlmMLE()を適用する関数を定義。
fn_mod1 <- function(parm) {
V_estim <- parm[1]
W_estim <- parm[2]
m.0 <- 10 # 状態の初期平均を仮設定
C.0 <- 50 # 状態の初期標準偏差を仮設定
# return( dlm( m0=m.0, C0=C.0, FF=1, V=V_estim, GG=1, W=W_estim ) ) # dlm()を使うとエラーが出る?!?
return( dlmModPoly( order=1, m0=m.0, C0=C.0, dV=exp(V_estim), dW=exp(W_estim) ) )
}
dlm()を使うと推定の際になぜかエラーが出る。原因は不明。なので、dlmModPoly()で定義しておく。
parm_start <- c( 10*runif(1), 10*runif(1) ) # 推定するパラメータ（V、W）の初期値を適当に設定
res2 <- dlmMLE(y=df1$y[-1], parm=parm_start, build=fn_mod1)
V_est <- exp(res2$par[1])
7.681681
W_est <- exp(res2$par[2])
2.406207
V=7.681681、W=2.406207と推定されたが、推定結果は真の値であるV=3、W=6と異なる。
データ数が20で少ないのかと思い、データ数を増やしたりしていろいろとシミュレーションしてみた。数百データがあるとだいたい真の値に近づく傾向があるようである。
次回以降、さまざまな時系列モデルを状態空間表現で定義し、この流れを見て行く。
続く。
R 最適化計算 - optimize(), optim()
optimize(), optim() の使い方に関するメモ。
（１）optimize()による1パラメータ推定
ベクトルXの平均値を（ややこしく）求めてみる。
X <- c(1, 2, 3, 5, 8)
mean(X) # 3.8
optimize()を適用し、パラメータの最適値を求める対象の関数をまず定義する。
find_mean <- function(x){
return(function(m){
return( sum( (x-m)^2 ) )
})
}
xを外側関数、mを内側関数が受け取る2重関数。
find_mean(x)(m)と実行すると計算結果が返ってくる。
ベクトルxを渡し、その平均mを最適化計算の関数（optimize、optim）で返す処理を行う場合に定義する。
sum( (x-m)^2 )は、xをベクトルとしたとき、x-mの二乗和を最小にするmを求めることは、xの平均値を求めることに等しくなる。
1-10の間でfind_mean(X)関数を最小とするmの解（Xの平均）を探す。
optimize(find_mean(X), interval=c(1,10), tol=0.0001, maximum=FALSE)
$minimum [1] 3.8 $objective [1] 30.8
minimunは推定したパラメータ値。
objectiveが、そのパラメータを代入したときの関数の値（最小値）。
（２）optim()による2パラメータ推定
最尤法により、正規分布の平均と分散を
Y <- rnorm(50)
mean(Y) # -0.02269269
sd(Y) # 1.032003
ベクトルx（正規分布を仮定）の平均mu、標準偏差sigの2つのパラメータを最尤推定する処理を行う場合に定義する。
find_norm_para <- function(x){
return(function(para){
mu <- para[1]
sig <- para[2]
# 対数尤度を返す
return( - length(x)/2 * log(sig) - 1/2 * sum( (x - mu)^2 ) / sig )
})
}
para[1]=mu=-1、para[2]=sig=3を初期値（適当に設定）として、それらの真の解（mean(Y)、sd(Y)）を探す。
求めるパラメータが2つ以上の場合、optim(初期値, 関数. ...)を使用。
optim(par=c(-1,3), fn=find_norm_para(Y), control=list(fnscale=-1))
parが適当に設定する初期値。ベクトルで与える。
control=list(fnscale = -1) で最大化を指示できる。デフォルトは最小化。
$par
[1] -0.02262334 1.04377987
$value
[1] -26.07003
$counts
function gradient
67 NA
$convergence
[1] 0
$message
NULL
parがパラメータの推定結果。
valueが推定したパラメータの場合の式の値（最大対数尤度）。
countsのfunctionが計算繰り返し回数。
convergenceには、収束したなら0が返る。
（３）optim()を使用したその他の例
data1
score1 score2 score3 ans 1 0 0 1 0 2 1 1 1 1 3 2 1 0 0 4 1 0 0 0 5 0 0 0 0 6 0 1 2 1 7 1 0 0 0 8 2 0 0 0 9 0 1 2 0 10 0 0 2 1 11 0 0 0 0
このようなdata1があったとする。
各スコアは0,1,2を取る、よってscore1+2+3の最大は6となる場合がある。
scoreの重み付き線形和とansの差の二乗和が最小となる重みを求めたいとする。
重みの範囲は[0,1]とする。
要するに、p1、p2、p3が推定するパラメータだったとき（max_score=6）、solとdata1のans列の差の2乗和を最小にする問題。
sol = (p1*score1 + p2*score2 + p3*score3) / max_score
sum( (sol-ans)^2 ) ⇒ 最小化
これまで見てきたように、optim_fn(データ)(パラメータ)形式の関数と定義する。
max_score <- 6
optim_fn <- function(data){
score1 <- data[,1]
score2 <- data[,2]
score3 <- data[,3]
ans <- data[,4]
return(function(para){
p1 <- para[1]
p2 <- para[2]
p3 <- para[3]
sol <- (p1*score1 + p2*score2 + p3*score3) / max_score
return( sum((sol-ans)^2) ) # want to minimize this
})
}
optim(par=para_list, fn=optim_fn(data1), method="L-BFGS-B", lower = 0, upper = 1)
methodに"L-BFGS-B"を指定すると、探索値の範囲指定（lower=0、upper=1）が行える。
$par
[1] 0.1818182 1.0000000 1.0000000
$value
[1] 1.434343
$counts
function gradient
6 6
$convergence
[1] 0
$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"
最適解がいくつあるか分からない問題（これはシンプルだが...）の場合、ランダムに複数回初期値を選んで計算させてみる。
10回実施。
n_try <- 10
n_para <- 3
res_list <- list()
count <- 1
for(i in 1:n_try){
para_list <- runif(n_para)
res <- optim(par=para_list, fn=optim_fn(data1), method="L-BFGS-B", lower = 0, upper = 1)
res_listcount <- c(para_list, res$par, res$value)
count <- count+1
}
結果をデータフレームにまとめる。
mtx1 <- matrix( unlist(res_list), ncol=7, byrow=TRUE )
res_data <- data.frame(st_para1=mtx1[,1], st_para2=mtx1[,2], st_para2=mtx1[,3],
end_para1=mtx1[,4], end_para2=mtx1[,5], end_para3=mtx1[,6], value=mtx1[,7])
res_data
st_para1 st_para2 st_para2.1 end_para1 end_para2 end_para3 value 1 0.15277128 0.14019332 0.3100025 0.1818182 1 1 1.434343 2 0.96345322 0.29611349 0.6449935 0.1818182 1 1 1.434343 3 0.57237208 0.52925936 0.3130487 0.1818182 1 1 1.434343 4 0.52376898 0.34242409 0.1968113 0.1818182 1 1 1.434343 5 0.13490861 0.70259612 0.9709726 0.1818182 1 1 1.434343 6 0.42852674 0.74335816 0.6959915 0.1818182 1 1 1.434343 7 0.55659350 0.42611310 0.5567565 0.1818182 1 1 1.434343 8 0.63883635 0.99873597 0.5545664 0.1818182 1 1 1.434343 9 0.11357835 0.05827815 0.1648482 0.1818182 1 1 1.434343 10 0.08554884 0.15794897 0.7886483 0.1818182 1 1 1.434343
end_para1、end_para2、end_para3が推定した結果となる。今回は10回とも同じ。
以上。
R データフレームのFactor型変数の水準順序の入れ替え方法
スマートなやり方かどうかは分からないが、データフレームのFactor型変数の水準順序の入れ替え方法に関するメモ。
季節気温の疑似データを作成。
Season=c( rep("春",5), rep("夏",3), rep("秋",4), rep("冬",3) )
Temp=c( c(20,21,15,15,21), c(32,33,35), c(16,10,20,9), c(0,7,9) )
( TempData <- data.frame(Season, Temp) )
Season Temp 1 春 20 2 春 21 3 春 15 4 春 15 5 春 21 6 夏 32 7 夏 33 8 夏 35 9 秋 16 10 秋 10 11 秋 20 12 秋 9 13 冬 0 14 冬 7 15 冬 9
季節ごとに平均気温を集計。
tapply(TempData$Temp, TempData$Season, mean)
夏 秋 春 冬 33.333333 13.750000 18.400000 5.333333
このように、Season変数は文字コード順に従い、"夏","秋","春","冬"順で分析結果が出力される。
当然、"春","夏","秋","冬"順で分析結果を得たい。
意図的に水準の順序を入れたデータフレームを作成するかたちで対処。
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/16.html
を参考に、factor関数でlevels引数を指定。
Season <- factor(TempData$Season, levels=c("春","夏","秋","冬"))
TempData2として新しいデータフレームを作成。
TempData2 <- data.frame(Season, Temp=TempData$Temp)
tapply(TempData2$Temp, TempData2$Season, mean)
春 夏 秋 冬 18.400000 33.333333 13.750000 5.333333
R データ確認（一変量集計）のための関数1
一変量集計や分布、度数のプロットが一度にできないとめんどくさいので、それ用の関数をとりあえず作成。
色々なデータではまだ試していません。
dataCheck1 <- function(Data, ngr=3){
print( "データのヘッド部分, head(Data)" )
print( head(Data) )
print( "データの構造, str(Data)" )
print( str(Data) )
print("データの要約統計量, summary(Data)")
print( summary(Data) )
# 変数の数を取得
r <- length(Data)
# 変数名を取得
vname <- names(Data)
index1 <- ngr
par( mfrow=c(ceiling(r/index1), index1) )
on.exit(par( mfrow=c(1,1) )) #関数が終了するときの処理
for(i in 1:r){
if( is.factor(Data[,i]) == TRUE ){
barplot( table(Data[,i]), main=vname[i] )
}else{
hist( Data[,i], xlab="", main=vname[i] )
}
}
}
dataCheck1(データ, ngr=プロットを横に何列並べるか) で実行。
dataCheck1(iris)
ngrは指定しなければ、デフォルトで3になる。
[1] "データのヘッド部分, head(Data)" Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa [1] "データの構造, str(Data)" 'data.frame': 150 obs. of 5 variables: $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ... $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ... $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ... $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ... $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ... NULL [1] "データの要約統計量, summary(Data)" Sepal.Length Sepal.Width Petal.Length Petal.Width Species Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100 setosa :50 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300 versicolor:50 Median :5.800 Median :3.000 Median :4.350 Median :1.300 virginica :50 Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800 Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
