Introduction

カルマンフィルタを使用した状態空間モデルをRで実装する方法についてのメモ．
今回は外生変数を組み込んだな時変係数モデルを実装する．
参考図書：時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装

ライブラリインポート

library(KFAS)  # モデル作成
library(ggplot2)  # グラフ描画
library(patchwork)  # グラフ配置

シミュレーションデータの作成

2つの外生変数と目的変数からなるデータを作成．

set.seed(1)
N=100

# 各種標準偏差
sigma_v = 10
sigma_w = 5
sigma_t1 = 0.6
sigma_t2 = 0.2

# 説明変数
x1 = rnorm(n=N, mean=5, sd=10)
x2 = rnorm(n=N, mean=10, sd=5)

# 切片，回帰係数，観測値
mu = cumsum(rnorm(n=N, mean=0, sd=sigma_w))
b1 = cumsum(rnorm(n=N, mean=0, sd=sigma_t1)) + 10
b2 = cumsum(rnorm(n=N, mean=0, sd=sigma_t1)) + 5
y = mu + b1*x1 + b2*x2 + rnorm(n=N, mean=0, sd=sigma_v)

dataset = data.frame(y=y, x1=x1, x2=x2)
coefset = data.frame(b1=b1, b2=b2, mu=mu)

データの確認

# 外生変数と目的変数
par(mfrow=c(3,1))
plot(x1)
plot(x2)
plot(y)

Out：外生変数と目的変数

# 回帰係数と切片
par(mfrow=c(3,1))
plot(b1)
plot(b2)
plot(mu)

Out：回帰係数と切片

テストデータと訓練データに分割

p = 80
train = dataset[1:p,]
test = dataset[81:N,]

train_coef = coefset[1:p,]
test_coef = coefset[81:N,]

モデリング

2つの外生変数をもつ時変係数モデルは以下のように定式化される．

$$
\begin{eqnarray}
\beta_{1,t} &=& \beta_{1,t-1} + \tau_{1,t}, \ \ \ \ &\tau_{1,t}&～Normal(0,\sigma^2_{\tau_1}) \\
\beta_{2,t} &=& \beta_{2,t-1} + \tau_{2,t}, \ \ &\tau_{2,t}&～Normal(0,\sigma^2_{\tau_2}) \\
\mu_t &=& \mu_{t-1} + w_t, \ \ &w_t&～Normal(0,\sigma^2_w) \\
y_t &=& \mu_t + \beta_{1,t}x_{1,t} + \beta_{2,t}x_{2,t} + v_t, \ \ \ \ &v_t&～Normal(0,\sigma^2_v)
\end{eqnarray}
$$

$\beta_t$：回帰係数
$\mu_t$：切片
$x_t$：外生変数
$y_t$：目的変数
$\tau_t$：回帰係数誤差
$w_t$：過程誤差
$v_t$：観測誤差

# Step1:モデルの構造を決める
build_kfas = SSModel(H = NA, train$y ~ SSMtrend(degree = 1, Q = NA) +
             SSMregression( ~ train$x1 + train$x2 , Q = diag(rep(NA, 2))))

# Step2:パラメータ推定
fit_kfas = fitSSM(build_kfas, inits = c(1, 1, 1, 1))

# Step3:フィルタリング
result_kfas = KFS(fit_kfas$model, filtering=c('state', 'mean'), smoothing=c('state', 'mean'))

モデル評価

訓練データの結果を図示．

# グラフ用データフレーム作成
df_kfas = data.frame(time=1:p, y=train$y, y_pred=result_kfas$m[1:length(train$y)])
df_kfas$y_pred[1] = NA

ggplot(data=df_kfas, aes(x=time, y=y))+
  geom_point(alpha=0.5)+
  geom_line(aes(y=y_pred),size=1.2, col = 'red')

Out：黒点が観測値，赤線が予測値

切片および回帰係数の推定結果

df_kfas_coef = data.frame(time=1:p, b1=train_coef$b1, b2=train_coef$b2, mu=train_coef$mu,
                          b1_pred = coef(result_kfas)[,1],
                          b2_pred = coef(result_kfas)[,2],
                          mu_pred = coef(result_kfas)[,3])

