R/ggplot.R

Defines functions gghistogram geom_forecast blendHex forecast2plotdf fortify.ts autoplot.msts autoplot.mts autoplot.ts autolayer.mforecast autolayer.forecast autolayer.ts autolayer.msts autolayer.mts autoplot.seas autoplot.StructTS autoplot.stl autoplot.splineforecast ggseasonplot ggsubseriesplot ggmonthplot gglagchull gglagplot ggtsdisplay autoplot.mforecast autoplot.forecast autoplot.bats autoplot.tbats autoplot.ets autoplot.decomposed.ts autoplot.ar autoplot.Arima ggtaperedpacf ggtaperedacf autoplot.mpacf ggCcf ggPacf ggAcf autoplot.acf ggtsbreaks ggAddExtras autolayer

Documented in autolayer autolayer.forecast autolayer.mforecast autolayer.msts autolayer.mts autolayer.ts autoplot.acf autoplot.ar autoplot.Arima autoplot.bats autoplot.decomposed.ts autoplot.ets autoplot.forecast autoplot.mforecast autoplot.mpacf autoplot.msts autoplot.mts autoplot.seas autoplot.splineforecast autoplot.stl autoplot.StructTS autoplot.tbats autoplot.ts fortify.ts geom_forecast ggAcf ggCcf gghistogram gglagchull gglagplot ggmonthplot ggPacf ggseasonplot ggsubseriesplot ggtaperedacf ggtaperedpacf ggtsdisplay

globalVariables(".data")

#' @inherit ggplot2::autolayer
#' @export
autolayer <- function(object, ...){
  UseMethod("autolayer")
}

#' @importFrom ggplot2 autoplot
#' @export
ggplot2::autoplot

ggAddExtras <- function(xlab=NA, ylab=NA, main=NA) {
  dots <- eval.parent(quote(list(...)))
  extras <- list()
  if ("xlab" %in% names(dots) || is.null(xlab) || any(!is.na(xlab))) {
    if ("xlab" %in% names(dots)) {
      extras[[length(extras) + 1]] <- ggplot2::xlab(dots$xlab)
    }
    else {
      extras[[length(extras) + 1]] <- ggplot2::xlab(paste0(xlab[!is.na(xlab)], collapse = " "))
    }
  }
  if ("ylab" %in% names(dots) || is.null(ylab) || any(!is.na(ylab))) {
    if ("ylab" %in% names(dots)) {
      extras[[length(extras) + 1]] <- ggplot2::ylab(dots$ylab)
    }
    else {
      extras[[length(extras) + 1]] <- ggplot2::ylab(paste0(ylab[!is.na(ylab)], collapse = " "))
    }
  }
  if ("main" %in% names(dots) || is.null(main) || any(!is.na(main))) {
    if ("main" %in% names(dots)) {
      extras[[length(extras) + 1]] <- ggplot2::ggtitle(dots$main)
    }
    else {
      extras[[length(extras) + 1]] <- ggplot2::ggtitle(paste0(main[!is.na(main)], collapse = " "))
    }
  }
  if ("xlim" %in% names(dots)) {
    extras[[length(extras) + 1]] <- ggplot2::xlim(dots$xlim)
  }
  if ("ylim" %in% names(dots)) {
    extras[[length(extras) + 1]] <- ggplot2::ylim(dots$ylim)
  }
  return(extras)
}

ggtsbreaks <- function(x) {
  # Make x axis contain only whole numbers (e.g., years)
  return(unique(round(pretty(floor(x[1]):ceiling(x[2])))))
}


#' ggplot (Partial) Autocorrelation and Cross-Correlation Function Estimation
#' and Plotting
#'
#' Produces a ggplot object of their equivalent Acf, Pacf, Ccf, taperedacf and
#' taperedpacf functions.
#'
#' If \code{autoplot} is given an \code{acf} or \code{mpacf} object, then an
#' appropriate ggplot object will be created.
#'
#' ggtaperedpacf
#' @param object Object of class \dQuote{\code{acf}}.
#' @param x a univariate or multivariate (not Ccf) numeric time series object
#' or a numeric vector or matrix.
#' @param y a univariate numeric time series object or a numeric vector.
#' @param ci coverage probability for confidence interval. Plotting of the
#' confidence interval is suppressed if ci is zero or negative.
#' @param lag.max maximum lag at which to calculate the acf.
#' @param type character string giving the type of acf to be computed. Allowed
#' values are "\code{correlation}" (the default), \dQuote{\code{covariance}} or
#' \dQuote{\code{partial}}.
#' @param plot logical. If \code{TRUE} (the default) the resulting ACF, PACF or
#' CCF is plotted.
#' @param na.action function to handle missing values. Default is
#' \code{\link[stats]{na.contiguous}}.  Useful alternatives are
#' \code{\link[stats]{na.pass}} and \code{\link{na.interp}}.
#' @param demean Should covariances be about the sample means?
#' @param calc.ci If \code{TRUE}, confidence intervals for the ACF/PACF
#' estimates are calculated.
#' @param level Percentage level used for the confidence intervals.
#' @param nsim The number of bootstrap samples used in estimating the
#' confidence intervals.
#' @param ... Other plotting parameters to affect the plot.
#' @return A ggplot object.
#' @author Mitchell O'Hara-Wild
#' @seealso \code{\link[stats]{plot.acf}}, \code{\link{Acf}},
#' \code{\link[stats]{acf}}, \code{\link{taperedacf}}
#' @examples
#'
#' library(ggplot2)
#' ggAcf(wineind)
#' wineind %>% Acf(plot=FALSE) %>% autoplot
#' \dontrun{
#' wineind %>% taperedacf(plot=FALSE) %>% autoplot
#' ggtaperedacf(wineind)
#' ggtaperedpacf(wineind)}
#' ggCcf(mdeaths, fdeaths)
#'
#' @export
autoplot.acf <- function(object, ci=0.95, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!inherits(object, "acf")) {
      stop("autoplot.acf requires a acf object, use object=object")
    }

    acf <- `dimnames<-`(object$acf, list(NULL, object$snames, object$snames))
    lag <- `dimnames<-`(object$lag, list(NULL, object$snames, object$snames))

    data <- as.data.frame.table(acf)[-1]
    data$lag <- as.numeric(lag)

    if (object$type == "correlation" & is.null(object$ccf)) {
      data <- data[data$lag != 0, ]
    }

    # Initialise ggplot object
    p <- ggplot2::ggplot(
      ggplot2::aes(x = .data[["lag"]], xend = .data[["lag"]], y = 0, yend = .data[["Freq"]]),
      data = data
    )
    p <- p + ggplot2::geom_hline(yintercept = 0)

    # Add data
    p <- p + ggplot2::geom_segment(lineend = "butt", ...)

    # Add ci lines (assuming white noise input)
    ci <- qnorm((1 + ci) / 2) / sqrt(object$n.used)
    p <- p + ggplot2::geom_hline(yintercept = c(-ci, ci), colour = "blue", linetype = "dashed")

    # Add facets if needed
    if(any(dim(object$acf)[2:3] != c(1,1))){
      p <- p + ggplot2::facet_grid(
        as.formula(paste0(colnames(data)[1:2], collapse = "~"))
      )
    }

    # Prepare graph labels
    if (!is.null(object$ccf)) {
      ylab <- "CCF"
      ticktype <- "ccf"
      main <- paste("Series:", object$snames)
      nlags <- round(dim(object$lag)[1] / 2)
    }
    else if (object$type == "partial") {
      ylab <- "PACF"
      ticktype <- "acf"
      main <- paste("Series:", object$series)
      nlags <- dim(object$lag)[1]
    }
    else if (object$type == "correlation") {
      ylab <- "ACF"
      ticktype <- "acf"
      main <- paste("Series:", object$series)
      nlags <- dim(object$lag)[1]
    }
    else {
      ylab <- NULL
    }

    # Add seasonal x-axis
    # Change ticks to be seasonal and prepare default title
    if (!is.null(object$tsp)) {
      freq <- object$tsp[3]
    } else {
      freq <- 1
    }
    if (!is.null(object$periods)) {
      periods <- object$periods
      periods <- periods[periods != freq]
      minorbreaks <- periods * seq(-20:20)
    }
    else {
      minorbreaks <- NULL
    }
    p <- p + ggplot2::scale_x_continuous(breaks = seasonalaxis(
      freq,
      nlags, type = ticktype, plot = FALSE
    ), minor_breaks = minorbreaks)
    p <- p + ggAddExtras(ylab = ylab, xlab = "Lag", main = main)
    return(p)
  }
}

#' @rdname autoplot.acf
#' @export
ggAcf <- function(x, lag.max = NULL,
                  type = c("correlation", "covariance", "partial"),
                  plot = TRUE, na.action = na.contiguous, demean=TRUE, ...) {
  object <- Acf(x, lag.max = lag.max, type = type, na.action = na.action, demean = demean, plot = FALSE)
  object$tsp <- tsp(x)
  object$periods <- attributes(x)$msts
  object$series <- deparse(substitute(x))
  if (plot) {
    return(autoplot(object, ...))
  }
  else {
    return(object)
  }
}

#' @rdname autoplot.acf
#' @export
ggPacf <- function(x, lag.max = NULL,
                   plot = TRUE, na.action = na.contiguous, demean=TRUE, ...) {
  object <- Acf(x, lag.max = lag.max, type = "partial", na.action = na.action, demean = demean, plot = FALSE)
  object$tsp <- tsp(x)
  object$periods <- attributes(x)$msts
  object$series <- deparse(substitute(x))
  if (plot) {
    return(autoplot(object, ...))
  } else {
    return(object)
  }
}

#' @rdname autoplot.acf
#' @export
ggCcf <- function(x, y, lag.max=NULL, type=c("correlation", "covariance"),
                  plot=TRUE, na.action=na.contiguous, ...) {
  object <- Ccf(x, y, lag.max = lag.max, type = type, na.action = na.action, plot = FALSE)
  object$snames <- paste(deparse(substitute(x)), "&", deparse(substitute(y)))
  object$ccf <- TRUE
  if (plot) {
    return(autoplot(object, ...))
  }
  else {
    return(object)
  }
}

#' @rdname autoplot.acf
#' @export
autoplot.mpacf <- function(object, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!inherits(object, "mpacf")) {
      stop("autoplot.mpacf requires a mpacf object, use object=object")
    }
    if (!is.null(object$lower)) {
      data <- data.frame(Lag = 1:object$lag, z = object$z, sig = (object$lower < 0 & object$upper > 0))
      cidata <- data.frame(Lag = rep(1:object$lag, each = 2) + c(-0.5, 0.5), z = rep(object$z, each = 2), upper = rep(object$upper, each = 2), lower = rep(object$lower, each = 2))
      plotpi <- TRUE
    }
    else {
      data <- data.frame(Lag = 1:object$lag, z = object$z)
      plotpi <- FALSE
    }
    # Initialise ggplot object
    p <- ggplot2::ggplot()
    p <- p + ggplot2::geom_hline(ggplot2::aes(yintercept = 0), linewidth = 0.2)

    # Add data
    if (plotpi) {
      p <- p + ggplot2::geom_ribbon(ggplot2::aes(x = .data[["Lag"]], ymin = .data[["lower"]], ymax = .data[["upper"]]), data = cidata, fill = "grey50")
    }
    p <- p + ggplot2::geom_line(ggplot2::aes(x = .data[["Lag"]], y = .data[["z"]]), data = data)
    if (plotpi) {
      p <- p + ggplot2::geom_point(ggplot2::aes(x = .data[["Lag"]], y = .data[["z"]], colour = .data[["sig"]]), data = data)
    }

    # Change ticks to be seasonal
    freq <- frequency(object$x)
    msts <- is.element("msts", class(object$x))

    # Add seasonal x-axis
    if (msts) {
      periods <- attributes(object$x)$msts
      periods <- periods[periods != freq]
      minorbreaks <- periods * seq(-20:20)
    }
    else {
      minorbreaks <- NULL
    }

    p <- p + ggplot2::scale_x_continuous(
      breaks = seasonalaxis(frequency(object$x), length(data$Lag), type = "acf", plot = FALSE),
      minor_breaks = minorbreaks
    )

    if (object$type == "partial") {
      ylab <- "PACF"
    }
    else if (object$type == "correlation") {
      ylab <- "ACF"
    }

    p <- p + ggAddExtras(ylab = ylab)

    return(p)
  }
}

