R/utilityfcns.R

Defines functions getnames vcov.arfima BIC.arfima AIC.arfima coef.arfima logLik.arfima InfoMatrix print.summary.arfima summary.arfima resid.arfima fitted.arfima print.arfima print.predarfima BIC removeMode

Documented in AIC.arfima BIC BIC.arfima coef.arfima fitted.arfima logLik.arfima print.arfima print.predarfima print.summary.arfima removeMode resid.arfima summary.arfima vcov.arfima

#' Removes a mode from an \code{arfima} fit.
#'
#'
#' This function is useful if one suspects a mode is spurious and does not want
#' to call the weed function.
#'
#' @param object An object of class "arfima".
#' @param num The number of the mode as in the printed value of the object.
#' @return The original object with the mode removed.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima}}
#' @keywords ts
#' @examples
#'\donttest{
#' set.seed(8765)
#' sim <- arfima.sim(1000, model = list(phi = 0.4, theta = 0.9, dfrac = 0.4))
#' fit <- arfima(sim, order = c(1, 0, 1), back=TRUE)
#' fit
#' fit <- removeMode(fit, 3)
#' fit
#'}
#' @export removeMode
removeMode <- function(object, num) {
    if (!object$weeded)
        stop("Please weed the object first before manually removing a mode.")
    m <- length(object$modes)
    if (num > m)
        stop("incorrect mode specification")
    modes <- object$modes
    modes[[num]] <- NULL
    object$modes <- modes
    object <- weed(object, eps2 = 0)
    object
}

BIC <- function(object, ...) {
    UseMethod("BIC")
}



#' Prints predictions and prediction intervals
#'
#' Prints the output of \code{predict} on an \code{arfima} object
#'
#' Prints all the relavent output of the prediction function of the
#' \code{arfima} package
#'
#' @param x An object of class "predarfima"
#' @param digits The number of digits to print
#' @param \dots Currently not used
#' @return \code{x} is returned invisibly
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima}}, \code{\link{predict.arfima}},
#' \code{\link{predict}}, \code{\link{plot.predarfima}}
#' @keywords ts
#' @examples
#'
#' set.seed(82365)
#' sim <- arfima.sim(1000, model = list(dfrac = 0.4, theta=0.9, dint = 1))
#' fit <- arfima(sim, order = c(0, 1, 1), back=TRUE)
#' fit
#' pred <- predict(fit, n.ahead = 5)
#' pred
#' plot(pred)
#'
print.predarfima <- function(x, digits = max(6, getOption("digits") - 3), ...) {

    n <- length(x$z)

    m <- length(x) - arfima_const_predvarnum

    seed <- x$seed
    limiting <- x$limiting

    predint <- x$predint
    m <- x$m
    ret <- vector("list", m)
    nam <- paste("Mode", 1:m)
    for (i in 1:m) {
        exact <- (length(x[[i]]$exactSD)>0)

        ans <- rbind(as.numeric(x[[i]]$Forecast), x[[i]]$exactSD, x[[i]]$limitSD)
        n.ahead <- dim(ans)[2]
        rownames(ans) <- c("Forecasts",
                           if (exact) "Exact SD" else NULL,
                           if (limiting) "Limiting SD" else NULL)
        colnames(ans) <- 1:n.ahead

        ret[[i]] <- list(`Forecasts and SDs` = ans)
    }
    names(ret) <- nam

    print(ret, digits = digits, ...)

    invisible(x)
}



