読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

状態空間モデル - R dlm (5)

今回は回帰モデルを扱う。これを用いることにより、説明変数を時系列データのモデリングに含めることができるようになる。

 

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")

f:id:High_School_Student:20141221145326j:plain

このような時系列データ(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)                     状態

f:id:High_School_Student:20141221144609j:plain

 

例。

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:id:High_School_Student:20141221144855j:plain

 

これまで見てきたのと同じように、要素の分解を行うことができる。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)

f:id:High_School_Student:20141221150401j:plain

これまでの流れの通り、関数を定義し、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)

f:id:High_School_Student:20141221151456j:plain

が推定値となる。

 

前に季節要素モデルを取り上げたときと同様。季節要素alphaはだいたい綺麗に推定できるようである。

 

β0とβ2の推定された分散が真の分散より小さかったので、実際より滑らかになっている。ただし、傾向はそれなりに捉えているように思う。