#' @rdname autoplot.acf
#' @export
ggtaperedacf <- function(x, lag.max=NULL, type=c("correlation", "partial"),
                         plot=TRUE, calc.ci=TRUE, level=95, nsim=100, ...) {
  cl <- match.call()
  if (plot) {
    cl$plot <- FALSE
  }
  cl[[1]] <- quote(taperedacf)
  object <- eval.parent(cl)
  if (plot) {
    return(autoplot(object, ...))
  }
  else {
    return(object)
  }
}

#' @rdname autoplot.acf
#' @export
ggtaperedpacf <- function(x, ...) {
  ggtaperedacf(x, type = "partial", ...)
}

#' @rdname plot.Arima
#' @export
autoplot.Arima <- function(object, type = c("both", "ar", "ma"), ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (is.Arima(object)) {
      # Detect type
      type <- match.arg(type)
      q <- p <- 0
      if (length(object$model$phi) > 0) {
        test <- abs(object$model$phi) > 1e-09
        if (any(test)) {
          p <- max(which(test))
        }
      }
      if (length(object$model$theta) > 0) {
        test <- abs(object$model$theta) > 1e-09
        if (any(test)) {
          q <- max(which(test))
        }
      }

      if (type == "both") {
        type <- c("ar", "ma")
      }
    }
    else if (inherits(object, "ar")) {
      type <- "ar"
      p <- length(arroots(object)$roots)
      q <- 0
    }
    else {
      stop("autoplot.Arima requires an Arima object")
    }

    # Remove NULL type
    type <- intersect(type, c("ar", "ma")[c(p > 0, q > 0)])

    # Prepare data
    arData <- maData <- NULL
    allRoots <- data.frame(roots = numeric(0), type = character(0))
    if ("ar" %in% type && p > 0) {
      arData <- arroots(object)
      allRoots <- rbind(allRoots, data.frame(roots = arData$roots, type = arData$type))
    }
    if ("ma" %in% type && q > 0) {
      maData <- maroots(object)
      allRoots <- rbind(allRoots, data.frame(roots = maData$roots, type = maData$type))
    }
    allRoots$Real <- Re(1 / allRoots$roots)
    allRoots$Imaginary <- Im(1 / allRoots$roots)
    allRoots$UnitCircle <- factor(ifelse((abs(allRoots$roots) > 1), "Within", "Outside"))

    # Initialise general ggplot object
    p <- ggplot2::ggplot(ggplot2::aes(x = .data[["Real"]], y = .data[["Imaginary"]], colour = .data[["UnitCircle"]]), data = allRoots)
    p <- p + ggplot2::coord_fixed(ratio = 1)
    p <- p + ggplot2::annotate(
      "path", x = cos(seq(0, 2 * pi, length.out = 100)),
      y = sin(seq(0, 2 * pi, length.out = 100))
    )
    p <- p + ggplot2::geom_vline(xintercept = 0)
    p <- p + ggplot2::geom_hline(yintercept = 0)
    p <- p + ggAddExtras(xlab = "Real", ylab = "Imaginary")

    if (NROW(allRoots) == 0) {
      return(p + ggAddExtras(main = "No AR or MA roots"))
    }

    p <- p + ggplot2::geom_point(size = 3)

    if (length(type) == 1) {
      p <- p + ggAddExtras(main = paste("Inverse", toupper(type), "roots"))
    }
    else {
      p <- p + ggplot2::facet_wrap(~ type, labeller = function(labels) lapply(labels, function(x) paste("Inverse", as.character(x), "roots")))
    }
  }
  return(p)
}

#' @rdname plot.Arima
#' @export
autoplot.ar <- function(object, ...) {
  autoplot.Arima(object, ...)
}

#' @rdname autoplot.seas
#' @export
autoplot.decomposed.ts <- function(object, labels=NULL, range.bars = NULL, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!inherits(object, "decomposed.ts")) {
      stop("autoplot.decomposed.ts requires a decomposed.ts object")
    }

    if (is.null(labels)) {
      labels <- c("trend", "seasonal", "remainder")
    }

    cn <- c("data", labels)

    data <- data.frame(
      datetime = rep(time(object$x), 4),
      y = c(object$x, object$trend, object$seasonal, object$random),
      parts = factor(rep(cn, each = NROW(object$x)), levels = cn)
    )

    # Initialise ggplot object
    p <- ggplot2::ggplot(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]), data = data)

    # Add data
    int <- as.numeric(object$type == "multiplicative")
    p <- p + ggplot2::geom_line(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]), data = subset(data, data$parts != cn[4]), na.rm = TRUE)
    p <- p + ggplot2::geom_segment(
      ggplot2::aes(x = .data[["datetime"]], xend = .data[["datetime"]], y = int, yend = .data[["y"]]),
      data = subset(data, data$parts == cn[4]), lineend = "butt", na.rm = TRUE
    )
    p <- p + ggplot2::facet_grid("parts ~ .", scales = "free_y", switch = "y")
    p <- p + ggplot2::geom_hline(ggplot2::aes(yintercept = .data[["y"]]),
                                 data = data.frame(y = int, parts = factor(cn[4], levels = cn)))

    if (is.null(range.bars)) {
      range.bars <- object$type == "additive"
    }
    if (range.bars) {
      yranges <- vapply(split(data$y, data$parts), function(x) range(x, na.rm = TRUE), numeric(2))
      xranges <- range(data$datetime)
      barmid <- apply(yranges, 2, mean)
      barlength <- min(apply(yranges, 2, diff))
      barwidth <- (1 / 64) * diff(xranges)
      barpos <- data.frame(
        left = xranges[2] + barwidth, right = xranges[2] + barwidth * 2,
        top = barmid + barlength / 2, bottom = barmid - barlength / 2,
        parts = factor(colnames(yranges), levels = cn),
        datetime = xranges[2], y = barmid
      )
      p <- p + ggplot2::geom_rect(
        ggplot2::aes(
          xmin = .data[["left"]], xmax = .data[["right"]],
          ymax = .data[["top"]], ymin = .data[["bottom"]]
        ), data = barpos, fill = "gray75", colour = "black", size = 1 / 3)
    }

    # Add axis labels
    p <- p + ggAddExtras(
      main = paste("Decomposition of", object$type, "time series"), xlab = "Time",
      ylab = ""
    )

    # Make x axis contain only whole numbers (e.g., years)
    p <- p + ggplot2::scale_x_continuous(breaks = unique(round(pretty(data$datetime))))

    return(p)
  }
}

#' @rdname plot.ets
#' @export
autoplot.ets <- function(object, range.bars = NULL, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!is.ets(object)) {
      stop("autoplot.ets requires an ets object, use object=object")
    }

    names <- c(y = "observed", l = "level", b = "slope", s1 = "season")
    data <- cbind(object$x, object$states[, colnames(object$states) %in% names(names)])
    cn <- c("y", c(colnames(object$states)))
    colnames(data) <- cn <- names[stats::na.exclude(match(cn, names(names)))]

    # Convert to longform
    data <- data.frame(
      datetime = rep(time(data), NCOL(data)), y = c(data),
      parts = factor(rep(cn, each = NROW(data)), levels = cn)
    )

    # Initialise ggplot object
    p <- ggplot2::ggplot(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]), data = data, ylab = "")

    # Add data
    p <- p + ggplot2::geom_line(na.rm = TRUE)
    p <- p + ggplot2::facet_grid(parts ~ ., scales = "free_y", switch = "y")
    if (is.null(range.bars)) {
      range.bars <- is.null(object$lambda)
    }
    if (range.bars) {
      yranges <- vapply(split(data$y, data$parts), function(x) range(x, na.rm = TRUE), numeric(2))
      xranges <- range(data$datetime)
      barmid <- apply(yranges, 2, mean)
      barlength <- min(apply(yranges, 2, diff))
      barwidth <- (1 / 64) * diff(xranges)
      barpos <- data.frame(
        left = xranges[2] + barwidth, right = xranges[2] + barwidth * 2,
        top = barmid + barlength / 2, bottom = barmid - barlength / 2,
        parts = factor(colnames(yranges), levels = cn),
        datetime = xranges[2], y = barmid
      )
      p <- p + ggplot2::geom_rect(
        ggplot2::aes(
          xmin = .data[["left"]], xmax = .data[["right"]],
          ymax = .data[["top"]], ymin = .data[["bottom"]]
        ), data = barpos, fill = "gray75", colour = "black", size = 1 / 3)
    }

    p <- p + ggAddExtras(xlab = NULL, ylab = "", main = paste("Components of", object$method, "method"))
    return(p)
  }
}

#' @rdname plot.bats
#' @export
autoplot.tbats <- function(object, range.bars = FALSE, ...) {
  cl <- match.call()
  cl[[1]] <- quote(autoplot.bats)
  eval.parent(cl)
}

#' @rdname plot.bats
#' @export
autoplot.bats <- function(object, range.bars = FALSE, ...) {
  data <- tbats.components(object)

  cn <- colnames(data)
  # Convert to longform
  data <- data.frame(
    datetime = rep(time(data), NCOL(data)), y = c(data),
    parts = factor(rep(cn, each = NROW(data)), levels = cn)
  )

  # Initialise ggplot object
  p <- ggplot2::ggplot(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]), data = data, ylab = "")

  # Add data
  p <- p + ggplot2::geom_line(na.rm = TRUE)
  p <- p + ggplot2::facet_grid(parts ~ ., scales = "free_y", switch = "y")

  if (range.bars) {
    yranges <- vapply(split(data$y, data$parts), function(x) range(x, na.rm = TRUE), numeric(2))
    xranges <- range(data$datetime)
    barmid <- apply(yranges, 2, mean)
    barlength <- min(apply(yranges, 2, diff))
    barwidth <- (1 / 64) * diff(xranges)
    barpos <- data.frame(
      left = xranges[2] + barwidth, right = xranges[2] + barwidth * 2,
      top = barmid + barlength / 2, bottom = barmid - barlength / 2,
      parts = factor(colnames(yranges), levels = cn),
      datetime = xranges[2], y = barmid
    )
    p <- p + ggplot2::geom_rect(
      ggplot2::aes(
        xmin = .data[["left"]], xmax = .data[["right"]],
        ymax = .data[["top"]], ymin = .data[["bottom"]]
      ), data = barpos, fill = "gray75", colour = "black", size = 1 / 3)
  }

  p <- p + ggAddExtras(xlab = NULL, ylab = "", main = paste("Components of", object$method, "method"))
  return(p)
}

