演習の概要

注意事項

  • 慶應義塾大学SFCで開講している「統計解析」の授業履修者向けの演習用ページです。
  • 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.6.0を使っています。
  • Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

分析の準備

パッケージのインストール

install.packages("datasets")
install.packages("ggplot2")
install.packages("dlm")
install.packages("rstan")
install.packages("forecast")
install.packages("tibble")
library(datasets)
library(ggplot2)
library(dlm)
library(rstan)
library(forecast)
library(tibble)

オプションの設定

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
needs::prioritize(magrittr) 

データの準備

  • datasetsパッケージに含まれるNileデータを用います。このデータは、1871-1970年の100年間のナイル川での年間河川流量(\(10^8 m^3\))に関する時系列データです。1898年のアスワンダム建設で、1989年と1899年の間に流量が変化しています。
data(Nile)

モデル推定

最尤推定

  • dlmパッケージを用いて、ローカルレベルモデル、フィルタリング、スムージングなどを行えます。
  • まず、分析用dataframeを作成します。
Nile.df <- data.frame(
  "flow" = c(Nile),
  "year" = as.numeric(c(1871:1970), units="years")
)
  • 次に、作成したNile.dfからナイル川の流量を可視化します。
Nile.df %>% 
  ggplot(., aes(x=year)) +
  geom_line(aes(y=flow), col="blue") +
  labs(y="Annual Flow", x ="Year") 

ローカルレベルモデル

  • まず、ローカルレベルモデルを推定します。そのために、以下のような関数を作成します。
buildModPoly1 <- function(theta){
  dlmModPoly(order=1,           # ローカルレベルモデル。2はローカルトレンドモデル
             dV=exp(theta[1]),  # 観測方程式のノイズ
             dW=exp(theta[2]),  # 状態の分散行列の対角成分
             m0=theta[3])       # 初期値mu0
}
  • dlmMLE()関数を使ってモデルを推定します。
Nile.local <- dlmMLE(Nile.df$flow, parm=c(1,1,Nile.df$flow[1]), buildModPoly1)
Nile.local
$par
[1]    9.622363    7.292343 1119.999910

$value
[1] 549.63

$counts
function gradient 
      29       29 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
  • 推定したパラメータを使ってモデルを組み直します。
mod.Nile <- buildModPoly1(Nile.local$par)
mod.Nile
$FF
     [,1]
[1,]    1

$V
         [,1]
[1,] 15098.68

$GG
     [,1]
[1,]    1

$W
         [,1]
[1,] 1469.009

$m0
[1] 1120

$C0
      [,1]
[1,] 1e+07

フィルタリング

  • 次に、カルマンフィルターによるフィルタリングを行います。
Nile.filt <- dlmFilter(Nile.df$flow, mod.Nile)
Nile.filt$m
  [1] 1119.9999 1120.0000 1140.9141 1072.8136 1117.3106 1129.9719 1138.4573 1048.8600 1098.0306 1171.3000 1162.9020
 [12] 1117.9508 1069.0283 1079.9775 1057.0087 1047.1237 1023.8553 1065.5552  994.3711  984.6583 1026.1416 1045.8651
 [23] 1089.6964 1105.8001 1144.3077 1175.2026 1187.1655 1145.1955 1133.1263 1037.2240  984.5566  955.0332  885.3260
 [34]  899.9263  882.0541  833.7048  855.6812  811.9712  867.5239  916.2530  930.3387  903.8110  856.3277  749.4229
 [45]  769.3382  751.3560  849.7999  916.6142  894.0185  859.2979  849.0707  827.4213  832.1156  840.6301  846.3368
 [56]  806.7244  816.9457  797.4660  797.0745  861.9462  834.4551  820.1803  832.1491  835.5808  864.5334  896.4362
 [67]  896.5868  876.6689  912.2741  874.5478  821.5269  775.4547  794.2934  799.0218  783.7945  788.3891  855.5801
 [78]  856.7604  861.3641  857.7953  866.3954  833.7105  811.0891  818.2755  880.1560  890.2620  915.8282  884.0959
 [89]  894.4850  915.9860  889.0183  923.9961  919.1904  914.3328  982.6070  963.7519  905.6027  909.1803  858.1269
[100]  819.6388  798.3718
  • フィルタリングの結果を可視化します(オレンジ色)。
Nile.df %>% 
  ggplot(., aes(x=year)) +
  geom_line(aes(y=flow), col = "blue") +
  geom_line(aes(y=Nile.filt$m[2:101]), col="orange") +
  labs(y="Annual Flow", x="Year")