5
votes
$\begingroup$

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"
$\endgroup$

locked by whuber Nov 17 '16 at 20:32

This question exists because it has historical significance, but it is not considered a good, on-topic question for this site so please do not use it as evidence that you can ask similar questions here. This question and its answers are frozen and cannot be changed. See the help center for guidance on writing a good question.

Read more about locked posts here.

  • $\begingroup$ Which version of R is this? 3.0.2? $\endgroup$ – Glen_b Nov 16 '13 at 23:39
  • $\begingroup$ See edit. @Glen_b $\endgroup$ – Graeme Walsh Nov 17 '13 at 0:25
  • $\begingroup$ Feel free to edit and trim down the details if you like. Wasn't sure how much detail you wanted. :) $\endgroup$ – Graeme Walsh Nov 17 '13 at 0:49
  • $\begingroup$ Well, for myself I was just after the version number in the interest of checking the code for arima (the machine I was on at the time didn't reproduce the suggested behavior; it had an older version of R), but since this is intended to be a permanent repository and future versions of R will continue to change the behavior, it's hard for me to be sure what parts to take out. Someone more knowledgeable than me might be better at guessing what won't ever be relevant. $\endgroup$ – Glen_b Nov 17 '13 at 15:04
  • 2
    $\begingroup$ @Glen_b I think I understand the problem with the arima output that you encountered. The forecast package had been loaded on my machine and this has the effect of changing the output of the arima() function! Not ideal behaviour and I only noticed it today, myself. Not loading the forecast package and calling arima() returns only the aic in the output - not the BIC and AICc. Compare the outputs when the forecast package is loaded and not loaded to see what I mean. Apologies for not pointing this out earlier - I've now added a note to my answer to highlight the issue. $\endgroup$ – Graeme Walsh Nov 17 '13 at 15:20
11
votes
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Ah! Excellent! I had tried the AIC() and BIC() functions! I wasn't aware that the AIC() function could provide the BIC. Thanks very much, @Stat! This method saves a whole lot of bother. $\endgroup$ – Graeme Walsh Nov 17 '13 at 1:20
7
votes
$\begingroup$

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))))
$\endgroup$
  • $\begingroup$ This answer, rather than the ones that merely call built-in black-box functions, is the one that actually is informative and of interest on this site. $\endgroup$ – whuber Jun 10 '16 at 16:05
0
votes
$\begingroup$

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)
$\endgroup$
0
votes
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Are you suggesting there is a typo in his code which does show them after sue of arima? $\endgroup$ – mdewey Jun 10 '16 at 16:21
0
votes
$\begingroup$

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.. :)

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.