#' @rdname plot.forecast
#' @export
autoplot.forecast <- function(object, include, PI=TRUE, shadecols=c("#596DD5", "#D5DBFF"), fcol="#0000AA", flwd=0.5, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!is.forecast(object)) {
      stop("autoplot.forecast requires a forecast object, use object=object")
    }
    if (is.null(object$lower) || is.null(object$upper) || is.null(object$level)) {
      PI <- FALSE
    }
    else if (!is.finite(max(object$upper))) {
      PI <- FALSE
    }

    if (!is.null(object$model$terms) && !is.null(object$model$model)) {
      # Initialise original dataset
      mt <- object$model$terms
      if (!is.null(object$series)) {
        yvar <- object$series
      } else {
        yvar <- deparse(mt[[2]])
      } # Perhaps a better way to do this
      xvar <- attr(mt, "term.labels")
      vars <- c(yvar = yvar, xvar = xvar)
      data <- object$model$model
      colnames(data) <- names(vars)[match(colnames(data), vars)]
      if (!is.null(object$model$lambda)) {
        data$yvar <- InvBoxCox(data$yvar, object$model$lambda)
      }
    }
    else {
      if (!is.null(object$x)) {
        data <- data.frame(yvar = c(object$x))
      }
      else if (!is.null(object$residuals) && !is.null(object$fitted)) {
        data <- data.frame(yvar = c(object$residuals + object$fitted))
      }
      else {
        stop("Could not find data")
      }
      if (!is.null(object$series)) {
        vars <- c(yvar = object$series)
      }
      else if (!is.null(object$model$call)) {
        vars <- c(yvar = deparse(object$model$call$y))
        if (vars == "object") {
          vars <- c(yvar = "y")
        }
      }
      else {
        vars <- c(yvar = "y")
      }
    }

    # Initialise ggplot object
    p <- ggplot2::ggplot()

    # Cross sectional forecasts
    if (!is.element("ts", class(object$mean))) {
      if (length(xvar) > 1) {
        stop("Forecast plot for regression models only available for a single predictor")
      }
      if (NCOL(object$newdata) == 1) { # Make sure column has correct name
        colnames(object$newdata) <- xvar
      }
      flwd <- 2 * flwd # Scale for points

      # Data points
      p <- p + ggplot2::geom_point(ggplot2::aes(x = .data[["xvar"]], y = .data[["yvar"]]), data = data)
      p <- p + ggplot2::labs(y = vars["yvar"], x = vars["xvar"])

      # Forecasted intervals
      if (PI) {
        levels <- NROW(object$level)
        interval <- data.frame(xpred = rep(object$newdata[[1]], levels), lower = c(object$lower), upper = c(object$upper), level = rep(object$level, each = NROW(object$newdata[[1]])))
        interval <- interval[order(interval$level, decreasing = TRUE), ] # Must be ordered for gg z-index
        p <- p + ggplot2::geom_linerange(ggplot2::aes(x = .data[["xpred"]], ymin = .data[["lower"]], ymax = .data[["upper"]], colour = .data[["level"]]), data = interval, linewidth = flwd)
        if (length(object$level) <= 5) {
          p <- p + ggplot2::scale_colour_gradientn(breaks = object$level, colours = shadecols, guide = "legend")
        }
        else {
          p <- p + ggplot2::scale_colour_gradientn(colours = shadecols, guide = "colourbar")
        }
      }

      # Forecasted points
      predicted <- data.frame(object$newdata, object$mean)
      colnames(predicted) <- c("xpred", "ypred")
      p <- p + ggplot2::geom_point(ggplot2::aes(x = .data[["xpred"]], y = .data[["ypred"]]), data = predicted, color = fcol, size = flwd)

      # Line of best fit
      coef <- data.frame(int = 0, m = 0)
      i <- match("(Intercept)", names(object$model$coefficients))
      if (i != 0) {
        coef$int <- object$model$coefficients[i]
        if (NROW(object$model$coefficients) == 2) {
          coef$m <- object$model$coefficients[-i]
        }
      }
      else {
        if (NROW(object$model$coefficients) == 1) {
          coef$m <- object$model$coefficients
        }
      }
      p <- p + ggplot2::geom_abline(intercept = coef$int, slope = coef$m)
    }
    else {
      # Time series objects (assumed)
      if(!missing(shadecols)){
        warning(
"The `schadecols` argument is deprecated for time series forecasts.
Interval shading is now done automatically based on the level and `fcol`.",
          call. = FALSE)
      }
      # Data points
      if (!is.null(time(object$x))) {
        timex <- time(object$x)
      }
      else if (!is.null(time(object$model$residuals))) {
        timex <- time(object$model$residuals)
      }
      data <- data.frame(yvar = as.numeric(data$yvar), datetime = as.numeric(timex))
      if (!missing(include)) {
        data <- tail(data, include)
      }
      p <- p + ggplot2::scale_x_continuous()
      p <- p + ggplot2::geom_line(ggplot2::aes(x = .data[["datetime"]], y = .data[["yvar"]]), data = data) +
        ggplot2::labs(y = vars["yvar"], x = "Time")

      # Forecasted intervals
      p <- p + autolayer(object, PI = PI, colour = fcol, size = flwd)

      # predicted <- data.frame(xvar = time(object$mean), yvar = object$mean)
      # colnames(predicted) <- c("datetime", "ypred")
      # if (PI) {
      #   levels <- NROW(object$level)
      #   interval <- data.frame(datetime = rep(predicted$datetime, levels), lower = c(object$lower), upper = c(object$upper), level = rep(object$level, each = NROW(object$mean)))
      #   interval <- interval[order(interval$level, decreasing = TRUE), ] # Must be ordered for gg z-index
      #   p <- p + ggplot2::geom_ribbon(ggplot2::aes_(x = ~datetime, ymin = ~lower, ymax = ~upper, group = ~-level, fill = ~level), data = interval)
      #   if (min(object$level) < 50) {
      #     scalelimit <- c(1, 99)
      #   }
      #   else {
      #     scalelimit <- c(50, 99)
      #   }
      #   if (length(object$level) <= 5) {
      #     p <- p + ggplot2::scale_fill_gradientn(breaks = object$level, colours = shadecols, limit = scalelimit, guide = "legend")
      #   }
      #   else {
      #     p <- p + ggplot2::scale_fill_gradientn(colours = shadecols, limit = scalelimit)
      #   }
      #   # Negative group is a work around for missing z-index
      # }

      # # Forecasted points
      # p <- p + ggplot2::geom_line(ggplot2::aes_(x = ~datetime, y = ~ypred), data = predicted, color = fcol, size = flwd)
    }

    p <- p + ggAddExtras(main = paste("Forecasts from ", object$method, sep = ""))
    return(p)
  }
}

#' @rdname plot.mforecast
#' @export
autoplot.mforecast <- function(object, PI = TRUE, facets = TRUE, colour = FALSE, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!is.mforecast(object)) {
      stop("autoplot.mforecast requires a mforecast object, use object=object")
    }
    if (is.ts(object$forecast[[1]]$mean)) {
      # ts forecasts
      p <- autoplot(getResponse(object), facets = facets, colour = colour) + autolayer(object, ...)
      if (facets) {
        p <- p + ggplot2::facet_wrap(
          ~ series,
          labeller = function(labels) {
            if (!is.null(object$method)) {
              lapply(labels, function(x) paste0(as.character(x), "\n", object$method[as.character(x)]))
            }
            else {
              lapply(labels, function(x) paste0(as.character(x)))
            }
          },
          ncol = 1,
          scales = "free_y"
        )
      }
      p <- p + ggAddExtras(ylab = NULL)
      return(p)
    }
    else {
      # lm forecasts
      if (!requireNamespace("grid")) {
        stop("grid is needed for this function to work. Install it via install.packages(\"grid\")", call. = FALSE)
      }

      K <- length(object$forecast)
      if (K < 2) {
        warning("Expected at least two plots but forecast required less.")
      }

      # Set up vector arguments
      if (missing(PI)) {
        PI <- rep(TRUE, K)
      }

      # Set up grid
      # ncol: Number of columns of plots
      # nrow: Number of rows needed, calculated from # of cols
      gridlayout <- matrix(seq(1, K), ncol = 1, nrow = K)

      grid::grid.newpage()
      grid::pushViewport(grid::viewport(layout = grid::grid.layout(nrow(gridlayout), ncol(gridlayout))))

      for (i in 1:K) {
        partialfcast <- object$forecast[[i]]
        partialfcast$model <- mlmsplit(object$model, index = i)
        matchidx <- as.data.frame(which(gridlayout == i, arr.ind = TRUE))
        print(
          autoplot(
            structure(partialfcast, class = "forecast"),
            PI = PI[i], ...
          ) + ggAddExtras(ylab = names(object$forecast)[i]),
          vp = grid::viewport(
            layout.pos.row = matchidx$row,
            layout.pos.col = matchidx$col
          )
        )
      }
    }
  }
}

