# bbeed: Bayesian Estimation of Extremal Dependence In ExtremalDep: Extremal Dependence Models

## Description

Nonparametric Bayesian estimation of the extremal dependence function (angular measure and Pickands dependence function) for two-dimensional data on the basis of Bernstein polynomials representations.

## Usage

 ```1 2 3``` ```bbeed(data, pm0, param0, k0, hyperparam, nsim, nk = 70, prior.k = c('pois','nbinom'), prior.pm = c('unif','beta', 'null'), lik = TRUE) ```

## Arguments

 `data` (n x 2) matrix of component-wise maxima (`n` is the number of replications). `pm0` list; the initial state of the Markov chain of point masses p0 and p1. Randomly generated if missing. `param0` real vector in (0,1) of length `k0`+1; the initial state of the Markov chain of η coefficients (see Details). Randomly generated if missing. `k0` postive integer indicating the initial state of the Markov chain of the Bernstein polynomial's degree. Randomly generated if missing. `hyperparam` list containing the hyperparameters (see Details). `nsim` postive integer indicating the number of the iterations of the chain. `nk` positive integer, maximum value of k that can be reached (`nk = 70` by default). `prior.k` string, indicating the prior distribution on the degree of the Bernstein polynomial, k. If missing, `prior.k = "pois"` by default. `prior.pm` string, indicating the prior distribution on the point masses p0 and p1. If missing, `prior.pm = "unif"` by default. `lik` logical; if `FALSE`, an approximation of the likelihood, by means of the angular measure in Bernstein form is used (`TRUE` by default).

## Details

The hyperparameters and starting values may not be specified and in this case a warning message will be printed and default values are set. In particular if `hyperparam` is missing we have

``` hyperparam <- NULL hyperparam\$a.unif <- 0 hyperparam\$b.unif <- .5 hyperparam\$a.beta <- c(.8,.8) hyperparam\$b.beta <- c(5,5) mu.pois <- hyperparam\$mu.pois <- 4 mu.nbinom <- hyperparam\$mu.nbinom <- 4 var.nbinom <- 8 pnb <- hyperparam\$pnb <- mu.nbinom/var.nbinom rnb <- hyperparam\$rnb <- mu.nbinom^2 / (var.nbinom-mu.nbinom) ```

The routine returns an estimate of the Pickands dependence function using the Bernstein polynomials approximation proposed in Marcon et al. (2016). The method is based on a preliminary empirical estimate of the Pickands dependence function. If you do not provide such an estimate, this is computed by the routine. In this case, you can select one of the empirical methods available. `est = "ht"` refers to the Hall-Tajvidi estimator (Hall and Tajvidi 2000). With `est = "cfg"` the method proposed by Caperaa et al. (1997) is considered. Note that in the multivariate case the adjusted version of Gudendorf and Segers (2011) is used. Finally, with `est = "md"` the estimate is based on the madogram defined in Marcon et al. (2016).

For each row of the (m x d) design matrix `x` is a point in the unit `d`-dimensional simplex,

S_d := { (w_1,..., w_d) in [0,1]^{d}: ∑_{i=1}^{d} w_i = 1}.

With this "regularization" method, the final estimate satisfies the neccessary conditions in order to be a Pickands dependence function.

A(w) = ∑_{α in Γ_k} β_α b_α (w;k).

The estimates are obtained by solving an optimization quadratic problem subject to the constraints. The latter are represented by the following conditions: A(e_i)=1; \max(w_i)<=A(w)<=1; for all i=1,...,d; (convexity).

The order of polynomial `k` controls the smoothness of the estimate. The higher `k` is, the smoother the final estimate is. Higher values are better with strong dependence (e. g. `k=23`), whereas small values (e.g. `k=6` or `k=10`) are enough with mild or weak dependence.

An empirical transformation of the marginals is performed when `margin="emp"`. A max-likelihood fitting of the GEV distributions is implemented when `margin="est"`. Otherwise it refers to marginal parametric GEV theorethical distributions (`margin = "exp", "frechet", "gumbel"`).

If `bounds = TRUE`, a modification can be implemented to satisfy \max(w,1-w) <= A_n(w) <= 1 for all 0 <= w <= 1:

A_n(w) = min(1, max{A_n(w), w, 1-w}).

## Value

return(list(pm=spm, eta=seta, k=sk, accepted=accepted/nsim, prior = list(hyperparam=hyperparam, k=prior.k, pm=prior.pm)))

 `beta` vector of polynomial coefficients `A` numeric vector of the estimated Pickands dependence function A `Anonconvex` preliminary non-convex function `extind` extremal index

## Note

The number of coefficients depends on both the order of polynomial `k` and the dimension `d`. The number of parameters is choose(k+d-1,d-1).

The size of the vector `beta` must be compatible with the polynomial order `k` chosen.

With the estimated polynomial coefficients, the extremal coefficient, i.e. d*A(1/d,...,1/d) is computed.

## Author(s)

Simone Padoan, [email protected], http://faculty.unibocconi.it/simonepadoan; Boris Beranger, [email protected] http://www.borisberanger.com; Giulia Marcon, [email protected]

## References

Marcon G., Padoan, S.A. and Antoniano-Villalobos I. (2016) Bayesian Inference for the Extremal Dependence. Electronic Journal of Statistics, 10(2), 3310-3337.

`summary.bbeed`
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19``` ```# This reproduces some of the results showed in Fig. 1 (Marcon, 2016). set.seed(1890) data <- evd::rbvevd(n=100, dep=.6, asy=c(0.8,0.3), model="alog", mar1=c(1,1,1)) nsim = 500000 burn = 400000 mu.nbinom = 3.2 var.nbinom = 4.48 hyperparam <- list(a.unif=0, b.unif=.5, mu.nbinom=mu.nbinom, var.nbinom=var.nbinom) k0 = 5 pm0 = list(p0=0.06573614, p1=0.3752118) eta0 = ExtremalDep:::rcoef(k0, pm0) mcmc <- bbeed(data, pm0, eta0, k0, hyperparam, nsim, prior.k = "nbinom", prior.pm = "unif") w <- seq(0.001, .999, length=100) summary.mcmc <- summary.bbeed(w, mcmc, burn, nsim, plot=TRUE) ```