R: Non-parametric extremal dependence estimation

fExtDep.np {ExtremalDep}

R Documentation

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

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

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.

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.

Author(s)

Simone Padoan, simone.padoan@unibocconi.it, https://faculty.unibocconi.it/simonepadoan/; Boris Beranger, borisberanger@gmail.com 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.

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

# 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

## Not run: 
  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!                     

## End(Not run)

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

[Package ExtremalDep version 0.0.4-1 Index]