ExtQ {ExtremalDep}

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)

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".

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, the 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}.

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 returs 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.

Author(s)

Simone Padoan, [e-mail], https://faculty.unibocconi.it/simonepadoan/; Boris Beranger, [e-mail] https://www.borisberanger.com

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 ###
##################################################

## Not run: 

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 1:length(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)
}


## End(Not run)
  
##########################################################
### Example - Simulated data from Frechet distirbution ###
##########################################################
  
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 1:length(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)

}