#' @rdname tsdisplay
#'
#' @examples
#' library(ggplot2)
#' ggtsdisplay(USAccDeaths, plot.type="scatter", theme=theme_bw())
#'
#' @export
ggtsdisplay <- function(x, plot.type=c("partial", "histogram", "scatter", "spectrum"),
                        points=TRUE, smooth=FALSE,
                        lag.max, na.action=na.contiguous, theme=NULL, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else if (!requireNamespace("grid", quietly = TRUE)) {
    stop("grid is needed for this function to work. Install it via install.packages(\"grid\")", call. = FALSE)
  }
  else {
    if (NCOL(x) > 1) {
      stop("ggtsdisplay is only for univariate time series")
    }
    plot.type <- match.arg(plot.type)
    main <- deparse(substitute(x))

    if (!is.ts(x)) {
      x <- ts(x)
    }
    if (missing(lag.max)) {
      lag.max <- round(min(max(10 * log10(length(x)), 3 * frequency(x)), length(x) / 3))
    }

    dots <- list(...)
    if (is.null(dots$xlab)) {
      dots$xlab <- ""
    }
    if (is.null(dots$ylab)) {
      dots$ylab <- ""
    }
    labs <- match(c("xlab", "ylab", "main"), names(dots), nomatch = 0)

    # Set up grid for plots
    gridlayout <- matrix(c(1, 2, 1, 3), nrow = 2)
    grid::grid.newpage()
    grid::pushViewport(grid::viewport(layout = grid::grid.layout(nrow(gridlayout), ncol(gridlayout))))

    # Add ts plot with points
    matchidx <- as.data.frame(which(gridlayout == 1, arr.ind = TRUE))
    tsplot <- do.call(ggplot2::autoplot, c(object = quote(x), dots[labs]))
    if (points) {
      tsplot <- tsplot + ggplot2::geom_point(size = 0.5)
    }
    if (smooth) {
      tsplot <- tsplot + ggplot2::geom_smooth(method = "loess", se = FALSE)
    }
    if (is.null(tsplot$labels$title)) { # Add title if missing
      tsplot <- tsplot + ggplot2::ggtitle(main)
    }
    if (!is.null(theme)) {
      tsplot <- tsplot + theme
    }
    print(
      tsplot,
      vp = grid::viewport(
        layout.pos.row = matchidx$row,
        layout.pos.col = matchidx$col
      )
    )

    # Prepare Acf plot
    acfplot <- do.call(ggAcf, c(x = quote(x), lag.max = lag.max, na.action = na.action, dots[-labs])) + ggplot2::ggtitle(NULL)
    if (!is.null(theme)) {
      acfplot <- acfplot + theme
    }

    # Prepare last plot (variable)
    if (plot.type == "partial") {
      lastplot <- ggPacf(x, lag.max = lag.max, na.action = na.action) + ggplot2::ggtitle(NULL)
      # Match y-axis
      acfplotrange <- ggplot2::layer_scales(acfplot)$y$range$range
      pacfplotrange <- ggplot2::layer_scales(lastplot)$y$range$range
      yrange <- range(c(acfplotrange, pacfplotrange))
      acfplot <- acfplot + ggplot2::ylim(yrange)
      lastplot <- lastplot + ggplot2::ylim(yrange)
    }
    else if (plot.type == "histogram") {
      lastplot <- gghistogram(x, add.normal = TRUE, add.rug = TRUE) + ggplot2::xlab(main)
    }
    else if (plot.type == "scatter") {
      scatterData <- data.frame(y = x[2:NROW(x)], x = x[1:NROW(x) - 1])
      lastplot <- ggplot2::ggplot(ggplot2::aes(y = .data[["y"]], x = .data[["x"]]), data = scatterData) +
        ggplot2::geom_point() + ggplot2::labs(x = expression(Y[t - 1]), y = expression(Y[t]))
    }
    else if (plot.type == "spectrum") {
      specData <- spec.ar(x, plot = FALSE)
      specData <- data.frame(spectrum = specData$spec, frequency = specData$freq)
      lastplot <- ggplot2::ggplot(ggplot2::aes(y = .data[["spectrum"]], x = .data[["frequency"]]), data = specData) +
        ggplot2::geom_line() + ggplot2::scale_y_log10()
    }
    if (!is.null(theme)) {
      lastplot <- lastplot + theme
    }

    # Add ACF plot
    matchidx <- as.data.frame(which(gridlayout == 2, arr.ind = TRUE))
    print(
      acfplot,
      vp = grid::viewport(
        layout.pos.row = matchidx$row,
        layout.pos.col = matchidx$col
      )
    )

    # Add last plot
    matchidx <- as.data.frame(which(gridlayout == 3, arr.ind = TRUE))
    print(
      lastplot,
      vp = grid::viewport(
        layout.pos.row = matchidx$row,
        layout.pos.col = matchidx$col
      )
    )
  }
}


#' Time series lag ggplots
#'
#' Plots a lag plot using ggplot.
#'
#' \dQuote{gglagplot} will plot time series against lagged versions of
#' themselves. Helps visualising 'auto-dependence' even when auto-correlations
#' vanish.
#'
#' \dQuote{gglagchull} will layer convex hulls of the lags, layered on a single
#' plot. This helps visualise the change in 'auto-dependence' as lags increase.
#'
#' @param x a time series object (type \code{ts}).
#' @param lags number of lag plots desired, see arg set.lags.
#' @param set.lags vector of positive integers specifying which lags to use.
#' @param diag logical indicating if the x=y diagonal should be drawn.
#' @param diag.col color to be used for the diagonal if(diag).
#' @param do.lines if TRUE, lines will be drawn, otherwise points will be
#' drawn.
#' @param colour logical indicating if lines should be coloured.
#' @param continuous Should the colour scheme for years be continuous or
#' discrete?
#' @param labels logical indicating if labels should be used.
#' @param seasonal Should the line colour be based on seasonal characteristics
#' (TRUE), or sequential (FALSE).
#' @param \dots Not used (for consistency with lag.plot)
#' @return None.
#' @author Mitchell O'Hara-Wild
#' @seealso \code{\link[stats]{lag.plot}}
#' @examples
#'
#' gglagplot(woolyrnq)
#' gglagplot(woolyrnq,seasonal=FALSE)
#'
#' lungDeaths <- cbind(mdeaths, fdeaths)
#' gglagplot(lungDeaths, lags=2)
#' gglagchull(lungDeaths, lags=6)
#'
#' @export
gglagplot <- function(x, lags=ifelse(frequency(x) > 9, 16, 9),
                      set.lags = 1:lags, diag=TRUE, diag.col="gray", do.lines = TRUE, colour = TRUE,
                      continuous = frequency(x) > 12, labels = FALSE, seasonal = TRUE, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    freq <- frequency(x)
    if (freq > 1) {
      linecol <- cycle(x)
      if (freq > 24) {
        continuous <- TRUE
      }
    }
    else {
      seasonal <- FALSE
      continuous <- TRUE
    }
    if (!seasonal) {
      continuous <- TRUE
    }

    # Make sure lags is evaluated
    tmp <- lags
    x <- as.matrix(x)

    # Prepare data for plotting
    n <- NROW(x)
    data <- data.frame()
    for (i in 1:NCOL(x)) {
      for (lagi in set.lags) {
        sname <- colnames(x)[i]
        if (is.null(sname)) {
          sname <- deparse(match.call()$x)
        }
        data <- rbind(
          data,
          data.frame(
            lagnum = 1:(n - lagi),
            freqcur = if(seasonal) linecol[(lagi + 1):n] else (lagi + 1):n,
            orig = x[(lagi + 1):n, i],
            lagged = x[1:(n - lagi), i],
            lagVal = rep(lagi, n - lagi),
            series = factor(rep(sname, n - lagi))
          )
        )
      }
    }
    if (!continuous) {
      data$freqcur <- factor(data$freqcur)
    }

    # Initialise ggplot object
    p <- ggplot2::ggplot(ggplot2::aes(x = .data[["lagged"]], y = .data[["orig"]]), data = data)

    if (diag) {
      p <- p + ggplot2::geom_abline(colour = diag.col, linetype = "dashed")
    }
    if (labels) {
      linesize <- 0.25 * (2 - do.lines)
    }
    else {
      linesize <- 0.5 * (2 - do.lines)
    }
    plottype <- if (do.lines) {
      function(...) ggplot2::geom_path(..., linewidth = linesize)
    }
    else {
      function(...) ggplot2::geom_point(..., size = linesize)
    }
    if (colour) {
      p <- p + plottype(ggplot2::aes(colour = .data[["freqcur"]]))
    }
    else {
      p <- p + plottype()
    }

    if (labels) {
      p <- p + ggplot2::geom_text(ggplot2::aes(label = .data[["lagnum"]]))
    }
    # Ensure all facets are of size size (if extreme values are excluded in lag specification)
    if (max(set.lags) > NROW(x) / 2) {
      axissize <- rbind(aggregate(orig ~ series, data = data, min), aggregate(orig~ series, data = data, max))
      axissize <- data.frame(series = rep(axissize$series, length(set.lags)), orig = rep(axissize$orig, length(set.lags)), lagVal = rep(set.lags, each = NCOL(x)))
      p <- p + ggplot2::geom_blank(ggplot2::aes(x = .data[["orig"]], y = .data[["orig"]]), data = axissize)
    }

    # Facet
    labellerFn <- function(labels) {
      if (!is.null(labels$series)) {
        # Multivariate labels
        labels$series <- as.character(labels$series)
      }
      labels$lagVal <- paste("lag", labels$lagVal)
      return(labels)
    }
    if (NCOL(x) > 1) {
      p <- p + ggplot2::facet_wrap(~series + lagVal, scales = "free", labeller = labellerFn)
    }
    else {
      p <- p + ggplot2::facet_wrap(~lagVal, labeller = labellerFn)
    }
    p <- p + ggplot2::theme(aspect.ratio = 1)
    if (colour) {
      if (seasonal) {
        if (freq == 4L) {
          title <- "Quarter"
        } else if (freq == 12L) {
          title <- "Month"
        } else if (freq == 7L) {
          title <- "Day"
        } else if (freq == 24L) {
          title <- "Hour"
        } else {
          title <- "Season"
        }
      }
      else {
        title <- "Time"
      }
      if (continuous) {
        p <- p + ggplot2::guides(colour = ggplot2::guide_colourbar(title = title))
      }
      else {
        p <- p + ggplot2::guides(colour = ggplot2::guide_legend(title = title))
      }
    }

    p <- p + ggAddExtras(ylab = NULL, xlab = NULL)

    return(p)
  }
}

#' @rdname gglagplot
#'
#' @examples
#' gglagchull(woolyrnq)
#'
#' @export
gglagchull <- function(x,
                       lags=ifelse(frequency(x) > 1, min(12, frequency(x)), 4),
                       set.lags = 1:lags, diag=TRUE, diag.col="gray", ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    # Make sure lags is evaluated
    tmp <- lags
    x <- as.matrix(x)

    # Prepare data for plotting
    n <- NROW(x)
    data <- data.frame()
    for (i in 1:NCOL(x)) {
      for (lag in set.lags) {
        sname <- colnames(x)[i]
        if (is.null(sname)) {
          sname <- deparse(substitute(x))
        }
        data <- rbind(data, data.frame(orig = x[(lag + 1):n, i], lagged = x[1:(n - lag), i], lag = rep(lag, n - lag), series = rep(sname, n - lag))[grDevices::chull(x[(lag + 1):n, i], x[1:(n - lag), i]), ])
      }
    }

    # Initialise ggplot object
    p <- ggplot2::ggplot(ggplot2::aes(x = .data[["orig"]], y = .data[["lagged"]]), data = data)
    if (diag) {
      p <- p + ggplot2::geom_abline(colour = diag.col, linetype = "dashed")
    }
    p <- p + ggplot2::geom_polygon(ggplot2::aes(group = .data[["lag"]], colour = .data[["lag"]], fill = .data[["lag"]]), alpha = 1 / length(set.lags))
    p <- p + ggplot2::guides(colour = ggplot2::guide_colourbar(title = "lag"))
    p <- p + ggplot2::theme(aspect.ratio = 1)

    # Facet
    if (NCOL(x) > 1) {
      p <- p + ggplot2::facet_wrap(~series, scales = "free")
    }

    p <- p + ggAddExtras(ylab = "lagged", xlab = "original")

    return(p)
  }
}


#' Create a seasonal subseries ggplot
#'
#' Plots a subseries plot using ggplot. Each season is plotted as a separate
#' mini time series. The blue lines represent the mean of the observations
#' within each season.
#'
#' The \code{ggmonthplot} function is simply a wrapper for
#' \code{ggsubseriesplot} as a convenience for users familiar with
#' \code{\link[stats]{monthplot}}.
#'
#' @param x a time series object (type \code{ts}).
#' @param labels A vector of labels to use for each 'season'
#' @param times A vector of times for each observation
#' @param phase A vector of seasonal components
#' @param \dots Not used (for consistency with monthplot)
#' @return Returns an object of class \code{ggplot}.
#' @author Mitchell O'Hara-Wild
#' @seealso \code{\link[stats]{monthplot}}
#' @examples
#'
#' ggsubseriesplot(AirPassengers)
#' ggsubseriesplot(woolyrnq)
#'
#' @export
ggmonthplot <- function(x, labels = NULL, times = time(x), phase = cycle(x), ...) {
  ggsubseriesplot(x, labels, times, phase, ...)
}

