Extremal index estimators based on interexceedance time and gap of exceedances

The function implements the maximum likelihood estimator and iteratively reweighted least square estimators of Suveges (2007) as well as the intervals estimator. The implementation differs from the presentation of the paper in that an iteration limit is enforced to make sure the iterative procedure terminates. Multiple thresholds can be supplied.

Usage

ext.index(
  xdat,
  q = 0.95,
  method = c("wls", "mle", "intervals"),
  plot = FALSE,
  warn = FALSE
)

Arguments

xdat: numeric vector of observations
q: a vector of quantile levels in (0,1). Defaults to 0.95
method: a string specifying the chosen method. Must be either wls for weighted least squares, mle for maximum likelihood estimation or intervals for the intervals estimator of Ferro and Segers (2003). Partial match is allowed.
plot: logical; if TRUE, plot the extremal index as a function of q
warn: logical; if TRUE, receive a warning when the sample size is too small

Value

a vector or matrix of estimated extremal index of dimension length(method) by length(q).

Details

The iteratively reweighted least square is a procedure based on the gaps of exceedances \(S_n=T_n-1\) The model is first fitted to non-zero gaps, which are rescaled to have unit exponential scale. The slope between the theoretical quantiles and the normalized gap of exceedances is \(b=1/\theta\), with intercept \(a=\log(\theta)/\theta\). As such, the estimate of the extremal index is based on \(\hat{\theta}=\exp(\hat{a}/\hat{b})\). The weights are chosen in such a way as to reduce the influence of the smallest values. The estimator exploits the dual role of \(\theta\) as the parameter of the mean for the interexceedance time as well as the mixture proportion for the non-zero component.

The maximum likelihood is based on an independence likelihood for the rescaled gap of exceedances, namely \(\bar{F}(u_n)S(u_n)\). The score equation is equivalent to a quadratic equation in \(\theta\) and the maximum likelihood estimate is available in closed form. Its validity requires however condition \(D^{(2)}(u_n)\) to apply; this should be checked by the user beforehand.

A warning is emitted if the effective sample size is less than 50 observations.

References

Ferro and Segers (2003). Inference for clusters of extreme values, JRSS: Series B, 65(2), 545-556.

Suveges (2007) Likelihood estimation of the extremal index. Extremes, 10(1), 41-55.

Suveges and Davison (2010), Model misspecification in peaks over threshold analysis. Annals of Applied Statistics, 4(1), 203-221.

Fukutome, Liniger and Suveges (2015), Automatic threshold and run parameter selection: a climatology for extreme hourly precipitation in Switzerland. Theoretical and Applied Climatology, 120(3), 403-416.

Author

Leo Belzile

Examples

set.seed(234)
#Moving maxima model with theta=0.5
a <- 1; theta <-  1/(1+a)
sim <- rgev(10001, loc=1/(1+a),scale=1/(1+a),shape=1)
x <- pmax(sim[-length(sim)]*a,sim[-1])
q <- seq(0.9,0.99,by=0.01)
ext.index(xdat=x,q=q,method=c('wls','mle'))
#>           90%      91%       92%       93%       94%       95%       96%
#> wls 0.6314522 0.633461 0.6381835 0.6410253 0.6501045 0.6677165 0.6685871
#> mle 0.4801204 0.480535 0.4832988 0.4884225 0.4864203 0.4868270 0.4931644
#>           97%       98%       99%
#> wls 0.6453696 0.6821534 0.6503444
#> mle 0.4909889 0.4989578 0.5003508