R/xyplot.Arima.R

In tactile: New and Extended Plots, Methods, and Panel Functions for 'lattice'

#' Diagnostic Plots for ARIMA Models
#'
#' Diagnostic plots modelled after [stats::tsdiag()] with some modifications
#' and corrections of p-values in the Box--Ljung test.
#'
#' @param x A fitted time-series model of class `Arima`.
#' @param data Ignored
#' @param which A sequence of integers between 1 and 4, identifying the
#'   plots to be shown.
#' @param na.action Treatment of `NA`s.
#' @param main Optional titles for the plots. Can also be `TRUE`, in which
#'   case a default list of titles will be added.
#' @param lag.max Number of lags to compute ACF for.
#' @param qq.aspect Aspect of Q-Q plot (see [lattice::qqmath()]).
#' @param gof.lag The maximum number of lags for the Ljung--Box test.
#' @param layout Either a numeric vector with (columns, rows) to use in the call
#'   to [gridExtra::grid.arrange()], or a layout matrix which will then be
#'   passed as the `layout_matrix` in `grid.arrange()`.
#' @param ... Parameters to pass to [lattice::xyplot()].
#'
#' @seealso [stats::tsdiag()], [stats::arima()], [lattice::xyplot()],
#'   [gridExtra::grid.arrange()], [stats::Box.test()].
#' @return Plots a lattice plot and returns a `trellis` object.
#' @export
#'
#' @examples
#' fit <- arima(lh, order = c(1, 1, 0))
#' xyplot(fit, layout = c(2, 2))
#' xyplot(fit, which = c(1:2, 4), layout = rbind(c(1, 1), c(2, 3)))
xyplot.Arima <- function(
  x,
  data = NULL,
  which = 1:4,
  lag.max = NULL,
  gof.lag = NULL,
  qq.aspect = "iso",
  na.action = na.pass,
  main = NULL,
  layout = NULL,
  ...
) {
  # Assertions
  stopifnot(
    length(which) > 0,
    length(which) < 5,
    all(which > 0),
    all(which < 5)
  )

  # Warnings
  if (!is.null(data)) {
    warning("'data' is ignored.")
  }

  show <- rep.int(FALSE, 4L)
  show[which] <- TRUE

  titles <- rep.int("", 4L)

  if (is.character(main)) {
    stopifnot(length(main) == sum(show))
    titles[which] <- main
  } else if (isTRUE(main)) {
    titles <- c(
      "Standardized residuals over time",
      "Q-Q Plot of Standardized residuals",
      "ACF of residuals",
      "Ljung-Box tests"
    )
  }

  plots <- vector("list", 4L)
  dots <- list(...)
  add.line <- trellis.par.get("add.line")

  na.action <- get_fun(na.action)

  r <- na.action(x$residuals)
  rstd <- r / sqrt(x$sigma2)

  # Standardized residuals
  if (show[1L]) {
    plots[[1L]] <- do.call(
      xyplot,
      updateList(
        list(
          rstd,
          main = main[1L],
          ylab = "Standardized residuals",
          panel = function(...) {
            panel.abline(h = 0L, reference = TRUE)
            panel.xyplot(...)
          }
        ),
        dots
      )
    )
  }

  # Q-Q-diagram of standardized residuals
  if (show[2L]) {
    plots[[2L]] <- do.call(
      qqmath,
      updateList(
        list(
          ~rstd,
          xlab = "Theoretical quantiles",
          ylab = "Standardized residuals",
          main = main[2L],
          aspect = qq.aspect,
          prepanel = prepanel.qqmathline,
          panel = function(...) {
            panel.qqmathci(...)
            panel.qqmathline(...)
            panel.qqmath(...)
          }
        ),
        dots
      )
    )
  }

  if (any(show[3L:4L])) {
    df <- sum(x$arma[c(1, 3, 4, 7)], na.rm = TRUE)
    period <- x$arma[5]
  }

  # ACF of residuals
  if (show[3L]) {
    if (is.null(lag.max)) {
      lag.max <- if (period < 6) 20 else 3 * period
      if (lag.max <= df + 8) {
        lag.max <- df + 8
      }
    }
    autocor <- acf(r, na.action = na.action, lag.max = lag.max, plot = FALSE)

    plots[[3L]] <- do.call(
      xyplot,
      updateList(
        list(
          x = autocor,
          main = main[3L]
        ),
        dots
      )
    )
  }

  # Box-Ljung p.tests
  if (show[4L]) {
    if (is.null(gof.lag)) {
      gof.lag <- if (period < 6) 20 else 3 * period
      if (gof.lag <= df + 8) {
        gof.lag <- df + 8
      }
    } else if (gof.lag < df) {
      stop("'gof.lag' cannot be < df (p + q + P + Q)")
    }

    pval <- double(gof.lag)
    pval[] <- NA
    for (i in (df + 1):gof.lag) {
      pval[i] <- Box.test(r, i, "Ljung-Box", df)$p.value
    }

    plots[[4L]] <- do.call(
      xyplot,
      updateList(
        list(
          x = pval ~ seq_along(pval),
          xlab = "Lag",
          ylab = "Ljung-Box p-values",
          main = main[4L],
          ylim = range(c(0, max(pval, na.rm = TRUE) * 1.08, 0.1)),
          panel = function(x, y, ...) {
            panel.abline(h = 0.05, lty = 2)
            panel.xyplot(x, y, ...)
          }
        ),
        dots
      )
    )
  }

  grid_wrap(plots, layout = layout)
}