#' @rdname ggmonthplot
#' @export
ggsubseriesplot <- function(x, labels = NULL, times = time(x), phase = cycle(x), ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!inherits(x, "ts")) {
      stop("ggsubseriesplot requires a ts object, use x=object")
    }

    if (round(frequency(x)) <= 1) {
      stop("Data are not seasonal")
    }

    if("1" %in% dimnames(table(table(phase)))[[1]]){
      stop(paste("Each season requires at least 2 observations.",
                 ifelse(frequency(x)%%1 == 0,
                        "Your series length may be too short for this graphic.",
                        "This may be caused from specifying a time-series with non-integer frequency.")
      )
      )
    }

    data <- data.frame(y = as.numeric(x), year = trunc(time(x)), season = as.numeric(phase))
    seasonwidth <- (max(data$year) - min(data$year)) * 1.05
    data$time <- data$season + 0.025 + (data$year - min(data$year)) / seasonwidth
    avgLines <- stats::aggregate(data$y, by = list(data$season), FUN = mean)
    colnames(avgLines) <- c("season", "avg")
    data <- merge(data, avgLines, by = "season")

    # Initialise ggplot object
    # p <- ggplot2::ggplot(ggplot2::aes_(x=~interaction(year, season), y=~y, group=~season), data=data, na.rm=TRUE)
    p <- ggplot2::ggplot(
      ggplot2::aes(x = .data[["time"]], y = .data[["y"]], group = .data[["season"]]),
      data = data, na.rm = TRUE
    )

    # Remove vertical break lines
    p <- p + ggplot2::theme(panel.grid.major.x = ggplot2::element_blank())

    # Add data
    p <- p + ggplot2::geom_line()

    # Add average lines
    p <- p + ggplot2::geom_line(ggplot2::aes(y = .data[["avg"]]), col = "#0000AA")

    # Create x-axis labels
    xfreq <- frequency(x)
    if (!is.null(labels)) {
      if (xfreq != length(labels)) {
        stop("The number of labels supplied is not the same as the number of seasons.")
      } else {
        xbreaks <- labels
      }
    }
    else if (xfreq == 4) {
      xbreaks <- c("Q1", "Q2", "Q3", "Q4")
      xlab <- "Quarter"
    }
    else if (xfreq == 7) {
      xbreaks <- c(
        "Sunday", "Monday", "Tuesday", "Wednesday", "Thursday",
        "Friday", "Saturday"
      )
      xlab <- "Day"
    }
    else if (xfreq == 12) {
      xbreaks <- month.abb
      xlab <- "Month"
    }
    else {
      xbreaks <- 1:frequency(x)
      xlab <- "Season"
    }

    # X-axis
    p <- p + ggplot2::scale_x_continuous(breaks = 0.5 + (1:xfreq), labels = xbreaks)

    # Graph labels
    p <- p + ggAddExtras(ylab = deparse(substitute(x)), xlab = xlab)
    return(p)
  }
}

#' @rdname seasonplot
#'
#' @param continuous Should the colour scheme for years be continuous or
#' discrete?
#' @param polar Plot the graph on seasonal coordinates
#'
#' @examples
#' ggseasonplot(AirPassengers, col=rainbow(12), year.labels=TRUE)
#' ggseasonplot(AirPassengers, year.labels=TRUE, continuous=TRUE)
#'
#' @export
ggseasonplot <- function(x, season.labels=NULL, year.labels=FALSE, year.labels.left=FALSE, type=NULL, col=NULL, continuous=FALSE, polar=FALSE, labelgap=0.04, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  if (!inherits(x, "ts")) {
    stop("autoplot.seasonplot requires a ts object, use x=object")
  }
  if (!is.null(type)) {
    message("Plot types are not yet supported for seasonplot()")
  }

  # Check data are seasonal and convert to integer seasonality
  s <- round(frequency(x))
  if (s <= 1) {
    stop("Data are not seasonal")
  }

  # Grab name for plot title
  xname <- deparse(substitute(x))

  tspx <- tsp(x)
  x <- ts(x, start = tspx[1], frequency = s)

  data <- data.frame(
    y = as.numeric(x),
    year = trunc(round(time(x), 8)),
    cycle = as.numeric(cycle(x)),
    time = as.numeric((cycle(x) - 1) / s)
  )
  data$year <- if (continuous) {
    as.numeric(data$year)
  }
  else {
    as.factor(data$year)
  }
  if (polar) {
    startValues <- data[data$cycle == 1, ]
    if (data$cycle[1] == 1) {
      startValues <- startValues[-1, ]
    }
    startValues$time <- 1 - .Machine$double.eps
    levels(startValues$year) <- as.numeric(levels(startValues$year)) - 1
    data <- rbind(data, startValues)
  }
  # Initialise ggplot object
  p <- ggplot2::ggplot(ggplot2::aes(x = .data[["time"]], y = .data[["y"]], group = .data[["year"]], colour = .data[["year"]]), data = data, na.rm = TRUE)
  # p <- p + ggplot2::scale_x_continuous()

  # Add data
  p <- p + ggplot2::geom_line()

  if (!is.null(col)) {
    if(is.numeric(col)){
      col <- palette()[(col-1)%%(length(palette())) + 1]
    }
    if (continuous) {
      p <- p + ggplot2::scale_color_gradientn(colours = col)
    }
    else {
      ncol <- length(unique(data$year))
      if (length(col) == 1) {
        p <- p + ggplot2::scale_color_manual(guide = "none", values = rep(col, ncol))
      }
      else {
        p <- p + ggplot2::scale_color_manual(values = rep(col, ceiling(ncol / length(col)))[1:ncol])
      }
    }
  }

  if (year.labels) {
    yrlab <- stats::aggregate(time ~ year, data = data, FUN = max)
    yrlab <- cbind(yrlab, offset = labelgap)
  }
  if (year.labels.left) {
    yrlabL <- stats::aggregate(time ~ year, data = data, FUN = min)
    yrlabL <- cbind(yrlabL, offset = -labelgap)
    if (year.labels) {
      yrlab <- rbind(yrlab, yrlabL)
    }
  }
  if (year.labels || year.labels.left) {
    yrlab <- merge(yrlab, data)
    yrlab$time <- yrlab$time + yrlab$offset
    p <- p + ggplot2::guides(colour = "none")
    p <- p + ggplot2::geom_text(ggplot2::aes(x = .data[["time"]], y = .data[["y"]], label = .data[["year"]]), data = yrlab)
  }

  # Add seasonal labels
  if (s == 12) {
    labs <- month.abb
    xLab <- "Month"
  }
  else if (s == 4) {
    labs <- paste("Q", 1:4, sep = "")
    xLab <- "Quarter"
  }
  else if (s == 7) {
    labs <- c("Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat")
    xLab <- "Day"
  }
  else if (s == 52) {
    labs <- 1:s
    xLab <- "Week"
  }
  else if (s == 24) {
    labs <- 0:(s - 1)
    xLab <- "Hour"
  }
  else if (s == 48) {
    labs <- seq(0, 23.5, by = 0.5)
    xLab <- "Half-hour"
  }
  else {
    labs <- 1:s
    xLab <- "Season"
  }

  if (!is.null(season.labels)) {
    if (length(season.labels) != length(labs)) {
      warning(paste0("Provided season.labels have length ", length(season.labels), ", but ", length(labs), " are required. Ignoring season.labels."))
    }
    else {
      labs <- season.labels
    }
  }

  breaks <- sort(unique(data$time))
  if (polar) {
    breaks <- head(breaks, -1)
    p <- p + ggplot2::coord_polar()
  }

  p <- p + ggplot2::scale_x_continuous(breaks = breaks, minor_breaks = NULL, labels = labs)

  # Graph title and axes
  p <- p + ggAddExtras(main = paste("Seasonal plot:", xname), xlab = xLab, ylab = NULL)
  return(p)
}

#' @rdname plot.forecast
#' @export
autoplot.splineforecast <- function(object, PI=TRUE, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    p <- autoplot(object$x) + autolayer(object)
    p <- p + ggplot2::geom_point(size = 2)
    fit <- data.frame(datetime = as.numeric(time(object$fitted)), y = as.numeric(object$fitted))
    p <- p + ggplot2::geom_line(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]), colour = "red", data = fit)
    p <- p + ggAddExtras(ylab = deparse(object$model$call$x))
    if (!is.null(object$series)) {
      p <- p + ggplot2::ylab(object$series)
    }
    return(p)
  }
}

#' @rdname autoplot.seas
#' @export
autoplot.stl <- function(object, labels = NULL, range.bars = TRUE, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!inherits(object, "stl")) {
      stop("autoplot.stl requires a stl object, use x=object")
    }
    # Re-order series as trend, seasonal, remainder
    object$time.series <- object$time.series[, c("trend", "seasonal", "remainder")]
    if (is.null(labels)) {
      labels <- colnames(object$time.series)
    }

    data <- object$time.series
    cn <- c("data", labels)
    data <- data.frame(
      datetime = rep(time(data), NCOL(data) + 1), y = c(rowSums(data), data),
      parts = factor(rep(cn, each = NROW(data)), levels = cn)
    )

    # Initialise ggplot object
    p <- ggplot2::ggplot(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]), data = data)

    # Add data
    # Time series lines
    p <- p + ggplot2::geom_line(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]), data = subset(data, data$parts != cn[4]), na.rm = TRUE)
    p <- p + ggplot2::geom_segment(
      ggplot2::aes(x = .data[["datetime"]], xend = .data[["datetime"]], y = 0, yend = .data[["y"]]),
      data = subset(data, data$parts == cn[4]), lineend = "butt"
    )

    # Rangebars
    if (range.bars) {
      yranges <- vapply(split(data$y, data$parts), function(x) range(x, na.rm = TRUE), numeric(2))
      xranges <- range(data$datetime)
      barmid <- apply(yranges, 2, mean)
      barlength <- min(apply(yranges, 2, diff))
      barwidth <- (1 / 64) * diff(xranges)
      barpos <- data.frame(
        left = xranges[2] + barwidth, right = xranges[2] + barwidth * 2,
        top = barmid + barlength / 2, bottom = barmid - barlength / 2,
        parts = factor(colnames(yranges), levels = cn),
        datetime = xranges[2], y = barmid
      )
      p <- p + ggplot2::geom_rect(ggplot2::aes(xmin = .data[["left"]], xmax = .data[["right"]], ymax = .data[["top"]], ymin = .data[["bottom"]]), data = barpos, fill = "gray75", colour = "black", size = 1 / 3)
    }

    # Remainder
    p <- p + ggplot2::facet_grid("parts ~ .", scales = "free_y", switch = "y")
    p <- p + ggplot2::geom_hline(ggplot2::aes(yintercept = .data[["y"]]),
                                 data = data.frame(y = 0, parts = factor(cn[4], levels = cn)))

    # Add axis labels
    p <- p + ggAddExtras(xlab = "Time", ylab = "")

    # Make x axis contain only whole numbers (e.g., years)
    p <- p + ggplot2::scale_x_continuous(breaks = unique(round(pretty(data$datetime))))
    # ^^ Remove rightmost x axis gap with `expand=c(0.05, 0, 0, 0)` argument when assymetric `expand` feature is supported
    # issue: tidyverse/ggplot2#1669

    return(p)
  }
}

