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

Posted 2023-02-19

【中文标题】ARIMA 的 R 和 Stata 之间的主要差异

【英文标题】：Major discrepancies between R and Stata for ARIMA

【发布时间】：2014-04-22 00:15:22

【问题描述】：

使用历史 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

【问题讨论】：

我不使用Stata，但也许从R中提取每个模型的对数似然和每个模型的参数数量并自己计算AIC。然后检查您的 AIC 值是否与 R 报告的值匹配。这可能是第一步。

感谢马克的建议。 AIC 计算正确。确实，R 和 Stata 之间的参数估计不同，因此导致 AIC 不同。我将 AIC 用于表格，因为它更容易一目了然地注意到某些 p,q 的回归结果中的巨大差异

@tbenst 您介意添加参数估计比较吗？我目前无法访问 Stata。

请使示例可重现。我在尝试复制时收到Error in models[[p + 1]][[q + 1]] <- model (from #4) : object 'models' not found。

@user12202013：添加了参数估计的链接。在每个链接中使用特定的 aic 执行 control-f 以查看比较。

【参考方案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 的估计一致，我不会感到惊讶。

【讨论】：

这非常有用。最好（1）可视化似然面和（2）找出哪个答案实际上最接近正确 - 即BFGS或BHHH（或其他东西）是否给出最佳答案？ [在同一平台/包中比较可能性/AIC 值很容易，但比较起来可能会很棘手......] 通过arma 调试并提取@987654324 可能会更深入地解决这个问题@函数进一步探索...

有见地，谢谢。编辑 2 中的新结果。 (4,3) 和 (3,3) 得到了很大改进。其他值仍然存在很大差异。 (3,2) 例如曾经匹配和现在非常不同。

当您更改方法时，某些估计值收敛，这表明优化存在数值问题。有些价值观仍然不一致，这很有趣。我怀疑这仍然是由于该数据的似然函数的数值不稳定性。如果函数是平坦的或数值不稳定，则结果可能很大程度上取决于起始值和其他参数。我对 STATA 的了解还不够，无法了解它们与 R 的默认值相比如何，但我敢打赌这是问题所在。