tidy時系列データに対してのVARモデル fable
これからのRにおいてVARモデルで予測したい場合には、varsよりもfableを使えば良さそうというもの。
データ
VARモデルについて詳しく解説してくださっているものと同じCanadaデータを用いることにします。
tjo.hatenablog.com
> library(tidyverse) > library(tsibble) > library(lubridate) > library(fable) > data(Canada, package = "vars")
fableではtsibbleクラスを用いるため、tsクラスから変換します。
> Canada_tsble <- Canada %>% + as_tsibble(pivot_longer = FALSE) > Canada_tsble # A tsibble: 84 x 5 [1Q] index e prod rw U <qtr> <dbl> <dbl> <dbl> <dbl> 1 1980 Q1 930. 405. 386. 7.53 2 1980 Q2 930. 405. 388. 7.70 3 1980 Q3 930. 404. 391. 7.47 4 1980 Q4 931. 404. 394. 7.27 5 1981 Q1 933. 405. 397. 7.37 6 1981 Q2 934. 404. 400. 7.13 7 1981 Q3 934. 403. 401. 7.4 8 1981 Q4 933. 402. 406. 8.33 9 1982 Q1 932. 402. 409. 8.83 10 1982 Q2 931. 401. 411. 10.4 # … with 74 more rows
モデルの推定
2年分のデータを残し、VARモデルを当てはめます。
> Canada_tsble_train <- Canada_tsble %>% + filter(year(index) < 1999) > var_fit <- Canada_tsble_train %>% + model( + VAR = fable::VAR(vars(e, prod, rw, U) ~ AR(p = 3:5)), + VARc = fable::VAR(vars(e, prod, rw, U) ~ AR(p = 3:5) + 1) + )
AR(p)のpを複数指定することで、情報量基準で最適なモデルを選択するようになります。
それぞれのモデルの対数尤度と各情報量基準を確認します。
> var_fit %>% + glance() # A tibble: 2 x 6 .model sigma2 log_lik AIC AICc BIC <chr> <list> <dbl> <dbl> <dbl> <dbl> 1 VAR <dbl [4 × 4… -131. 423. -1017. 605. 2 VARc <dbl [4 × 4… -121. 410. -689. 601.
定数項を含めたモデルの方が良さそうなので、こちらの詳細を確認します。
> var_fit %>% + select(VARc) %>% + report() Series: e, prod, rw, U Model: VAR(4) w/ mean Coefficients for e: lag(e,1) lag(prod,1) lag(rw,1) lag(U,1) lag(e,2) lag(prod,2) lag(rw,2) lag(U,2) 1.6044 0.2200 -0.1167 0.3128 -0.7695 -0.0354 -0.0692 0.1302 s.e. 0.1821 0.0643 0.0552 0.2203 0.2868 0.1008 0.0722 0.2593 lag(e,3) lag(prod,3) lag(rw,3) lag(U,3) lag(e,4) lag(prod,4) lag(rw,4) lag(U,4) 0.3964 -0.0184 -0.0425 0.4448 0.3274 -0.0306 0.0125 0.2742 s.e. 0.2761 0.0959 0.0713 0.2506 0.2067 0.0660 0.0559 0.2285 constant -497.7130 s.e. 117.3299 Coefficients for prod: lag(e,1) lag(prod,1) lag(rw,1) lag(U,1) lag(e,2) lag(prod,2) lag(rw,2) lag(U,2) -0.1439 1.0897 0.0858 -0.7909 -0.3212 -0.1565 -0.2165 0.5142 s.e. 0.3875 0.1369 0.1176 0.4688 0.6103 0.2145 0.1535 0.5518 lag(e,3) lag(prod,3) lag(rw,3) lag(U,3) lag(e,4) lag(prod,4) lag(rw,4) lag(U,4) 0.4989 0.1124 0.0918 0.5520 -0.1199 -0.0916 0.0709 -0.3805 s.e. 0.5876 0.2041 0.1517 0.5332 0.4398 0.1405 0.1189 0.4863 constant 87.2793 s.e. 249.6723 Coefficients for rw: lag(e,1) lag(prod,1) lag(rw,1) lag(U,1) lag(e,2) lag(prod,2) lag(rw,2) lag(U,2) -0.4102 -0.0663 0.8890 0.2429 0.6148 -0.2489 -0.1361 -0.5850 s.e. 0.4459 0.1575 0.1353 0.5394 0.7022 0.2468 0.1767 0.6349 lag(e,3) lag(prod,3) lag(rw,3) lag(U,3) lag(e,4) lag(prod,4) lag(rw,4) lag(U,4) -0.3956 0.3220 0.0304 -0.2078 0.1726 -0.0992 0.1980 0.1803 s.e. 0.6761 0.2348 0.1745 0.6134 0.5060 0.1617 0.1368 0.5595 constant 68.4229 s.e. 287.2648 Coefficients for U: lag(e,1) lag(prod,1) lag(rw,1) lag(U,1) lag(e,2) lag(prod,2) lag(rw,2) lag(U,2) -0.6387 -0.1610 0.0440 0.3868 0.3479 0.0787 0.0873 -0.1734 s.e. 0.1551 0.0548 0.0471 0.1877 0.2443 0.0859 0.0615 0.2209 lag(e,3) lag(prod,3) lag(rw,3) lag(U,3) lag(e,4) lag(prod,4) lag(rw,4) lag(U,4) -0.0729 0.0160 -0.0093 -0.0105 0.0001 0.0174 0.0058 -0.0048 s.e. 0.2352 0.0817 0.0607 0.2134 0.1760 0.0563 0.0476 0.1946 constant 314.5203 s.e. 99.9397 Residual covariance matrix: e prod rw U e 0.1019 -0.0152 -0.0578 -0.0610 prod -0.0152 0.4615 0.0343 0.0022 rw -0.0578 0.0343 0.6110 0.0528 U -0.0610 0.0022 0.0528 0.0739 log likelihood = -120.94 AIC = 409.87 AICc = -688.59 BIC = 601.11
予測と評価
残しておいたデータが含まれる2年先までの予測を行います。
> var_forecast <- var_fit %>% + forecast(h = "2 years") > > var_forecast %>% + autoplot(Canada_tsble) + + theme(legend.position = "bottom")
最後に予測精度の評価指標を確認します。
> var_forecast %>% + accuracy(Canada_tsble) # A tibble: 8 x 10 .model .response .type ME RMSE MAE MPE MAPE MASE ACF1 <chr> <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 VAR e Test 0.965 1.37 1.04 0.100 0.108 0.483 0.643 2 VAR prod Test 3.65 4.12 3.65 0.876 0.876 2.26 0.722 3 VAR rw Test -1.09 1.32 1.09 -0.233 0.233 0.255 0.361 4 VAR U Test -0.299 0.664 0.612 -4.66 8.66 0.618 0.690 5 VARc e Test -1.11 1.21 1.11 -0.116 0.116 0.519 0.316 6 VARc prod Test 3.96 4.53 3.96 0.951 0.951 2.46 0.715 7 VARc rw Test -0.970 1.50 0.977 -0.207 0.208 0.228 0.432 8 VARc U Test 1.08 1.10 1.08 15.1 15.1 1.09 0.0429
まとめ
普通のVARモデル（SVARなどは未実装）で予測したいならば、これからはfableを使えば良い。tidyverseな世界でデータを扱えるので、一癖ある時系列クラスのデータと格闘しなくてよくなる。複数のモデルも一度に推定から予測、評価まで行うことができる。
