plot_ExtDep.np {ExtremalDep} | R Documentation |
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, ...)
Arguments
out |
An output of the |
type |
A character string indicating the type of graphical summary to be plotted. Takes values |
summary.mcmc |
The output of the |
burn |
The burn-in period. Only required when |
y |
A 2-column matrix of unobserved thresholds at which the returns are calculated. Required when |
probs |
The probability of joint exceedances, the output of the |
A_true |
A vector representing the true pickands dependence function evaluated at the grid points on the simplex given by |
h_true |
A vector representing the true angular density function evaluated at the grid points on the simplex given by |
est.out |
A list containing:
Note that a posterior summary is made of its mean and Only required when using a Bayesian estimation method ( |
mar1 , mar2 |
Vectors of marginal GEV parameters. Required when |
dep |
A logical value; if |
QatCov1 , QatCov2 |
Matrices representing the value of the covariates at which extreme quantile regions should be computed. Required when |
P |
A vector indicating the probabilities associated with the quantiles to be computed. Required when |
labels |
A bivariate vector of character strings providing labels for extreme quantile regions. Required when |
CEX |
Label and axis sizes. |
xlim , ylim |
Limits of the x and y axis when computing extreme quantile regions. Required when |
col.data , col.Qfull , col.Qfade |
Colors for data, estimate of extreme quantile regions and its credible interval (when applicable). Required when |
data |
A 2-column matrix providing the original data to be plotted when |
... |
Additional graphical parameters |
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 is drawn.
If
type="k"
, a plot of the prior against the posterior for the polynomial degree 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).
Value
a graph depending on argument type
.
Author(s)
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
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.
See Also
Examples
###########################################################
### 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 ###
####################################################
## Not run:
### 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"))
## End(Not run)