状態空間モデル - R dlm （５）
今回は回帰モデルを扱う。これを用いることにより、説明変数を時系列データのモデリングに含めることができるようになる。
x1 <- rnorm(50)
x2 <- 0.3*x1 + rnorm(50)
y1 <- 10*x1 + 5*x2 + rnorm(50)
d1 <- data.frame(time=1:50, x1, x2, y1)
par(mfrow = c(3, 1))
plot(x=d1$time, y=d1$x1, type="o")
plot(x=d1$time, y=d1$x2, type="o")
plot(x=d1$time, y=d1$y1, type="o")
このような時系列データ（x1、x2、y1）があった場合、x1とx2によって、y1を説明したい場合に有用である。
まず簡単な、回帰モデル（説明変数2つ）で説明。上のデータの場合に対応。
以下が、説明変数2つの回帰モデル。
y[t] = β0[t] +β1[t] x1[t] +β2[t] x2[t] + v[t], v[t]～N(0, σ2)
βi[t] = βi[t-1] + wi[t], wi～N(0, w_i) i=0,1,2
x1、x2は、通常の回帰モデルと同様、確率変数でない。通常の回帰モデルと、この時系列の回帰モデルの違いは、係数β1、β1、β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) 状態
例。
dlmModReg(X=d1[2:3], dV=3, dW=c(1,4,5))
$FF
[,1] [,2] [,3]
[1,] 1 1 1
$V
[,1]
[1,] 3
$GG
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
$W
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 4 0
[3,] 0 0 5
$JFF
[,1] [,2] [,3]
[1,] 0 1 2
$X
x1 x2
[1,] -0.7554 -0.5918
[2,] -1.32 1.552
[3,] ...
$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
では、ローカルレベルモデル + 回帰モデル（説明変数2つ）の場合は。
回帰モデルのβ0が、ローカルレベルモデルに対応するので実際には意味がないが、動的線形モデル表現の練習のために書いてみる。
y[t] = μ[t] + β0[t] +β1[t] x1[t] +β2[t] x2[t] + v[t], v[t]～N(0, σ2)
μ[t] = μ[t] + w[t], w[t] ～N(0, w_l)
βi[t] = βi[t-1] + wi[t], wi[t]～N(0, w_i) i=0,1,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) 状態
これまで見てきたのと同じように、要素の分解を行うことができる。Fの一番左の要素、GとWは一番左上の要素がローカルレベルモデル部分に対応する。
例。
dlmModPoly(order=1, dV=0.5, dW=0.75) + dlmModReg(X=d1[2:3], dV=3, dW=c(1,4,5))
$FF
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
$V
[,1]
[1,] 3.5
$GG
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 0 1 0
[4,] 0 0 0 1
$W
[,1] [,2] [,3] [,4]
[1,] 0.75 0 0 0
[2,] 0.00 1 0 0
[3,] 0.00 0 4 0
[4,] 0.00 0 0 5
$JFF
[,1] [,2] [,3] [,4]
[1,] 0 0 1 2
$X
x1 x2
[1,] -0.7554 -0.5918
[2,] -1.32 1.552
[3,] ...
$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
$Fの場合は左から一つ目部分、$GGや$Wは一番左上の要素がローカルモデル部分である。
Rの出力では、$FFは、回帰の説明変数になっている部分は1となっている。回帰を含めると$JFFというのが出てくるが、詳細は不明。
では、季節要素モデル + 回帰モデル（説明変数2つ）でサンプルデータを作成してみる。
以下のモデルから出力されたサンプルデータを作成とする。
dlmModSeas(frequency=4, dV=0.5, dW=c(0.75,0,0)) + dlmModReg(X=d1[2:3], dV=1.5, dW=c(1,1.6,2.3))
サンプル数は50。
alpha1 <- 3
alpha2 <- 9
alpha3 <- -4
alpha4 <- -alpha1 -alpha2 -alpha3 # -8
seas <- c(alpha1, alpha2, alpha3, alpha4)
V <- 1.5+0.5
w_seas <- 0.75
w_reg <- c(1,1.6,2.3)
tm <- c(0)
alpha <- c(NA)
beta0 <- c(0)
beta1 <- c(0)
beta2 <- c(0)
y <- c(NA)
# 50サンプル
for(i in 1:50){
tm <- append(tm, i)
alpha_next <- seas[i%%4+1] + rnorm(1, mean=0, sd=sqrt(w_seas))
alpha <- append(alpha, alpha_next)
beta0_next <- beta0[i] + rnorm(1, mean=0, sd=sqrt(w_reg[1]))
beta0 <- append(beta0, beta0_next)
beta1_next <- beta1[i] + rnorm(1, mean=0, sd=sqrt(w_reg[2]))
beta1 <- append(beta1, beta1_next)
beta2_next <- beta2[i] + rnorm(1, mean=0, sd=sqrt(w_reg[3]))
beta2 <- append(beta2, beta2_next)
y_next <- alpha_next + beta0_next + beta1_next*x1[i] + beta2_next*x2[i] + rnorm(1, mean=0, sd=sqrt(V))
y <- append(y, y_next)
}
x1a <- c(NA,x1)
x2a <- c(NA,x2)
df1 <- data.frame(time=tm, alpha=alpha, beta0=beta0, x1=x1a, beta1=beta1, x2=x2a, beta2=beta2, y=y)
プロット。
# 4*2分割でプロット
par(mfrow = c(4, 2))
plot(x=df1$time, y=df1$alpha, type="o")
plot(x=df1$time, y=df1$x1, type="o")
plot(x=df1$time, y=df1$beta0, type="o")
abline(h=0, lty=2)
plot(x=df1$time, y=df1$x2, type="o")
plot(x=df1$time, y=df1$beta1, type="o")
abline(h=0, lty=2)
plot(x=df1$time, y=df1$y, type="o")
plot(x=df1$time, y=df1$beta2, type="o")
abline(h=0, lty=2)
これまでの流れの通り、関数を定義し、MLEでパラメータ推定。
# fn定義
fn_mod1 <- function(parm) {
V <- parm[1]
w_seas <- parm[2]
w_0 <- parm[3]
w_1 <- parm[4]
w_2 <- parm[5]
return( dlmModSeas(frequency=4, dV=exp(V), dW=c(exp(w_seas),0,0)) + dlmModReg(X=d1[2:3], dV=0, dW=c(exp(w_0),exp(w_1),exp(w_2))) )
}
# parm = c( log(V), log(w_seas), log(w_0), log(w_1), log(w_2) )
# X: n*2
# 推定
parm_start <- c(10*runif(1), 10*runif(1), 10*runif(1), 10*runif(1), 10*runif(1))
res_fn_mod1 <- dlmMLE(y=df1$y[-1], parm=parm_start, build=fn_mod1)
( V_est <- exp(res_fn_mod1$par[1]) ) # ans: 1.5+0.5=2
1.355866
( w_seas_est <- exp(res_fn_mod1$par[2]) ) # ans: 0.75
0.4222892
( w_0_est <- exp(res_fn_mod1$par[3]) ) # ans: 1
0.1048544
( w_1_est <- exp(res_fn_mod1$par[4]) ) # ans: 1.6
2.159517
( w_2_est <- exp(res_fn_mod1$par[5]) ) # ans: 2.3
0.3573823
今回はエラーも出ることなく値も安定しており、また、推定値も真の値から大きく異なることない結果が返ってくる。β0とβ2がやや小さいようである。
推定したパラメータによるカルマンスムーサーの実行。
# 推定したパラメータを関数に代入
mod1_est <- fn_mod1( res_fn_mod1$par )
# フィルタリング
fit_mod1_est <- dlmSmooth(dropFirst(df1$y), mod1_est)
str( fit_mod1_est, 1 )
List of 3 $ s : num [1:51, 1:6] 4.27 10.29 -5.92 -8.64 4.45 ... $ U.S:List of 51 $ D.S: num [1:51, 1:6] 2.81 2.44 2.21 1.85 1.64 ...
必要なのは$sの1、4、5、6列。
alpha_est <- fit_mod1_est$s[,1]
beta0_est <- fit_mod1_est$s[,4]
beta1_est <- fit_mod1_est$s[,5]
beta2_est <- fit_mod1_est$s[,6]
yの推定値を計算。
y_est <- alpha_est + beta0_est + beta1_est*x1a + beta2_est*x2a
データフレームにまとめる。
df2 <- data.frame( df1[1:2], alpha_est, df1[3], beta0_est, df1[4:5], beta1_est, df1[6:7], beta2_est, df1[8], y_est )
head(df2)
time alpha alpha_est beta0 beta0_est x1 beta1 1 0 NA 4.272656 0.000000 0.4634815 NA 0.0000000 2 1 8.138018 10.286522 -1.224692 0.4634815 0.59047979 0.2163145 3 2 -3.478060 -5.923984 -1.498450 0.4184211 -0.08259835 0.6751069 4 3 -7.006637 -8.635194 -1.214040 0.4127776 0.04707254 -0.1822476 5 4 2.311228 4.454132 1.374272 0.4162269 1.45996620 -0.8441979 6 5 8.872011 9.764824 2.155021 0.4518083 0.16616192 -0.4990168 beta1_est x2 beta2 beta2_est y y_est 1 0.2835804 NA 0.0000000 1.388990 NA NA 2 0.2835805 -0.6922756 -2.4077708 1.388990 10.538561 9.955888 3 -0.2644056 0.8080324 -0.6270543 1.495312 -4.785160 -4.275463 4 -0.8794455 0.2187214 0.8822229 1.710190 -8.007337 -7.889760 5 -1.4856702 0.2016472 2.2524059 1.931846 2.675382 3.090882 6 -1.1257253 0.1435495 3.2286180 2.175587 10.438714 10.341884
プロット。
# 4*2分割でプロット
par(mfrow = c(4, 2))
yl = c(min(na.omit(df2$alpha))-1, max(na.omit(df2$alpha))+1)
plot(x=df2$time, y=df2$alpha, ylim=yl, type="o")
par(new=TRUE)
plot(x=df2$time, y=df2$alpha_est, ylim=yl, type="o", col="blue", ylab="")
plot(x=df2$time, y=df2$x1, type="o")
yl = c(min(na.omit(df2$beta0))-1, max(na.omit(df2$beta0))+1)
plot(x=df2$time, y=df2$beta0, ylim=yl, type="o")
par(new=TRUE)
plot(x=df2$time, y=df2$beta0_est, ylim=yl, type="o", col="blue", ylab="")
abline(h=0, lty=2)
plot(x=df2$time, y=df2$x2, type="o")
yl = c(min(na.omit(df2$beta1))-1, max(na.omit(df2$beta1))+1)
plot(x=df2$time, y=df2$beta1, ylim=yl, type="o")
par(new=TRUE)
plot(x=df2$time, y=df2$beta1_est, ylim=yl, type="o", col="blue", ylab="")
abline(h=0, lty=2)
yl = c(min(na.omit(df2$y))-1, max(na.omit(df2$y))+1)
plot(x=df2$time, y=df2$y, ylim=yl, type="o")
par(new=TRUE)
plot(x=df2$time, y=df2$y_est, ylim=yl, type="o", col="blue", ylab="")
yl = c(min(na.omit(df2$beta2))-1, max(na.omit(df2$beta2))+1)
plot(x=df2$time, y=df2$beta2, ylim=yl, type="o")
par(new=TRUE)
plot(x=df2$time, y=df2$beta2_est, ylim=yl, type="o", col="blue", ylab="")
abline(h=0, lty=2)
青が推定値となる。
前に季節要素モデルを取り上げたときと同様。季節要素alphaはだいたい綺麗に推定できるようである。
β0とβ2の推定された分散が真の分散より小さかったので、実際より滑らかになっている。ただし、傾向はそれなりに捉えているように思う。
