KFASパッケージを使ったガウス状態空間モデルによる変化点検出の試み

##1.はじめに
時系列データの分析に状態空間モデルが使われています。
ガウス状態空間モデルの説明は、ここで説明するよりも有名な次のサイトをご覧ください。
[logic of blues  なぜ状態空間モデルを使うのか](https://logics-of-blue.com/%e3%81%aa%e3%81%9c%e7%8a%b6%e6%85%8b%e7%a9%ba%e9%96%93%e3%83%a2%e3%83%87%e3%83%ab%e3%82%92%e4%bd%bf%e3%81%86%e3%81%ae%e3%81%8b/)

Rでは状態空間モデルは、KFAS,dlm,ベイズ推論（Rstan)等を用いることが多いと思いますが、ここではKFASパッケージを使ってみたいと思います。

異常値検出というのは、時系列問題でよくテーマになるかと思いますが、ここではナイル川の流量データを用いたいと思います。

課題は、アスワンダムが建設された年を知ることです。

前提知識として、1900年頃に建設されたらしいということは知っていますが、正確な年代はわかりません。

そこで、流出量のデータから変化点（ダムの完成による流出量の減少）を読み取ることを目的とします。

###2.パッケージの有効化

```{r}
library(tidyverse)
library(KFAS)
```
##3.データの概要
ナイル川の流量データです。

```{r}
d = Nile
ts.plot(d)
```
![01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/ca443d3c-abb9-4dc1-c27b-c4460dba7486.png)

1900年前後から流量が減った傾向があるように見えます。

###4.ローカルレベル :カルマンフィルタの適用

ローカルレベルについても、ここで説明するよりも有名な次のサイトをご覧ください。
[logic of blues ローカルレベルモデル](https://logics-of-blue.com/%E3%83%AD%E3%83%BC%E3%82%AB%E3%83%AB%E3%83%AC%E3%83%99%E3%83%AB%E3%83%A2%E3%83%87%E3%83%AB/)
まず、最初にローカルレベルで平滑化をしてみます。

##4.1Step1:modelの構造を決める

```{r}
build_kfas = SSModel(H=NA,d ~ SSMtrend(degree=1,Q=NA))
```

H:観測誤差の分散　未知数なのでNAとしています。
Q:システム誤差の分散　　未知数なのでNAとしています。
SSM:トレンドモデルの指定　degreeはローカルレベルである1を指定します。

##4.2Step2:パラメータの推定
```{r}
fit_kfas = fitSSM(build_kfas, inits = c(1,1))
```
パラメータを最尤法で推定します。
初期値はシステム誤差、観測誤差ともに１にしています。
（１に意味はありません）

パラメータの推定結果です。
・システム誤差

```{r}
fit_kfas$model$Q

$> , , 1
$>
$>         [,1]
$>[1,] 1396.146
```
・観測誤差

```{r}
fit_kfas$model$H

$> , , 1
$> 
$>         [,1]
$>[1,] 15298.12
```
##4.3Step3:スムージング(平滑化)

```{r}
#KFS関数はカルマンフィルタリング、それに続くスムージングを同時に実行する
result_kfas = KFS(fit_kfas$model,
                  filtering = c("state","mean"),
                  smoothing = c("state","mean"))

#predict関数でモデルからスムージングの結果を取り出す
smooth_conf = predict(fit_kfas$model,
                      interval = "confidence",level = 0.95)

#data.frame型に変更、観測値、年代を追加
smooth_conf = data.frame(smooth_conf,
                         y=as.numeric(d))

smooth_conf$year = 1871:1970
```
##4.4Step4:図示化

```{r}
ggplot(smooth_conf,aes(x=year,y=y)) +
  geom_point(alpha=0.5) +
  geom_line(aes(y=fit),size=1.2) +
  geom_ribbon(aes(ymin=lwr,ymax=upr),alpha=0.3)
```
![02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/531b8d14-881b-fac3-1548-8701605c67fe.png)

状態の推移から、1900年前後に流量の変化があるようですが、この図からは年代の特定は難しいです。
##5　外生変数の追加

ダムの建設をyear_damとする外生変数をモデルに追加します。
year_damは{0,1}の値をとり、1がダムが存在することを意味します。
建設年が不明のため、1900年から前後10年間をスタートとしてyear_damに当てはめてみます。

```{r}
reg=0
t_max　=　length(d)
for(i in 1:20){
year_dam = rep(0,t_max)
year_dam[which(1890+i <=time(d))] =1
#SSMregressionでyear_damを説明変数とする回帰分析を追加（外生変数）
build_reg = SSModel(H=NA,d ~ SSMtrend(degree=1,Q=NA) +
                      SSMregression(~year_dam,Q=NA))

fit_reg = fitSSM(build_reg,inits = c(1,1,1))
result_reg = KFS(
  fit_reg$model,
  filtering = c("state","mean"),
  smoothing =  c("state","mean")
)
interval_conf = predict(fit_reg$model,states = "regression",
                       interval = "confidence",level=0.95)
interval_conf=data.frame(interval_conf)

#説明変数の係数が最も低い値を記録する
reg=c(reg,min(interval_conf$fit))
}
```
##6　外生変数の比較による年代の推定
各年毎の外生変数の係数が最も低い値を比較する。

```{r}
#最初に0を入れてるので取り除く
reg = reg[-1]
df_reg=data.frame(reg=reg,time=1891:(1891+19))
df_reg

$> reg        time
$><dbl>      <int>
$>0.00000     1891		
$>0.00000     1892		
$>0.00000     1893		
$>0.00000     1894		
$>-140.84011  1895		
$>-234.85192  1896		
$>-244.64355  1897		
$>-244.48592  1898		
$>-248.03319  1899		
$>-235.77715  1900
$>-227.25009  1901		
$>-221.15336  1902		
$>-140.21783  1903		
$>-145.71143  1904		
$>-52.75103   1905		
$>-75.26610   1906		
$>0.00000     1907		
$>0.00000     1908		
$>0.00000     1909		
$>-93.64160   1910	
```
・考察
年間流量は1899年に-248で最も低くなったことから、ダムの完成した年代は1899年と推定するのが妥当である。


##参考
馬場真哉,「時系列分析と状態空間モデルの基礎：RとStanで学ぶ理論と実装」,プレアデス出版,2018.
[logic of blues](https://logics-of-blue.com/)