R機械学習/カルマンフィルタによる状態空間モデル2

# Introduction
カルマンフィルタを使用した状態空間モデルをRで実装する方法についてのメモ．
今回は外生変数を組み込んだな時変係数モデルを実装する．
参考図書：[時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装](https://www.kinokuniya.co.jp/f/dsg-01-9784903814872)

# ライブラリインポート

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

# シミュレーションデータの作成
2つの外生変数と目的変数からなるデータを作成．
```r
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)
```

データの確認

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

Out：外生変数と目的変数
![1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/8adaf149-2165-8d65-0ceb-2c5d68f3171d.png)


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

Out：回帰係数と切片
![2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/3fa7f544-4f40-3006-a654-ab75d371a816.png)


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

```r
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$：観測誤差

```r
# 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'))
```
# モデル評価
訓練データの結果を図示．

```r
# グラフ用データフレーム作成
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：黒点が観測値，赤線が予測値
![3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/b687b510-f16c-b121-e100-6c7d76f6f19e.png)

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

```r
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：黒点が真値，赤線が推定値
![4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/26232154-2523-7da8-5abe-a6b5bdb39d24.png)

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

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

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

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

```r
# 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))
```

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

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

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

```r
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：紫垂直線の右がテストデータ
![5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/36617add-3614-16c8-a77a-29ebbe20cc42.png)


切片と回帰係数を図示

```r
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：紫垂直線の右がテストデータ
![6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/2aff4498-aaff-731f-8aa9-81af7ca995e0.png)

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

-------------
-------------

# Code

```r
# ライブラリインポート
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)
```