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))))
}
Run the code above in your browser using DataLab