266 lines
8.8 KiB
R
266 lines
8.8 KiB
R
### PR#14682 : https://bugs.r-project.org/show_bug.cgi?id=14682
|
|
## ========
|
|
## Subject: getQ0() returns a non-positive covariance matrix
|
|
## Date: Tue, 20 Sep 2011 12:06:16 -0400
|
|
## ReportedBy: raphaelrossignol@...
|
|
|
|
## ...........
|
|
|
|
## 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
|
|
## (eq. (3.3.8), (3.3.9), p.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, tol=.Machine$double.eps) {
|
|
## 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,r-q-1))
|
|
|
|
A1 <- matrix(0,r,p)
|
|
C <- (col(A1)+row(A1)-1)
|
|
B <- (C <= p) ## == (col(A1)+row(A1) <= p+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 <- C2 <- matrix(0,r2,r2)
|
|
F <- row(C1)-col(C1)+1
|
|
E <- (1 <= F) & (F <= p+1)
|
|
C1[E] <- tphi[F[E]]
|
|
|
|
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, tol=tol)
|
|
## --------- --
|
|
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 <- 1 <= F & 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)
|
|
} else rXZ <- numeric()
|
|
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)))
|
|
}
|
|
A %*% S %*% t(A)
|
|
## == 2 x 2 Block matrix product; A = [A1 | A2 ]
|
|
## == A1 SX A1' + A1 SXZ A2' + (A1 SXZ A2')' + A2 A2'
|
|
|
|
}## {Q0bis}
|
|
|
|
## 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 <- V <- matrix(0,r,r)
|
|
if (p) T[1:p,1] <- phi
|
|
if (r >= 2) T[1:(r-1),2:r] <- diag(r-1)
|
|
ttheta <- c(1,theta)
|
|
V[1:(q+1),1:(q+1)] <- ttheta %x% t(ttheta)
|
|
S <- diag(r*r) - T %x% T
|
|
Q0 <- solve(S, c(V))
|
|
matrix(Q0, ncol=r)
|
|
}
|
|
|
|
Q0.orig <- function(phi,theta) .Call(stats:::C_getQ0, phi, theta)
|
|
|
|
Q0bisC <- function(phi,theta, tol=.Machine$double.eps)
|
|
.Call(stats:::C_getQ0bis, phi, theta, tol=tol)
|
|
|
|
##' The k smallest eigenvalues of m
|
|
EV.k <- function(m, k = 2) {
|
|
ev <- eigen(m, only.values=TRUE)$values
|
|
m <- length(ev)
|
|
ev[m:(m-k+1)]
|
|
}
|
|
|
|
chkQ0 <- function(phi,theta, tol=.Machine$double.eps^0.5,
|
|
tolC=1e-15, strict=TRUE, doEigen=FALSE)
|
|
{
|
|
Q0 <- Q0.orig(phi, theta)
|
|
Q0bis <- Q0bis(phi, theta)
|
|
Q0ter <- Q0ter(phi, theta)
|
|
|
|
eig <- if(doEigen) rbind("0" = EV.k(Q0), bis = EV.k(Q0bis), ter = EV.k(Q0ter))
|
|
## else NULL
|
|
|
|
a.eq <- list(cRC = all.equal(Q0bis,Q0bisC(phi,theta), tol= tolC),
|
|
c12 = all.equal(Q0, Q0bis, tol=tol),
|
|
c13 = all.equal(Q0, Q0ter, tol=tol),
|
|
c23 = all.equal(Q0bis,Q0ter, tol=tol))
|
|
if(strict) do.call(stopifnot, a.eq)
|
|
invisible(list(Q0 = Q0, Q0bis = Q0bis, Q0ter = Q0ter,
|
|
all.eq = a.eq, eigen = eig))
|
|
}
|
|
|
|
##' @title AR-phi corresponding to AR(1) + Seasonality(s)
|
|
##' @param s: seasonality
|
|
##' @param phi1, phis: phi[1], phi[s] .. defaults: close to non-stationarity
|
|
mkPhi <- function(s, phi1 = 0.0001, phis = 0.99) {
|
|
stopifnot(length(s) > 0, s == as.integer(s), s >= 2,
|
|
length(phi1) == 1, is.numeric(phi1), length(phis) == 1)
|
|
c(phi1, rep(0, s-2), phis, -phi1*phis)
|
|
}
|
|
|
|
##--{end of function defs}-------------------------------------------------------
|
|
|
|
## cases with p=0, q=0 :
|
|
chkQ0(numeric(), numeric())
|
|
chkQ0( .5, numeric())
|
|
chkQ0(numeric(), .7)
|
|
chkQ0(numeric(), c(.7, .2))
|
|
|
|
chkQ <- function(s, theta) chkQ0(mkPhi(s=s), theta=theta, tol = 0, strict=FALSE)
|
|
all.eq2num <- function(ae) as.numeric(sub(".* difference: ", '', ae))
|
|
getN12 <- function(r) all.eq2num(r$all.eq$c12)
|
|
ss <- setNames(,2:20)
|
|
chk0 <- lapply(ss, chkQ, theta= numeric())
|
|
chk1 <- lapply(ss, chkQ, theta= 0.75)
|
|
chk2 <- lapply(ss, chkQ, theta= c(0.75, -0.5))
|
|
chks <- list(q0 = chk0, q1 = chk1, q2 = chk2)
|
|
## Quite platform dependent, in F19, 32 bit looks slightly better than 64:
|
|
(re <- sapply(chks, function(C) sapply(C, getN12)))
|
|
matplot(ss, re, type = "b", log="y", pch = paste(0:2))
|
|
stopifnot(re[paste(2:7),] < 1e-7, # max(.) seen 9.626e-9
|
|
re < 0.9) # max(.) seen 0.395
|
|
|
|
## The smallest few eigen values:
|
|
round(t(sapply(lapply(chk1, `[[`, "Q0"), EV.k, k=3)), 3)
|
|
ev3.0 <- lapply(chks, function(ck) t(sapply(lapply(ck, `[[`, "Q0"), EV.k, k=3)))
|
|
lapply(ev3.0, round, digits=3) ## problem for q >= 1 (none for q=0)
|
|
ev3.bis <- lapply(chks, function(ck) t(sapply(lapply(ck, `[[`, "Q0bis"), EV.k, k=3)))
|
|
lapply(ev3.bis[-1], round, digits=3) ## all fine
|
|
e1.bis <- sapply(ev3.bis, function(m) m[,1])
|
|
min(e1.bis) # -7.1e-15 , -7.5e-15
|
|
stopifnot(e1.bis > -1e-12)
|
|
|
|
|
|
## Now Rossignol's example
|
|
phi <- mkPhi(s = 12)
|
|
theta <- 0.7
|
|
true.cf <- c(ar1=phi[1], ma1=theta, sar1=phi[12])
|
|
tt <- chkQ0(phi,theta, tol=0.50, doEigen=TRUE)
|
|
tt$eigen
|
|
|
|
out.0 <- makeARIMA(phi, theta, NULL)
|
|
out.R <- makeARIMA(phi, theta, NULL, SSinit="Rossignol")
|
|
|
|
set.seed(7)
|
|
x <- arima.sim(1000,model=list(ar=phi,ma=theta))
|
|
str(k0 <- KalmanLike(x, mod=out.0))
|
|
str(kS <- KalmanLike(x, mod=out.R))
|
|
stopifnot(sapply(kS, is.finite))
|
|
|
|
ini.ph <- true.cf
|
|
## Default method = "CSS-ML" works fine
|
|
fm1 <- arima(x, order= c(1,0,1), seasonal= list(period=12, order=c(1,0,0)),
|
|
include.mean=FALSE, init=ini.ph)
|
|
stopifnot(all.equal(true.cf, coef(fm1), tol = 0.05))
|
|
|
|
## Using 'ML' seems "harder" :
|
|
e1 <- try(
|
|
arima(x, order= c(1,0,1), seasonal= list(period=12, order=c(1,0,0)),
|
|
include.mean=FALSE, init=ini.ph, method='ML')
|
|
)
|
|
## Error: NAs in 'phi'
|
|
e2 <- try(
|
|
arima(x, order= c(1,0,1), seasonal= list(period=12, order=c(1,0,0)),
|
|
include.mean=FALSE, init=ini.ph, method='ML', transform.pars=FALSE)
|
|
)
|
|
## Error in optim(init[mask], armafn, ..): initial value in 'vmmin' is not finite
|
|
|
|
## MM: The new Q0 does *not* help here, really:
|
|
e3 <- try(
|
|
arima(x, order= c(1,0,1), seasonal= list(period=12, order=c(1,0,0)),
|
|
include.mean=FALSE, init=ini.ph, method='ML', SSinit = "Rossi")
|
|
)
|
|
## actually fails still, but *not* transforming parameters works :
|
|
fm2 <-
|
|
arima(x, order= c(1,0,1), seasonal= list(period=12, order=c(1,0,0)),
|
|
include.mean=FALSE, init=ini.ph, method='ML', SSinit = "Rossi", transform.p=FALSE)
|
|
|
|
stopifnot(all.equal(confint(fm1),
|
|
confint(fm2), tol = 4e-4))
|
|
|
|
###---------- PR#16278 --------------------------------------
|
|
|
|
## xreg *and* differentiation order d >= 1 :
|
|
set.seed(0)
|
|
n <- 5
|
|
x <- cumsum(rnorm(n, sd=0.01))
|
|
Vr <- var(diff(x)) # 6.186e-5 : REML
|
|
V. <- var(diff(x)) * (n-2) / (n-1) # 4.640e-5 : ML
|
|
|
|
f00 <- arima0(x, c(0,1,0), method="ML", xreg=1:n)
|
|
(fit1 <- arima (x, c(0,1,0), method="ML", xreg=1:n))
|
|
stopifnot(all.equal(fit1$sigma2, V.), fit1$nobs == n-1,
|
|
all.equal(fit1$loglik, 14.28, tol=4e-4),
|
|
all.equal(f00$sigma2, fit1$sigma2),
|
|
all.equal(f00$loglik, fit1$loglik))
|
|
|
|
(fit2 <- arima (x, c(0,2,0), method="ML", xreg=(1:n)^2))
|
|
stopifnot(all.equal(fit2$sigma2, 0.000109952342),
|
|
all.equal(fit2$loglik, 9.4163797), fit2$nobs == n-2)
|
|
|
|
## "well"-fitting higher order model {optim failed in R <= 3.0.1)
|
|
n <- length(x. <- c(1:4,3:-2,2*(0:3),4:5,5:-4)/32)
|
|
xr <- poly(x., 3)
|
|
x. <- cumsum(cumsum(cumsum(x.))) + xr %*% 10^(0:2)
|
|
(fit3 <- arima (x., c(0,3,0), method="ML", xreg = xr))
|
|
stopifnot(fit3$ nobs == n-3,
|
|
all.equal(fit3$ sigma2, 0.00859843, tol = 1e-6),
|
|
all.equal(fit3$ loglik, 22.06043, tol = 1e-6),
|
|
all.equal(unname(coef(fit3)),
|
|
c(0.70517, 9.9415, 100.106), tol = 1e-5))
|
|
|
|
x.[5:6] <- NA
|
|
(fit3N <- arima (x., c(0,3,0), method="ML", xreg = xr))
|
|
stopifnot(fit3N$ nobs == n-3-2, # == #{obs} - d - #{NA}
|
|
all.equal(fit3N$ sigma2, 0.009297345, tol = 1e-6),
|
|
all.equal(fit3N$ loglik, 16.73918, tol = 1e-6),
|
|
all.equal(unname(coef(fit3N)),
|
|
c(0.64904, 9.92660, 100.126), tol = 1e-5))
|
|
|