R/arma_Q0dotdotstats.R

In GeoBosh/sarima: Simulation and Prediction with Seasonal ARIMA Models

## Do not edit this file manually.
## It has been automatically generated from *.org sources.

## arma_Q0Gardner <- function(phi, theta, tol = .Machine$double.eps){
##     ## tol is not used here
##     .Call(stats:::C_getQ0, phi, theta)
## }

arma_Q0naive <- function(phi, theta, tol = .Machine$double.eps){
    p <- length(phi)
    q <- length(theta)
    r <- max(p, q+1)

    T <- cbind(numeric(r), rbind(diag(r - 1), numeric(r - 1)))
    if(p > 0)
        T[1:p, 1] <- phi

    Sigma <- matrix(0, r, r)
    Sigma[1:(q+1), 1:(q + 1)] <- c(1, theta) %o% c(1, theta)

    A <- T %x% T
    v <- solve(diag(nrow(A)) - A, as.vector(Sigma))
    matrix(v, r)
}

arma_Q0gnbR <- function(phi, theta, tol = .Machine$double.eps){
    p <- length(phi)
    q <- length(theta)
    r <- max(p, q+1)

    if(r == 1) # q = 0;  p = 0 or 1 here
        return( matrix(if(p == 1) 1/(1-phi^2) else 1.0, 1, 1) )

    rphi <- if(p == r) phi  else c(phi, numeric(r - p))

    Sigma <- matrix(0, r, r)
    Sigma[1:(q+1), 1:(q + 1)] <- c(1, theta) %o% c(1, theta)

    A <- diag(r)
    b <- numeric(r)
    for(k in seq_len(r)){
        b[k] <- Sigma[1, k]
        A[k,2] <- A[k,2] - rphi[k]
            # ak1 <- rphi[r - k + 1] * rphi[k + (r - k + 1) - 1]
        ak1 <- rphi[r - k + 1] * rphi[r]
        if(k < r){
            for(i in 1:(r - k)){
                b[k] <- b[k] + Sigma[1 + i, k + i]
                A[k, k + i] <- A[k, k + i] - rphi[i]
                if(2 + i <= r)
                    A[k, 2 + i] <- A[k, 2 + i] - rphi[k + i]
            
                ak1 <- ak1 + rphi[i] * rphi[k + i - 1]
            }
        }
        A[k, 1] <- A[k, 1] - ak1
    }

    g1 <- solve(A, b)

    res <- matrix(0.0, r, r)
    res[1, ] <- g1
    res[r, r] <- rphi[r]^2 * g1[1] + Sigma[r, r]
    if(r > 2){
        for(i in (r-1):2){
            res[i, r] <- rphi[i] * rphi[r] * g1[1] + rphi[r] * g1[i + 1] + Sigma[i, r]
            for(j in (r-1):i){
                res[i, j] <- rphi[i] * rphi[j] * g1[1] +
                    rphi[i] * g1[j + 1] + rphi[j] * g1[i + 1] +
                    Sigma[i,j] + res[i + 1, j + 1]
            }
        }
    }
    ## fill the elements under the diagonal, symmetrically;
    ## r > 1 here, so we don't check
    res[2,1] <- res[1,2]
    if(r > 2){
        for(i in (1:(r - 1)))
            for(j in (i + 1) : r)
                res[j,i] <- res[i, j]
    }
 
    res
}