Fix Krupskii et Joe method (still numerically unstable)

R/rmev.R

@@ -114,7 +114,7 @@ rmev <- function(
     grid = FALSE,
     dist = NULL,
     ...) {
-
+ model <- match.arg(model)
   # Create gridded values if specification is for random field discretization
   ellips <- list(...)
   if(is.null(coord) && !is.null(ellips$loc)){
@@ -239,6 +239,7 @@ rmevspec <- function(n, d, param, sigma, model = c("log", "neglog", "bilog", "ne
     weights = NULL, vario = NULL, coord = NULL, grid = FALSE, dist = NULL, ...) {
   # Dummy algorithm argument for internal checks
   alg <- "sm"
+   model <- match.arg(model)
   if(model %in% c("alog","aneglog", "maxlik")){
     stop("Invalid model: cannot simulate from angular distribution of asymmetric models.")
   }

R/taildep2.R

@@ -3,7 +3,7 @@
 #' Estimators proposed by Krupskii and Joe under second order expansion
 #' for the coefficient of tail dependence \eqn{\eta} and the
 #' joint tail orthant probability
-#' 
+#'
 #' @param data a matrix of observations
 #' @param q vector of quantile levels
 #' @param ptail tail probability smaller than \code{q}. Default to \code{NULL}
@@ -12,17 +12,25 @@
 #' @param type integer indicating the estimator type
 #' @return a list with elements
 #' \itemize{
-#' \item \code{p} probability
+#' \item \code{p} quantile level for estimation
 #' \item \code{eta} estimated coefficient of tail dependence \eqn{\eta}
 #' \item \code{eta_sd} estimated standard error of \eqn{\eta}
+#' \item \code{k1} parameter of the tail expansion
+#' \item \code{pat} proportion of observations above the threshold
+#' \item \code{lambda} tail dependence coefficient (sic)
+#' \item \code{tailprob} tail probability, if \code{ptail} is provided
 #' }
+#'
+#' @note The numerical optimization of the likelihood surface is difficult, as the function is ill-behaved. VIsual inspection of estimates is necessary to check for non-convergence.
 #' @examples
-#' rho <- runif(1, -1, 1)
 #' d <- 2
+#' rho <- 0.9
+#' Sigma <- matrix(rho, d, d) + diag(1 - rho, d)
+#' eta_true <- 1/sum(Sigma)
 #' data <- mev::mvrnorm(
-#'    n = 1e5,
+#'    n = 1e4,
 #'    mu = rep(0, d),
-#'   Sigma = matrix(rho, d, d) + diag(1 - rho, d))
+#'   Sigma = Sigma)
 #' q <- seq(0.95, 0.995, by = 0.005)
 #' taildep <- kjtail(data = data, q = q)
 #' with(taildep,
@@ -66,30 +74,30 @@ kjtail <- function(
   unif <- mev::spunif(x = data, thresh = qulvl)
   }
   # Consider minimum of tail
-  Ud <- apply(unif, 1, function(x){max(1-x)})
+  Ud <- 1-apply(unif, 1, min)
   q <- sort(q, decreasing = TRUE)
   etahat <- etasd <- k1hat <- nabove <- jtprob <- plev <- rep(NA_real_, length(q))
   neglik_fn <- function(par, xdat, Tp){
     k1 <- par[1]
     eta <- par[2]
-    - sum(log(pmax(0, k1+(1-k1*Tp)/eta*exp((1/eta-1)*log(xdat) - log(Tp)/eta))))
+    return(- sum(log(pmax(1e-20, k1 + (1 - k1 * Tp) / eta * exp((1/eta-1) * log(xdat) - log(Tp)/eta)))))
   }
   for(j in seq_along(q)){
     plev[j] <- Tp <- quantile(Ud, q[j])
     sUd <- Ud[Ud < Tp]
     nabove[j] <- length(sUd)
     if(j == 1L){
-      start <- c(1/Tp, 1/(d-0.5))
+      start <- c(1, 1/(d-0.5))
     } else{
       start <- est$pars
     }
-  est <- Rsolnp::solnp(
-     pars = c(1/Tp, 1/(d-0.5)),
+    est <- Rsolnp::solnp(
+     pars = c(1, 1/(d-0.5)),
      fun = neglik_fn,
      ineqfun = function(par, xdat, Tp){
        # Probability must be positive
        # k1 must be positive?
-      par[1]*range(xdat)+(1-par[1]*Tp)*exp((log(range(xdat)) - log(Tp)/par[2]))
+      range(par[1]*xdat+(1-par[1]*Tp)*exp((log(xdat) - log(Tp))/par[2]))
     },
     ineqLB = rep(0, 2),
     ineqUB = rep(1, 2),
@@ -98,11 +106,13 @@ kjtail <- function(
     Tp = Tp,
     xdat = sUd,
     control = list(trace = 0))
-  est$hessian <- numDeriv::hessian(func = neglik_fn, x = est$pars, Tp = Tp, xdat = sUd)
+  # browser()
+  est$hessian <- est$hessian[3:4,3:4]
+  # numDeriv::hessian(func = neglik_fn, x = est$pars, Tp = Tp, xdat = sUd)
     if(isTRUE(est$convergence == 0)){
      k1hat[j] <- est$pars[1]
      etahat[j] <- est$pars[2]
-     etasd_j <- try(silent = TRUE, suppressWarnings(sqrt(solve(est$hessian[1:2,1:2])[2,2])))
+     etasd_j <- try(silent = TRUE, suppressWarnings(sqrt(solve(est$hessian)[2,2])))
      if(!inherits(etasd_j, "try-error")){
        etasd[j] <- etasd_j
      }
@@ -113,6 +123,7 @@ kjtail <- function(
     p = 1-as.numeric(plev),
     eta = as.numeric(etahat),
     eta_sd = etasd,
+    k1 = k1hat,
     pat = nabove/n,
     lambda = as.numeric(lambdahat))
   if(!is.null(ptail)){