R/plot.predARFIMA.R

In arfima: Fractional ARIMA (and Other Long Memory) Time Series Modeling

#' Plots the original time series, the predictions, and the prediction
#' intervals for a \code{predarfima} object.
#' 
#' This function takes a \code{predarfima} object generated by
#' \code{predict.arfima} and plots all of the information contained in it.  The
#' colour code is as follows:
#'
#' grey: exact prediction
#' red: exact prediction intervals (PIs)
#' orange: limiting PIs
#'
#' See \code{\link{predict.arfima}}.
#'
#'
#' @param x A \code{predarfima} object
#' @param xlab Optional
#' @param ylab Optional
#' @param main Optional
#' @param ylim Optional
#' @param numback The number of last values of the original series to plot
#' defined by the user. The default is five
#' @param xlim Optional
#' @param \dots Currently not used
#' @return None. Generates a plot
#' @author JQ (Justin) Veenstra
#' @seealso \code{\link{predict.arfima}}, \code{\link{print.predarfima}}
#' @references Veenstra, J.Q. Persistence and Antipersistence:  Theory and
#' Software (PhD Thesis)
#' @keywords ts
#' @examples
#'\donttest{
#' 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)
#' #Let's look at more context
#' plot(pred, numback = 50)
#'}
"plot.predarfima" <- function(x, xlab = NULL, ylab = NULL, main = NULL, ylim = NULL, numback = 5,
    xlim = NULL, ...) {

    op <- par(no.readonly = TRUE)
    predint <- x$predint
    z <- x$z
    n <- length(z)
    if (length(xlab) == 0)
        xlab <- "Time"
    if (length(ylab) == 0)
        ylab <- "Mode "
    if (length(main) == 0)
        main <- paste("Time Series and Predictions of ", x$name, sep = "")
    m <- x$m
    n.ahead <- length(as.vector(x[[1]]$Forecast))
    seed <- x$seed
    limiting <- x$limiting


    numSD <- qnorm(1 - (1 - predint)/2, sd = 1)

    if (numback > 0) {
        ys <- z[(n - numback + 1):n]
        xs <- (n - numback + 1):n
    } else {
        cat("setting numback to 1\n")
        ys <- z[n]
        xs <- n
    }
    islim <- !is.null(x[[1]]$limitSD)

    minn <- if (numback > 0)
        min(ys) else Inf
    maxx <- if (numback > 0)
        max(ys) else -Inf
    for (i in 1:m) {
        if(length(x[[i]]$exactSD)>0) {
          exactPIU <- x[[i]]$Forecast + numSD * x[[i]]$exactSD
          exactPIL <- x[[i]]$Forecast - numSD * x[[i]]$exactSD
          exact <- TRUE
        }
        else {
          exactPIU <- exactPIL <- NULL
          exact <- FALSE
        }

        limitPIU <- if (islim)
            x[[i]]$Forecast + numSD * x[[i]]$limitSD else NULL
        limitPIL <- if (islim)
            x[[i]]$Forecast - numSD * x[[i]]$limitSD else NULL
        maxx <- max(maxx, exactPIU, limitPIU, x[[i]]$uppernp)
        minn <- min(minn, exactPIL, limitPIL, x[[i]]$lowernp)
        x[[i]]$exactPIU <- exactPIU
        x[[i]]$exactPIL <- exactPIL
        x[[i]]$limitPIU <- limitPIU
        x[[i]]$limitPIL <- limitPIL
    }
    ranger <- maxx - minn
    if (length(ylim) == 0)
        ylim <- c(minn - ranger/20, maxx + ranger/20)

    leg <- c("Exact prediction")
    ltt <- c(1)
    coll <- c("gray")

    if(exact) {
      ltt <- c(ltt, 2)
      coll <- c(coll, "red")
      leg <- c(leg, paste("Exact", predint * 100, "% PI"))
    }

    if (islim) {
        ltt <- c(ltt, 2)
        coll <- c(coll, "orange")
        leg <- c(leg, paste("Limiting", predint * 100, "% PI"))
    }


    nf <- layout(matrix(c(1:(m + 1)), m + 1, 1), widths = rep(1, m + 1), heights = c(0.5,
        rep(1, m)), FALSE)
    par(mar = c(0, 0, 0, 0))
    plot.default(0, 0, axes = FALSE, type = "n")
    legend("bottom", legend = leg, col = coll, lty = ltt, bty = "n", ncol = ceiling(length(leg)/2))
    text(0, 0.5, main, cex = 2)
    for (i in 1:m) {
        par(mar = c(3, 4.2, 0, 0))
        xx <- x[[i]]
        xx$z <- ys
        xx$xs <- c(xs, xs[length(xs)] + 1:(n.ahead + 1))
        ylab1 <- paste(ylab, i, sep = "")
        plotpredARFIMA(xx, xlab = xlab, ylab = ylab1, ylim = ylim)
    }
    par(op)
}


"plotpredARFIMA" <- function(xy, xlab = NULL, ylab = NULL, ylim = NULL, ...) {
    nn <- length(xy$z)
    plot(x = xy$xs[1:nn], y = xy$z, ylab = ylab, xlab = xlab, xlim = c(xy$xs[1], xy$xs[length(xy$xs)]),
        ylim = ylim, type = "l", cex.lab = 1.4)
    zlast <- xy$z[nn]
    lines(xy$xs[nn:(length(xy$xs) - 1)], c(zlast, xy$Forecast), col = "gray")
    lines(xy$xs[(nn + 1):(length(xy$xs) - 1)], xy$exactPIU, col = "red", lty = 2)
    lines(xy$xs[(nn + 1):(length(xy$xs) - 1)], xy$exactPIL, col = "red", lty = 2)

    if (!is.null(xy$limitPIU)) {
        lines(xy$xs[(nn + 1):(length(xy$xs) - 1)], xy$limitPIU, col = "orange", lty = 2)
        lines(xy$xs[(nn + 1):(length(xy$xs) - 1)], xy$limitPIL, col = "orange", lty = 2)
    }

}