#' @rdname autoplot.seas
#' @export
autoplot.StructTS <- function(object, labels = NULL, range.bars = TRUE, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!inherits(object, "StructTS")) {
      stop("autoplot.StructTS requires a StructTS object.")
    }

    if (is.null(labels)) {
      labels <- colnames(object$fitted)
    }

    data <- object$fitted
    cn <- c("data", labels)
    data <- data.frame(
      datetime = rep(time(data), NCOL(data) + 1), y = c(object$data, data),
      parts = factor(rep(cn, each = NROW(data)), levels = cn)
    )

    # Initialise ggplot object
    p <- ggplot2::ggplot(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]), data = data)

    # Add data
    p <- p + ggplot2::geom_line(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]), na.rm = TRUE)
    p <- p + ggplot2::facet_grid("parts ~ .", scales = "free_y", switch = "y")

    if (range.bars) {
      yranges <- vapply(split(data$y, data$parts), function(x) range(x, na.rm = TRUE), numeric(2))
      xranges <- range(data$datetime)
      barmid <- apply(yranges, 2, mean)
      barlength <- min(apply(yranges, 2, diff))
      barwidth <- (1 / 64) * diff(xranges)
      barpos <- data.frame(
        left = xranges[2] + barwidth, right = xranges[2] + barwidth * 2,
        top = barmid + barlength / 2, bottom = barmid - barlength / 2,
        parts = factor(colnames(yranges), levels = cn),
        datetime = xranges[2], y = barmid
      )
      p <- p + ggplot2::geom_rect(ggplot2::aes(xmin = .data[["left"]], xmax = .data[["right"]], ymax = .data[["top"]], ymin = .data[["bottom"]]), data = barpos, fill = "gray75", colour = "black", size = 1 / 3)
    }

    # Add axis labels
    p <- p + ggAddExtras(xlab = "Time", ylab = "")

    # Make x axis contain only whole numbers (e.g., years)
    p <- p + ggplot2::scale_x_continuous(breaks = unique(round(pretty(data$datetime))))

    return(p)
  }
}


#' Plot time series decomposition components using ggplot
#'
#' Produces a ggplot object of seasonally decomposed time series for objects of
#' class \dQuote{\code{stl}} (created with \code{\link[stats]{stl}}), class
#' \dQuote{\code{seas}} (created with \code{\link[seasonal]{seas}}), or class
#' \dQuote{\code{decomposed.ts}} (created with \code{\link[stats]{decompose}}).
#'
#' @param object Object of class \dQuote{\code{seas}}, \dQuote{\code{stl}}, or
#' \dQuote{\code{decomposed.ts}}.
#' @param labels Labels to replace \dQuote{seasonal}, \dQuote{trend}, and
#' \dQuote{remainder}.
#' @param range.bars Logical indicating if each plot should have a bar at its
#' right side representing relative size. If \code{NULL}, automatic selection
#' takes place.
#' @param ... Other plotting parameters to affect the plot.
#' @return Returns an object of class \code{ggplot}.
#' @author Mitchell O'Hara-Wild
#' @seealso \code{\link[seasonal]{seas}}, \code{\link[stats]{stl}},
#' \code{\link[stats]{decompose}}, \code{\link[stats]{StructTS}},
#' \code{\link[stats]{plot.stl}}.
#' @examples
#'
#' library(ggplot2)
#' co2 %>%
#'   decompose() %>%
#'   autoplot()
#' nottem %>%
#'   stl(s.window = "periodic") %>%
#'   autoplot()
#' \dontrun{
#' library(seasonal)
#' seas(USAccDeaths) %>% autoplot()
#' }
#'
#' @export
autoplot.seas <- function(object, labels = NULL, range.bars = NULL, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  if (!inherits(object, "seas")) {
    stop("autoplot.seas requires a seas object")
  }
  if (is.null(labels)) {
    if ("seasonal" %in% colnames(object$data)) {
      labels <- c("trend", "seasonal", "irregular")
    } else {
      labels <- c("trend", "irregular")
    }
  }
  data <- cbind(object$x, object$data[, labels])
  colnames(data) <- cn <- c("data", labels)

  data <- data.frame(
    datetime = rep(time(data), NCOL(data)), y = c(data),
    parts = factor(rep(cn, each = NROW(data)), levels = cn)
  )

  # Is it additive or multiplicative?
  freq <- frequency(object$data)
  sum_first_year <- try(sum(seasonal(object)[seq(freq)]), silent=TRUE)
  if(!inherits(sum_first_year, "try-error")) {
    int <- as.integer(sum_first_year > 0.5) # Closer to 1 than 0.
  } else {
    int <- 0
  }

  # Initialise ggplot object
  p <- ggplot2::ggplot(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]), data = data)
  # Add data
  p <- p + ggplot2::geom_line(ggplot2::aes(x = .data[["datetime"]], y = .data[["y"]]),
                              data = subset(data, data$parts != tail(cn,1)),
                              na.rm = TRUE)
  p <- p + ggplot2::geom_segment(
    ggplot2::aes(x = .data[["datetime"]], xend = .data[["datetime"]], y = int, yend = .data[["y"]]),
    data = subset(data, data$parts == tail(cn,1)), lineend = "butt"
  )
  p <- p + ggplot2::facet_grid("parts ~ .", scales = "free_y", switch = "y")
  p <- p + ggplot2::geom_hline(ggplot2::aes(yintercept = .data[["y"]]),
    data = data.frame(y = int, parts = factor(tail(cn,1), levels = cn))
  )

  # Rangebars
  if (is.null(range.bars)) {
    range.bars <- object$spc$transform$`function` == "none"
  }
  if (range.bars) {
    yranges <- vapply(split(data$y, data$parts), function(x) range(x, na.rm = TRUE), numeric(2))
    xranges <- range(data$datetime)
    barmid <- apply(yranges, 2, mean)
    barlength <- min(apply(yranges, 2, diff))
    barwidth <- (1 / 64) * diff(xranges)
    barpos <- data.frame(
      left = xranges[2] + barwidth, right = xranges[2] + barwidth * 2,
      top = barmid + barlength / 2, bottom = barmid - barlength / 2,
      parts = factor(colnames(yranges), levels = cn),
      datetime = xranges[2], y = barmid
    )
    p <- p + ggplot2::geom_rect(ggplot2::aes(xmin = .data[["left"]], xmax = .data[["right"]], ymax = .data[["top"]], ymin = .data[["bottom"]]), data = barpos, fill = "gray75", colour = "black", size = 1 / 3)
  }

  # Add axis labels
  p <- p + ggAddExtras(xlab = "Time", ylab = "")

  # Make x axis contain only whole numbers (e.g., years)
  p <- p + ggplot2::scale_x_continuous(breaks = unique(round(pretty(data$datetime))))

  return(p)
}

#' @rdname autoplot.ts
#' @export
autolayer.mts <- function(object, colour = TRUE, series = NULL, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    cl <- match.call()
    cl[[1]] <- quote(autolayer)
    cl$object <- quote(object[, i])
    if (length(series) != NCOL(object)) {
      if (colour) {
        message("For a multivariate time series, specify a seriesname for each time series. Defaulting to column names.")
      }
      series <- colnames(object)
    }
    out <- list()
    for (i in 1:NCOL(object)) {
      cl$series <- series[i]
      out[[i]] <- eval(cl)
    }
    return(out)
  }
}

#' @rdname autoplot.ts
#' @export
autolayer.msts <- function(object, series = NULL, ...) {
  if (NCOL(object) > 1) {
    class(object) <- c("mts", "ts", "matrix")
  }
  else {
    if (is.null(series)) {
      series <- deparse(substitute(series))
    }
    class(object) <- c("ts")
  }
  attr(object, "msts") <- NULL
  autolayer(object, series = series, ...)
}

#' @rdname autoplot.ts
#' @export
autolayer.ts <- function(object, colour=TRUE, series=NULL, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    tsdata <- data.frame(
      timeVal = as.numeric(time(object)),
      series = ifelse(is.null(series), deparse(substitute(object)), series),
      seriesVal = as.numeric(object)
    )
    if (colour) {
      ggplot2::geom_line(ggplot2::aes(x = .data[["timeVal"]], y = .data[["seriesVal"]], group = .data[["series"]], colour = .data[["series"]]), data = tsdata, ..., inherit.aes = FALSE)
    }
    else {
      ggplot2::geom_line(ggplot2::aes(x = .data[["timeVal"]], y = .data[["seriesVal"]], group = .data[["series"]]), data = tsdata, ..., inherit.aes = FALSE)
    }
  }
}

#' @rdname plot.forecast
#' @export
autolayer.forecast <- function(object, series = NULL, PI = TRUE, showgap = TRUE, ...) {
  PI <- PI & !is.null(object$level)
  data <- forecast2plotdf(object, PI = PI, showgap = showgap)
  mapping <- ggplot2::aes(x = .data[["x"]], y = .data[["y"]])
  if (!is.null(object$series)) {
    data[["series"]] <- object$series
  }
  if (!is.null(series)) {
    data[["series"]] <- series
    mapping$colour <- quote(series)
  }
  if (PI) {
    mapping$level <- quote(level)
    mapping$ymin <- quote(ymin)
    mapping$ymax <- quote(ymax)
  }
  geom_forecast(mapping = mapping, data = data, stat = "identity", ..., inherit.aes = FALSE)
}

#' @rdname plot.mforecast
#' @export
autolayer.mforecast <- function(object, series = NULL, PI = TRUE, ...) {
  cl <- match.call()
  cl[[1]] <- quote(autolayer)
  cl$object <- quote(object$forecast[[i]])
  if (!is.null(series)) {
    if (length(series) != length(object$forecast)) {
      series <- names(object$forecast)
    }
  }
  out <- list()
  for (i in 1:length(object$forecast)) {
    cl$series <- series[i]
    out[[i]] <- eval(cl)
  }
  return(out)
}


#' Automatically create a ggplot for time series objects
#'
#' \code{autoplot} takes an object of type \code{ts} or \code{mts} and creates
#' a ggplot object suitable for usage with \code{stat_forecast}.
#'
#' \code{fortify.ts} takes a \code{ts} object and converts it into a data frame
#' (for usage with ggplot2).
#'
#' @param object Object of class \dQuote{\code{ts}} or \dQuote{\code{mts}}.
#' @param series Identifies the time series with a colour, which integrates well
#' with the functionality of \link{geom_forecast}.
#' @param facets If TRUE, multiple time series will be faceted (and unless
#' specified, colour is set to FALSE). If FALSE, each series will be assigned a
#' colour.
#' @param colour If TRUE, the time series will be assigned a colour aesthetic
#' @param model Object of class \dQuote{\code{ts}} to be converted to
#' \dQuote{\code{data.frame}}.
#' @param data Not used (required for \code{\link[ggplot2]{fortify}} method)
#' @param ... Other plotting parameters to affect the plot.
#' @inheritParams plot.forecast
#' @return None. Function produces a ggplot graph.
#' @author Mitchell O'Hara-Wild
#' @seealso \code{\link[stats]{plot.ts}}, \code{\link[ggplot2]{fortify}}
#' @examples
#'
#' library(ggplot2)
#' autoplot(USAccDeaths)
#'
#' lungDeaths <- cbind(mdeaths, fdeaths)
#' autoplot(lungDeaths)
#' autoplot(lungDeaths, facets=TRUE)
#'
#' @export
autoplot.ts <- function(object, series=NULL, xlab = "Time", ylab = deparse(substitute(object)),
                        main = NULL,  ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!is.ts(object)) {
      stop("autoplot.ts requires a ts object, use object=object")
    }

    # Create data frame with time as a column labelled x
    # and time series as a column labelled y.
    data <- data.frame(y = as.numeric(object), x = as.numeric(time(object)))
    if (!is.null(series)) {
      data <- transform(data, series = series)
    }

    # Initialise ggplot object
    p <- ggplot2::ggplot(ggplot2::aes(y = .data[["y"]], x = .data[["x"]]), data = data)

    # Add data
    if (!is.null(series)) {
      p <- p + ggplot2::geom_line(ggplot2::aes(group = .data[["series"]], colour = .data[["series"]]), na.rm = TRUE, ...)
    }
    else {
      p <- p + ggplot2::geom_line(na.rm = TRUE, ...)
    }

    # Add labels
    p <- p + ggAddExtras(xlab = xlab, ylab = ylab, main = main)

    # Make x axis contain only whole numbers (e.g., years)
    p <- p + ggplot2::scale_x_continuous(breaks = ggtsbreaks)
    return(p)
  }
}

