extremo: Pairwise extremogram for max-risk functional

In mev: Modelling of Extreme Values

Pairwise extremogram for max-risk functional

Description

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

## Not run: 
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))))

## End(Not run)

mev documentation built on Sept. 11, 2024, 8:14 p.m.