一昨日ぐらいにカルマンフィルターを使った時系列データのトレンド変化の要因分析をかけたみたわけですが、ある系列のトレンド（状態変数）を、観察可能な外生変数ではなく、別の系列のトレンド（状態変数）で回帰するケースを考えてみたいと思います。
データ生成
2系統の状態変数ベクトル
と
に
と言う関係があると設定して、データを生成します。
は時点、
は誤差項です。
状態変数ベクトル
は観察できず、それにノイズが乗った観測値
のみが分析に使えます。
set.seed(505)
n <- 100
Gt <- matrix(c(5, 0, 0, 0.5), 2, 2)
y <- x <- matrix(NA, 2, n)
x[1, 1] <- 0
x[2, ] <- 1 + sin(2*pi/n*(1:n))
y[, 1] <- 0
beta <- 0.75
for(t in 2:n){
x[1, t] <- x[1, t - 1] + beta*x[2, t - 1] + rnorm(1)
y[, t] <- x[, t] + Gt %*% rnorm(2)
}
y <- t(y)
2系列同時カルマンフィルター
本題に入る前の準備として、2系列を同時に同じ行列で処理してローカル線形トレンドを取り出せることを確認します。
library(KFAS)
Zt <- matrix(c(1, 0, 0, 0,
0, 0, 1, 0), 2, 4, byrow = TRUE)
Tt <- matrix(c(1, 1, 0, 0,
0, 1, 0, 0,
0, 0, 1, 1,
0, 0, 0, 1), 4, 4, byrow = TRUE)
Ht <- diag(c(NA, NA))
Qt <- diag(c(0, NA, 0, NA))
model_IND <- SSModel(y ~ -1 + SSMcustom(
Z = Zt, T = Tt, R = diag(4), Q = Qt,
a1 = matrix(c(0, 0, 0, 0), 4, 1),
P1 = diag(4) * 10
),
H = Ht
)
m <- sum(is.na(Qt)) + sum(is.na(Ht))
fit_IND <- fitSSM(model_IND, inits = rep(1, m), method = "L-BFGS-B")
r_IND <- KFS(fit_IND$model, filtering = c("state", "signal"), smoothing = c("state", "signal"))
plot_order <- c(3, 1, 5, 4, 2, 6)
plotted_series <- cbind(y, t(x), r_IND$att[, c(1, 3)])[, plot_order]
plot_cnames <- c("観測値y₁", "観測値y₂", "観察不能な状態変数x₁", "観察不能な状態変数x₂", "推定されたαₜ₁", "推定されたαₜ₃")[plot_order]
colnames(plotted_series) <- plot_cnames
plot(plotted_series, main = "2系列に同時にカルマンフィルターによるローカル線形トレンド抽出")
プロットの添付は省略します。
時系列トレンドを時系列トレンドで回帰
さて本題です。前節の基本的な処理をちょっといじります。遷移行列
の成分に
を入れます。
カルマンフィルターによっては逐次推定できませんが、KFASパッケージは最尤法カルマンフィルターなので、やや変則処理になりますが推定できます。
Tt_CAU <- Tt
Tt_CAU[1, 3] <- NA
model_CAU <- SSModel(y ~ -1 + SSMcustom(
Z = Zt, T = Tt_CAU, R = diag(4), Q = Qt,
a1 = matrix(c(0, 0, 0, 0), 4, 1),
P1 = diag(4) * 10
),
H = Ht
)
m <- sum(is.na(Qt)) + sum(is.na(Ht)) + sum(is.na(Tt_CAU))
fit_CAU <- fitSSM(model_CAU, inits = rep(1, m), method = "L-BFGS-B", updatefn = function(pars, model, ...){
model$Q[is.na(model$Q)] <- exp(pars[1:(m-3)])
model$H[is.na(model$H)] <- exp(pars[(m-2):(m-1)])
model$T[is.na(model$T)] <- pars[m]
return(model)
},
hessian = TRUE, control = list(trace = 3, maxit = 200, lmm = 10, factr = 1e1))
beta <- fit_CAU$optim.out$par
vcov <- -solve(fit_CAU$optim.out$hessian)
se <- sqrt(abs(diag(vcov)))
a <- 0.05
conf <- qnorm(c(a/2, 1 - a/2), mean = beta[m], sd = se[m])
print(sprintf("β̂ : %f (%.0f%% CI %f %f)", beta[m], 100*(1-a), conf[1], conf[2]))
r_CAU <- KFS(fit_CAU$model, filtering = c("state", "signal"), smoothing = c("state", "signal"))
plot_order <- c(3, 1, 5, 4, 2, 6)
plotted_series <- cbind(y, t(x), r_CAU$att[, c(1, 3)])[, plot_order]
colnames(plotted_series) <- plot_cnames
plot(plotted_series, main = "状態変数と状態変数の関係を推定した場合のトレンド")
トレンド自体は上手く抽出できていると言えるでしょう。観測不能な状態変数
の推定量が
に、
の推定量が
になりますが、だいたい似ています。
β̂ : 0.910916 (95% CI 0.433817 1.388015)
の推定量がそれらしくなるには、サンプルサイズがもう10倍ぐらいは必要そうです。収束自体はしているのですが、ヘッシアンの値が何やらおかしいことになっている*1ので、もやもや感をもって利用してください('-' )＼(--；)BAKI
疑似相関のケース
疑似相関のケースでの推定量がどうなるか確認しておきましょう。
DGPを以下のように書き換えて、
set.seed(505)
n <- 100
Gt <- matrix(c(5, 0, 0, 0.5), 2, 2)
y <- x <- matrix(NA, 2, n)
x[1, 1] <- 0
x[2, ] <- 1 + sqrt(1:n)
y[, 1] <- 0
beta <- 0.75
for(t in 2:n){
x[1, t] <- 1 + x[1, t - 1] + rnorm(1)
y[, t] <- x[, t] + Gt %*% rnorm(2)
}
y <- t(y)
（2系列同時カルマンフィルターの節のコードの変数宣言をしたあと）時系列トレンドを時系列トレンドで回帰するコードを走らせてみます。
とに関係があるように誤認してしまいそうですが、
β̂ : -0.030803 (95% CI -0.117191 0.055585)
と、疑似相関に騙されることなく、ほぼゼロの係数を推定できました。
まとめ
先に二つ試しましたが、これが一番、お題の解に近い気がしています。
実用上は、
- カルマンフィルターと言う良く知られた既存手法の応用
- 人気のRのパッケージKFASで計算できる
- 「トレンドをトレンドで回帰しました」と口頭で説明しやすい
- プロットが出るので可視化されている
- プロットに対応した推定量が出る
- 疑似相関にそれなり*2頑強
と言う利点があります。昔流行ったデータ可視化マンになって使ってみてください*3。
