Skip to contents

Pairwise extremogram for max-risk functional

Source: R/extremogram.R
extremo.Rd

The function computes the pairwise \(chi\) estimates and plots them as a function of the distance between sites.

Usage

extremo(dat, margp, coord, scale = 1, rho = 0, plot = FALSE, ...)

Arguments

dat

data matrix

margp

marginal probability above which to threshold observations

coord

matrix of coordinates (one site per row)

scale

geometric anisotropy scale parameter

rho

geometric anisotropy angle parameter

plot

logical; should a graph of the pairwise estimates against distance? Default to FALSE

...

additional arguments passed to plot

Value

an invisible matrix with pairwise estimates of chi along with distance (unsorted)

Examples

if (FALSE) {
lon <- seq(650, 720, length = 10)
lat <- seq(215, 290, length = 10)
# Create a grid
grid <- expand.grid(lon,lat)
coord <- as.matrix(grid)
dianiso <- distg(coord, 1.5, 0.5)
sgrid <- scale(grid, scale = FALSE)
# Specify marginal parameters `loc` and `scale` over grid
eta <- 26 + 0.05*sgrid[,1] - 0.16*sgrid[,2]
tau <- 9 + 0.05*sgrid[,1] - 0.04*sgrid[,2]
# Parameter matrix of Huesler--Reiss
# associated to power variogram
Lambda <- ((dianiso/30)^0.7)/4
# Regular Euclidean distance between sites
di <- distg(coord, 1, 0)
# Simulate generalized max-Pareto field
set.seed(345)
simu1 <- rgparp(n = 1000, thresh = 50, shape = 0.1, riskf = "max",
                scale = tau, loc = eta, sigma = Lambda, model = "hr")
extdat <- extremo(dat = simu1, margp = 0.98, coord = coord,
                  scale = 1.5, rho = 0.5, plot = TRUE)

# Constrained optimization
# Minimize distance between extremal coefficient from fitted variogram
mindistpvario <- function(par, emp, coord){
alpha <- par[1]; if(!isTRUE(all(alpha > 0, alpha < 2))){return(1e10)}
scale <- par[2]; if(scale <= 0){return(1e10)}
a <- par[3]; if(a<1){return(1e10)}
rho <- par[4]; if(abs(rho) >= pi/2){return(1e10)}
semivariomat <- power.vario(distg(coord, a, rho), alpha = alpha, scale = scale)
  sum((2*(1-pnorm(sqrt(semivariomat[lower.tri(semivariomat)]/2))) - emp)^2)
}

hin <- function(par, ...){
  c(1.99-par[1], -1e-5 + par[1],
    -1e-5 + par[2],
    par[3]-1,
    pi/2 - par[4],
    par[4]+pi/2)
  }
opt <- alabama::auglag(par = c(0.7, 30, 1, 0),
                       hin = hin,
                        fn = function(par){
                          mindistpvario(par, emp = extdat[,'prob'], coord = coord)})
stopifnot(opt$kkt1, opt$kkt2)
# Plotting the extremogram in the deformed space
distfa <- distg(loc = coord, opt$par[3], opt$par[4])
plot(
 x = c(distfa[lower.tri(distfa)]), 
 y = extdat[,2], 
 pch = 20,
 yaxs = "i", 
 xaxs = "i", 
 bty = 'l',
 xlab = "distance", 
 ylab= "cond. prob. of exceedance", 
 ylim = c(0,1))
lines(
  x = (distvec <- seq(0,200, length = 1000)), 
  col = 2, lwd = 2,
  y = 2*(1-pnorm(sqrt(power.vario(distvec, alpha = opt$par[1],
                               scale = opt$par[2])/2))))
}