基礎から分かる時系列分析 3.カルマンフィルタ(ローカルレベルモデルの実装)

[前回：基礎から分かる時系列分析 2.カルマンフィルタ (数理的な導出)](https://qiita.com/_R2_/items/76b078d6af7c1715b48b)

## 3.1モデルの定義

ローカルレベルモデルとは、次のようなランダムウォークを表す状態空間モデルのこと。遷移行列が単位行列。

```math 
\begin{cases}
{\bf x}_t&={\bf x}_{t-1}+{\bf x}_t\ \ \ \ {\bf w}_t\sim\mathcal{N}({\bf 0},W)\\
{\bf y}_t&={\bf x}_t+{\bf v}_t\ \ \ \ \ {\bf v}_t\sim\mathcal{N}({\bf 0},V)
\end{cases}
```

```R 
#---
#　モデルの定義
#---

mod = dlmModPoly(order = 1)
```

dlmライブラリのdlmModPoly()関数でモデルを定義した。この関数は多項式モデルを設定するもので、order=1とすることで次数1の多項式モデル(ローカルレベルモデル)となる。他にもシステムノイズや観測ノイズ、事前分布も設定できる。システムノイズと観測ノイズはモデルのパラメータであるから、

## 3.2 パラメータの特定

システムノイズと観測ノイズを最尤推定する。

```R
#---
#　パラメータの特定
#---

build_dlm = function(par){
  # モデルのパラメータを設定する関数
  
  mod$W[1,1] = exp(par[1]) # 負の領域を探索しないよう指数変換を行う
  mod$V[1,1] = exp(par[2])
  
  return(mod)
}

fit_dlm = dlmMLE(y = Nile, parm = c(0,0), build = build_dlm)
mod     = build_dlm(fit_dlm$par)
```

モデルを定義、構築する関数build_dlmを作り、dlmMLEでパラメータの最尤推定値を探索させ、推定結果をもとにモデルを構築している。fit_dlm$parが返す値は指数変換前の値なのでそのまま代入して良い。

## 3.3 フィルタリング

フィルタリングを実施する。

```R
#---
#　フィルタリングの実施
#---

dlmFiltered_obj = dlmFilter(y = Nile, mod = mod)

# フィルタ分布の平均、標準偏差の導出
m               = dropFirst(dlmFiltered_obj$m) # 先頭の事前分布の平均をドロップ
m_sdev          = sqrt(
  dropFirst( # 先頭を事前分布の分散をドロップ
    as.numeric(
      dlmSvd2var(dlmFiltered_obj$U.C, dlmFiltered_obj$D.C) # 特異値を分散に戻す
    ))
  )
m_quant         = list(m + qnorm(0.025, sd = m_sdev), m + qnorm(0.975, sd = m_sdev)) # 95%信頼区間

# 結果のプロット
oldpar = par(no.readonly = T) 
par(family = "HiraginoSans-W4")
ts.plot(cbind(Nile, m, do.call('cbind', m_quant)),
        col = c('lightgray', 'black', 'black', 'black'),
        lty = c('solid', 'solid', 'dashed', 'dashed'))
legend(legend = c('観測値', '平均(フィルタリング分布)', '95%区間'),
       lty    = c('solid', 'solid', 'dashed', 'dashed'),
       col    = c('lightgray', 'black', 'black', 'black'),
       x      = 'topright', text.width = 32, cex = 0.6)
par(oldpar)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/686685/67866472-09cc-cff3-2f97-05be0e1925e4.png)

dlmFilter()関数でフィルタリングを行う。この関数は観測値yとモデルmodの他に

- m：フィルタ分布の平均
- U.C.,D.C.：フィルタ分布の共分散行列の特異値分解。dlmSvd2varで共分散行列に戻る。

などを返す。

## 3.4予測

10期先までの予測を行う。

```R
#---
#　予測
#---

dlmForecasted_obj = dlmForecast(mod = dlmFiltered_obj, nAhead = 10)

# 予測分布の平均、標準偏差の導出
a                 = ts(data = dlmForecasted_obj$a, start = c(1971,1))
a_sdev            = sqrt( as.numeric( dlmForecasted_obj$R ) )
a_quant = list(a + qnorm(0.025, sd = a_sdev), a + qnorm(0.975, sd = a_sdev))

# 結果のプロット
oldpar = par(no.readonly = T) 
par(family = "HiraginoSans-W4")
ts.plot(cbind(Nile, a, do.call('cbind', a_quant)),
        col = c('lightgray', 'black', 'black', 'black'),
        lty = c('solid', 'solid', 'dashed', 'dashed'))
legend(legend = c('観測値', '平均(予測分布)', '95%区間'),
       lty    = c('solid', 'solid', 'dashed', 'dashed'),
       col    = c('lightgray', 'black', 'black', 'black'),
       x      = 'topright', text.width = 26, cex = 0.6)
par(oldpar)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/686685/5a51a085-d824-a4c3-9e50-84d0a2231f56.png)


dlmForecast()関数で予測を行う。この関数は

- a：予測分布の平均ベクトル
- R：予測分布の共分散行列

などを返す。

[次回：基礎から分かる時系列分析 4.カルマンフィルタ(ローカルトレンドモデルの実装)]()