ExtremalDep (version 1.0.0)

ExtQ: Univariate Extreme Quantile

Description

Computes the extreme-quantiles of a univariate random variable corresponding to some exceedance probabilities.

Usage

ExtQ(
  P = NULL,
  method = "Frequentist",
  pU = NULL,
  cov = NULL,
  param = NULL,
  param_post = NULL
)

Value

When method == "frequentist", the function returns a vector of length length(P) if ncol(cov) = 1 (constant mean) or a (length(P) x nrow(cov)) matrix if ncol(cov) > 1.

When method == "bayesian", the function returns a

(length(param_post) x length(P)) matrix if ncol(cov) = 1

or a list of ncol(cov) elements each taking a

(length(param_post) x length(P)) matrix if ncol(cov) > 1.

Arguments

P

A vector with values in \([0,1]\) indicating the probabilities of the quantiles to be computed.

method

A character string indicating the estimation method. Takes value "bayesian" or "frequentist".

pU

A value in \([0,1]\) indicating the probability of exceeding a high threshold. In the estimation procedure, observations below the threshold are censored.

cov

A \(q \times c\) matrix indicating \(q\) observations of \(c-1\) covariates for the location parameter.

param

A \((c + 2)\) vector indicating the estimated parameters. Required when method="Frequentist".

param_post

A \(n \times (c + 2)\) matrix indicating the posterior sample for the parameters, where \(n\) is the number of MCMC replicates after removal of the burn-in period. Required when method="Bayesian".

Author

Simone Padoan, , https://faculty.unibocconi.it/simonepadoan/; Boris Beranger, https://www.borisberanger.com

Details

The first column of cov is a vector of 1s corresponding to the intercept.

When pU is NULL (default), then it is assumed that a block maxima approach was taken and quantiles are computed using the qGEV function. When pU is provided, it is assumed that a threshold exceedances approach is taken and the quantiles are computed as $$\mu + \sigma \left(\left(\frac{pU}{P}\right)^\xi - 1\right) \frac{1}{\xi}$$

References

Beranger, B., Padoan, S. A. and Sisson, S. A. (2021). Estimation and uncertainty quantification for extreme quantile regions. Extremes, 24, 349-375.

See Also

fGEV, qGEV

Examples

##################################################
### Example - Pollution levels in Milan, Italy ###
##################################################

if (FALSE) {

data(MilanPollution)

# Frequentist estimation
fit <- fGEV(Milan.winter$PM10)
fit$est

q1 <- ExtQ(
  P = 1 / c(600, 1200, 2400),
  method = "Frequentist",
  param = fit$est
)
q1

# Bayesian estimation with high threshold
cov <- cbind(
  rep(1, nrow(Milan.winter)),
  Milan.winter$MaxTemp,
  Milan.winter$MaxTemp^2
)
u <- quantile(Milan.winter$PM10, prob = 0.9, type = 3, na.rm = TRUE)

fit2 <- fGEV(
  data = Milan.winter$PM10,
  par.start = c(50, 0, 0, 20, 1),
  method = "Bayesian",
  u = u,
  cov = cov,
  sig0 = 0.1,
  nsim = 5e+4
)

r <- range(Milan.winter$MaxTemp, na.rm = TRUE)
t <- seq(from = r[1], to = r[2], length = 50)
pU <- mean(Milan.winter$PM10 > u, na.rm = TRUE)

q2 <- ExtQ(
  P = 1 / c(600, 1200, 2400),
  method = "Bayesian",
  pU = pU,
  cov = cbind(rep(1, 50), t, t^2),
  param_post = fit2$param_post[-c(1:3e+4), ]
)

R <- c(min(unlist(q2)), 800)
qseq <- seq(from = R[1], to = R[2], length = 512)
Xl <- "Max Temperature"
Yl <- expression(PM[10])

for (i in seq_along(q2)) {
  K_q2 <- apply(q2[[i]], 2, function(x) density(x, from = R[1], to = R[2])$y)
  D <- cbind(expand.grid(t, qseq), as.vector(t(K_q2)))
  colnames(D) <- c("x", "y", "z")
  fields::image.plot(
    x = t, y = qseq, z = matrix(D$z, 50, 512),
    xlim = r, ylim = R, xlab = Xl, ylab = Yl
  )
}

}

##########################################################
### Example - Simulated data from Frechet distribution ###
##########################################################

if (interactive()) {

set.seed(999)
data <- extraDistr::rfrechet(n = 1500, mu = 3, sigma = 1, lambda = 1/3)

u <- quantile(data, probs = 0.9, type = 3)
fit3 <- fGEV(
  data = data,
  par.start = c(1, 2, 1),
  method = "Bayesian",
  u = u,
  sig0 = 1,
  nsim = 5e+4
)

pU <- mean(data > u)
P <- 1 / c(750, 1500, 3000)

q3 <- ExtQ(
  P = P,
  method = "Bayesian",
  pU = pU,
  param_post = fit3$param_post[-c(1:3e+4), ]
)

### Illustration

# Tail index estimation
ti_true <- 3
ti_ps <- fit3$param_post[-c(1:3e+4), 3]

K_ti <- density(ti_ps) # KDE of the tail index
H_ti <- hist(
  ti_ps, prob = TRUE, col = "lightgrey",
  ylim = range(K_ti$y), main = "", xlab = "Tail Index",
  cex.lab = 1.8, cex.axis = 1.8, lwd = 2
)
ti_ic <- quantile(ti_ps, probs = c(0.025, 0.975))

points(x = ti_ic, y = c(0, 0), pch = 4, lwd = 4)
lines(K_ti, lwd = 2, col = "dimgrey")
abline(v = ti_true, lwd = 2)
abline(v = mean(ti_ps), lwd = 2, lty = 2)

# Quantile estimation
q3_true <- extraDistr::qfrechet(
  p = P, mu = 3, sigma = 1, lambda = 1/3, lower.tail = FALSE
)

ci <- apply(log(q3), 2, function(x) quantile(x, probs = c(0.025, 0.975)))
K_q3 <- apply(log(q3), 2, density)

R <- range(log(c(q3_true, q3, data)))
Xlim <- c(log(quantile(data, 0.95)), R[2])
Ylim <- c(0, max(K_q3[[1]]$y, K_q3[[2]]$y, K_q3[[3]]$y))

plot(
  0, main = "", xlim = Xlim, ylim = Ylim,
  xlab = expression(log(x)), ylab = "Density",
  cex.lab = 1.8, cex.axis = 1.8, lwd = 2
)

cval <- c(211, 169, 105)
for (j in seq_along(P)) {
  col <- rgb(cval[j], cval[j], cval[j], 0.8 * 255, maxColorValue = 255)
  col2 <- rgb(cval[j], cval[j], cval[j], 255, maxColorValue = 255)
  polygon(K_q3[[j]], col = col, border = col2, lwd = 4)
}

points(log(data), rep(0, n), pch = 16)
# add posterior means
abline(v = apply(log(q3), 2, mean), lwd = 2, col = 2:4)
# add credible intervals
abline(v = ci[1, ], lwd = 2, lty = 3, col = 2:4)
abline(v = ci[2, ], lwd = 2, lty = 3, col = 2:4)

}