[R-SIG-Finance] The C function getQ0 returns a non-positive covariance matrix and causes errors in arima()

[Raphael Rossignol]

Hi,

I sent this message a few days ago to the r-devel list, but the only  
answer I got back was to advise me to post it here. So here it is.

The function makeARIMA(), designed to construct some state space
representation of an ARIMA model, uses a C function called getQ0,
which can be found at the end of arima.c in R source files (library
stats). getQ0 takes two arguments, phi and theta, and returns the
covariance matrix of the state prediction error at time zero. The
reference for getQ0 (cited by help(arima)) is:
Gardner, G, Harvey, A. C. and Phillips, G. D. A. (1980) Algorithm
       AS154. An algorithm for exact maximum likelihood estimation of
       autoregressive-moving average models by means of Kalman filtering.
       _Applied Statistics_ *29*, 311-322.
where it is called subroutine STARMA (and coded in fortran 77).

My problem is that getQ0 returns incorrect covariance matrices in
certain cases. Indeed, below is an example of a
SARIMA(1,0,1)x(1,0,0)_12 where getQ0 gives a covariance matrix which
possess negative eigenvalues ! Below, I obtain getQ0 results through
makeARIMA().
Example:
> s <- 12
> phis <- 0.95
> phi1 <- 0.0001
> phi <- c(phi1,rep(0,s-2),phis,-phi1*phis)
> theta <- 0.7
> out <- makeARIMA(phi,theta,NULL)
> min(eigen(out$Pn)$value)
[1] -19.15890

There are consequences of this "bug" in the functions KalmanLike() and
arima(). Indeed, arima() in its default behaviour uses first CSS
method to get the initial value for an MLE search through optim. To
compute the likelihood, it uses getQ0 at the initialization of the
Kalman Filter. It may happen that getQ0 returns a covariance matrix
which possesses negative eigenvalues and that this gives a negative
gain in the Kalman filter, which in turn gives a likelihood equal to
NaN. Here is a reproducible example where I forced the use of 'ML'.
> set.seed(1)
> x <- arima.sim(100,model=list(ar=phi,ma=theta))
> KalmanLike(x,mod=out,fast=FALSE)
$Lik
ssq
NaN

$s2
       ssq
0.971436

