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

東京に棲む日々

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

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

 

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)                     状態

 

f:id:High_School_Student:20141130144533j:plain

 

この場合、状態式は、

θ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)                     状態

 

f:id:High_School_Student:20141130145009j:plain

 最後のWであるが、

f:id:High_School_Student:20141130145036j:plain

と計算できる。

このように、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)

f:id:High_School_Student:20141130145650j:plain

黒 : 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)

f:id:High_School_Student:20141130150138j:plain

理論と異なり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)

 f:id:High_School_Student:20141130150322j:plain

理論と異なり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データを切り落とした。

f:id:High_School_Student:20141130151023j:plain

 

モデルを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つ。

状態は以下。

f:id:High_School_Student:20141130151216j:plain

ちなみに、

f:id:High_School_Student:20141130151228j:plain

 

モデル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つ。

f:id:High_School_Student:20141130151410j:plain

 

モデルC、ARMA項はなし。

dlmModPoly( order=1, dV=V, dW=w_loc ) + dlmModSeas(frequency=4, dW=c(w_seas,0,0))

パラメータは3つ。

f:id:High_School_Student:20141130151459j:plain

 

また、モデル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)f:id:High_School_Student:20141130153102j:plain

:モデル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