R/chiplot.R

In eFCM: Exponential Factor Copula Model

#' Chi Plot for Fitted eFCM Model
#'
#' Plots the eFCM conditional exceedance probability \eqn{\chi_h(u)}.
#'
#' @method chiplot fcm
#' @export
#'
#' @param object An object of class `"fcm"` returned by [fcm()].
#' @param h A positive numeric distance in kilometers. If `NULL` and `emp = TRUE`, inferred from `coord[which, ]`.
#' @param method Character. Method for computing confidence intervals. One of `"default"`, `"hessian"`, or `"boot"`.
#' @param ci Confidence level for interval estimation.
#' @param emp Logical. Whether to add empirical chi estimates.
#' @param which Integer vector of length 2. Locations to compute empirical chi.
#' @param ... Further arguments passed to base plotting functions (e.g., `main`, `xlab`, `ylab`, etc.).
#'
#' @return A (invisible) list containing chi estimates and confidence bounds:
#' \describe{
#'   \item{chi.u}{Estimated chi values.}
#'   \item{chi.lower}{Lower confidence bounds (if applicable).}
#'   \item{chi.upper}{Upper confidence bounds (if applicable).}
#'   \item{chi.emp}{Empirical chi curve (if `emp = TRUE`).}
#'   \item{h}{Distance used.}
#' }
#'
#' @seealso [chi()], [Echi()]
#'
#' @references
#' Castro-Camilo, D. and Huser, R. (2020).
#' Local likelihood estimation of complex tail dependence structures,
#' applied to US precipitation extremes.
#' *Journal of the American Statistical Association*, 115(531), 1037–1054.
#'
#' @srrstats {G1.5} *Include all code necessary to reproduce performance claims.*
#' @srrstats {G1.6} *Include code to compare with alternative implementations.*
#' @srrstats {G2.1, G2.2, G2.3a, G2.10}
#' @rdname chiplot
chiplot.fcm <- function(object, h = NULL,
                        method = c("default", "hessian", "boot"),
                        ci = 0.95, emp = TRUE, which = c(1, 2),...) {

  stopifnot(inherits(object, "fcm"))
  stopifnot(length(which) == 2)

  # Default distance from coordinates
  if (emp && is.null(h)) {
    h <- rdist.earth(object$coord[which, ], miles = FALSE)[1, 2]
  }
  if (is.null(h) || !is.numeric(h) || length(h) != 1 || h <= 0) {
    stop("`h` must be a positive numeric scalar.")
  }

  method <- match.arg(method)
  u <- seq(0.5, 0.999, 0.001)
  out <- list(h = h)

  if (method == "hessian") {
    # 1) 拟合点和参数提取
    theta <- coef(object, ci = ci, method = method)[,1]
    lambda  <- theta[1]
    delta <- theta[2]

    # 2) 计算协方差矩阵 Cov_theta
    H <- object$hessian
    ev <- eigen(H, symmetric = TRUE, only.values = TRUE)$values
    # 如果 Hessian 是对数似然的 Hessian（多数为负），取 -H
    if (all(ev < 0)) {
      Cov_theta <- solve(-H)
    } else {
      Cov_theta <- solve( H)
    }

    grad_mat <- t(sapply(u, function(ui) {
      numDeriv::grad(function(p) {
        chi_fcm(ui, h = h,
                lambda  = p[1],
                delta = p[2],
                nu    = object$arg$nu)
      }, x = c(lambda, delta))
    }))    # 结果是 length(u) × 2 矩阵

    # 4) 计算标准误
    se <- apply(grad_mat, 1, function(g) {
      sqrt(as.numeric( t(g) %*% Cov_theta %*% g ))
    })

    # 5) 计算点估计 & CI
    chi.u <- sapply(u, function(ui) {
      chi_fcm(ui, h = h,
              lambda  = lambda,
              delta = delta,
              nu    = object$arg$nu)
    })
    zcrit     <- qnorm(1 - (1 - ci)/2)
    chi.lower <- chi.u - zcrit * se
    chi.upper <- chi.u + zcrit * se

    plot(u, chi.u, type = "l", ylab = expression(chi(u)), xlab = "u", ...)
    lines(u, chi.upper, col = "red", lty = 2)
    lines(u, chi.lower, col = "red", lty = 2)

    out$chi.u <- chi.u
    out$chi.lower <- chi.lower
    out$chi.upper <- chi.upper

  } else if (method == "boot") {
    theta <- coef(object, method = method)[, 1]
    chiv_mat <- chiv(u, object$boot, h, object$arg$nu)
    chi.u <- sapply(u, function(ui) {
      chi_fcm(ui, h = h, lambda  = theta[1], delta = theta[2], nu = object$arg$nu)
    })
    chi.lower <- apply(chiv_mat, 1, quantile, (1 - ci) / 2, na.rm = TRUE)
    chi.upper <- apply(chiv_mat, 1, quantile, 1 - (1 - ci) / 2, na.rm = TRUE)

    plot(u, chi.u, type = "l", ylab = expression(chi(u)), xlab = "u", ...)
    lines(u, chi.upper, col = "red", lty = 2)
    lines(u, chi.lower, col = "red", lty = 2)

    out$chi.u <- chi.u
    out$chi.lower <- chi.lower
    out$chi.upper <- chi.upper

  } else {
    theta <- coef(object)
    lambda <- theta[1]
    delta <- theta[2]
    chi.u <- sapply(u, function(x) chi_fcm(x, h = h, lambda = lambda, delta = delta, nu = object$arg$nu))
    plot(u, chi.u, type = "l", ylab = expression(chi(u)), xlab = "u", ...)
    out$chi.u <- chi.u
  }

  # Empirical estimator
  if (emp) {
    chi.emp <- sapply(u, Echi, object = object, which = which)
    lines(u, chi.emp, col = "blue", lty = 3)
    out$chi.emp <- chi.emp
  }

  invisible(out)
}
#' @rdname chiplot
#' @export
chiplot <- function(object, ...) {
  UseMethod("chiplot")
}