fExtDep.np: Non-parametric extremal dependence estimation

Description

This function estimates the bivariate extremal dependence structure using a non-parametric approach based on Bernstein Polynomials.

Usage

fExtDep.np(method, data, cov1=NULL, cov2=NULL, u, mar.fit=TRUE, 
           mar.prelim=TRUE, par10, par20, sig10, sig20, param0=NULL, 
           k0=NULL, pm0=NULL, prior.k="nbinom", prior.pm="unif", 
           nk=70, lik=TRUE, 
           hyperparam = list(mu.nbinom=3.2, var.nbinom=4.48), 
           nsim, warn=FALSE, type="rawdata")

Value

Outputs take the form of list including:

method:: The argument.
type:: whether it is "maxima" or "rawdata" (in the broader sense that a threshold exceedance model was taken).

If method="Bayesian" the list also includes:

mar.fit:: The argument.
pm:: The posterior sample of probability masses.
eta:: The posterior sample for the coeficients of the Bernstein polynomial.
k:: The posterior sample for the Bernstein polynoial order.
accepted:: A binary vector indicating if the proposal was accepted.
acc.vec:: A vector containing the acceptance probabilities for the dependence parameters at each iteration.
prior:: A list containing hyperparam, prior.pm and prior.k.
nsim:: The argument.
data:: The argument.

In addition if the marginal parameters are estimated (mar.fit=TRUE):

cov1, cov2:: The arguments.
accepted.mar, accepted.mar2:: Binary vectors indicating if the marginal proposals were accepted.
straight.reject1, straight.reject2:: Binary vectors indicating if the marginal proposals were rejected straight away as not respecting existence conditions (proposal is multivariate normal).
acc.vec1, acc.vec2:: Vectors containing the acceptance probabilities for the marginal parameters at each iteration.
sig1.vec, sig2.vec:: Vectors containing the values of the scaling parameter in the marginal proposal distributions.

Finally, if the argument u is provided, the list also contains:

threshold:: A bivariate vector indicating the threshold level for both margins.
kn:: The empirical estimate of the probability of being greater than the threshold.

When method="Frequentist", the list includes:

When method="Empirical", the list includes:

hhat:: An estimate of the angular density.
Hhat:: An estimate of the angular measure.
p0, p1:: The estimates of the probability mass at 0 and 1.
Ahat:: An estimate of the PIckands dependence function.
w:: A sequence of value on the bivariate unit simplex.
q:: A real in $[0,1]$ indicating the quantile associated with the threshold u. Takes value NULL if type="maxima".
data:: The data on the unit Frechet scale (empirical transformation) if type="rawdata" and mar.fit=TRUE. Data on the original scale otherwise.

Arguments

method: A character string indicating the estimation method inlcuding "Bayesian", "Frequentist" and "Empirical".
data: A matrix containing the data.
cov1, cov2: A covariate vector/matrix for linear model on the location parameter of the marginal distributions. length(cov1)/nrow(cov1) needs to match nrow(data). If NULL it is assumed to be constant. Required when method="Bayesian".
u: When method="Bayesian": a bivariate indicating the threshold for the exceedance approach. If logical TRUE the threshold are calculated by default as the 90% quantiles. When missing, a block maxima approach is considered. When method="Frequentist": the associated quantile is used to select observations with the largest radial components. If missing, the 90% quantile is used.
mar.fit: A logical value indicated whether the marginal distributions should be fitted. When method="Frequentist", data are empirically transformed to unit Frechet scale if mar.fit=TRUE. Not required when method="Empirical".
rawdata: A character string specifying if the data is "rawdata" or "maxima". Required when method="Frequentist".
mar.prelim: A logical value indicated whether a preliminary fit of marginal distributions should be done prior to estimating the margins and dependence. Required when method="Bayesian".
par10, par20: Vectors of starting values for the marginal parameter estimation. Required when method="Bayesian" and mar.fit=TRUE
sig10, sig20: Positive reals representing the initial value for the scaling parameter of the multivariate normal proposal distribution for both margins. Required when method="Bayesian" and mar.fit=TRUE
param0: A vector of initial value for the Bernstein polynomial coefficients. It should be a list with elements $eta and $beta which can be generated by the internal function rcoef if missing. Required when method="Bayesian".
k0: An integer indicating the initial value of the polynomial order. Required when method="Bayesian".
pm0: A list of initial values for the probability masses at the boundaries of the simplex. It should be a list with two elements $p0 and $p1. Randomly generated if missing. Required when method="Bayesian".
prior.k: A character string indicating the prior distribution on the polynomial order. By default prior.k="nbinom" (negative binomial) but can also be "pois" (Poisson). Required when method="Bayesian".
prior.pm: A character string indicating the prior on the probability masses at the endpoints of the simplex. By default prior.pm="unif" (uniform) but can also be "beta" (beta). Required when method="Bayesian".
nk: An integer indicating the maximum polynomial order. Required when method="Bayesian".
lik: A logical value; if FALSE, an approximation of the likelihood, by means of the angular measure in Bernstein form is used (TRUE by default). Required when method="Bayesian".
hyperparam: A list of the hyper-parameters depending on the choice of prior.k and prior.pm. See details. Required when method="Bayesian".
nsim: An integer indicating the number of iterations in the Metropolis-Hastings algorithm. Required when method="Bayesian".
warn: A logical value. If TRUE warnings are printed when some specifics (e.g., param0, k0, pm0 and hyperparam) are not provided and set by default. Required when method="Bayesian".
type: A character string indicating whther the data are "rawdata" or "maxima". Required when method="Bayesian".

