Extract BIC and AICc from arima() object

Question

Problem: I would like to extract the BIC and AICc from an arima() object in R.

Background: The arima() function produces an output of results, which includes the estimated coefficients, standard errors, AIC, BIC, and AICc. Let's run some sample code to see what this looks like:

# 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 

The output of results for the model appears like this:

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

On the bottom line, we can see values for AIC, BIC, and AICc. (Note: this is the output shown by arima() when the forecast package has been loaded, i.e. library(forecast))

Accessing the AIC value is quite easy. One can simply type:

> model$aic
[1] 23563.39

Access to the AIC value in this manner is made possible due to the fact that it's stored as one of the model's attributes. The following code and output will make this clear:

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

$class
[1] "Arima"

Notice, however, that bic and aicc are not model attributes, so the following code is no use to us:

> model$bic
NULL
> model$aicc
NULL

The BIC and AICc values are, indeed, calculated by the arima() function, but the object that it returns does not give us direct access to their values. This is inconvenient and I've come across others who've raised the issue. Unfortunately, I've not found a solution to the problem.

Can anyone out there help? Which method can I use to access the BIC and AICc from the Arima class of object.

Note: I've suggested an answer below, but would like to hear improvements and suggestions.

Edit (Version details as requested):

> 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"

Answer 1

For the BIC and AIC, you can simply use AIC function as follow:

> 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

The function AIC can provide both AIC and BIC. Look at ?AIC.

Answer 2

Answer: One possible solution, although no claim to be the best, is as follows; it's a hack that I've come up with after looking at some source code.

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)

Now that the bic and aicc have been stored as objects - using, solely, output from the arima() function - we can now set them as attributes to the model object.

# 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

Pass on these attributes to a new object (we don't want to overwrite model).

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

We can now access the BIC and AICc in a similar manner as to how we accessed the AIC value. The following code should make this clear:

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

Edit: Based on the very useful information provided by @Stat about the AIC() function, the following code may be useful as alternative ways of getting the AIC, BIC, AICc, and HQC. Attach them as attributes to the model object and work away.

# 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))))

Answer 3

Here is a function my TA in my time series analysis course at UC Davis wrote to extract the AICc

Function aicc() computes the AICc of a given ARIMA model.

INPUT: an ARIMA model object produced by arima()

OUTPUT: AICc value for the given model object

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

Example:

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

Answer 4

Once you have loaded forecast package, you must use Arima() function for AIC, AICc and BIC. Notice upper "A" in Arima() function. If you use arima() function with lower "a", then R will use the function that comes with base R.

Answer 5

what it seems like is you are using incorrect arima function to get the values. Note that, arima() is not part of forecast library, you will have to use Arima() instead.

once the library(forecast) is imported use below function to extract the values:

model=Arima(grow, order=c(2,0,0))

attributes(model)

$names

[1] "coef" "sigma2" "var.coef" "mask" "loglik" "aic"
[7] "arma" "residuals" "call" "series" "code" "n.cond"
[13] "nobs" "model" "aicc" "bic" "x"

$class [1] "ARIMA" "Arima"

model$bic

[1] 1069.786

model=arima(grow, order=c(2,0,0))

attributes(model)

$names

[1] "coef" "sigma2" "var.coef" "mask" "loglik" "aic"
[7] "arma" "residuals" "call" "series" "code" "n.cond"
[13] "nobs" "model"

$class [1] "Arima"

Hope this helps.. :)