The R Project SVN
R
(root)/branches/R-ucrt/src/library/stats/tests/arimaML.R – Rev 81322
Rev 80719 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
### 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 Xr2 <- max(p+q, p+1)tphi <- c(1,-phi)C1 <- C2 <- matrix(0,r2,r2)F <- row(C1)-col(C1)+1E <- (1 <= F) & (F <= p+1)C1[E] <- tphi[F[E]]F <- col(C2)+row(C2)-1E <- (F <= p+1) & col(C2) >= 2C2[E] <- tphi[F[E]]Gam <- C1 + C2g <- matrix(0,r2,1)g[1] <- 1rU <- solve(Gam, g, tol=tol)## --------- --SU <- toeplitz(rU[1:(p+q),1])## End of the difference equations method## Then, compute correlation matrix of XA2 <- matrix(0,p,p+q)C <- col(A2)-row(A2)+1B <- (1 <= C) & (C <= q+1)A2[B] <- ttheta[C[B]]SX <- A2 %*% SU %*% t(A2)## Now, compute correlation matrix between X and ZC1 <- matrix(0,q,q)F <- row(C1)-col(C1)+1E <- 1 <= F & F <= p+1C1[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 >= 1SXZ[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] <- phiif (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% TQ0 <- 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 mEV.k <- function(m, k = 2) {ev <- eigen(m, only.values=TRUE)$valuesm <- 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 NULLa.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-stationaritymkPhi <- 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-9re < 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 finee1.bis <- sapply(ev3.bis, function(m) m[,1])min(e1.bis) # -7.1e-15 , -7.5e-15stopifnot(e1.bis > -1e-12)## Now Rossignol's examplephi <- mkPhi(s = 12)theta <- 0.7true.cf <- c(ar1=phi[1], ma1=theta, sar1=phi[12])tt <- chkQ0(phi,theta, tol=0.50, doEigen=TRUE)tt$eigenout.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 finefm1 <- 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 <- 5x <- cumsum(rnorm(n, sd=0.01))Vr <- var(diff(x)) # 6.186e-5 : REMLV. <- var(diff(x)) * (n-2) / (n-1) # 4.640e-5 : MLf00 <- 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))