### PR#14682 : https://bugs.r-project.org/bugzilla3/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))