1.はじめに
時系列データの分析に状態空間モデルが使われています。
ガウス状態空間モデルの説明は、ここで説明するよりも有名な次のサイトをご覧ください。
logic of blues なぜ状態空間モデルを使うのか
Rでは状態空間モデルは、KFAS,dlm,ベイズ推論（Rstan)等を用いることが多いと思いますが、ここではKFASパッケージを使ってみたいと思います。
異常値検出というのは、時系列問題でよくテーマになるかと思いますが、ここではナイル川の流量データを用いたいと思います。
課題は、アスワンダムが建設された年を知ることです。
前提知識として、1900年頃に建設されたらしいということは知っていますが、正確な年代はわかりません。
そこで、流出量のデータから変化点（ダムの完成による流出量の減少）を読み取ることを目的とします。
2.パッケージの有効化
library(tidyverse)
library(KFAS)
3.データの概要
ナイル川の流量データです。
d = Nile
ts.plot(d)
1900年前後から流量が減った傾向があるように見えます。
4.ローカルレベル :カルマンフィルタの適用
ローカルレベルについても、ここで説明するよりも有名な次のサイトをご覧ください。
logic of blues ローカルレベルモデル
まず、最初にローカルレベルで平滑化をしてみます。
4.1Step1:modelの構造を決める
build_kfas = SSModel(H=NA,d ~ SSMtrend(degree=1,Q=NA))
H:観測誤差の分散 未知数なのでNAとしています。
Q:システム誤差の分散 未知数なのでNAとしています。
SSM:トレンドモデルの指定 degreeはローカルレベルである1を指定します。
4.2Step2:パラメータの推定
fit_kfas = fitSSM(build_kfas, inits = c(1,1))
パラメータを最尤法で推定します。
初期値はシステム誤差、観測誤差ともに１にしています。
（１に意味はありません）
パラメータの推定結果です。
・システム誤差
fit_kfas$model$Q
$> , , 1
$>
$> [,1]
$>[1,] 1396.146
・観測誤差
fit_kfas$model$H
$> , , 1
$>
$> [,1]
$>[1,] 15298.12
4.3Step3:スムージング(平滑化)
#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:図示化
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)
状態の推移から、1900年前後に流量の変化があるようですが、この図からは年代の特定は難しいです。
5 外生変数の追加
ダムの建設をyear_damとする外生変数をモデルに追加します。
year_damは{0,1}の値をとり、1がダムが存在することを意味します。
建設年が不明のため、1900年から前後10年間をスタートとしてyear_damに当てはめてみます。
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 外生変数の比較による年代の推定
各年毎の外生変数の係数が最も低い値を比較する。
#最初に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