CODEWITHSUNDEEP
rtime-series
javlacalle
javlacalle

Reputation: 1049

accuracy of C_ARIMA_Gradtrans in R stats::arima

The R function stats::arima calls to C_ARIMA_transPars to transform the autoregressive coefficients and keep them in the region of stationarity.

The derivatives of this transformation are needed to obtain the standard errors of parameter estimates. They are obtained by means of C_ARIMA_Gradtrans.

In the case of an AR(2) model I was working with, the derivatives returned by the latter function do not look accurate. As shown below, they differ from the values obtained by numerical methods (based on package numDeriv) and the values obtained based on the analytical expressions that are easy to get for the AR(2) model.

Example for the AR(2) model:

Generate and transform the AR coefficients:

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

Derivatives of the transformation with respect to the parameters:

# 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]

The second and third methods match each other but differ with the result obtained with 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"

Looking at the source code, it is not straightforward to me to find out which approach or operations are involved (numerical or analytical derivatives).

It is strange that cosh (which is part of the derivative of tanh employed in the transformation) is not used. Also, a variable eps defined as 1e-3 is used, but this does not show up in the analytical derivatives. Maybe cosh is computed by other (less accurate) means?

Is there any reason explaining the less accurate results shown above?

Upvotes: 1

Views: 184

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226732

(Should be a comment but written as an answer for formatting ...)

It looks like this block is indeed computing finite differences with an epsilon of 1e-3 (defined just above):

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 */
}

I wonder if the simple analytical expressions that you can use for AR(2) break down/get too unwieldy for higher-order AR models?

I would suggest that the next debugging step (not easy, unless you're already set up for it) would be to compile a standalone version of ARIMA_Gradtrans with a smaller value of eps to see if that helps in this case ...

Or use numDeriv::grad() with method="simple", method.args=list(eps=1e-3) and see if it matches the ARIMA_Gradtrans result.

Following your code above:

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

The next question is whether a relative difference of 5e-4 is actually important for your application ...

Upvotes: 1

Related Questions