plot_ExtDep.np: Graphical summaries of non-parametric representations of extremal dependence.

Description

This function displays several summaries of extremal dependence represented through Bernstein polynomials.

Usage

plot_ExtDep.np(out, type, summary.mcmc, burn, y, probs, 
               A_true, h_true, est.out, mar1, mar2, dep, 
               QatCov1=NULL, QatCov2=QatCov1, P, 
               labels=c(expression(y[1]),expression(y[2])), 
               CEX=1.5, xlim, ylim, col.data, 
               col.Qfull, col.Qfade, data=NULL, ...)

Value

a graph depending on argument type.

Arguments

out

An output of the fExtDep.np function.

type

A character string indicating the type of graphical summary to be plotted. Takes values "summary", "returns", "A", "h", "pm", "k" or "Qsets".

summary.mcmc

The output of the summary_ExtDep function. Only required when out is obtained using a Bayesian estimation method (out$est=="Bayesian").

burn

The burn-in period. Only required when out is obtained using a Bayesian estimation method (out$est=="Bayesian").

y

A 2-column matrix of unobserved thresholds at which the returns are calculated. Required when type="returns".

probs

The probability of joint exceedances, the output of the return function.

A_true

A vector representing the true pickands dependence function evaluated at the grid points on the simplex given by summary.mcmc$w.

h_true

A vector representing the true angular density function evaluated at the grid points on the simplex given by summary.mcmc$w.

est.out

A list containing:

ghat:: a 3-row matrix giving the posterior summary for the estimate of the bound;

Shat and Shat_post:

a matrix of the posterior and a 3-row matrix giving the posterior summary for the estimate of the basic set $S$;

nuShat and nuShat_post:

a matrix of the posterior and a 3-row matrix giving the posterior summary for the estimate of the measure $\nu(S)$;

Note that a posterior summary is made of its mean and $90\%$credibility interval.

Only required when using a Bayesian estimation method (out$est=="Bayesian") and type="Qsets".

mar1, mar2

Vectors of marginal GEV parameters. Required when type="Qsets" and either out$method=="Bayesian" if the marginal parameter weren't fitted or "empirical".

dep

A logical value; if TRUE the estimate of the dependence is plotted when computing extreme quantile regions (type="Qsets").

QatCov1, QatCov2

Matrices representing the value of the covariates at which extreme quantile regions should be computed. Required when type="Qsets". See arguments cov1 and cov2 from fExtDep.

A vector indicating the probabilities associated with the quantiles to be computed. Required when type="Qsets".

labels

A bivariate vector of character strings providing labels for extreme quantile regions. Required when type="Qsets".

CEX

Label and axis sizes.

xlim, ylim

Limits of the x and y axis when computing extreme quantile regions. Required when type="Qsets".

col.data, col.Qfull, col.Qfade

Colors for data, estimate of extreme quantile regions and its credible interval (when applicable). Required when type="Qsets".

data

A 2-column matrix providing the original data to be plotted when type="Qsets".

...

Additional graphical parameters

Author

Simone Padoan, simone.padoan@unibocconi.it, https://faculty.unibocconi.it/simonepadoan/; Boris Beranger, borisberanger@gmail.com https://www.borisberanger.com; Giulia Marcon, giuliamarcongm@gmail.com

Details

If type="returns", a (contour) plot of the probabilities of exceedances for some threshold is returned. This corresponds to the output of the returns function.
If type="A", a plot of the estimated Pickands dependence function is drawn. If A_true is specified the plot includes the true Pickands dependence function and a functional boxplot for the estimated function. If type="h", a plot of the estimated angular density function is drawn. If h_true is specified the plot includes the true angular density and a functional boxplot for the estimated function. If type="pm", a plot of the prior against the posterior for the mass at $\{0\}$ is drawn. If type="k", a plot of the prior against the posterior for the polynomial degree $k$ is drawn. If type="summary", when the estimation was performed in a Bayesian framework then a 2 by 2 plot with types "A", "h", "pm" and "k" is returned. Otherwise a 1 by 2 plot with types "A" and "h" is returned. If type="Qsets", extreme quantile regions are computed according to the methodology developped in Beranger et al. (2021).

References

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

Marcon, G., Padoan, S.A., Naveau, P., Muliere, P., Segers, J. (2017) Multivariate Nonparametric Estimation of the Pickands Dependence Function using Bernstein Polynomials. Journal of Statistical Planning and Inference, 183, 1-17.

Examples

Run this code


###########################################################
### Example 1 - Wind Speed and Differential of pressure ###
###########################################################

data(WindSpeedGust)

