ARIMA 的 R 和 Stata 之间的主要差异

Question

使用历史 Lynx Pelt 数据 (https://www.dropbox.com/s/v0h9oywa4pdjblu/Lynxpelt.csv)，这里有两个 AIC 值表，来自 R 和 Stata，用于 ARIMA(p,q) 模型，0

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 参数估计：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 参数估计：http://pastie.org/8942232

以下是在R中创建AIC表的代码。注意我强制使用最大似然，不转换参数，并增加了最大迭代次数。

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

这是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

---

编辑：添加到参数估计的链接并在 R 代码中添加行以修复“未找到模型”错误

编辑 2：根据 iacobus 的建议，手动强制 Stata 使用 BFGS 作为优化方法。 (4,3) & (3,3) 有了很大的改进。其他值仍然存在很大差异。 (3,2) 的例子曾经匹配和现在非常不同。

STATA results with technique(bfgs):
           c1         c2         c3         c4         c5
r1  145.25614  100.20123   87.45929  77.570744  85.863777
r2  101.54848  84.916921   82.11809  86.444131  74.263937
r3  63.411671  49.424167  44.149023  40.966325  42.760294
r4  52.260723  49.196628  40.442078  43.498413  43.622292
r5  46.196192  48.195322  42.396986  42.289595          0

R results from above for easy comparison:

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

Answer 1

我认为您的数据会产生数值不稳定的似然函数，尤其是对于高阶模型。 R（至少对我而言）在某些高阶模型上向我发出警告，并且您在使用 Stata 使用不受限制的 MLE 来拟合它们时遇到问题，这一事实表明可能存在一些数值问题。 SAS 还警告我左右收敛。

如果可能性存在数值问题，这可能会影响优化步骤。默认情况下，Stata 似乎使用 Berndt-Hall-Hall-Hausman 算法使用 5 步，然后使用 BFGS 算法使用 10 步，根据需要重复组合直到收敛。另一方面，R 默认使用 BFGS。您可以使用 optim.method 参数更改它，但 R 并不像 Stata 那样轻松支持使用 BHHH 或在 BHHH 和 BFGS 之间移动。

在 R 中使用各种不同的优化器处理您的数据表明，通过在优化器之间进行更改，产生的 AIC 会发生相当大的变化。我怀疑这是导致 Stata 和 R 的估计值不同的原因。

我建议去 Stata 并设置最大化选项 BFGS（有关如何执行此操作的详细信息，请参阅http://www.stata.com/help.cgi?arima#maximize_options）。如果 Stata 的估计在做出改变后与 R 的估计一致，我不会感到惊讶。