> arima(x,order=c(1,0,1),seasonal=list(period=12,order=c(1,0,0)),include.mean=FALSE,init=c(phi1,theta,phis),method='ML')
Erreur dans optim(init[mask], armafn, method = optim.method, hessian =
TRUE,  :
    valeur non-finie fournie par optim

If needed, I can send a more natural example in which the above
behaviour is obtained. This error message ("Error in optim ...
non-finite finite-difference value") was already noted in the
following message, which remained without answer:
https://stat.ethz.ch/pipermail/r-devel/2009-February/052009.html

I could not figure out whether there is a real bug in getQ0 or if this
is due to some numerical instability. However, I tried to replace
getQ0 in two ways. The first one is to compute first the covariance
matrix of (X_{t-1},...,X_{t-p},Z_t,...,Z_{t-q}) and this is achieved
through the method of difference equations (page 93 of Brockwell and
Davis). This way was apparently suggested by a referee to Gardner et
al. paper (see page 314 of their paper).
Q0bis <- function(phi,theta){
    ## Computes the initial covariance matrix for the state space
representation of Gardner et al.
    p <- length(phi)
    q <- length(theta)
    r <- max(p,q+1)
    ttheta <- c(1,theta,rep(0,max(p,q+1)-q-1))
    A1 <- matrix(0,r,p)
    B <- (col(A1)+row(A1)<=p+1)
    C <- (col(A1)+row(A1)-1)
    A1[B] <- phi[C[B]]
    A2 <- matrix(0,r,q+1)
    C <- (col(A2)+row(A2)-1)
    B <- (C<=q+1)
    A2[B] <- ttheta[C[B]]
    A <- cbind(A1,A2)
    if (p==0) {
      S <- diag(q+1)
    }
    else {
      ## Compute the autocovariance function of U, the AR part of X
      r2 <- max(p+q,p+1)
      tphi <- c(1,-phi)
      C1 <- matrix(0,r2,r2)
      F <- row(C1)-col(C1)+1
      E <- (1<=F)&(F<=p+1)
      C1[E] <- tphi[F[E]]
      C2 <- matrix(0,r2,r2)
      F <- col(C2)+row(C2)-1
      E <- (F<=p+1) & col(C2)>=2
      C2[E] <- tphi[F[E]]
      Gam <- (C1+C2)
      g <- matrix(0,r2,1)
      g[1] <- 1
      rU <- solve(Gam,g)
      SU <- toeplitz(rU[1:(p+q),1])
      ## End of the difference equations method
      ## Then, compute correlation matrix of X
      A2 <- matrix(0,p,p+q)
      C <- col(A2)-row(A2)+1
      B <- (1<=C)&(C<=q+1)
      A2[B] <- ttheta[C[B]]
      SX <- A2%*%SU%*%t(A2)
      ## Now, compute correlation matrix between X and Z
      C1 <- matrix(0,q,q)
      F <- row(C1)-col(C1)+1
      E <- (F>=1) & (F<=p+1)
      C1[E] <- tphi[F[E]]
      g <- matrix(0,q,1)
      if (q) g[1:q,1] <- ttheta[1:q]
      rXZ <- forwardsolve(C1,g)
      SXZ <- matrix(0, p, q+1)
      F <- col(SXZ)-row(SXZ)
      E <- F>=1
      SXZ[E] <- rXZ[F[E]]
      S <- rbind(cbind(SX,SXZ),cbind(t(SXZ),diag(q+1)))
    }
    Q0 <- A%*%S%*%t(A)
}

The second way is to resolve brutally the equation of Gardner et al.
in the form (12), page 314 of their paper.

Q0ter <- function(phi,theta){
    p <- length(phi)
    q <- length(theta)
    r <- max(p,q+1)
    T <- matrix(0,r,r)
    if (p) T[1:p,1] <- phi
    if (r) T[1:(r-1),2:r] <- diag(r-1)
    V <- matrix(0,r,r)
    ttheta <- c(1,theta)
    V[1:(q+1),1:(q+1)] <- ttheta%x%t(ttheta)
    V <- matrix(V,ncol=1)
    S <- diag(r*r)-T%x%T
    Q0 <- solve(S,V)
    Q0 <- matrix(Q0,ncol=r)
}

My conclusion for now is that these two other ways give the same
answers (as returned by all.equal) and sometimes different answers
than getQ0. It may happen that they give a Q0 with negative
eigenvalues, but they are very very small, and then, the likelihood
computed by KalmanLike is a number (and not NaN). Here is a function
allowing to compare the three methods
test <- function(phi,theta){
    out <- makeARIMA(phi,theta,NULL)
    Q0bis <- Q0bis(phi,theta)
    Q0ter <- Q0ter(phi,theta)
    mod <- out
    modbis <- out
    modter <- out
    modbis$Pn <- Q0bis
    modter$Pn <- Q0ter
    set.seed(1)
    x <- arima.sim(100,model=list(ar=phi,ma=theta))
    s <- KalmanLike(x,mod=mod,fast=FALSE)
    sbis <- KalmanLike(x,modbis)
    ster <- KalmanLike(x,modter)
    test12 <- all.equal(out$Pn,Q0bis)
    test13 <- all.equal(out$Pn,Q0ter)
    test23 <- all.equal(Q0bis,Q0ter)

list(eigen=min(eigen(out$Pn)$value),eigenbis=min(eigen(Q0bis)$value),eigenter=min(eigen(Q0ter)$value),test12=test12,test13=test13,test23=test23,s=s,sbis=sbis,ster=ster)
}

And the test on the values of phi and theta above:
> test(phi,theta)
$eigen
[1] -19.15890

$eigenbis
[1] -9.680598e-23

$eigenter
[1] 1.859864e-23

$test12
[1] "Mean relative difference: 0.4255719"

$test13
[1] "Mean relative difference: 0.4255724"

$test23
[1] TRUE

$s
$s$Lik
ssq
NaN

$s$s2
       ssq
0.971436


$sbis
$sbis$Lik
        ssq
0.1322859

$sbis$s2
        ssq
0.9789775


$ster
$ster$Lik
        ssq
0.1322859

$ster$s2
        ssq
0.9789775



Here are my questions:
1) Does someone understand this strange behaviour of Q0 ?
2) Should I report this as a bug ?

By the way, Q0bis is only twice slower than makeARIMA() (but it is not
optimised), and Q0ter is 50 times slower than Q0bis.

For information:
> sessionInfo()

R version 2.10.1 (2009-12-14)
x86_64-pc-linux-gnu

locale:
   [1] LC_CTYPE=fr_FR.utf8       LC_NUMERIC=C
   [3] LC_TIME=fr_FR.utf8        LC_COLLATE=fr_FR.utf8
   [5] LC_MONETARY=C             LC_MESSAGES=fr_FR.utf8
   [7] LC_PAPER=fr_FR.utf8       LC_NAME=C
   [9] LC_ADDRESS=C              LC_TELEPHONE=C
[11] LC_MEASUREMENT=fr_FR.utf8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base


Best,

Raphael Rossignol

--
Assistant Professor of Mathematics
Univ. Paris-Sud 11