R/spectral_likelihood.R

In mvPot: Multivariate Peaks-over-Threshold Modelling for Spatial Extreme Events

#' Spectral log-likelihood function
#'
#' Compute the negative spectral log-likelihood function for Brown--Resnick model with peaks-over-threshold.
#'
#' The function compute the negative log-likelihood function based on the spectral representation developed
#' by Engelke et al. (2015). This simplified expression is obtained by conditioning on the event
#' `\code{sum(x)} exceeds a high threshold \code{u > 1}'. Margins must have been standardized.
#'
#' @param obs List of observations vectors for which \code{sum(x)} exceeds a high threshold.
#' @param loc Matrix of coordinates as given by \code{expand.grid()}.
#' @param vario Semi-variogram function taking a vector of coordinates as input.
#' @param nCores Number of cores used for the computation
#' @param cl Cluster instance as created by \code{makeCluster} of the \code{parallel} package.
#' @return Negative spectral log-likelihood function evaluated at the set of observations \code{obs} with semi-variogram \code{vario}.
#' @examples
#' #Define semi-variogram function
#' vario <- function(h){
#'    1 / 2 * norm(h,type = "2")^1.5
#' }
#'
#' #Define locations
#' loc <- expand.grid(1:4, 1:4)
#'
#' #Simulate data
#' obs <- simulPareto(1000, loc, vario)
#'
#' #Evaluate risk functional
#' sums <- sapply(obs, sum)
#'
#' #Select exceedances
#' exceedances <- obs[sums > quantile(sums, 0.9)]
#'
#' #Evaluate the spectral function
#' spectralLikelihood(exceedances, loc, vario)
#' @export
#' @references Engelke, S. et al. (2015). Estimation of Huesler-Reiss distributions and Brown-Resnick processes. Journal of the Royal Statistical Society: Series B, 77(1), 239-265

spectralLikelihood <- function(obs, loc, vario, nCores = 1L, cl = NULL){
  if(is.matrix(obs)){ #Not converted to list
    obs <- split(obs, row(obs)) #create list
  }
  if(!inherits(obs, "list") || length(obs) < 1 || !inherits(obs[[1]], c("numeric","integer"))){
    stop('obs must be a list of vectors')
  }
  if(!inherits(loc, c("matrix", "data.frame"))) {
    stop('`loc` must be a data frame of coordinates as generated by `expand.grid()` or a matrix of locations (one site per row)')
  }
  if(!is.numeric(nCores) || nCores < 1) {
    stop('`nCores` must a positive number of cores to use for parallel computing.')
  }
  if(nCores > 1 && !inherits(cl, "cluster")) {
    stop('For parallel computation, `cl` must an cluster created by `makeCluster` of the package `parallel`.')
  }

  n <- length(obs)
  dim <- nrow(loc)

  if(dim != length(obs[[1]])){
    stop('The size of the vectors of observations does not match the grid size.')
  }

  gamma <- tryCatch({
    dists <- lapply(1:ncol(loc), function(i) {
      outer(loc[,i],loc[,i], "-")
    })

    computeVarMat <- sapply(1:length(dists[[1]]), function(i){
      h <- rep(0,ncol(loc))
      for(j in 1:ncol(loc)){
        h[j] = dists[[j]][i]
      }
      vario(h)
    })
    matrix(computeVarMat, dim, dim)
  }, warning = function(war) {
    war
  }, error = function(err) {
    stop('The semi-variogram is not valid for the locations provided.')
  })

  psi <- (outer(gamma[-1,1], gamma[-1,1], "+") - (gamma[-1,-1]))
  invPsi <- MASS::ginv(psi)

  productCovMatrix = function(i){
    omega <- log(obs[[i]][-1]/obs[[i]][1]) + gamma[-1,1]
    res <- t(omega) %*% invPsi %*% omega
  }

  if(nCores > 1){
    likelihood <- parallel::parLapply(cl,1:n, productCovMatrix)
  } else {
    likelihood <- lapply(1:n, productCovMatrix)
  }

  logdetA = determinant(psi, logarithm = TRUE)$modulus

  (n/2 * logdetA + 1/2 * sum(unlist(likelihood)))
}