ExtremalDep (version 1.0.0)

pFailure: Probability of Falling into a Failure Region

Description

Compute the empirical probability of falling into two types of failure regions, based on exceedances simulated from a fitted extreme-value model.

Usage

pFailure(n, beta, u1, u2, mar1, mar2, type, plot, xlab, ylab, nlevels = 10)

Value

A list with elements AND and/or OR, depending on type. Each element is a matrix of size length(u1) x length(u2).

Arguments

n

Integer. Number of points generated to compute the empirical probability.

beta

Numeric vector. Bernstein polynomial coefficients.

u1, u2

Numeric vectors of positive thresholds.

mar1, mar2

Numeric vectors. Marginal GEV parameters.

type

Character. Type of failure region:

• "or" - at least one exceedance,

• "and" - both exceedances,

• "both" - compute both probabilities.

plot

Logical. If TRUE, display contour plots of the probabilities.

xlab, ylab

Character strings for axis labels in plots.

nlevels

Integer. Number of contour levels for plots.

Author

Simone Padoan (https://faculty.unibocconi.it/simonepadoan/) \ Boris Beranger (https://www.borisberanger.com)

Details

The two failure regions are: $$A_u = \{ (v_1, v_2): v_1 > u_1 \ \textrm{or}\ v_2 > u_2 \}$$ and $$B_u = \{ (v_1, v_2): v_1 > u_1 \ \textrm{and}\ v_2 > u_2 \}$$ for \((v_1, v_2) \in \mathbb{R}_+^2\), with thresholds \(u_1,u_2 > 0\).

Exceedance samples are generated following Algorithm 3 of Marcon et al. (2017). The empirical estimates are $$\hat{P}_{A_u} = \frac{1}{n}\sum_{i=1}^n I(y_{i1} > u_1^* \ \textrm{or}\ y_{i2} > u_2^*)$$ and $$\hat{P}_{B_u} = \frac{1}{n}\sum_{i=1}^n I(y_{i1} > u_1^* \ \textrm{and}\ y_{i2} > u_2^*)$$.

References

Marcon, G., Naveau, P. and Padoan, S.A. (2017). A semi-parametric stochastic generator for bivariate extreme events. Stat, 6, 184--201.

See Also

dExtDep, rExtDep, fExtDep, fExtDep.np

Examples

# Example: Wind speed and gust data
data(WindSpeedGust)
years <- format(ParcayMeslay$time, format = "%Y")
attach(ParcayMeslay[years %in% 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(x = extdata, method = "Frequentist", k0 = 10, type = "maxima")

# \donttest{
pF <- pFailure(
  n = 50000, beta = SP_mle$Ahat$beta,
  u1 = seq(19, 28, length = 200), mar1 = pars.WS,
  u2 = seq(40, 60, length = 200), mar2 = pars.DP,
  type = "both", plot = TRUE,
  xlab = "Daily-maximum Wind Speed (m/s)",
  ylab = "Differential of Pressure (mbar)",
  nlevels = 15
)
# }