使用するデータ

以前にARIMAモデルの記事で使用した横浜市の月別火災発生件数を題材とする。
【統計】時系列分析（ARIMA、SARIMA、SARIMAXをRで実行する） - こちにぃるの日記

ソース：
横浜市のオープンデータ
 気象庁ホームページ

2000年1月～2019年12月の20年間240カ月分を用いる。
ARIMAモデルの記事と同様に、17年分（2000～2016年）を学習期間、３年（2017～2019年）をテスト期間する（ホールドアウト法）。
また、火災と関連がありそうなものとして、横浜市の月別平均気温を回帰成分（外生変数）として用いる。

#-- import
FIRE <- read_csv("./yokohama_kasai.csv")
TEMP <- read.csv("./yokohama_kion.csv")

#-- データ加工（時系列、回帰成分）
colnames(FIRE) <- c("YEARJ","YEAR","MONTH","FIRE")
ADTS <- data.frame(FIRE, TEMP=TEMP$平均気温)
ADTS$TESTFL <- ifelse(ADTS$YEAR >= 2017, 1, 0)  #-- 2017年以降をテストデータ
ADTS$CNT    <- ifelse(ADTS$TESTFL == 0, ADTS$FIRE, NA)  #-- テストデータはNA

状態空間モデルの実装

今回は、ローカル線形トレンドモデルを使って以下の３モデルを比較してみる。
・ローカル線形トレンドモデル
・ローカル線形トレンドモデル＋季節成分
・ローカル線形トレンドモデル＋季節成分＋回帰成分（平均気温）

#-- ローカル線形トレンド
mod1 <- SSModel(ADTS$CNT ~ SSMtrend(2, Q=list(NA,NA)), H=NA)
fit1 <- fitSSM(mod1, inits=numeric(3), method="BFGS")
kfs1 <- KFS(fit1$model, smoothing=c("state","mean","disturbance"))

#-- ローカル線形トレンド + 季節成分
mod2 <- SSModel(ADTS$CNT ~ SSMtrend(2, Q=list(NA,NA)) + SSMseasonal(12, Q=NA), H=NA)
fit2 <- fitSSM(mod2, inits=numeric(4), method="BFGS")
kfs2 <- KFS(fit2$model, smoothing=c("state","mean","disturbance"))

#-- ローカル線形トレンド + 季節成分 + 回帰成分
mod3 <- SSModel(ADTS$CNT ~ SSMtrend(2, Q=list(NA,NA)) + SSMseasonal(12, Q=NA) + SSMregression(~ ADTS$TEMP, Q=NA), H=NA)
fit3 <- fitSSM(mod3, inits=numeric(5), method="BFGS")
kfs3 <- KFS(fit3$model, smoothing=c("state","mean","disturbance"))

実行結果は以下の通り（学習データ：2000～2016年、予測：2017～2019年）。

モデル評価

残差分析は今回割愛する。本来は残差を見て外生変数の取捨選択や経時的な構造変化の有無などを検討する。
例えば、大きな災害や技術革新によって構造変化があれば、変化時点＝１、それ以外＝０としたダミー変数を設けるなどしてモデルの改善を目指す。

以下はテストデータに対するAICとRMSE、MAE、MAPEによる評価を記載。
長期予測の評価においては、AICは不向きなようである（AICは１期先）。
RMSE、MAE、MAPEなどの評価指標については、データの特性に応じて用いる指標を決める。
例えばスパイクの多いデータではRMSEは不向きかもしれない（スパイクに過剰に反応してしまうかもしれない）。0件が発生するデータではMAPEは不向きかもしれない（分母ゼロで計算できない）。

> #-- AIC
> logLik.1 <- kfs1$logLik - sum(kfs1$Finf>0) * log(2*pi)/2
> logLik.2 <- kfs2$logLik - sum(kfs2$Finf>0) * log(2*pi)/2
> logLik.3 <- kfs3$logLik - sum(kfs3$Finf>0) * log(2*pi)/2
> AIC1 <- -2*logLik.1 + 2*(3 + 1)           #-- 攪乱項(NA数)+トレンド成分
> AIC2 <- -2*logLik.2 + 2*(4 + 1 + 11)      #-- 攪乱項(NA数)+トレンド成分+季節成分
> AIC3 <- -2*logLik.3 + 2*(5 + 1 + 11 + 1)  #-- 攪乱項(NA数)+トレンド成分+季節成分+回帰成分 ←回帰成分の数え方はこれで良い？
> print(paste0(round(c(AIC1,AIC2,AIC3),4)))
[1] "1787.9992" "1640.0171" "1643.6485"

> #-- RMSE, MAE, MAPE
> ADTS2 <- ADTS %>% mutate(KFAS1=c(kfs1$muhat), KFAS2=c(kfs2$muhat), KFAS3=c(kfs3$muhat))
> ADTS3 <- ADTS2[ADTS2$TESTFL==1,] #-- テストデータのみにする
> 
> fnc1 <- function(KFAS){
+   RMSE <- sqrt(1/length(ADTS3$FIRE) * sum(KFAS - ADTS3$FIRE)^2)
+   MAE  <- 1/length(ADTS3$FIRE) * sum(abs(KFAS - ADTS3$FIRE))
+   MAPE <- 1/length(ADTS3$FIRE) * sum(abs(KFAS - ADTS3$FIRE) / ADTS3$FIRE) * 100
+   print(paste0(round(c(RMSE,MAE,MAPE),4)))
+ }
> fnc1(ADTS3$KFAS1)
[1] "133.343" "23.4793" "43.7605"
> fnc1(ADTS3$KFAS2)
[1] "2.1789"  "10.7435" "16.9743"
> fnc1(ADTS3$KFAS3)
[1] "0.0949"  "10.7321" "16.8962"