g1 = ggplot(data=df_kfas_coef, aes(x=time, y=b1)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=b1_pred), size=1, col='red')

g2 = ggplot(data=df_kfas_coef, aes(x=time, y=b2)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=b2_pred), size=1, col='red')

g3 = ggplot(data=df_kfas_coef, aes(x=time, y=mu)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=mu_pred), size=1, col='red')

g1+g2+g3+plot_layout(ncol=1)

Out：黒点が真値，赤線が推定値

逐次更新

訓練データで推定したパラメータを使用し，テストデータにおいても状態の更新をしたい．
少し強引な方法ではあるがKFASパッケージを使って逐次更新する方法について記す．

データ準備
テストデータのサイズ分訓練データ内の目的変数にNAを代入する．

train2 = dataset
train2[81:N,1] = NA

KFASパッケージでパラメタの推定まではこれまでと同様に行う．

# Step1:モデルの構造を決める
build_kfas2 = SSModel(H = NA, train2$y ~ SSMtrend(degree = 1, Q = NA) +
              SSMregression( ~ train2$x1 + train2$x2 , Q = diag(rep(NA, 2))))

# Step2:パラメータ推定
fit_kfas2 = fitSSM(build_kfas2, inits = c(1, 1, 1, 1))

パラメータ推定結果を格納したオブジェクトに全データを格納しフィルタリングを行う．

fit_kfas2[["model"]][["y"]][1:N] = dataset$y
# Step3:フィルタリング
result_kfas2 = KFS(fit_kfas2$model, filtering=c('state', 'mean'), smoothing=c('state', 'mean'))

目的変数と外生変数を図示

df_kfas2 = data.frame(time=1:N, y=dataset$y, y_pred=result_kfas2$m[1:N],
                       x1=dataset$x1, x2=dataset$x2)
df_kfas2$y_pred[1] = NA


g1 = ggplot(data=df_kfas2, aes(x=time, y=y))+
      geom_point(alpha=0.5)+
      geom_line(aes(y=y_pred),size=1.2, col = 'red') +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g2 = ggplot(data=df_kfas2, aes(x=time, y=x1))+
      geom_point(alpha=0.5) +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g3 = ggplot(data=df_kfas2, aes(x=time, y=x2))+
      geom_point(alpha=0.5) +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g1+g2+g3+plot_layout(ncol=1)

Out：紫垂直線の右がテストデータ

切片と回帰係数を図示

df_kfas_coef = data.frame(time=1:N, b1=coefset$b1, b2=coefset$b2, mu=coefset$mu,
                          b1_pred = coef(result_kfas2)[,1],
                          b2_pred = coef(result_kfas2)[,2],
                          mu_pred = coef(result_kfas2)[,3])

g1 = ggplot(data=df_kfas_coef, aes(x=time, y=b1)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=b1_pred), size=1, col='red') +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g2 = ggplot(data=df_kfas_coef, aes(x=time, y=b2)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=b2_pred), size=1, col='red') +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g3 = ggplot(data=df_kfas_coef, aes(x=time, y=mu)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=mu_pred), size=1, col='red') +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g1+g2+g3+plot_layout(ncol=1)

Out：紫垂直線の右がテストデータ

Conclusion

今回は時変係数モデルの構築，予測，逐次更新を実装した．
状態空間モデルは時系列データを柔軟に表現できるモデルとされているため，今後も学習を続けていこうと思う．

Code

# ライブラリインポート
library(KFAS)  # モデル作成
library(ggplot2)  # グラフ描画
library(patchwork)  # グラフ配置

# シミュレーションデータの作成
set.seed(1)
N=100

sigma_v = 10
sigma_w = 5
sigma_t1 = 0.6
sigma_t2 = 0.2

x1 = rnorm(n=N, mean=5, sd=10)
x2 = rnorm(n=N, mean=10, sd=5)

mu = cumsum(rnorm(n=N, mean=0, sd=sigma_w))
b1 = cumsum(rnorm(n=N, mean=0, sd=sigma_t1)) + 10
b2 = cumsum(rnorm(n=N, mean=0, sd=sigma_t1)) + 5
y = mu + b1*x1 + b2*x2 + rnorm(n=N, mean=0, sd=sigma_v)