#' Prints a Fitted Object
#'
#' Prints a fitted \code{arfima} object's relevant details
#'
#'
#' @param x A fitted \code{arfima} object
#' @param digits The number of digits to print
#' @param \dots Optional arguments.  See \code{\link{print}}.
#' @return The object is returned invisibly
#' @author JQ (Justin) Veenstra
#' @references Veenstra, J.Q. Persistence and Antipersistence:  Theory and
#' Software (PhD Thesis)
#' @keywords ts
print.arfima <- function(x, digits = max(6, getOption("digits") - 3), ...) {
    dmean <- x$dmean
    itmean <- x$itmean

    if (!x$weeded) {
        warning("The object has not been weeded:  there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
    }
    ans <- x
    cat("Number of modes:", length(ans$modes), "\n")
    coeffs <- NULL
    nm <- namer <- NULL
    r <- ans$r
    s <- ans$s
    b <- ans$b


    period <- ans$period
    for (i in 1:length(ans$modes)) {
        phi <- ans$modes[[i]]$phi
        theta <- ans$modes[[i]]$theta
        phiseas <- ans$modes[[i]]$phiseas
        thetaseas <- ans$modes[[i]]$thetaseas
        phip <- ans$modes[[i]]$phip
        thetap <- ans$modes[[i]]$thetap
        phiseasp <- ans$modes[[i]]$phiseasp
        thetaseasp <- ans$modes[[i]]$thetaseasp

        dfrac <- ans$modes[[i]]$dfrac
        dfs <- ans$modes[[i]]$dfs
        H <- ans$modes[[i]]$H
        Hs <- ans$modes[[i]]$Hs
        alpha <- ans$modes[[i]]$alpha
        alphas <- ans$modes[[i]]$alphas

        delta <- ans$modes[[i]]$delta
        omega <- ans$modes[[i]]$omega
        logl <- ans$modes[[i]]$loglik
        sigma2 <- ans$modes[[i]]$sigma2
        muHat <- ans$modes[[i]]$muHat

        se <- ans$modes[[i]]$se
        if (i == 1) {
            nm <- getnames(ans)
            nm <- c(nm, paste("logl", sep = ""), paste("sigma^2", sep = ""))
            counter <- 0
            if (length(phi) > 1) {
                nm <- c(nm, paste("phi_p(", 1:length(phi), ")", sep = ""))
                counter <- counter + length(phi)
            }
            if (length(theta) > 1) {
                nm <- c(nm, paste("theta_p(", 1:length(theta), ")", sep = ""))
                counter <- counter + length(theta)
            }
            if (length(phiseas) > 1) {
                nm <- c(nm, paste("phi_p.", period, "(", 1:length(phiseas), ")", sep = ""))
                counter <- counter + length(phiseas)
            }
            if (length(thetaseas) > 1) {
                nm <- c(nm, paste("theta_p.", period, "(", 1:length(thetaseas), ")", sep = ""))
                counter <- counter + length(thetaseas)
            }
        }


        coeff <- c(phi, theta, phiseas, thetaseas, dfrac, H, alpha, dfs, Hs, alphas,
                   omega, delta, muHat, logl, sigma2)

        if (length(phi) > 1)
            coeff <- c(coeff, phip)
        if (length(theta) > 1)
            coeff <- c(coeff, thetap)
        if (length(phiseas) > 1)
            coeff <- c(coeff, phiseasp)
        if (length(thetaseas) > 1)
            coeff <- c(coeff, thetaseasp)
        if (x$getHess) {
            ses <- c(se, -Inf, -Inf, if (counter > 0) rep(-Inf, counter) else numeric(0))

            coeff <- rbind(coeff, ses)
        } else coeff <- matrix(coeff, nrow = 1)

        coeff <- signif(coeff, digits = digits)

        coeffs <- rbind(coeffs, coeff)
        star <- if (i %in% ans$index)
            "*" else ""
        nam <- if (x$getHess)
            c(paste("Coef.", i, star, ":", sep = ""), paste("SE.", i, star, ":", sep = "")) else paste("Coef.", i, star, ":", sep = "")
        namer <- c(namer, nam)
    }
    cat("\nCall:", deparse(x$call, width.cutoff = 75L), "", sep = "\n")
    cat("Coefficients for fits:\n", sep = "")

    colnames(coeffs) <- nm
    rownames(coeffs) <- namer
    coeffs <- t(coeffs)
    indss <- which(!is.finite(coeffs) & !is.na(coeffs), arr.ind = TRUE)
    inds2 <- which(is.na(coeffs), arr.ind = TRUE)
    coeffs2 <- coeffs
    pos <- which(is.finite(coeffs) & !is.na(coeffs) & (coeffs >= 0), arr.ind = TRUE)
    coeffs2[indss] <- ""
    coeffs2[inds2] <- "NA"
    coeffs2[pos] <- paste(" ", coeffs[pos], sep = "")
    print.default(coeffs2, print.gap = 2, quote = FALSE)
    cat("Starred fits are close to invertibility/stationarity boundaries\n")
    if (ans$numcut > 0)
        cat("Wall detection on: ", ans$numcut, " mode(s) were believed to be spurious and were cut\n",
            sep = "")
    if (any(is.na(coeffs)))
        cat("NAs come from singular Hessians\n")
    if(x$differencing)
      cat("Note that the fitted mean/zbar is of the differenced series.")
    invisible(x)
}



#' Extract Model Fitted Values
#'
#' Extract fitted values from an \code{arfima} object.
#'
#'
#' @param object A \code{arfima} object.
#' @param \dots Optional parameters.  Currently not used.
#' @return A list of vectors of fitted values, one for each mode.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima}}, \code{\link{resid.arfima}}
#' @references Veenstra, J.Q. Persistence and Antipersistence:  Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#' 
#'\donttest{
#' set.seed(8564)
#' sim <- arfima.sim(1000, model = list(phi = c(0.2, 0.1), dfrac = 0.4, theta = 0.9))
#' fit <- arfima(sim, order = c(2, 0, 1), back=TRUE)
#'
#' fit
#'
#' resid <- resid(fit)
#' par(mfrow = c(1, 3))
#' fitted <- fitted(fit)
#' plot(fitted[[1]], resid[[1]])
#' plot(fitted[[2]], resid[[2]])
#' plot(fitted[[3]], resid[[3]])
#' par(mfrow = c(1, 1))
#'}
fitted.arfima <- function(object, ...) {

    if (!object$weeded) {
        warning("The object has not been weeded:  there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
    }
    ans <- object


    m <- length(ans$modes)

    ret <- vector("list", m)

    for (i in 1:m) {
        ret[[i]] <- ans$modes[[i]]$fitted
    }

    names(ret) <- paste("Mode", 1:m, sep = "")
    ret
}

#' Extract the Residuals of a Fitted Object
#'
#' Extracts the residuals or regression residuals from a fitted \code{arfima}
#' object
#'
#'
#' @aliases resid.arfima
#' @param object A fitted \code{arfima} object
#' @param reg Whether to extract the regression residuals instead.  If
#' \code{TRUE}, throws an error if no regression was done.
#' @param \dots Optional parameters.  Currently not used.
#' @return A list of vectors of residuals, one for each mode.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima}}, \code{\link{fitted.arfima}}
#' @references Veenstra, J.Q. Persistence and Antipersistence:  Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#'
#'\donttest{
#' set.seed(8564)
#' sim <- arfima.sim(1000, model = list(phi = c(0.2, 0.1), dfrac = 0.4, theta = 0.9))
#' fit <- arfima(sim, order = c(2, 0, 1), back=TRUE)
#'
#' fit
#'
#' resid <- resid(fit)
#' par(mfrow = c(1, 3))
#' plot(resid[[1]])
#' plot(resid[[2]])
#' plot(resid[[3]])
#' fitted <- fitted(fit)
#' plot(fitted[[1]], resid[[1]])
#' plot(fitted[[2]], resid[[2]])
#' plot(fitted[[3]], resid[[3]])
#' par(mfrow = c(1, 1))
#'}
#'
residuals.arfima <- resid.arfima <- function(object, reg = FALSE, ...) {

    if (!object$weeded) {
        warning("The object has not been weeded:  there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
    }
    ans <- object


    xreg <- object$xreg
    if (reg && is.null(xreg))
        stop("Regression residuals requested, but no regression performed.")


    m <- length(ans$modes)

    ret <- vector("list", m)

    for (i in 1:m) {
        cmode <- ans$modes[[i]]
        if (!reg)
            ret[[i]] <- cmode$residuals else ret[[i]] <- cmode$regResiduals
    }

    names(ret) <- paste("Mode", 1:m, sep = "")
    ret
}




#' Extensive Summary of an Object
#'
#' Provides a very comprehensive summary of a fitted \code{arfima} object.
#' Includes correlation and covariance matrices (observed and expected), the
#' Fisher Information matrix of those parameters for which it is defined, and
#' more, for each mode.
#'
#'
#' @param object A fitted \code{arfima} object
#' @param digits The number of digits to print
#' @param \dots Optional arguments, currently not used.
#' @return A list of lists (one for each mode) of all relevant information
#' about the fit that can be passed to \code{print.summary.arfima}.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima}}, \code{\link{iARFIMA}},
#' \code{\link{vcov.arfima}}
#' @references Veenstra, J.Q. Persistence and Antipersistence:  Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#'
#' data(tmpyr)
#'
#' fit <- arfima(tmpyr, order = c(1, 0, 1), back=TRUE)
#' fit
#'
#' summary(fit)
#'
summary.arfima <- function(object, digits = max(4, getOption("digits") - 3), ...) {
    if (!object$weeded) {
        warning("The object has not been weeded:  there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
    }
    ans <- object
    cov <- TRUE
    cor <- TRUE
    FInfo <- TRUE
    cf <- coef(ans)
    vcovs <- vcov(ans, type = "b", summ = TRUE)
    m <- length(ans$modes)
    dmean <- ans$dmean
    cff <- vector("list", m)
    if (FInfo)
        finfo <- vector("list", m) else finfo <- NULL
    if (cov && cor)
        corrs <- vcov(ans, cor = cor, type = "b", summ = TRUE) else if (cov)
        corrs <- vcovs else corrs <- NULL
    logl <- logLik(ans)
    aics <- AIC(logl)
    bics <- AIC(logl, k=log(length(object$z)))
    ans1 <- list()
    sigmas <- rep(0, m)
    tseflag <- TRUE
    for (i in 1:m) {
        ses <- sqrt(diag(vcovs[[i]]$observed))
        ses[ans$indfixx] <- NA
        if ((!is.logical(dmean)) || (!dmean))
            ses <- c(ses, NA)
        if (!vcovs$warnH && !vcovs$warnX)
            tses <- c(sqrt(diag(vcovs[[i]]$expected)), NA) else {
            tses <- NULL
            tseflag <- FALSE
        }
        if (!is.null(finfo))
            finfo[[i]] <- solve(vcovs[[i]]$expected * ans$n)
        zval <- cf[i, ]/ses
        ps <- 2 * pnorm(-abs(zval))
        cff[[i]] <- cbind(Estimate = cf[i, ], `Std. Error` = ses, `Th. Std. Err.` = tses,
            `z-value` = zval, `Pr(>|z|)` = ps)
        sigmas[i] <- ans$modes[[i]]$sigma2
    }
    ans1$tseflag <- tseflag
    if (cor)
        ans1$corrs <- corrs else ans1$covs <- corrs
    ans1$sigma2 <- sigmas
    ans1$logl <- logl
    ans1$coef <- cff
    ans1$aics <- aics
    ans1$bics <- bics
    ans1$finfo <- finfo
    ans1$m <- m
    ans1$call <- ans$call
    class(ans1) <- "summary.arfima"
    ans1
}



#' Prints the output of a call to \code{summary} on an \code{arfima} object
#'
#' Prints the output of a call to \code{summary} on an \code{arfima} object
#'
#'
#' @param x A \code{summary.arfima} object
#' @param digits The number of digits to print
#' @param signif.stars Whether to print stars on significant output
#' @param \dots Currently not used
#' @return Returns the object \code{x} invisibly
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{arfima}}, \code{\link{print.arfima}},
#' \code{\link{summary.arfima}}, \code{\link{print}}
#' @references Veenstra, J.Q. Persistence and Antipersistence:  Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#'
#' set.seed(54678)
#' sim <- arfima.sim(1000, model = list(phi = 0.9, H = 0.3))
#' fit <- arfima(sim, order = c(1, 0, 0), lmodel = "g", back=TRUE)
#' summary(fit)
#'
print.summary.arfima <- function(x, digits = max(6, getOption("digits") - 3), signif.stars = getOption("show.signif.stars"),
    ...) {
    cat("\nCall:\n ", deparse(x$call, width.cutoff = 75L), "", sep = "\n")

    for (i in 1:x$m) {
        cat("\nMode", i, "Coefficients:\n")
        printCoefmat(x$coef[[i]], digits = digits, signif.stars = signif.stars, ...)
        cat(paste("sigma^2 estimated as ", format(x$sigma2[i], digits = digits), "; Log-likelihood = ",
            format(x$logl[i], digits = digits), "; AIC = ", format(x$aics[i], digits = digits),
            "; BIC = ", format(x$bics[i], digits = digits), "\n", sep = ""))
        if (!is.null(x$corrs)) {
            cat("\nNumerical Correlations of Coefficients:\n")
            correl <- format(round(x$corrs[[i]]$observed, 2), nsmall = 2, digits = digits)

            print(correl, quote = FALSE)
            cat("\nTheoretical Correlations of Coefficients:\n")
            correl <- format(round(x$corrs[[i]]$expected, 2), nsmall = 2, digits = digits)
            print(correl, quote = FALSE)
        } else if (!is.null(x$covs)) {
            cat("\nNumerical Covariances of Coefficients:\n")
            correl <- format(round(x$covs[[i]]$observed, 2), nsmall = 2, digits = digits)

            print(correl, quote = FALSE)
            cat("\nTheoretical Covariances of Coefficients:\n")
            correl <- format(round(x$covs[[i]]$expected, 2), nsmall = 2, digits = digits)
            print(correl, quote = FALSE)
        }
        if (!is.null(x$finfo)) {
            ff <- format(round(x$finfo[[i]], 2), nsmall = 2, digits = digits)
            cat("\nExpected Fisher Information Matrix of Coefficients:\n")
            print(ff, quote = FALSE)
        }
    }
    if (!x$tseflag)
        cat("\nTheoretical standard errors not reported since there is either FGN or (dynamic) regression in the model\n")
    invisible(x)
}

InfoMatrix <- function(object, digits = max(6, getOption("digits") - 3), tapprox = FALSE,
    ...) {
    if (!inherits(object, "arfima"))
        stop("only arfima objects can be passed to this function")
    if (!object$weeded) {
        warning("The object has not been weeded:  there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
    }
    ans <- object

    m <- length(ans$modes)

    ret <- vector("list", m)
    ifdf <- length(ans$modes[[1]]$dfrac) > 0
    ifdfs <- length(ans$modes[[1]]$dfs) > 0
    period <- ans$period
    for (i in 1:m) {
        cmode <- ans$modes[[i]]
        if (length(cmode$H) > 0 || length(cmode$Hs) > 0 || length(cmode$alpha) > 0 || length(cmode$alphas) >
            0)
            warning("Please note FGN and HD do not have known information matrix components: \n\n\t\treporting matrix without those components")
        ret[[i]] <- iARFIMA(phi = cmode$phi, theta = cmode$theta, phiseas = cmode$phiseas,
            thetaseas = cmode$thetaseas, dfrac = ifdf, dfs = ifdfs, period = period, exact = !tapprox)
    }
    names(ret) <- paste("Mode", 1:m)
    ret
}





#' Extract Log-Likelihood Values
#'
#' Extracts log-likelihood values from a \code{arfima} fit.
#'
#' Uses the function \code{DLLoglikelihood} from the package
#' \code{\link{ltsa}}.  The log-likelihoods returned are exact up to an
#' additive constant.
#'
#' @param object A fitted \code{arfima} object
#' @param \dots Optional arguments not currently used.
#' @return A vector of log-likelihoods, one for each mode, is returned, along
#' with the degrees of freedom.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{AIC.arfima}}
#' @references Veenstra, J.Q. Persistence and Antipersistence:  Theory and
#' Software (PhD Thesis)
#' @keywords ts
logLik.arfima <- function(object, ...) {
    if (!object$weeded) {
        warning("The object has not been weeded:  there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
    }
    ans <- object
    logl <- rep(0, length(ans$modes))
    for (i in 1:length(ans$modes)) logl[i] <- ans$modes[[i]]$loglik

    len <- ncol(coef(ans))

    if (is.numeric(ans$dmean) && !ans$itmean)
        len <- len - 1

    # sigma^2 is the +1 in the below.
    attr(logl, "df") <- len + 1
    attr(logl, "nobs") <- attr(logl, "nall") <- ans$n
    class(logl) <- "logLik"
    logl
}

#' Extract Model Coefficients
#'
#' Extracts the coefficients from a \code{arfima} fit.
#'
#'
#' @param object A fitted \code{arfima} object.
#' @param tpacf If \code{TRUE}, the (ARMA) coefficients are in the transformed
#' PACF space.
#' @param digits The number of digits to print
#' @param \dots Other optional arguments.  Currently not used.
#' @return A matrix of coefficients.  The rows are for the modes, and the
#' columns are for the model variables.
#' @author JQ (Justin) Veenstra
#' @keywords ts
#' @examples
#'
#'\donttest{
#' set.seed(8564)
#' sim <- arfima.sim(1000, model = list(phi = c(0.2, 0.1), dfrac = 0.4, theta = 0.9))
#' fit <- arfima(sim, order = c(2, 0, 1), back=TRUE)
#'
#' fit
#' coef(fit)
#'}
coef.arfima <- function(object, tpacf = FALSE, digits = max(4, getOption("digits") - 3),
    ...) {
    if (!object$weeded) {
        warning("The object has not been weeded:  there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
    }
    ans <- object

    dmean <- ans$dmean
    itmean <- ans$itmean
    m <- length(ans$modes)
    s <- object$s
    r <- object$r
    b <- object$b
    coeffs <- NULL
    period <- ans$period
    for (i in 1:m) {
        phi <- ans$modes[[i]]$phi
        theta <- ans$modes[[i]]$theta
        phiseas <- ans$modes[[i]]$phiseas
        thetaseas <- ans$modes[[i]]$thetaseas
        phip <- ans$modes[[i]]$phip
        thetap <- ans$modes[[i]]$thetap
        phiseasp <- ans$modes[[i]]$phiseasp
        thetaseasp <- ans$modes[[i]]$thetaseasp

        dfrac <- ans$modes[[i]]$dfrac
        dfs <- ans$modes[[i]]$dfs
        H <- ans$modes[[i]]$H
        Hs <- ans$modes[[i]]$Hs
        alpha <- ans$modes[[i]]$alpha
        alphas <- ans$modes[[i]]$alphas

        delta <- ans$modes[[i]]$delta
        omega <- ans$modes[[i]]$omega

        muHat <- ans$modes[[i]]$muHat

        xreglist <- ans$modes[[i]]$xreglist
        if (i == 1) {
            nm <- NULL

            if (tpacf) {
                if (length(phip) > 0)
                  nm <- c(nm, paste("phi_p(", 1:length(phi), ")", sep = ""))

                if (length(thetap) > 0)
                  nm <- c(nm, paste("theta_p(", 1:length(theta), ")", sep = ""))

                if (length(phiseasp) > 0)
                  nm <- c(nm, paste("phi_p.", period, "(", 1:length(phiseas), ")", sep = ""))

                if (length(thetaseasp) > 0)
                  nm <- c(nm, paste("theta_p.", period, "(", 1:length(thetaseas), ")", sep = ""))
            } else {
                if (length(phi) > 0)
                  nm <- c(nm, paste("phi(", 1:length(phi), ")", sep = ""))

                if (length(theta) > 0)
                  nm <- c(nm, paste("theta(", 1:length(theta), ")", sep = ""))

                if (length(phiseas) > 0)
                  nm <- c(nm, paste("phi.", period, "(", 1:length(phiseas), ")", sep = ""))

                if (length(thetaseas) > 0)
                  nm <- c(nm, paste("theta.", period, "(", 1:length(thetaseas), ")", sep = ""))
            }

            if (length(dfrac) > 0)
                nm <- c(nm, paste("d.f", sep = ""))
            if (length(H) > 0)
                nm <- c(nm, paste("H", sep = ""))
            if (length(alpha) > 0)
                nm <- c(nm, paste("alpha", sep = ""))

            if (length(dfs) > 0)
                nm <- c(nm, paste("d.f.", period, sep = ""))
            if (length(Hs) > 0)
                nm <- c(nm, paste("H.", period, sep = ""))
            if (length(alphas) > 0)
                nm <- c(nm, paste("alpha.", period, sep = ""))


            if(ans$strReg){
              nm <- c(nm, names(xreglist$omega))
            }
            else {
              if (length(omega) > 0) {
                  for (j in 1:length(s)) nm <- c(nm, names(xreglist$omega[[j]]))
              }
              if (length(delta) > 0) {
                  for (j in 1:length(r)) nm <- c(nm, names(xreglist$delta[[j]]))
              }

            }
            if (itmean) {
              nm <- c(nm, "It. Fit. mean")
            } else  {
              if (is.logical(dmean) && !dmean)
                nm <- c(nm, "zbar") else if (is.logical(dmean) && dmean) {
                  if(ans$strReg && !ans$intindex)
                    nm <- c(nm, "Intercept")
                  else if(!ans$intindex)
                    nm <- c(nm, "Fitted mean")
                  else
                    nm <- c(nm, ans$intname)
                }
              else nm <- c(nm, "Set mean")
            }

        }

        if (tpacf)
            coeff <- c(phip, thetap, phiseasp, thetaseasp, dfrac, H, alpha, dfs, Hs, alphas,
                       omega, delta, muHat)
        else
            coeff <- c(phi, theta, phiseas, thetaseas, dfrac, H, alpha, dfs, Hs, alphas,
                       omega, delta, muHat)

        coeffs <- rbind(coeffs, coeff)
    }
    colnames(coeffs) <- nm
    rownames(coeffs) <- paste("Mode", 1:m)
    coeffs
}

#' Information criteria for \code{arfima} objects
#'
#' Computes information criteria for \code{arfima} objects.  See
#' \code{\link{AIC}} for more details.
#'
#'
#' @aliases AIC.arfima BIC BIC.arfima
#' @param object An object of class "arfima".  Note these functions can only be
#' called on one object at a time because of possible multimodality.
#' @param \dots Other models fit to data for which to extract the AIC/BIC.  Not
#' recommended, as an \code{arfima} object can be multimodal.
#' @param k The penalty term to be used.  See \code{\link{AIC}}.
#' @return The information criteria for each mode in a vector.
#' @author JQ (Justin) Veenstra
#' @keywords ts
#' @examples
#'
#'\donttest{
#' set.seed(34577)
#' sim <- arfima.sim(500, model = list(theta = 0.9, phi = 0.5, dfrac = 0.4))
#' fit1 <- arfima(sim, order = c(1, 0, 1), cpus = 2, back=TRUE)
#' fit2 <- arfima(sim, order = c(1, 0, 1), cpus = 2, lmodel = "g", back=TRUE)
#' fit3 <- arfima(sim, order = c(1, 0, 1), cpus = 2, lmodel = "h", back=TRUE)
#'
#' AIC(fit1)
#' AIC(fit2)
#' AIC(fit3)
#'}
AIC.arfima <- function(object, ..., k = 2) {
    if (!object$weeded) {
        warning("The object has not been weeded:  there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
    }
    ans <- object

    logl <- logLik(ans, ...)
    AIC(logl, k = k)
}

BIC.arfima <- function(object, ...) {
    if (!object$weeded) {
        warning("The object has not been weeded:  there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
    }
    ans <- object

    logl <- logLik(ans, ...)
    AIC(logl, k = log(length(object$z)))
}




#' Extracts the Variance-Covariance Matrix
#'
#' Extracts the variance-covariance matrices (one or two for each mode) from a
#' fitted \code{arfima} object.
#'
#'
#' @param object A fitted \code{arfima} object
#' @param type Which type of covariance matrix to return: "o" is the observed
#' matrix (from solving the Hessian), "e" is the expected matrix (from solving
#' the information matrix), and "b" is both.
#' @param cor Whether or not the correlation matrix should be returned instead.
#' @param digits The number of digits to print.
#' @param tapprox Whether or not to use an approximation to find the expected
#' matrix.  Highly recommended to be FALSE, as it takes much longer, and is an
#' approximation.
#' @param summ Whether the call is from the \code{summary.arfima} function.
#' Should not be used except internally.
#' @param \dots Optional arguments, currently not used.
#' @return A list of lists (one for each mode) with components \code{observed}
#' and/or \code{expected}.
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{summary.arfima}}, \code{\link{arfima}}
#' @references Veenstra, J.Q. Persistence and Antipersistence:  Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#'\donttest{
#' set.seed(1234)
#' sim <- arfima.sim(1000, model = list(dfrac = 0.4, phi = .8, theta = -0.5))
#' fit1 <- arfima(sim, order = c(1, 0, 1), back=TRUE)
#' fit2 <- arfima(sim, order = c(1, 0, 1), lmodel = "g", back=TRUE)
#' fit3 <- arfima(sim, order = c(1, 0, 1), lmodel = "h", back=TRUE)
#' fit1
#' fit2
#' fit3
#' vcov(fit1)
#' vcov(fit2)
#' vcov(fit2)
#'}
vcov.arfima <- function(object, type = c("b", "o", "e"), cor = FALSE, digits = max(4, getOption("digits") -
    3), tapprox = FALSE, summ = FALSE, ...) {
    if (!object$weeded) {
        warning("The object has not been weeded:  there may be duplicate modes, as well as excessive output. \n For more sensible output, please call weed() on the object")
    }
    ans <- object
    
    if (is.vector(type))
      type <- type[1]
    
    if (type %in% c("b", "o", "e") && !summ)
        cat("Returning both observed and theoretical covariance matrices: use type = \"o\" for just\n\t observed and \"e\" for just expected (from the expected Fisher information)\n")

    xreg <- ans$xreg
    dmean <- ans$dmean
    itmean <- ans$itmean

    m <- length(ans$modes)

    ret <- vector("list", m)

    nm <- getnames(ans)

    if ((is.logical(ans$dmean) && !ans$dmean) || is.double(ans$dmean))
        nm <- nm[-length(nm)]
    period <- ans$period
    warnH <- warnX <- FALSE
    if (type %in% c("b", "o")) {
        for (i in 1:m) {
            temp <- solve(-ans$modes[[i]]$hess)
            if(all(!ans$indfixx)) {
              ret[[i]]$observed <- temp
            }
            else {
              num <- length(ans$indfixx)
              ret[[i]]$observed <- matrix(0, nrow = num, ncol = num)
              ret[[i]]$observed[!ans$indfixx, !ans$indfixx] <- temp
            }

            if (cor) {
                temp <- ret[[i]]$observed

                for (j in 1:nrow(temp)) {
                  for (k in 1:ncol(temp)) {
                    temp[j, k] <- ret[[i]]$observed[j, k]/sqrt(ret[[i]]$observed[j, j] *
                      ret[[i]]$observed[k, k])
                  }
                }
                ret[[i]]$observed <- temp
                ret[[i]]$observed[ans$indfixx, ans$indfixx] <- 0
            }
            ret[[i]]$observed <- signif(ret[[i]]$observed, digits = digits)
            rownames(ret[[i]]$observed) <- colnames(ret[[i]]$observed) <- nm
        }
    }
    if (type %in% c("b", "e")) {
        if ((length(ans$modes[[1]]$H) > 0 || length(ans$modes[[1]]$Hs) > 0 || length(ans$modes[[1]]$alpha) >
            0 || length(ans$modes[[1]]$alphas) > 0)) {
            if (!summ)
                warning("The expected Fisher Information matrix has no known closed form for FGN or HD processes.\n\n\t\t\t\tReporting this matrix for the non-FGN/non-HD processes: however, results may be skewed")
            warnH <- TRUE
        }
        if (!is.null(xreg)) {
            if (!summ)
                warning("The expected Fisher Information matrix does not include regression parameters.\n\n\t\t\t\tReporting this matrix for the ARFIMA processes: however, results may be skewed")
            warnX <- TRUE
        }
        ifdf <- (length(ans$modes[[1]]$dfrac) > 0)
        ifdfs <- (length(ans$modes[[1]]$dfs) > 0)
        for (i in 1:m) {
            cmode <- ans$modes[[i]]

            ret[[i]]$expected <- solve(iARFIMA(phi = cmode$phi, theta = cmode$theta, phiseas = cmode$phiseas,
                thetaseas = cmode$thetaseas, dfrac = ifdf, dfs = ifdfs, period = period,
                exact = !tapprox))/ans$n
            if (cor) {
                temp <- ret[[i]]$expected
                for (j in 1:nrow(temp)) {
                  for (k in 1:ncol(temp)) {
                    temp[j, k] <- ret[[i]]$expected[j, k]/sqrt(ret[[i]]$expected[j, j] *
                      ret[[i]]$expected[k, k])
                  }
                }
                ret[[i]]$expected <- temp
            }
            ret[[i]]$expected <- signif(ret[[i]]$expected, digits = digits)
        }
    }
    names(ret) <- paste("Mode ", 1:m, sep = "")

    if (summ) {
        ret$warnH <- warnH
        ret$warnX <- warnX
    }

    ret
}


getnames <- function(ans) {

    if (!inherits(ans, "arfima"))
        stop("the getnames function only works for arfima objects")

    nm <- colnames(coef(weed(ans)))

    nm
}
JQVeenstra/arfima documentation built on Jan. 31, 2024, 12:11 p.m.