R/qqplot.R

In eFCM: Exponential Factor Copula Model

#' Q–Q Plot for Fitted Factor Copula Model
#'
#' Produce a Q–Q plot comparing empirical exceedances to the fitted eFCM tail,
#' with an optional GPD tail overlay for diagnostic comparison.
#'
#' The function first selects a threshold \code{u} as the empirical
#' \code{thres}-quantile of the selected station series \code{x}.
#' It then forms exceedances \eqn{Y = X - u \mid X>u}, fits the eFCM
#' (implicitly via \code{qfcm()} and a scalar \eqn{\lambda} estimate),
#' and plots empirical exceedances against theoretical eFCM quantiles in the tail.
#' If \code{gpd=TRUE}, a GPD is fitted to the exceedances (threshold 0) and
#' its theoretical tail quantiles are added for visual comparison.
#'
#' @param object An object of class \code{"fcm"} returned by [fcm()].
#' @param which Integer scalar. Station (column) index to plot.
#' @param gpd Logical; if \code{TRUE}, add a GPD tail fit to the Q–Q plot.
#' @param thres Numeric in \code{(0,1)}; the probability threshold used to pick
#'   the empirical quantile \code{u = quantile(x, thres)}. Defaults to \code{0.9}.
#' @param main,xlab,ylab Character. Graphical labels passed to [plot()].
#' @param ... Additional graphical arguments forwarded to [plot()].
#' @return A numeric vector of fitted eFCM theoretical tail quantiles, invisibly returned.
#' @export
qqplot <- function(object, ...) {
  UseMethod("qqplot")
}

#' @rdname qqplot
#' @method qqplot fcm
#' @importFrom ismev gpd.fit
#' @importFrom ismev gpdq
#' @export
qqplot.fcm <- function(object, which = 1, gpd = TRUE, thres = 0.9,
                       main = "Q-Q plot", xlab = "Theoretical quantiles (exceedances)",
                       ylab = "Empirical exceedances", ...) {
  stopifnot(inherits(object, "fcm"))
  stopifnot(is.numeric(which), length(which) == 1)

  if (which > ncol(object$data) || which < 1) {
    stop("Station index out of bounds.")
  }

  x <- sort(object$data[, which])
  emp_u <- u_np(x)
  u <- quantile(x, thres)
  y <- sort(x[x > u] - u)
  n <- length(y)

  # Fit eFCM
  lambda <- optimize(L1,
                   lower = 0.5 * min(1 / max(x), min(x)),
                   upper = 2 * min(1 / min(x), 2),
                   w = x)$minimum
  fitted_fcm <- sapply(emp_u, qfcm, lambda = lambda)

  plot(tail(fitted_fcm, n) - u, y, main = main, xlab = xlab, ylab = ylab,
       pch = 16, col = "black", ...)
  abline(0, 1, col = "gray", lty = 2)

  # Optional GPD fit
  if (gpd) {
    fitgpd <- tryCatch(
      gpd.fit(y, threshold = 0, show = FALSE),
      error = function(e) NULL
    )
    if (!is.null(fitgpd)) {
      qth <- gpdq(fitgpd$mle, 0, 1 - (1:n) / (n + 1))
      points(qth, y, col = "blue", pch = 16, lwd = 2)
      legend("topleft", legend = c("eFCM", "GPD"),
             col = c("black", "blue"), pch = c(16, 16), lwd = c(2, 2))
    }
  } else {
    legend("topleft", legend = "eFCM", col = "black", pch = 16)
  }
  invisible(fitted_fcm)
}