dataset = data.frame(y=y, x1=x1, x2=x2)
coefset = data.frame(b1=b1, b2=b2, mu=mu)


par(mfrow=c(3,1))
plot(x1)
plot(x2)
plot(y)

par(mfrow=c(3,1))
plot(b1)
plot(b2)
plot(mu)

p = 80
train = dataset[1:p,]
test = dataset[81:N,]

train_coef = coefset[1:p,]
test_coef = coefset[81:N,]

# モデリング
build_kfas = SSModel(H = NA, train$y ~ SSMtrend(degree = 1, Q = NA) +
             SSMregression( ~ train$x1 + train$x2 , Q = diag(rep(NA, 2))))

fit_kfas = fitSSM(build_kfas, inits = c(1, 1, 1, 1))

result_kfas = KFS(fit_kfas$model, filtering=c('state', 'mean'), smoothing=c('state', 'mean'))

# モデル評価
df_kfas = data.frame(time=1:p, y=train$y, y_pred=result_kfas$m[1:length(train$y)])
df_kfas$y_pred[1] = NA

ggplot(data=df_kfas, aes(x=time, y=y))+
  geom_point(alpha=0.5)+
  geom_line(aes(y=y_pred),size=1.2, col = 'red')

df_kfas_coef = data.frame(time=1:p, b1=train_coef$b1, b2=train_coef$b2, mu=train_coef$mu,
                          b1_pred = coef(result_kfas)[,1],
                          b2_pred = coef(result_kfas)[,2],
                          mu_pred = coef(result_kfas)[,3])

g1 = ggplot(data=df_kfas_coef, aes(x=time, y=b1)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=b1_pred), size=1, col='red')

g2 = ggplot(data=df_kfas_coef, aes(x=time, y=b2)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=b2_pred), size=1, col='red')

g3 = ggplot(data=df_kfas_coef, aes(x=time, y=mu)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=mu_pred), size=1, col='red')

g1+g2+g3+plot_layout(ncol=1)

# 逐次更新
train2 = dataset
train2[81:N,1] = NA

build_kfas2 = SSModel(H = NA, train2$y ~ SSMtrend(degree = 1, Q = NA) +
              SSMregression( ~ train2$x1 + train2$x2 , Q = diag(rep(NA, 2))))

fit_kfas2 = fitSSM(build_kfas2, inits = c(1, 1, 1, 1))

fit_kfas2[["model"]][["y"]][1:N] = dataset$y
result_kfas2 = KFS(fit_kfas2$model, filtering=c('state', 'mean'), smoothing=c('state', 'mean'))

df_kfas2 = data.frame(time=1:N, y=dataset$y, y_pred=result_kfas2$m[1:N],
                       x1=dataset$x1, x2=dataset$x2)
df_kfas2$y_pred[1] = NA


g1 = ggplot(data=df_kfas2, aes(x=time, y=y))+
      geom_point(alpha=0.5)+
      geom_line(aes(y=y_pred),size=1.2, col = 'red') +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g2 = ggplot(data=df_kfas2, aes(x=time, y=x1))+
      geom_point(alpha=0.5) +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g3 = ggplot(data=df_kfas2, aes(x=time, y=x2))+
      geom_point(alpha=0.5) +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g1+g2+g3+plot_layout(ncol=1)

df_kfas_coef = data.frame(time=1:N, b1=coefset$b1, b2=coefset$b2, mu=coefset$mu,
                          b1_pred = coef(result_kfas2)[,1],
                          b2_pred = coef(result_kfas2)[,2],
                          mu_pred = coef(result_kfas2)[,3])

g1 = ggplot(data=df_kfas_coef, aes(x=time, y=b1)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=b1_pred), size=1, col='red') +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g2 = ggplot(data=df_kfas_coef, aes(x=time, y=b2)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=b2_pred), size=1, col='red') +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g3 = ggplot(data=df_kfas_coef, aes(x=time, y=mu)) +
      geom_point(alpha=0.5) +
      geom_line(aes(y=mu_pred), size=1, col='red') +
      geom_vline(xintercept=p, col = 'purple', linetype='dashed')

g1+g2+g3+plot_layout(ncol=1)