C_ARIMA_Gradtrans 在 R stats::arima 中的准确性

Question

R 函数 stats::arima 调用 C_ARIMA_transPars 来变换自回归系数并使它们保持在平稳区域。

需要这个变换的导数来获得参数估计的标准误差。它们是通过 C_ARIMA_Gradtrans.

获得的

在我使用的 AR(2) 模型的情况下，后一个函数返回的导数看起来不准确。如下所示，它们不同于数值方法（基于包numDeriv）得到的值和基于AR（2）模型容易得到的解析表达式得到的值

AR(2) 模型示例：

生成并转换 AR 系数：

require("numDeriv")
# random values for the AR coefficients
ARcoefs <- rnorm(2)
# transformation of AR coefficients
tp1 <- .Call(stats:::C_ARIMA_transPars, ARcoefs, c(2L,0L,0L,0L,0L,0L,0L), TRUE)[[1]]
# for the AR(2) model, the above reduces to:
tp2 <- tanh(ARcoefs)
tp2[1] <- tp2[1] - tp2[2] * tp2[1]
all.equal(tp1, tp2)
# [1] TRUE

关于参数的变换导数：

# based on package "stats"
g1 <- .Call(stats:::C_ARIMA_Gradtrans, ARcoefs, c(2L,0L,0L,0L,0L,0L,0L))
# based on package "numDeriv" (numerical derivatives) 
fnc <- function(x, i) .Call(stats:::C_ARIMA_transPars, 
  x, c(2L,0L,0L,0L,0L,0L,0L), TRUE)[[1]][i]
g2 <- matrix(0, nrow=2, ncol=2)
g2[,1] <- numDeriv::grad(func=fnc, x=ARcoefs, i=1)
g2[,2] <- numDeriv::grad(func=fnc, x=ARcoefs, i=2)
# based on analytical derivatives
g3 <- matrix(0, nrow=2, ncol=2)
tmp <- 1/cosh(ARcoefs)^2
g3[1,1] <- tmp[1] - tanh(ARcoefs[2]) * tmp[1]
g3[2,2] <- tmp[2]
g3[2,1] <- -tanh(ARcoefs[1]) * tmp[2]

第二种和第三种方法相互匹配，但与C_ARIMA_Gradtrans获得的结果不同。

all.equal(g2, g3)
# [1] TRUE
# output for ARcoefs = c(0.7, -0.4)
all.equal(g2, g1)
# [1] "Mean relative difference: 0.0004672373"
# output for ARcoefs = c(0.7, -0.4)
all.equal(g3, g1)
# [1] "Mean relative difference: 0.0004672373"

查看 source code，要找出涉及的方法或操作（数值或分析导数）对我来说并不简单。

奇怪的是没有使用 cosh（它是转换中使用的 tanh 的导数的一部分）。此外，使用了定义为 1e-3 的变量 eps，但这不会出现在分析导数中。也许 cosh 是通过其他（不太准确）的方式计算的？

是否有任何理由可以解释上面显示的不太准确的结果？

Answer 1

（应该是评论，但写成格式的答案...）

看起来这个块确实在计算 epsilon 为 1e-3（刚刚在上面定义）的有限差分：

for (int i = 0; i < mp; i++) w1[i] = raw[i];
partrans(mp, w1, w2);           /* transform raw (w1) to w2 */
for (int i = 0; i < mp; i++) {
    w1[i] += eps;
    partrans(mp, w1, w3);       /* w3 = trans(w1+eps) */
    for (int j = 0; j < mp; j++)  /* calc derivative */
         A[i + j*n] = (w3[j] - w2[j])/eps;
    w1[i] -= eps;               /* restore value of w1 */
}

我想知道您可以用于 AR(2) 的简单解析表达式是否会打破 down/get 对于高阶 AR 模型来说太笨重？

我建议下一个调试步骤（不容易，除非您已经为此做好了准备）是编译 ARIMA_Gradtrans 的独立版本，其值较小 eps看看在这种情况下是否有帮助...

或者使用 numDeriv::grad() 和 method="simple", method.args=list(eps=1e-3) 看看它是否匹配 ARIMA_Gradtrans 结果。

按照上面的代码：

g2A <- t(jacobian(fnc,ARcoefs))
g2B <- t(jacobian(fnc,ARcoefs,
        method="simple",method.args=list(eps=1e-3))
all.equal(g1,g2A)   ## as before
## [1] "Mean relative difference: 0.000467239"
all.equal(g1,g2B)
## [1] TRUE
all.equal(g2B,g3)
## [1] "Mean relative difference: 0.000467239"
all.equal(g2A,g3)
## [1] TRUE

下一个问题是 5e-4 的相对差异对您的应用程序是否真的很重要...

C_ARIMA_Gradtrans 在 R stats::arima 中的准确性

accuracy of C_ARIMA_Gradtrans in R stats::arima

r

time-series