状態空間モデル - 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
