Using historical Lynx Pelt data (https://www.dropbox.com/s/v0h9oywa4pdjblu/Lynxpelt.csv), here are two tables of AIC values from R and Stata for ARIMA(p,q) models for 0<=p<=5 and 0<=q<=5. Note that for (p,q) = (0,1), (0,2), (0,3), (1,0), (1,1), (1,2), (2,0), (2,1), (2,2), (2,3), (3,0), (3,1), (3,2), (4,0) and (4,1) the values are identical to seven significant digits. However, the remaining cases are wildly different--just look at (4,2)! The coefficient estimates are also very different when the AICs do not match. Is this a bug in the core ARIMA function, or what is going on?

AIC calculations from R for ARIMA(p,q) q0 q1 q2 q3 q4 p0 145.25613 100.20123 87.45927 77.57073 85.86376 p1 101.54847 84.91691 82.11806 77.15318 74.26392 p2 63.41165 49.42414 44.14899 40.96787 44.33848 p3 52.26069 49.19660 52.00560 43.50156 45.17175 p4 46.19617 48.19530 49.50422 42.43198 45.71375

R parameter estimates: http://pastie.org/8942238

AIC ( Stata ) FOR LOG MODELS q p 0 1 2 3 4 0 100.2012 87.45929 77.57074 83.86378 1 101.5485 84.91692 82.11809 86.44413 74.26394 2 63.41167 49.42417 44.14902 40.96633 40.76029 3 52.26072 49.19663 52.00562 40.37268 42.20399 4 46.19619 48.19532 40.39699 43.12795 na

Stata parameter estimates: http://pastie.org/8942232

Below is the code for creating the AIC table in R. Note that I force the used of Maximum Likelihood, no transformation of parameters, and increased the maximum iterations.

pelts <- read.csv("Lynxpelt.csv") pelts$log <- log(pelts$W7) models <- array(list(),5) aic <- data.frame(q0=rep(NA,5), q1=rep(NA,5), q2=rep(NA,5), q3=rep(NA,5), q4=rep(NA,5), row.names=c("p0", "p1", "p2", "p3", "p4")) makeModel <- function(p,q) { arima(pelts$log, order=c(p,0,q), transform.pars=FALSE, method="ML", optim.control=list(maxit=1000)) } options(warn=1) for (p in 0:4) { for (q in 0:4) { model <- makeModel(p,q) models[[p+1]][[q+1]] <- model aic[p+1,q+1] <- model$aic print(cat("p=",p,", q=",q)) } } aic

And here's the code for Stata:

insheet using Lynxpelt.csv save Lynxpelt, replace tsset year tsline w7 gen logw7=log(w7) label var logw7 "logarithm of w7" mat A=J(5,5,0) /*This matrix is a 5*5 matrix with 0s*/ mat list A /*show the matrix A*/ forvalues i=0/4 { forvalues j=0/4 { set more off quietly arima logw7, arima(`i',0,`j') estat ic matrix list r(S) matrix s=r(S) scalar alpha=s[1,5] mat A[`i'+1,`j'+1]=alpha } } * ARMA(4,4) cannot be done since stata cannot choose an initial value - we give one manually * * I will use the estimates from ARMA(3,4) * * Let's run ARMA(3,4) again * quietly arima logw7, ar(1/3) ma(1/4) matrix list e(b) mat B=e(b) *Now, let's run ARMA(4,4) with initial values from ARMA(3,4) * quietly arima logw7, ar(1/4) ma(1/4) from(B) estat ic matrix s=r(S) scalar alpha=s[1,5] mat A[5,5]=alpha

Edit: added links to parameter estimates & added line to R code to fix "models not found" error

Edit 2: Upon iacobus's advice, manually forced Stata to use BFGS as the optimization method. The (4,3) & (3,3) are much improved. Other values still differ wildly. The (3,2) for example used to match and now is very different.