Author

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

Details

When method="Bayesian", the vector of hyper-parameters is provided in the argument hyperparam. It should include:

If prior.pm="unif": requires hyperparam$a.unif and hyperparam$b.unif.
If prior.pm="beta": requires hyperparam$a.beta and hyperparam$b.beta.
If prior.k="pois": requires hyperparam$mu.pois.
If prior.k="nbinom": requires hyperparam$mu.nbinom and hyperparam$var.nbinom or hyperparam$pnb and hyperparam$rnb. The relationship is pnb = mu.nbinom/var.nbinom and rnb = mu.nbinom^2 / (var.nbinom-mu.nbinom).

When u is specified Algorithm 1 of Beranger et al. (2021) is applied whereas when it is not specified Algorithm 3.5 of Marcon et al. (2016) is considered.

When method="Frequentist", if type="rawdata" then pseudo-polar coordinates are extracted and only observations with a radial component above some high threshold (the quantile equivalent of u for the raw data) are retained. The inferential approach proposed in Marcon et al. (2017) based on the approximate likelihood is applied.

When method="Empirical", the empirical estimation procedure presented in Einmahl et al. (2013) is applied.

References

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

Einmahl, J. H. J., de Haan, L. and Krajina, A. (2013). Estimating extreme bivariate quantile regions. Extremes, 16, 121-145.

Marcon, G., Padoan, S. A. and Antoniano-Villalobos, I. (2016). Bayesian inference for the extremal dependence. Electronic Journal of Statistics, 10, 3310-3337.

Marcon, G., Padoan, S.A., Naveau, P., Muliere, P. and 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 Bayesian estimation, 
# Threshold exceedances approach, threshold set by default
# Joint estimation margins + dependence
# Default uniform prior on pm
# Default negative binomial prior on polynomial order
# Quadratic relationship between location and max temperature

if (FALSE) {
  data(MilanPollution)
  data <- Milan.winter[,c("NO2", "SO2", "MaxTemp")]
  data <- data[complete.cases(data),]
  
  covar <- cbind(rep(1,nrow(data)), data[,3], data[,3]^2)
  hyperparam <- list(mu.binom=6, var.binom=8, a.unif=0, b.unif=0.2)
  pollut <- fExtDep.np(method="Bayesian", data = data[,-3], u=TRUE,
                       cov1 = covar, cov2 = covar, mar.prelim=FALSE,
                       par10 = c(100,0,0,35,1), par20 = c(20,0,0,20,1),
                       sig10 = 0.1, sig20 = 0.1, k0 = 5,
                       hyperparam = hyperparam, nsim = 5e+4)
  # Warning: This is slow!                     
}

# Example Frequentist estimation
# Data are maxima

data(WindSpeedGust)
years <- format(ParcayMeslay$time, format="%Y")
attach(ParcayMeslay[which(years %in% c(2004:2013)),])
  
WS_th <- quantile(WS,.9)
DP_th <- quantile(DP,.9)
  
pars.WS <- evd::fpot(WS, WS_th, model="pp")$estimate
pars.DP <- evd::fpot(DP, DP_th, model="pp")$estimate
  
data_uf <- trans2UFrechet(cbind(WS,DP), type="Empirical")
  
rdata <- rowSums(data_uf)
r0 <- quantile(rdata, probs=.90)
extdata <- data_uf[rdata>=r0,]
  
SP_mle <- fExtDep.np(method="Frequentist", data=extdata, k0=10, 
                     type="maxima")

Run the code above in your browser using DataLab