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

[Raphael Rossignol]

Hi,

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