ARIMA模型在R中的奇妙运用案例

Question

在使用R中的函数arma{tseries}和arima{stats}拟合ARMA模型时，我发现了一些奇怪的事情。

这两个函数采用的估计过程有根本性的不同，arima{stats}中采用卡尔曼滤波器，而arma{tseries}中采用的是最大似然估计。

考虑到两个函数之间估计过程的差异，如果我们使用相同的时间序列，不应该预期两个函数的结果会有根本性的不同。

但事实证明它们确实存在巨大的差异！

生成以下时间序列并添加2个异常值。

set.seed(1010)
ts.sim <- abs(arima.sim(list(order = c(1,0,0), ar = 0.7), n = 50))
ts.sim[8] <- ts.sim[12]*8
ts.sim[35] <- ts.sim[32]*8

使用这两个函数拟合ARMA模型。

# Works perfectly fine
arima(ts.sim, order = c(1,0,0))
# Works perfectly fine
arma(ts.sim, order = c(1,0))

将时间序列的级别乘以十亿倍

# Introduce a multiplicative shift
ts.sim.1 <- ts.sim*1000000000
options(scipen = 999)
summary(ts.sim.1)

使用以下两个函数拟合ARMA模型：

# Works perfectly fine
arma(ts.sim.1, order = c(1,0))

# Does not work
arima(ts.sim.1, order = c(1,0,0))

## Error in solve.default(res$hessian * n.used, A): system is 
 computationally singular: reciprocal condition number = 1.90892e-19

当我使用SAS软件成功运行proc x12程序进行季节性测试时，我发现在R上进行相同操作却出现了以上错误。这让我非常怀疑SAS的结果，但事实证明，可能只是与arima{stats}有关。

有人能否详细说明以上错误的原因，以及它如何限制我们使用arima{stats}来拟合模型？

- Shreyes

Answer 1

在计算系数的协方差矩阵时，stats::arima函数存在问题。由于数字过大，代码对比例影响不太健壮，并且在计算Hessian矩阵逆矩阵时在此行崩溃：

var <- crossprod(A, solve(res$hessian * n.used, A))

通过简单地对数据进行缩放就可以避免这个问题。例如

arima(ts.sim.1/100, order = c(1,0,0))

会工作。

tseries::arma函数并不完美地工作。它会返回一个警告信息：

In arma(ts.sim.1, order = c(1, 0)): Hessian negative-semidefinite

这可以通过缩放来避免：

arma(ts.sim.1/1000, order = c(1,0))

- Rob Hyndman

---

感谢您的评论，Rob。我说arma()的结果非常好的原因是，我的beta估计系数基本保持不变，而我的截距增加了十亿，这在直觉上是有意义的。您是正确的，缩放可以解决我们在arima()函数中看到的问题，但仍然令人惊讶的是，即使没有使用arima()，估计也没有发生。 - Shreyes

如果我将初始数据缩小10倍，那么我是否应该将Arima产生的残差放大10倍？ - Frank

这取决于您想如何处理残差。如果您只是检查它们是否为白噪声，则比例无关紧要。如果您想使用它们来计算一步预测方差（例如），那么您需要重新调整它们的比例。 - Rob Hyndman