#' @rdname autoplot.ts
#' @export
autoplot.mts <- function(object, colour=TRUE, facets=FALSE, xlab = "Time", ylab = deparse(substitute(object)),
                         main = NULL, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (!stats::is.mts(object)) {
      stop("autoplot.mts requires a mts object, use x=object")
    }
    if (NCOL(object) <= 1) {
      return(autoplot.ts(object, ...))
    }

    cn <- colnames(object)
    if (is.null(cn)) {
      cn <- paste("Series", seq_len(NCOL(object)))
    }

    data <- data.frame(
      y = as.numeric(c(object)), x = rep(as.numeric(time(object)), NCOL(object)),
      series = factor(rep(cn, each = NROW(object)), levels = cn)
    )

    # Initialise ggplot object
    mapping <- ggplot2::aes(y = .data[["y"]], x = .data[["x"]], group = .data[["series"]])
    if (colour && (!facets || !missing(colour))) {
      mapping$colour <- quote(series)
    }
    p <- ggplot2::ggplot(mapping, data = data)
    p <- p + ggplot2::geom_line(na.rm = TRUE, ...)
    if (facets) {
      p <- p + ggplot2::facet_grid(series~., scales = "free_y")
    }
    p <- p + ggAddExtras(xlab = xlab, ylab = ylab, main = main)
    return(p)
  }
}

#' @rdname autoplot.ts
#' @export
autoplot.msts <- function(object, ...) {
  sname <- deparse(substitute(object))
  if (NCOL(object) > 1) {
    class(object) <- c("mts", "ts", "matrix")
  }
  else {
    class(object) <- c("ts")
  }
  attr(object, "msts") <- NULL
  autoplot(object, ...) + ggAddExtras(ylab = sname)
}

#' @rdname autoplot.ts
#' @export
fortify.ts <- function(model, data, ...) {
  # Use ggfortify version if it is loaded
  # to prevent cran errors
  if (exists("ggfreqplot")) {
    tsp <- attr(model, which = "tsp")
    dtindex <- time(model)
    if (any(tsp[3] == c(4, 12))) {
      dtindex <- zoo::as.Date.yearmon(dtindex)
    }
    model <- data.frame(Index = dtindex, Data = as.numeric(model))
    return(ggplot2::fortify(model))
  }
  else {
    model <- cbind(x = as.numeric(time(model)), y = as.numeric(model))
    as.data.frame(model)
  }
}

forecast2plotdf <- function(model, data=as.data.frame(model), PI=TRUE, showgap=TRUE, ...) {
  # Time series forecasts
  if (is.element("ts", class(model$mean))) {
    xVals <- as.numeric(time(model$mean)) # x axis is time
  }
  # Cross-sectional forecasts
  else if (!is.null(model[["newdata"]])) {
    xVals <- as.numeric(model[["newdata"]][, 1]) # Only display the first column of newdata, should be generalised.
    if (NCOL(model[["newdata"]]) > 1) {
      message("Note: only extracting first column of data")
    }
  }
  else {
    stop("Could not find forecast x axis")
  }
  Hiloc <- grep("Hi ", names(data))
  Loloc <- grep("Lo ", names(data))
  if (PI && !is.null(model$level)) { # PI
    if (length(Hiloc) == length(Loloc)) {
      if (length(Hiloc) > 0) {
        out <- data.frame(
          x = rep(xVals, length(Hiloc) + 1),
          y = c(rep(NA, NROW(data) * (length(Hiloc))), data[, 1]),
          level = c(as.numeric(rep(gsub("Hi ", "", names(data)[Hiloc]), each = NROW(data))), rep(NA, NROW(data))),
          ymax = c(unlist(data[, Hiloc]), rep(NA, NROW(data))), ymin = c(unlist(data[, Loloc]), rep(NA, NROW(data)))
        )
        numInterval <- length(model$level)
      }
    }
    else {
      warning("missing intervals detected, plotting point predictions only")
      PI <- FALSE
    }
  }
  if (!PI) { # No PI
    out <- data.frame(x = xVals, y = as.numeric(model$mean), level = rep(NA, NROW(model$mean)), ymax = rep(NA, NROW(model$mean)), ymin = rep(NA, NROW(model$mean)))
    numInterval <- 0
  }
  if (!showgap) {
    if (is.null(model$x)) {
      warning("Removing the gap requires historical data, provide this via model$x. Defaulting showgap to TRUE.")
    }
    else {
      intervalGap <- data.frame(
        x = rep(time(model$x)[length(model$x)], numInterval + 1),
        y = c(model$x[length(model$x)], rep(NA, numInterval)),
        level = c(NA, model$level)[seq_along(1:(numInterval + 1))],
        ymax = c(NA, rep(model$x[length(model$x)], numInterval)),
        ymin = c(NA, rep(model$x[length(model$x)], numInterval))
      )
      out <- rbind(intervalGap, out)
    }
  }
  return(out)
}

#' @rdname geom_forecast
#' @export
StatForecast <- ggplot2::ggproto(
  "StatForecast", ggplot2::Stat,
  required_aes = c("x", "y"),

  compute_group = function(data, scales, params, PI=TRUE, showgap=TRUE, series=NULL,
                           h=NULL, level=c(80, 95), fan=FALSE, robust=FALSE, lambda=NULL,
                           find.frequency=FALSE, allow.multiplicative.trend=FALSE, ...) {
    ## TODO: Rewrite
    tspx <- recoverTSP(data$x)
    if (is.null(h)) {
      h <- ifelse(tspx[3] > 1, 2 * tspx[3], 10)
    }
    tsdat <- ts(data = data$y, start = tspx[1], frequency = tspx[3])
    fcast <- forecast(
      tsdat, h = h, level = level, fan = fan, robust = robust,
      lambda = lambda, find.frequency = find.frequency,
      allow.multiplicative.trend = allow.multiplicative.trend
    )

    fcast <- forecast2plotdf(fcast, PI = PI, showgap = showgap)

    # Add ggplot & series information
    extraInfo <- as.list(data[1, !colnames(data) %in% colnames(fcast)])
    extraInfo$`_data` <- quote(fcast)
    if (!is.null(series)) {
      if (data$group[1] > length(series)) {
        message("Recycling series argument, please provide a series name for each time series")
      }
      extraInfo[["series"]] <- series[(abs(data$group[1]) - 1) %% length(series) + 1]
    }
    do.call("transform", extraInfo)
  }
)

#' @rdname geom_forecast
#' @export
GeomForecast <- ggplot2::ggproto(
  "GeomForecast", ggplot2::Geom, # Produces both point forecasts and intervals on graph
  required_aes = c("x", "y"),
  optional_aes = c("ymin", "ymax", "level"),
  default_aes = ggplot2::aes(
    colour = "blue", fill = "grey60", size = .5,
    linetype = 1, weight = 1, alpha = 1, level = NA
  ),
  draw_key = function(data, params, size) {
    lwd <- min(data$size, min(size) / 4)

    # Calculate and set colour
    linecol <- blendHex(data$col, "gray30", 1)
    fillcol <- blendHex(data$col, "#CCCCCC", 0.8)

    grid::grobTree(
      grid::rectGrob(
        width = grid::unit(1, "npc") - grid::unit(lwd, "mm"),
        height = grid::unit(1, "npc") - grid::unit(lwd, "mm"),
        gp = grid::gpar(
          col = fillcol,
          fill = scales::alpha(fillcol, data$alpha),
          lty = data$linetype,
          lwd = lwd * ggplot2::.pt,
          linejoin = "mitre"
        )
      ),
      grid::linesGrob(
        x = c(0, 0.4, 0.6, 1),
        y = c(0.2, 0.6, 0.4, 0.9),
        gp = grid::gpar(
          col = linecol,
          fill = scales::alpha(linecol, data$alpha),
          lty = data$linetype,
          lwd = lwd * ggplot2::.pt,
          linejoin = "mitre"
        )
      )
    )
  },

  handle_na = function(self, data, params) { ## TODO: Consider removing/changing
    data
  },

  draw_group = function(data, panel_scales, coord) {
    data <- if (!is.null(data$level)) {
      split(data, !is.na(data$level))
    } else {
      list(data)
    }

    # Draw forecasted points and intervals
    if (length(data) == 1) { # PI=FALSE
      ggplot2:::ggname(
        "geom_forecast",
        GeomForecastPoint$draw_panel(data[[1]], panel_scales, coord)
      )
    }
    else { # PI=TRUE
      ggplot2:::ggname(
        "geom_forecast",
        grid::addGrob(
          GeomForecastInterval$draw_group(data[[2]], panel_scales, coord),
          GeomForecastPoint$draw_panel(data[[1]], panel_scales, coord)
        )
      )
    }
  }
)

GeomForecastPoint <- ggplot2::ggproto(
  "GeomForecastPoint", GeomForecast, ## Produces only point forecasts
  required_aes = c("x", "y"),

  setup_data = function(data, params) {
    data[!is.na(data$y), ] # Extract only forecast points
  },

  draw_group = function(data, panel_scales, coord) {
    linecol <- blendHex(data$colour[1], "gray30", 1)
    # Compute alpha transparency
    data$alpha <- grDevices::col2rgb(linecol, alpha = TRUE)[4, ] / 255 * data$alpha

    # Select appropriate Geom and set defaults
    if (NROW(data) == 0) { # Blank
      ggplot2::GeomBlank$draw_panel
    }
    else if (NROW(data) == 1) { # Point
      GeomForecastPointGeom <- ggplot2::GeomPoint$draw_panel
      pointpred <- transform(data, fill = NA, colour = linecol, size = 1, shape = 19, stroke = 0.5)
    }
    else { # Line
      GeomForecastPointGeom <- ggplot2::GeomLine$draw_panel
      pointpred <- transform(data, fill = NA, colour = linecol, linewidth = size, size = NULL)
    }

    # Draw forecast points
    ggplot2:::ggname(
      "geom_forecast_point",
      grid::grobTree(GeomForecastPointGeom(pointpred, panel_scales, coord))
    )
  }
)

