R/glarmaNegBinPearson.R

In glarma: Generalized Linear Autoregressive Moving Average Models

glarmaNegBinPearson <- function(y, X, offset = NULL, delta, phiLags,
                                thetaLags, method  = "FS") {
    r <- ncol(X)
    n <- length(y)
    p <- length(phiLags)
    q <- length(thetaLags)
    s <- r + p + q + 1
    beta <- delta[1:r]
    phi <- delta[(r + 1):(r + p)]
    theta <- delta[(r + p + 1):(r + p + q)]
    alpha <- delta[s]
    u <- c(rep(0, s - 1), 1)
    mpq <- 0

    if ((p + q) > 0) {
        mpq <- max(phiLags[p], thetaLags[q])
    }

    nmpq <- n + mpq
    e <- array(0, nmpq)
    Z <- array(0, nmpq)
    W <- array(0, nmpq)
    mu <- array(0, nmpq)
    e.d <- array(0, c(s, nmpq))
    Z.d <- array(0, c(s, nmpq))
    W.d <- array(0, c(s, nmpq))
    if (method  == "NR") {
        e.dd <- array(0, c(s, s, nmpq))
        Z.dd <- array(0, c(s, s, nmpq))
        W.dd <- array(0, c(s, s, nmpq))
    }
    if(is.null(offset)) eta<-X %*% beta else eta<- X %*% beta + offset
    ll <- 0
    ll.d <- matrix(0, ncol = 1, nrow = s)
    ll.dd <- matrix(0, ncol = s, nrow = s)

    for (time in 1:n) {
        tmpq <- time + mpq

        if (p > 0) {
            Z.d[(r + 1):(r + p), tmpq] <- Z[tmpq - phiLags] + e[tmpq - phiLags]
            if (method  == "NR") {
                Z.dd[(r + 1):(r + p), , tmpq] <-
                    t((Z.d + e.d)[, (tmpq - phiLags)])
                Z.dd[, (r + 1):(r + p), tmpq] <-
                    Z.dd[, (r + 1):(r + p), tmpq] +
                        (Z.d + e.d)[, (tmpq - phiLags)]
            }
            for (i in 1:p) {
                Z[tmpq] <- Z[tmpq] + phi[i] * (Z + e)[tmpq - phiLags[i]]
                Z.d[, tmpq] <- Z.d[, tmpq] + phi[i] *
                               (Z.d[, tmpq - phiLags[i]] +
                               e.d[, tmpq - phiLags[i]])
                if (method  == "NR") {
                  Z.dd[, , tmpq] <- Z.dd[, , tmpq] + phi[i] *
                                    (Z.dd[, , tmpq - phiLags[i]] +
                                    e.dd[, , tmpq - phiLags[i]])
                }
            }
        }

        if (q > 0) {

            Z.d[(r + p + 1):(r + p + q), tmpq] <- e[tmpq - thetaLags]
            if (method  == "NR") {
                Z.dd[(r + p + 1):(r + p + q), , tmpq] <-
                    Z.dd[(r + p + 1):(r + p + q), , tmpq] +
                        t(e.d[, tmpq - thetaLags])
                Z.dd[, (r + p + 1):(r + p + q), tmpq] <-
                    Z.dd[, (r + p + 1):(r + p + q), tmpq] +
                        e.d[, tmpq - thetaLags]
            }
            for (i in 1:q) {

                Z[tmpq] <- Z[tmpq] + theta[i] * e[tmpq - thetaLags[i]]
                Z.d[, tmpq] <- Z.d[, tmpq] +
                               theta[i] * e.d[, tmpq - thetaLags[i]]
                if (method  == "NR") {
                  Z.dd[, , tmpq] <- Z.dd[, , tmpq] + theta[i] *
                                    e.dd[, , tmpq - thetaLags[i]]
                }
            }
        }

        W[tmpq] <- eta[time] + Z[tmpq]
        W.d[, tmpq] <- matrix(c(X[time, ], rep(0, p + q + 1)), ncol = 1) +
                       Z.d[, tmpq]
        if (method  == "NR") {
            W.dd[, , tmpq] <- Z.dd[, , tmpq]
        }
        mu[tmpq] <- exp(W[tmpq])

        mut <- mu[tmpq]
        yt <- y[time]
        var <- mut + mut^2/alpha


        e.W <- (-mut * (alpha * yt + mut * (alpha + 2 * yt)))/
               (2 * alpha * var^(3/2))
        e.a <- mut^2 * (yt - mut)/(2 * alpha^2 * var^(3/2))
        if (method  == "NR") {
            e.WW <- (yt * alpha^2 + mut * alpha * (2 * yt - alpha) +
                     2 * mut^2 *
                    (2 * yt + alpha))/(4 * (mut + alpha)^2 * var^(1/2))
            e.aW <- mut^3 * (alpha * yt - mut * (2 * yt + 3 * alpha))/
                    (4 * alpha^3 * var^(5/2))
            e.aa <- mut^3 * (mut - yt) * (mut + 4 * alpha)/
                    (4 * alpha^4 * var^(5/2))
        }
        e[tmpq] <- (yt - mut)/var^0.5
        e.d[, tmpq] <- e.W * W.d[, tmpq] + e.a * u
        if (method  == "NR") {
            e.dd[, , tmpq] <- (e.W * W.dd[, , tmpq] + e.WW * W.d[, tmpq] %o%
                              W.d[, tmpq] + e.aW *
                              (W.d[, tmpq] %o% u + u %o% W.d[, tmpq]) +
                              e.aa * u %o% u)
        }

        ## update likelihood and derivatives.

    #    ll <- ll + log(gamma(alpha + yt)/(gamma(alpha) * gamma(yt + 1))) +
    #          alpha * log(alpha/(alpha + mut)) + yt * log(mut/(alpha + mut))
    # altered 1/8/14 by WD to avoid overflow in gamma function call for large counts

      ll<- ll+ lgamma(alpha + yt) - lgamma(alpha) - lgamma(yt + 1) + alpha * log(alpha) + 
                  yt * log(mut + (yt == 0)) - (alpha + yt) * log(alpha + mut)

        ll.W <- alpha * (yt - mut)/(mut + alpha)
        ll.a <- (digamma(alpha + yt) - digamma(alpha) +
                (mut - yt)/(mut + alpha) + log(alpha/(mut + alpha)))
        ll.d <- ll.d + ll.W * W.d[, tmpq] + ll.a * u
        if (method  == "FS") {
            ## empirical Fisher Scoring used here
            ll.dt <- (ll.W * W.d[, tmpq] + ll.a * u)
            ll.dd <- ll.dd - ll.dt %o% ll.dt
            ## OLD ll.dd <- ll.dd + (-(alpha*mut*(yt+alpha)/(mut+alpha)^2)*W.d[,
            ## tmpq]%o%W.d[, tmpq] OLD +(trigamma(alpha+yt)-trigamma(alpha) +
            ## (yt*alpha + mut^2)/(alpha*(alpha+mut)^2))*u%o%u )
        }

        if (method  == "NR") {
            ll.dd <- ll.dd +
                     (alpha * (yt - mut)/(mut + alpha) * W.dd[, , tmpq] -
                     (alpha * mut * (yt + alpha)/(mut + alpha)^2) *
                     W.d[, tmpq] %o% W.d[, tmpq] + ((yt - mut) *
                     mut/(mut + alpha)^2) * (W.d[, tmpq] %o% u + u %o%
                     W.d[, tmpq]) + (trigamma(alpha + yt) - trigamma(alpha) +
                     (yt * alpha + mut^2)/(alpha * (alpha + mut)^2)) * u %o% u)
        }
    }

    list(delta = delta, ll = ll, ll.d = ll.d, ll.dd = ll.dd, eta = eta,
         W = W[mpq + 1:n], e = e[mpq + 1:n], mu = mu[mpq + 1:n],
         fitted.values = mu[mpq + 1:n])
}