years <- format(ParcayMeslay$time, format="%Y")
attach(ParcayMeslay[which(years %in% c(2004:2013)),])

# Marginal quantiles
WS_th <- quantile(WS,.9)
DP_th <- quantile(DP,.9)

# Standardisation to unit Frechet (requires evd package)
pars.WS <- evd::fpot(WS, WS_th, model="pp")$estimate
pars.DP <- evd::fpot(DP, DP_th, model="pp")$estimate
  
# transform the marginal distribution to common unit Frechet:
data_uf <- trans2UFrechet(cbind(WS,DP), type="Empirical")

# compute exceedances 
rdata <- rowSums(data_uf)
r0 <- quantile(rdata, probs=.90)
extdata_WSDP <- data_uf[rdata>=r0,]

# Fit
SP_mle <- fExtDep.np(method="Frequentist", data=extdata_WSDP, k0=10, type="maxima")

# Plot
plot_ExtDep.np(out=SP_mle, type="summary")

####################################################
### Example 2 - Pollution levels in Milan, Italy ###
####################################################
	
if (FALSE) {

### Here we will only model the dependence structure	
data(MilanPollution)

data <- Milan.winter[,c("NO2","SO2")] 
data <- as.matrix(data[complete.cases(data),])

# Thereshold
u <- apply(data, 2, function(x) quantile(x, prob=0.9, type=3))

# Hyperparameters
hyperparam <- list(mu.nbinom = 6, var.nbinom = 8, a.unif=0, b.unif=0.2)

### Standardise data to univariate Frechet margins

f1 <- fGEV(data=data[,1], method="Bayesian", sig0 = 0.0001, nsim = 5e+4)
diagnostics(f1)
burn1 <- 1:30000
gev.pars1 <- apply(f1$param_post[-burn1,],2,mean)
sdata1 <- trans2UFrechet(data=data[,1], pars=gev.pars1, type="GEV")

f2 <- fGEV(data=data[,2], method="Bayesian", sig0 = 0.0001, nsim = 5e+4)
diagnostics(f2)
burn2 <- 1:30000
gev.pars2 <- apply(f2$param_post[-burn2,],2,mean)
sdata2 <- trans2UFrechet(data=data[,2], pars=gev.pars2, type="GEV")

sdata <- cbind(sdata1,sdata2)

### Bayesian estimation using Bernstein polynomials

pollut1 <- fExtDep.np(method="Bayesian", data=sdata, u=TRUE,
                      mar.fit=FALSE, k0=5, hyperparam = hyperparam, nsim=5e+4)

diagnostics(pollut1)
pollut1_sum <- summary_ExtDep(mcmc=pollut1, burn=3e+4, plot=TRUE)
pl1 <- plot_ExtDep.np(out=pollut1, type="Qsets", summary.mcmc=pollut1_sum, 
                      mar1=gev.pars1, mar2=gev.pars2, P = 1/c(600, 1200, 2400),
                      dep=TRUE, data=data, xlim=c(0,400), ylim=c(0,400))

pl1b <- plot_ExtDep.np(out=pollut1, type="Qsets", summary.mcmc=pollut1_sum, est.out=pl1$est.out, 
                       mar1=gev.pars1, mar2=gev.pars2, P = 1/c(1200),
                       dep=FALSE, data=data, xlim=c(0,400), ylim=c(0,400))

### Frequentist estimation using Bernstein polynomials

pollut2 <- fExtDep.np(method="Frequentist", data=sdata, mar.fit=FALSE, type="rawdata", k0=8)
plot_ExtDep.np(out=pollut2, type = c("summary"), CEX=1.5)

pl2 <- plot_ExtDep.np(out=pollut2, type="Qsets", mar1=gev.pars1, mar2=gev.pars2, 
                      P = 1/c(600, 1200, 2400),
                      dep=TRUE, data=data, xlim=c(0,400), ylim=c(0,400),
                      labels=c(expression(NO[2]),expression(SO[2])),
                      col.Qfull = c("red", "green", "blue"))
                      
### Frequentist estimation using EKdH estimator

pollut3 <- fExtDep.np(method="Empirical", data=data)
plot_ExtDep.np(out=pollut3, type = c("summary"), CEX=1.5)

pl3 <- plot_ExtDep.np(out=pollut3, type="Qsets", mar1=gev.pars1, mar2=gev.pars2,
                      P = 1/c(600, 1200, 2400),
                      dep=TRUE, data=data, xlim=c(0,400), ylim=c(0,400),
                      labels=c(expression(NO[2]),expression(SO[2])),
                      col.Qfull = c("red", "green", "blue"))

}

Run the code above in your browser using DataLab