ExtremalDep (version 0.0.3-3)

returns: Compute return values

Description

Predicts the probability of future simultaneous exceedances

Usage

returns(mcmc, summary.mcmc, y, plot=FALSE)

Arguments

mcmc

The output of the bbeed function.

summary.mcmc

The output of the summary.bbeed function.

y

A 2-column matrix of unobserved thresholds.

plot

If plot=TRUE, then the plot.bbeed function is used.

Value

Returns a vector whose length is equal to the number of rows of the input value y.

Details

Computes for a range of unobserved extremes (larger than those observed in a sample), the pointwise mean from the posterior predictive distribution of such predictive values. The probabilities are calculated through

\( P(Y_1 > y_1, Y_2 > y_2) = \frac{2}{k} \sum_{j=0}^{k-2} (\eta_{j+1} - \eta_j) \times \left( \frac{(j+1) B(y_1/(y_1+y_2)| j+2, k-j-1)}{y_1} - \frac{(k-j-1) B(y_2/(y_1+y_2)| k-j, j+1)}{y_2} \right), \)

where \(B(x|a,b)\) denotes the cumulative distribution function of a Beta random variable with shape \(a,b>0\). See Marcon et al. (2016, p.3323) for details.

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.

Examples

# NOT RUN {
	# This reproduces some of the results shown 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)

	plot.bbeed(type = "A", x=w, mcmc=mcmc, summary.mcmc, nsim=nsim, burn=burn)
	plot.bbeed(type = "h", x=w, mcmc=mcmc, summary.mcmc, nsim=nsim, burn=burn)
	plot.bbeed(type = "pm", x=w, mcmc=mcmc, summary.mcmc, nsim=nsim, burn=burn)
	plot.bbeed(type = "k", x=w, mcmc=mcmc, summary.mcmc, nsim=nsim, burn=burn)

	y <- seq(10,100,2)
	y <- as.matrix(expand.grid(y,y))
	probs <- returns(mcmc = mcmc, summary.mcmc = summary.mcmc, y = y, plot = TRUE)

# }