从 arima() 对象中提取 BIC 和 AICc

机器算法验证 r 时间序列 有马
2022-04-18 19:11:11

问题：我想从 R 中的 arima() 对象中提取 BIC 和 AICc。

背景： arima() 函数产生结果输出，其中包括估计系数、标准误差、AIC、BIC 和 AICc。让我们运行一些示例代码来看看它是什么样子的：

# Load the sunspots dataset
data(sunspots)
# Build an ARIMA(2,0,2) model and store as an object
model <- arima(x=sunspots, order=c(2,0,2), method="ML")
# Show a summary of the model
model

模型的结果输出如下所示：

Series: sunspots 
ARIMA(2,0,2) with non-zero mean 

Coefficients:
         ar1     ar2      ma1      ma2  intercept
      0.9822  0.0004  -0.3997  -0.1135    51.2652
s.e.  0.1221  0.1196   0.1206   0.0574     8.1441

sigma^2 estimated as 247.9:  log likelihood=-11775.69
AIC=23563.39   AICc=23563.42   BIC=23599.05

在底线上，我们可以看到 AIC、BIC 和 AICc 的值。（注意：这是 arima() 在加载预测包时显示的输出，即 library(forecast)）

访问 AIC 值非常容易。可以简单地键入：

> model$aic
[1] 23563.39

以这种方式访问​​ AIC 值是可能的，因为它被存储为模型的属性之一。以下代码和输出将清楚地说明这一点：

> attributes(model)
$names
 [1] "coef"      "sigma2"    "var.coef"  "mask"      "loglik"   
 [6] "aic"       "arma"      "residuals" "call"      "series"   
[11] "code"      "n.cond"    "model"    

$class
[1] "Arima"

但是请注意，bic 和 aicc 不是模型属性，所以下面的代码对我们没有用处：

> model$bic
NULL
> model$aicc
NULL

BIC 和 AICc 值确实是由 arima() 函数计算的，但它返回的对象并没有让我们直接访问它们的值。这很不方便，我遇到过其他人提出这个问题。不幸的是，我还没有找到解决问题的方法。

有人可以帮忙吗？我可以使用哪种方法从 Arima 对象类访问 BIC 和 AICc。

注意：我在下面提出了一个答案，但想听听改进和建议。

编辑（要求的版本详细信息）：

> R.Version()
$platform
[1] "i686-pc-linux-gnu"

$arch
[1] "i686"

$os
[1] "linux-gnu"

$system
[1] "i686, linux-gnu"

$status
[1] ""

$major
[1] "3"

$minor
[1] "0.2"

$year
[1] "2013"

$month
[1] "09"

$day
[1] "25"

$`svn rev`
[1] "63987"

$language
[1] "R"

$version.string
[1] "R version 3.0.2 (2013-09-25)"

$nickname
[1] "Frisbee Sailing"
4个回答

对于 BIC 和 AIC，您可以简单地使用 AIC 函数，如下所示： 

> model <- arima(x=sunspots, order=c(2,0,2), method="ML")
> AIC(model)
[1] 23563.39
> bic=AIC(model,k = log(length(sunspots)))
> bic
[1] 23599.05

AIC 功能可以同时提供AIC 和BIC。看?AIC。

答：一种可能的解决方案，虽然没有声称是最好的，但如下；这是我在查看一些源代码后想出的一个技巧。

npar <- length(model$coef) + 1
nstar <- length(model$residuals) - model$arma[6] - model$arma[7] * model$arma[5]

bic <- model$aic + npar * (log(nstar) - 2)
aicc <- model$aic + 2 * npar * (nstar/(nstar - npar - 1) - 1)

现在 bic 和 aicc 已经被存储为对象 - 仅使用 arima() 函数的输出 - 我们现在可以将它们设置为模型对象的属性。

# Give model attributes for bic and aicc
attr(model,"bic") <- bic
attr(model,"aicc") <- aicc

> attributes(model)
$names
 [1] "coef"      "sigma2"    "var.coef"  "mask"      "loglik"   
 [6] "aic"       "arma"      "residuals" "call"      "series"   
[11] "code"      "n.cond"    "model"    

$class
[1] "Arima"

$bic
[1] 23599.05

$aicc
[1] 23563.42

将这些属性传递给一个新对象（我们不想覆盖模型）。

# Create new object with these attributes
model_2 <- attributes(model)

我们现在可以以与访问 AIC 值类似的方式访问 BIC 和 AICc。下面的代码应该清楚地说明这一点：

> model$aic
[1] 23563.39
> model_2$bic
[1] 23599.05
> model_2$aicc
[1] 23563.42

编辑：根据@Stat 提供的关于 AIC() 函数的非常有用的信息，以下代码可能可用作获取 AIC、BIC、AICc 和 HQC 的替代方法。将它们作为属性附加到模型对象并开始工作。

# AIC
AIC(arima(x=sunspots, order=c(2,0,2), method="ML"))
# BIC
AIC(arima(x=sunspots, order=c(2,0,2), method="ML"),k=log(length(sunspots)))
# AICc    
AIC(arima(x=sunspots, order=c(2,0,2), method="ML")) + 2 * npar * (nstar/(nstar - npar - 1) - 1)
# HQC
AIC(arima(x=sunspots, order=c(2,0,2), method="ML"), k=2*log(log(length(sunspots))))

这是我的 TA 在加州大学戴维斯分校的时间序列分析课程中写的一个函数，用于提取 AICc

函数 aicc() 计算给定 ARIMA 模型的 AICc。

INPUT：由 arima() 生成的 ARIMA 模型对象

输出：给定模型对象的 AICc 值

aicc = function(model){
n = model$nobs
p = length(model$coef)
aicc = model$aic + 2*p*(p+1)/(n-p-1)
return(aicc)
}

例子：

x = arima.sim(100,model=list(ar=(0.3)))
mod = arima(x,order=c(1,0,0))
aicc(mod)

加载预测包后，您必须对 AIC、AICc 和 BIC 使用 Arima() 函数。注意 Arima() 函数中的上方“A”。如果您使用带有较低“a”的 arima() 函数，那么 R 将使用 base R 附带的函数。

其它你可能感兴趣的问题
上一篇如何给范围为单位区间的随机变量添加噪声？ 下一篇HMM有什么替代品吗？