R：用寓言、小语和地图预测多个时间序列

Question

我正在尝试使用 R 包 tsibble 和 fable 来拟合一些时间序列，它们仍在建设中，以替代可信赖的 Rob Hyndman 的 forecast 包。该系列全部合并为一个 tsibble，然后我将其与 ARIMA 配合使用，该功能可替代 forecast::auto.arima。

我使用map_at，首先迭代除Date 之外的所有元素，然后再次使用fablelite::components 从适合每个系列的模型中提取模型信息。 （很多fable 函数实际上都在fablelite 中）。

这失败了，显然是因为组件需要 mdl_df 类的对象，而我的模型对象有 mdl_defn 类

这是一个（几乎）重现错误的玩具示例：

library(tidyverse)
library(tsibble)
library(fable)
set.seed(1)
ar1  <-  arima.sim(model=list(ar=.6), n=10)
ma1 <- arima.sim(model=list(ma=0.4), n=10)
Date  <- c(ymd("2019-01-01"):ymd("2019-01-10"),  ymd("2019-01-01"):ymd("2019-01-10"))
tb <- tibble(Date, ar1, ma1)

# Fit the whole series
tb_all <- tb   %>% 
map_at(.at =  c("ar1", "ma1"), .f = ARIMA)
names(arima_all[2:3])<- c("ar1", "ma1")

# Extract model components
tb_components <- tb %>%  
  map_at(.at = c("ar1", "ma1"), 
         .f = fablelite::components)

请注意，在这个玩具中，和我的真实数据一样，时间是每周 5 天，缺少周末

在这个玩具示例中，错误消息显示组件函数拒绝列表元素，因为类 ts 没有方法。在我的真实案例中，它使用更长的系列和更多的系列，但在我看来其他方面是相同的，元素被拒绝，因为它们属于 mdl_defn 类。请注意，如果我用str( ) 检查tb_all 的第二个和第三个元素，它们也会显示为类'mdl_defn'、'R6' 不确定错误消息中的ts 来自哪里。

Answer 1

这里有一个例子，希望能做你想做的事情。

首先，你需要创建一个 tsibble：

library(tidyverse)
library(tsibble)
library(fable)
library(lubridate)
set.seed(1)
ar1  <-  arima.sim(model=list(ar=.6), n=30)
ma1 <- arima.sim(model=list(ma=0.4), n=30)
Date  <- ymd(paste0("2019-01-",1:30))
tb <- bind_cols(Date=Date, ar1=ar1, ma1=ma1) %>%
  gather("Series", "value", -Date) %>%
  as_tsibble(index=Date, key=Series)
tb
#> # A tsibble: 60 x 3 [1D]
#> # Key:       Series [2]
#>    Date       Series   value
#>    <date>     <chr>    <dbl>
#>  1 2019-01-01 ar1    -2.07  
#>  2 2019-01-02 ar1    -0.118 
#>  3 2019-01-03 ar1    -0.116 
#>  4 2019-01-04 ar1    -0.0856
#>  5 2019-01-05 ar1     0.892 
#>  6 2019-01-06 ar1     1.36  
#>  7 2019-01-07 ar1     1.41  
#>  8 2019-01-08 ar1     1.76  
#>  9 2019-01-09 ar1     1.84  
#> 10 2019-01-10 ar1     1.18  
#> # … with 50 more rows

这包含两个系列：ar1 和 ma1 在相同的 30 天内。

接下来，您可以通过一个简单的函数将 ARIMA 模型拟合到两个系列。

tb_all <- tb %>% model(arima = ARIMA(value))
tb_all
#> # A mable: 2 x 2
#> # Key:     Series [2]
#>   Series arima                 
#>   <chr>  <model>               
#> 1 ar1    <ARIMA(0,0,2)>        
#> 2 ma1    <ARIMA(0,0,0) w/ mean>

最后，尚不清楚您要使用components() 提取什么，但也许以下之一可以满足您的要求：

tidy(tb_all)
#> # A tibble: 3 x 7
#>   Series .model term     estimate std.error statistic  p.value
#>   <chr>  <chr>  <chr>       <dbl>     <dbl>     <dbl>    <dbl>
#> 1 ar1    arima  ma1         0.810     0.198      4.09 0.000332
#> 2 ar1    arima  ma2         0.340     0.181      1.88 0.0705  
#> 3 ma1    arima  constant    0.295     0.183      1.61 0.118
glance(tb_all)
#> # A tibble: 2 x 9
#>   Series .model sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
#>   <chr>  <chr>   <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
#> 1 ar1    arima   0.695   -36.4  78.9  79.8  83.1 <cpl [0]> <cpl [2]>
#> 2 ma1    arima   1.04    -42.7  89.4  89.8  92.2 <cpl [0]> <cpl [0]>
augment(tb_all)
#> # A tsibble: 60 x 6 [1D]
#> # Key:       Series, .model [2]
#>    Series .model Date         value .fitted  .resid
#>    <chr>  <chr>  <date>       <dbl>   <dbl>   <dbl>
#>  1 ar1    arima  2019-01-01 -2.07    -0.515 -1.56  
#>  2 ar1    arima  2019-01-02 -0.118   -1.21   1.09  
#>  3 ar1    arima  2019-01-03 -0.116    0.511 -0.627 
#>  4 ar1    arima  2019-01-04 -0.0856  -0.155  0.0690
#>  5 ar1    arima  2019-01-05  0.892   -0.154  1.05  
#>  6 ar1    arima  2019-01-06  1.36     0.871  0.486 
#>  7 ar1    arima  2019-01-07  1.41     0.749  0.659 
#>  8 ar1    arima  2019-01-08  1.76     0.699  1.06  
#>  9 ar1    arima  2019-01-09  1.84     1.09   0.754 
#> 10 ar1    arima  2019-01-10  1.18     0.973  0.206 
#> # … with 50 more rows

要以传统方式查看模型输出，请使用report()：

tb_all %>% filter(Series=='ar1') %>% report()
#> Series: value 
#> Model: ARIMA(0,0,2) 
#> 
#> Coefficients:
#>          ma1     ma2
#>       0.8102  0.3402
#> s.e.  0.1982  0.1809
#> 
#> sigma^2 estimated as 0.6952:  log likelihood=-36.43
#> AIC=78.86   AICc=79.78   BIC=83.06
tb_all %>% filter(Series=='ma1') %>% report()
#> Series: value 
#> Model: ARIMA(0,0,0) w/ mean 
#> 
#> Coefficients:
#>       constant
#>         0.2950
#> s.e.    0.1833
#> 
#> sigma^2 estimated as 1.042:  log likelihood=-42.68
#> AIC=89.36   AICc=89.81   BIC=92.17