如何从 R 中的 forecast::auto.arima 获取前 x 个模型的所有详细信息？
How can I get all details of the top x models from forecast::auto.arima in R?
我的目标是从 R 中的 auto.arima 函数准确地重新创建前 3 个模型。我的示例使用以下系列：
> data <- c(79, 73, 102, 158, 235, 326, 216, 153, 82, 87, 94, 163, 119, 81, 143,
179, 182, 247, 248, 105, 143, 64, 45, 163, 122, 95, 47, 117, 345, 454,
596, 440, 207, 94, 100, 187, 106, 98, 66, 81, 178, 196, 279, 128, 80,
178, 112, 175, 91, 52, 35, 14, 20, 332, 386, 121, 167, 53, 94, 224,
76, 68, 42, 61, 125, 273, 431, 1901, 86, 65, 120, 227, 131, 84, 48, 153, 157)
> time_series <- ts(data, freq=12, start=c(2013, 10))
我使用下面显示的 auto.arima 函数来捕获最佳模型的模型信息，并在文本中查看顶级模型。
> library(forecast)
> auto_mod <- auto.arima(time_series, trace=TRUE)
ARIMA(2,0,2)(1,0,1)[12] with non-zero mean : 1067.045
ARIMA(0,0,0) with non-zero mean : 1057.566
ARIMA(1,0,0)(1,0,0)[12] with non-zero mean : 1057.504
ARIMA(0,0,1)(0,0,1)[12] with non-zero mean : 1057.727
ARIMA(0,0,0) with zero mean : 1092.138
ARIMA(1,0,0) with non-zero mean : 1055.976
ARIMA(1,0,0)(0,0,1)[12] with non-zero mean : 1057.629
ARIMA(1,0,0)(1,0,1)[12] with non-zero mean : 1058.261
ARIMA(2,0,0) with non-zero mean : 1058.174
ARIMA(1,0,1) with non-zero mean : 1058.187
ARIMA(0,0,1) with non-zero mean : 1056.126
ARIMA(2,0,1) with non-zero mean : Inf
ARIMA(1,0,0) with zero mean : 1070.847
Best model: ARIMA(1,0,0) with non-zero mean
> summary(auto_mod)
Series: time_series
ARIMA(1,0,0) with non-zero mean
Coefficients:
ar1 mean
0.2170 176.2688
s.e. 0.1105 32.0040
sigma^2 estimated as 49991: log likelihood=-524.82
AIC=1055.65 AICc=1055.98 BIC=1062.68
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.3042569 220.664 104.5716 -76.32587 96.76548 1.051865
ACF1
Training set 0.006605514
我以为我可以使用文本输出来准确地重新创建除前 1 个模型之外的模型，但我意识到文本输出中对模型的描述很接近，但似乎并不完全确定每个模型都是独一无二的。我认为它提供了有关是否存在截距或是否存在漂移的信息，但不能同时提供两者。我在下面找到了一个这样的例子，其中两个不同的模型，一个有截距，一个没有，具有相同的标签。
> mod1 <- Arima(time_series, order=c(1,0,0), seasonal=c(1,0,0), include.drift=TRUE,
+ include.mean=TRUE)
> summary(mod1)
Series: time_series
ARIMA(1,0,0)(1,0,0)[12] with drift
Coefficients:
ar1 sar1 intercept drift
0.1722 0.1976 132.4150 1.2188
s.e. 0.1175 0.2059 68.9234 1.5233
sigma^2 estimated as 50174: log likelihood=-524.15
AIC=1058.31 AICc=1059.15 BIC=1070.03
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.05272209 218.099 99.62455 -76.12583 94.97947 1.002104
ACF1
Training set 0.004756765
> mod2 <- Arima(time_series, order=c(1,0,0), seasonal=c(1,0,0), include.drift=TRUE,
+ include.mean=FALSE)
> summary(mod2)
Series: time_series
ARIMA(1,0,0)(1,0,0)[12] with drift
Coefficients:
ar1 sar1 drift
0.2014 0.2862 3.6940
s.e. 0.1205 0.1960 0.9009
sigma^2 estimated as 51238: log likelihood=-525.77
AIC=1059.53 AICc=1060.09 BIC=1068.91
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 19.14777 221.9052 107.1779 -57.65942 96.53528 1.078082
ACF1
Training set -0.01092874
所以总而言之，我的问题是：有没有办法在 R 中使用自动拟合过程来获得顶级 x arima 模型？我对任何解决方案持开放态度，但我对解决方案的偏好顺序是：
- 一种从
auto.arima 中提取前 x 个模型的模型对象的方法
- 一种在文本中提取有关前 x 个模型的信息或更好地理解跟踪参数以便自己重新创建模型的方法。
- 一些其他包或功能
在此先感谢您提供的任何帮助！
来自 forecast::auto.arima 轨迹的信息唯一确定模型。
让我们针对第三个模型明确验证这一点（使用您的示例数据）：
ARIMA(1,0,0)(1,0,0)[12] with non-zero mean : 1057.504
我们使用 base R 的 arima 手动拟合相同的模型
fit3 <- arima(
time_series,
order = c(1, 0, 0),
seasonal = list(order = c(1, 0, 0), period = 12L),
include.mean = TRUE)
使用ML估计模型参数时，arima returns是AIC值；另一方面，forecast::auto.arima 的轨迹 return 是校正后的（对于小样本）AIC。感谢this post on Cross Validated，我们可以编写一个自定义函数，根据模型参数计算 AICc。
get_AICc <- function(fit) {
npar <- length(fit$coef) + 1
nstar <- length(fit$residuals) - fit$arma[6] - fit$arma[7] * fit$arma[5]
fit$aic + 2 * npar * (nstar/(nstar - npar - 1) - 1)
}
我们确认修正后的 AIC 与 forecast::auto.arima 中报告的相同：
get_AICc(fit3)
#[1] 1057.504
forecast::Arima 自动 returns AIC 和 AICc 值，拟合结果与 arima 一致（不出所料，因为 forecast::Arima 内部调用 arima) 和 forecast::auto.arima:
Arima(
time_series,
order = c(1, 0, 0),
seasonal = c(1, 0, 0),
include.mean = TRUE)
)
#Series: time_series
#ARIMA(1,0,0)(1,0,0)[12] with non-zero mean
#
#Coefficients:
# ar1 sar1 mean
#0.1849 0.1780 179.5156
#s.e. 0.1168 0.2077 36.2321
#
#sigma^2 estimated as 49966: log likelihood=-524.47
#AIC=1056.95 AICc=1057.5 BIC=1066.32
与基本 R 的 arima 相比，forecast::Arima 允许在包含常量方面具有更大的灵活性：而 arima 允许通过 include.mean、forecast::Arima 还允许使用 time-dependent（确定性）漂移参数。有关详细信息，请参阅 Constants and ARIMA models in R。
我不知道 already-available 从 forecast::auto.arima 中获取前 n 个模型拟合列表的方法。跟踪仅将参数 space 中搜索的中间结果打印到控制台。过去，我定义并编写了自己的 SARIMA 模型参数参数网格搜索，以根据其中一个信息标准确定最佳拟合模型；然而，forecast::auto.arima 更复杂并且还考虑（即排除）单位根接近单位圆的模型拟合。您最好的选择是定义您自己的网格搜索并使其 return 成为前 n 个模型的列表；或者将forecast::auto.arima的代码修改为return这样的列表；或者（可能是最不优雅的解决方案）捕获并解析 auto.arima(... trace = TRUE) 的输出以提取模型参数和 re-fit 使用 arima 或 forecast::Arima.
更新 1
我确实说过捕获 auto.arima 的输出可能是最不优雅的方法，但这就是我最终所做的；也许这对您有帮助作为起点：
让我们定义一个函数来解析从 auto.arima
生成的 ARIMA 模型字符串
parse_SARIMA_string <- function(x) {
include.mean <- str_detect(x, "with non-zero mean")
x %>%
str_remove_all(":.+$") %>%
str_replace_all("\D", " ") %>%
trimws() %>%
enframe() %>%
separate(
value,
c("p", "d", "q", "P", "D", "Q", "period"),
sep = "\s+",
fill = "right") %>%
mutate(include.mean = include.mean)
}
然后我们可以使用capture.output来捕获forecast::auto.arima的踪迹并解析模型字符串
library(tidyverse)
capture.output(auto.arima(time_series, trace = TRUE)) %>%
str_subset("^ ARIMA") %>%
parse_SARIMA_string()
## A tibble: 13 x 9
# name p d q P D Q period include.mean
# <int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <lgl>
# 1 1 2 0 2 1 0 1 12 TRUE
# 2 2 0 0 0 NA NA NA NA TRUE
# 3 3 1 0 0 1 0 0 12 TRUE
# 4 4 0 0 1 0 0 1 12 TRUE
# 5 5 0 0 0 NA NA NA NA FALSE
# 6 6 1 0 0 NA NA NA NA TRUE
# 7 7 1 0 0 0 0 1 12 TRUE
# 8 8 1 0 0 1 0 1 12 TRUE
# 9 9 2 0 0 NA NA NA NA TRUE
#10 10 1 0 1 NA NA NA NA TRUE
#11 11 0 0 1 NA NA NA NA TRUE
#12 12 2 0 1 NA NA NA NA TRUE
#13 13 1 0 0 NA NA NA NA FALSE
现在您有了带有模型参数的 data.frame/tibble，可以直接在 arima 或 forecast::Arima.
中使用它们
更新 2
关于基于auto.arima模型字符串的SARIMA模型的可识别性，我们来看一些具体的模型。
有差分
模型 1 是具有常数的 ARMA 模型
模型 2 是一个有漂移的 ARIMA 模型
我们可以将模型展开得到
模型 3 是具有漂移和常数的 ARIMA 模型
当我们扩展模型 3 时，我们注意到常数项由于差异而取消，我们最终得到相同的模型 2。
因此 d=1 具有漂移和漂移+常数的模型描述了相同的过程。在这种情况下，将模型报告为例如“ARIMA(1,1,0) with drift”指的是与“ARIMA(1,1,0) with drift and non-zero mean”相同的模型。
无差异
对于 d=0，默认使用 include.drift = TRUE 会导致漂移和常数偏移参数，参见例如Constants and ARIMA models in R。如果你指定 include.drift = TRUE 然后手动强制 include.mean = FALSE 你会得到一个不同的模型（如你在 mod1 和 mod2 中所示）。据我了解，这永远不会发生在 forecast::auto.arima 中，其中 include.drift = TRUE 隐含地 总是 include.mean = TRUE。因此来自 forecast::auto.arima 的 SARIMA 字符串 ARIMA(1,0,0)(1,0,0)[12] with drift 将转换为具有漂移 和 [=136= 的 SARIMA(1,0,0)(1,0,0)[12] 模型]偏移量。
我的目标是从 R 中的 auto.arima 函数准确地重新创建前 3 个模型。我的示例使用以下系列：
> data <- c(79, 73, 102, 158, 235, 326, 216, 153, 82, 87, 94, 163, 119, 81, 143,
179, 182, 247, 248, 105, 143, 64, 45, 163, 122, 95, 47, 117, 345, 454,
596, 440, 207, 94, 100, 187, 106, 98, 66, 81, 178, 196, 279, 128, 80,
178, 112, 175, 91, 52, 35, 14, 20, 332, 386, 121, 167, 53, 94, 224,
76, 68, 42, 61, 125, 273, 431, 1901, 86, 65, 120, 227, 131, 84, 48, 153, 157)
> time_series <- ts(data, freq=12, start=c(2013, 10))
我使用下面显示的 auto.arima 函数来捕获最佳模型的模型信息，并在文本中查看顶级模型。
> library(forecast)
> auto_mod <- auto.arima(time_series, trace=TRUE)
ARIMA(2,0,2)(1,0,1)[12] with non-zero mean : 1067.045
ARIMA(0,0,0) with non-zero mean : 1057.566
ARIMA(1,0,0)(1,0,0)[12] with non-zero mean : 1057.504
ARIMA(0,0,1)(0,0,1)[12] with non-zero mean : 1057.727
ARIMA(0,0,0) with zero mean : 1092.138
ARIMA(1,0,0) with non-zero mean : 1055.976
ARIMA(1,0,0)(0,0,1)[12] with non-zero mean : 1057.629
ARIMA(1,0,0)(1,0,1)[12] with non-zero mean : 1058.261
ARIMA(2,0,0) with non-zero mean : 1058.174
ARIMA(1,0,1) with non-zero mean : 1058.187
ARIMA(0,0,1) with non-zero mean : 1056.126
ARIMA(2,0,1) with non-zero mean : Inf
ARIMA(1,0,0) with zero mean : 1070.847
Best model: ARIMA(1,0,0) with non-zero mean
> summary(auto_mod)
Series: time_series
ARIMA(1,0,0) with non-zero mean
Coefficients:
ar1 mean
0.2170 176.2688
s.e. 0.1105 32.0040
sigma^2 estimated as 49991: log likelihood=-524.82
AIC=1055.65 AICc=1055.98 BIC=1062.68
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.3042569 220.664 104.5716 -76.32587 96.76548 1.051865
ACF1
Training set 0.006605514
我以为我可以使用文本输出来准确地重新创建除前 1 个模型之外的模型，但我意识到文本输出中对模型的描述很接近，但似乎并不完全确定每个模型都是独一无二的。我认为它提供了有关是否存在截距或是否存在漂移的信息，但不能同时提供两者。我在下面找到了一个这样的例子，其中两个不同的模型，一个有截距，一个没有，具有相同的标签。
> mod1 <- Arima(time_series, order=c(1,0,0), seasonal=c(1,0,0), include.drift=TRUE,
+ include.mean=TRUE)
> summary(mod1)
Series: time_series
ARIMA(1,0,0)(1,0,0)[12] with drift
Coefficients:
ar1 sar1 intercept drift
0.1722 0.1976 132.4150 1.2188
s.e. 0.1175 0.2059 68.9234 1.5233
sigma^2 estimated as 50174: log likelihood=-524.15
AIC=1058.31 AICc=1059.15 BIC=1070.03
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.05272209 218.099 99.62455 -76.12583 94.97947 1.002104
ACF1
Training set 0.004756765
> mod2 <- Arima(time_series, order=c(1,0,0), seasonal=c(1,0,0), include.drift=TRUE,
+ include.mean=FALSE)
> summary(mod2)
Series: time_series
ARIMA(1,0,0)(1,0,0)[12] with drift
Coefficients:
ar1 sar1 drift
0.2014 0.2862 3.6940
s.e. 0.1205 0.1960 0.9009
sigma^2 estimated as 51238: log likelihood=-525.77
AIC=1059.53 AICc=1060.09 BIC=1068.91
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 19.14777 221.9052 107.1779 -57.65942 96.53528 1.078082
ACF1
Training set -0.01092874
所以总而言之，我的问题是：有没有办法在 R 中使用自动拟合过程来获得顶级 x arima 模型？我对任何解决方案持开放态度，但我对解决方案的偏好顺序是：
- 一种从
auto.arima 中提取前 x 个模型的模型对象的方法
- 一种在文本中提取有关前 x 个模型的信息或更好地理解跟踪参数以便自己重新创建模型的方法。
- 一些其他包或功能
在此先感谢您提供的任何帮助！
来自 forecast::auto.arima 轨迹的信息唯一确定模型。
让我们针对第三个模型明确验证这一点（使用您的示例数据）：
ARIMA(1,0,0)(1,0,0)[12] with non-zero mean : 1057.504
我们使用 base R 的 arima 手动拟合相同的模型
fit3 <- arima(
time_series,
order = c(1, 0, 0),
seasonal = list(order = c(1, 0, 0), period = 12L),
include.mean = TRUE)
使用ML估计模型参数时，arima returns是AIC值；另一方面，forecast::auto.arima 的轨迹 return 是校正后的（对于小样本）AIC。感谢this post on Cross Validated，我们可以编写一个自定义函数，根据模型参数计算 AICc。
get_AICc <- function(fit) {
npar <- length(fit$coef) + 1
nstar <- length(fit$residuals) - fit$arma[6] - fit$arma[7] * fit$arma[5]
fit$aic + 2 * npar * (nstar/(nstar - npar - 1) - 1)
}
我们确认修正后的 AIC 与 forecast::auto.arima 中报告的相同：
get_AICc(fit3)
#[1] 1057.504
forecast::Arima 自动 returns AIC 和 AICc 值，拟合结果与 arima 一致（不出所料，因为 forecast::Arima 内部调用 arima) 和 forecast::auto.arima:
Arima(
time_series,
order = c(1, 0, 0),
seasonal = c(1, 0, 0),
include.mean = TRUE)
)
#Series: time_series
#ARIMA(1,0,0)(1,0,0)[12] with non-zero mean
#
#Coefficients:
# ar1 sar1 mean
#0.1849 0.1780 179.5156
#s.e. 0.1168 0.2077 36.2321
#
#sigma^2 estimated as 49966: log likelihood=-524.47
#AIC=1056.95 AICc=1057.5 BIC=1066.32
与基本 R 的 arima 相比，forecast::Arima 允许在包含常量方面具有更大的灵活性：而 arima 允许通过 include.mean、forecast::Arima 还允许使用 time-dependent（确定性）漂移参数。有关详细信息，请参阅 Constants and ARIMA models in R。
我不知道 already-available 从 forecast::auto.arima 中获取前 n 个模型拟合列表的方法。跟踪仅将参数 space 中搜索的中间结果打印到控制台。过去，我定义并编写了自己的 SARIMA 模型参数参数网格搜索，以根据其中一个信息标准确定最佳拟合模型；然而，forecast::auto.arima 更复杂并且还考虑（即排除）单位根接近单位圆的模型拟合。您最好的选择是定义您自己的网格搜索并使其 return 成为前 n 个模型的列表；或者将forecast::auto.arima的代码修改为return这样的列表；或者（可能是最不优雅的解决方案）捕获并解析 auto.arima(... trace = TRUE) 的输出以提取模型参数和 re-fit 使用 arima 或 forecast::Arima.
更新 1
我确实说过捕获 auto.arima 的输出可能是最不优雅的方法，但这就是我最终所做的；也许这对您有帮助作为起点：
让我们定义一个函数来解析从 auto.arima
parse_SARIMA_string <- function(x) {
include.mean <- str_detect(x, "with non-zero mean")
x %>%
str_remove_all(":.+$") %>%
str_replace_all("\D", " ") %>%
trimws() %>%
enframe() %>%
separate(
value,
c("p", "d", "q", "P", "D", "Q", "period"),
sep = "\s+",
fill = "right") %>%
mutate(include.mean = include.mean)
}
然后我们可以使用capture.output来捕获forecast::auto.arima的踪迹并解析模型字符串
library(tidyverse)
capture.output(auto.arima(time_series, trace = TRUE)) %>%
str_subset("^ ARIMA") %>%
parse_SARIMA_string()
## A tibble: 13 x 9
# name p d q P D Q period include.mean
# <int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <lgl>
# 1 1 2 0 2 1 0 1 12 TRUE
# 2 2 0 0 0 NA NA NA NA TRUE
# 3 3 1 0 0 1 0 0 12 TRUE
# 4 4 0 0 1 0 0 1 12 TRUE
# 5 5 0 0 0 NA NA NA NA FALSE
# 6 6 1 0 0 NA NA NA NA TRUE
# 7 7 1 0 0 0 0 1 12 TRUE
# 8 8 1 0 0 1 0 1 12 TRUE
# 9 9 2 0 0 NA NA NA NA TRUE
#10 10 1 0 1 NA NA NA NA TRUE
#11 11 0 0 1 NA NA NA NA TRUE
#12 12 2 0 1 NA NA NA NA TRUE
#13 13 1 0 0 NA NA NA NA FALSE
现在您有了带有模型参数的 data.frame/tibble，可以直接在 arima 或 forecast::Arima.
更新 2
关于基于auto.arima模型字符串的SARIMA模型的可识别性，我们来看一些具体的模型。
有差分
模型 1 是具有常数的 ARMA 模型
模型 2 是一个有漂移的 ARIMA 模型
我们可以将模型展开得到
模型 3 是具有漂移和常数的 ARIMA 模型
当我们扩展模型 3 时，我们注意到常数项由于差异而取消，我们最终得到相同的模型 2。
因此 d=1 具有漂移和漂移+常数的模型描述了相同的过程。在这种情况下，将模型报告为例如“ARIMA(1,1,0) with drift”指的是与“ARIMA(1,1,0) with drift and non-zero mean”相同的模型。
无差异
对于 d=0，默认使用 include.drift = TRUE 会导致漂移和常数偏移参数，参见例如Constants and ARIMA models in R。如果你指定 include.drift = TRUE 然后手动强制 include.mean = FALSE 你会得到一个不同的模型（如你在 mod1 和 mod2 中所示）。据我了解，这永远不会发生在 forecast::auto.arima 中，其中 include.drift = TRUE 隐含地 总是 include.mean = TRUE。因此来自 forecast::auto.arima 的 SARIMA 字符串 ARIMA(1,0,0)(1,0,0)[12] with drift 将转换为具有漂移 和 [=136= 的 SARIMA(1,0,0)(1,0,0)[12] 模型]偏移量。