blendHex <- function(mixcol, seqcol, alpha=1) {
  requireNamespace("methods")
  if (is.na(seqcol)) {
    return(mixcol)
  }

  # transform to hue/lightness/saturation colorspace
  seqcol <- grDevices::col2rgb(seqcol, alpha = TRUE)
  mixcol <- grDevices::col2rgb(mixcol, alpha = TRUE)
  seqcolHLS <- suppressWarnings(methods::coerce(colorspace::RGB(R = seqcol[1, ] / 255, G = seqcol[2, ] / 255, B = seqcol[3, ] / 255), structure(NULL, class = "HLS")))
  mixcolHLS <- suppressWarnings(methods::coerce(colorspace::RGB(R = mixcol[1, ] / 255, G = mixcol[2, ] / 255, B = mixcol[3, ] / 255), structure(NULL, class = "HLS")))

  # copy luminence
  mixcolHLS@coords[, "L"] <- seqcolHLS@coords[, "L"]
  mixcolHLS@coords[, "S"] <- alpha * mixcolHLS@coords[, "S"] + (1 - alpha) * seqcolHLS@coords[, "S"]
  mixcolHex <- suppressWarnings(methods::coerce(mixcolHLS, structure(NULL, class = "RGB")))
  mixcolHex <- colorspace::hex(mixcolHex)
  mixcolHex <- ggplot2::alpha(mixcolHex, mixcol[4, ] / 255)
  return(mixcolHex)
}

GeomForecastInterval <- ggplot2::ggproto(
  "GeomForecastInterval", GeomForecast, ## Produces only forecasts intervals on graph
  required_aes = c("x", "ymin", "ymax"),

  setup_data = function(data, params) {
    data[is.na(data$y), ] # Extract only forecast intervals
  },

  draw_group = function(data, panel_scales, coord) {
    # If level scale from fabletools is not loaded, convert to colour
    if(is.numeric(data$level)){
      leveldiff <- diff(range(data$level))
      if (leveldiff == 0) {
        leveldiff <- 1
      }
      shadeVal <- (data$level - min(data$level)) / leveldiff * 0.2 + 8 / 15
      data$level <- rgb(shadeVal, shadeVal, shadeVal)
    }

    intervalGrobList <- lapply(
      split(data, data$level),
      FUN = function(x) {
        # Calculate colour
        fillcol <- blendHex(x$colour[1], x$level[1], 0.7)
        # Compute alpha transparency
        x$alpha <- grDevices::col2rgb(fillcol, alpha = TRUE)[4, ] / 255 * x$alpha

        # Select appropriate Geom and set defaults
        if (NROW(x) == 0) { # Blank
          ggplot2::GeomBlank$draw_panel
        }
        else if (NROW(x) == 1) { # Linerange
          GeomForecastIntervalGeom <- ggplot2::GeomLinerange$draw_panel
          x <- transform(x, colour = fillcol, fill = NA, linewidth = 1)
        }
        else { # Ribbon
          GeomForecastIntervalGeom <- ggplot2::GeomRibbon$draw_group
          x <- transform(x, colour = NA, fill = fillcol, linewidth = size, size = NULL)
        }
        # Create grob
        return(GeomForecastIntervalGeom(x, panel_scales, coord)) ## Create list pair with average ymin/ymax to order layers
      }
    )

    # Draw forecast intervals
    ggplot2:::ggname("geom_forecast_interval", do.call(grid::grobTree, rev(intervalGrobList))) # TODO: Find reliable method to stacking them correctly
  }
)



#' Forecast plot
#'
#' Generates forecasts from \code{forecast.ts} and adds them to the plot.
#' Forecasts can be modified via sending forecast specific arguments above.
#'
#' Multivariate forecasting is supported by having each time series on a
#' different group.
#'
#' You can also pass \code{geom_forecast} a \code{forecast} object to add it to
#' the plot.
#'
#' The aesthetics required for the forecasting to work includes forecast
#' observations on the y axis, and the \code{time} of the observations on the x
#' axis. Refer to the examples below. To automatically set up aesthetics, use
#' \code{autoplot}.
#'
#' @inheritParams ggplot2::layer
#' @param data The data to be displayed in this layer. There are three options:
#'
#' If \code{NULL}, the default, the data is inherited from the plot data as
#' specified in the call to \code{\link[ggplot2]{ggplot}}.
#'
#' A \code{data.frame}, or other object, will override the plot data. All
#' objects will be fortified to produce a data frame. See \code{\link[ggplot2]{fortify}}
#' for which variables will be created.
#'
#' A \code{function} will be called with a single argument, the plot data. The
#' return value must be a \code{data.frame}, and will be used as the layer
#' data.
#' @param stat The stat object to use calculate the data.
#' @param position Position adjustment, either as a string, or the result of a
#' call to a position adjustment function.
#' @param na.rm If \code{FALSE} (the default), removes missing values with a
#' warning.  If \code{TRUE} silently removes missing values.
#' @param show.legend logical. Should this layer be included in the legends?
#' \code{NA}, the default, includes if any aesthetics are mapped. \code{FALSE}
#' never includes, and \code{TRUE} always includes.
#' @param inherit.aes If \code{FALSE}, overrides the default aesthetics, rather
#' than combining with them. This is most useful for helper functions that
#' define both data and aesthetics and shouldn't inherit behaviour from the
#' default plot specification, e.g. \code{\link[ggplot2]{borders}}.
#' @param PI If \code{FALSE}, confidence intervals will not be plotted, giving
#' only the forecast line.
#' @param showgap If \code{showgap=FALSE}, the gap between the historical
#' observations and the forecasts is removed.
#' @param series Matches an unidentified forecast layer with a coloured object
#' on the plot.
#' @param ... Additional arguments for \code{\link{forecast.ts}}, other
#' arguments are passed on to \code{\link[ggplot2]{layer}}. These are often aesthetics,
#' used to set an aesthetic to a fixed value, like \code{color = "red"} or
#' \code{alpha = .5}. They may also be parameters to the paired geom/stat.
#' @return A layer for a ggplot graph.
#' @author Mitchell O'Hara-Wild
#' @seealso \code{\link[generics]{forecast}}, \code{\link[ggplot2]{ggproto}}
#' @examples
#'
#' \dontrun{
#' library(ggplot2)
#' autoplot(USAccDeaths) + geom_forecast()
#'
#' lungDeaths <- cbind(mdeaths, fdeaths)
#' autoplot(lungDeaths) + geom_forecast()
#'
#' # Using fortify.ts
#' p <- ggplot(aes(x=x, y=y), data=USAccDeaths)
#' p <- p + geom_line()
#' p + geom_forecast()
#'
#' # Without fortify.ts
#' data <- data.frame(USAccDeaths=as.numeric(USAccDeaths), time=as.numeric(time(USAccDeaths)))
#' p <- ggplot(aes(x=time, y=USAccDeaths), data=data)
#' p <- p + geom_line()
#' p + geom_forecast()
#'
#' p + geom_forecast(h=60)
#' p <- ggplot(aes(x=time, y=USAccDeaths), data=data)
#' p + geom_forecast(level=c(70,98))
#' p + geom_forecast(level=c(70,98),colour="lightblue")
#'
#' #Add forecasts to multivariate series with colour groups
#' lungDeaths <- cbind(mdeaths, fdeaths)
#' autoplot(lungDeaths) + geom_forecast(forecast(mdeaths), series="mdeaths")
#' }
#'
#' @export
geom_forecast <- function(mapping = NULL, data = NULL, stat = "forecast",
                          position = "identity", na.rm = FALSE, show.legend = NA,
                          inherit.aes = TRUE, PI=TRUE, showgap=TRUE, series=NULL, ...) {
  if (is.forecast(mapping) || is.mforecast(mapping)) {
    warning("Use autolayer instead of geom_forecast to add a forecast layer to your ggplot object.")
    cl <- match.call()
    cl[[1]] <- quote(autolayer)
    names(cl)[names(cl) == "mapping"] <- "object"
    return(eval.parent(cl))
  }
  if (is.ts(mapping)) {
    data <- data.frame(y = as.numeric(mapping), x = as.numeric(time(mapping)))
    mapping <- ggplot2::aes(y = .data[["y"]], x = .data[["x"]])
  }
  if (stat == "forecast") {
    paramlist <- list(na.rm = na.rm, PI = PI, showgap = showgap, series = series, ...)
    if (!is.null(series)) {
      if (inherits(mapping, "uneval")) {
        mapping$colour <- quote(ggplot2::after_stat(series))
      }
      else {
        mapping <- ggplot2::aes(colour = ggplot2::after_stat(series))
      }
    }
  }
  else {
    paramlist <- list(na.rm = na.rm, ...)
  }
  ggplot2::layer(
    geom = GeomForecast, mapping = mapping, data = data, stat = stat,
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = paramlist
  )
}

# Produce nice histogram with appropriately chosen bin widths
# Designed to work with time series data without issuing warnings.


#' Histogram with optional normal and kernel density functions
#'
#' Plots a histogram and density estimates using ggplot.
#'
#'
#' @param x a numerical vector.
#' @param add.normal Add a normal density function for comparison
#' @param add.kde Add a kernel density estimate for comparison
#' @param add.rug Add a rug plot on the horizontal axis
#' @param bins The number of bins to use for the histogram. Selected by default
#' using the Friedman-Diaconis rule given by \code{\link[grDevices]{nclass.FD}}
#' @param boundary A boundary between two bins.
#' @return None.
#' @author Rob J Hyndman
#' @seealso \code{\link[graphics]{hist}}, \code{\link[ggplot2]{geom_histogram}}
#' @examples
#'
#' gghistogram(lynx, add.kde=TRUE)
#'
#' @export
gghistogram <- function(x, add.normal=FALSE, add.kde=FALSE, add.rug=TRUE, bins, boundary=0) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
  }
  else {
    if (missing(bins)) {
      bins <- min(500, grDevices::nclass.FD(na.exclude(x)))
    }
    data <- data.frame(x = as.numeric(c(x)))
    # Initialise ggplot object and plot histogram
    binwidth <- (max(x, na.rm = TRUE) - min(x, na.rm = TRUE)) / bins
    p <- ggplot2::ggplot() +
      ggplot2::geom_histogram(ggplot2::aes(x), data = data, binwidth = binwidth, boundary = boundary) +
      ggplot2::xlab(deparse(substitute(x)))
    # Add normal density estimate
    if (add.normal || add.kde) {
      xmin <- min(x, na.rm = TRUE)
      xmax <- max(x, na.rm = TRUE)
      if (add.kde) {
        h <- stats::bw.SJ(x)
        xmin <- xmin - 3 * h
        xmax <- xmax + 3 * h
      }
      if (add.normal) {
        xmean <- mean(x, na.rm = TRUE)
        xsd <- sd(x, na.rm = TRUE)
        xmin <- min(xmin, xmean - 3 * xsd)
        xmax <- max(xmax, xmean + 3 * xsd)
      }
      xgrid <- seq(xmin, xmax, length.out = 512)
      if (add.normal) {
        df <- data.frame(x = xgrid, y = length(x) * binwidth * stats::dnorm(xgrid, xmean, xsd))
        p <- p + ggplot2::geom_line(ggplot2::aes(df$x, df$y), col = "#ff8a62")
      }
      if (add.kde) {
        kde <- stats::density(x, bw = h, from = xgrid[1], to = xgrid[512], n = 512)
        p <- p + ggplot2::geom_line(ggplot2::aes(x = kde$x, y = length(x) * binwidth * kde$y), col = "#67a9ff")
      }
    }
    if (add.rug) {
      p <- p + ggplot2::geom_rug(ggplot2::aes(x))
    }
    return(p)
  }
}

Try the forecast package in your browser

Any scripts or data that you put into this service are public.

forecast documentation built on June 22, 2024, 9:20 a.m.