W.diag: Wadsworth's univariate and bivariate exponential threshold...

In mev: Modelling of Extreme Values

Wadsworth's univariate and bivariate exponential threshold diagnostics

Description

Function to produce diagnostic plots and test statistics for the threshold diagnostics exploiting structure of maximum likelihood estimators based on the non-homogeneous Poisson process likelihood

Usage

W.diag(
  xdat,
  model = c("nhpp", "exp", "invexp"),
  u = NULL,
  k,
  q1 = 0,
  q2 = 1,
  par = NULL,
  M = NULL,
  nbs = 1000,
  alpha = 0.05,
  plots = c("LRT", "WN", "PS"),
  UseQuantiles = FALSE,
  changepar = TRUE,
  ...
)

Arguments

xdat	a numeric vector of data to be fitted.
model	string specifying whether the univariate or bivariate diagnostic should be used. Either nhpp for the univariate model, exp (invexp) for the bivariate exponential model with rate (inverse rate) parametrization. See details.
u	optional; vector of candidate thresholds.
k	number of thresholds to consider (if u unspecified).
q1	lowest quantile for the threshold sequence.
q2	upper quantile limit for the threshold sequence (q2 itself is not used as a threshold, but rather the uppermost threshold will be at the (q_2-1/k) quantile).
par	parameters of the NHPP likelihood. If missing, the fit.pp routine will be run to obtain values
M	number of superpositions or 'blocks' / 'years' the process corresponds to (can affect the optimization)
nbs	number of simulations used to assess the null distribution of the LRT, and produce the p-value
alpha	significance level of the LRT
plots	vector of strings indicating which plots to produce; LRT= likelihood ratio test, WN = white noise, PS = parameter stability. Use NULL if you do not want plots to be produced
UseQuantiles	logical; use quantiles as the thresholds in the plot?
changepar	logical; if TRUE, the graphical parameters (via a call to par) are modified.
...	additional parameters passed to plot, overriding defaults including

Details

The function is a wrapper for the univariate (non-homogeneous Poisson process model) and bivariate exponential dependence model. For the latter, the user can select either the rate or inverse rate parameter (the inverse rate parametrization works better for uniformity of the p-value distribution under the LR test.

There are two options for the bivariate diagnostic: either provide pairwise minimum of marginally exponentially distributed margins or provide a n times 2 matrix with the original data, which is transformed to exponential margins using the empirical distribution function.

Value

plots of the requested diagnostics and an invisible list with components

• MLEmaximum likelihood estimates from all thresholds

• Covjoint asymptotic covariance matrix for \xi, \eta or \eta^{-1}.

• WNvalues of the white noise process

• LRTvalues of the likelihood ratio test statistic vs threshold

• pvalP-value of the likelihood ratio test

• kfinal number of thresholds used

• threshthreshold selected by the likelihood ratio procedure

• qthreshquantile level of threshold selected by the likelihood ratio procedure

• cthreshvector of candidate thresholds

• qcthreshquantile level of candidate thresholds

• mle.umaximum likelihood estimates for the selected threshold

• modelmodel fitted

Author(s)

Jennifer L. Wadsworth

References

Wadsworth, J.L. (2016). Exploiting Structure of Maximum Likelihood Estimators for Extreme Value Threshold Selection, Technometrics, 58(1), 116-126, http://dx.doi.org/10.1080/00401706.2014.998345.

Examples

## Not run: 
set.seed(123)
# Parameter stability only
W.diag(xdat = abs(rnorm(5000)), model = 'nhpp',
       k = 30, q1 = 0, plots = "PS")
W.diag(rexp(1000), model = 'nhpp', k = 20, q1 = 0)
xbvn <- mvrnorm(n = 6000,
                mu = rep(0, 2),
                Sigma = cbind(c(1, 0.7), c(0.7, 1)))
# Transform margins to exponential manually
xbvn.exp <- -log(1 - pnorm(xbvn))
#rate parametrization
W.diag(xdat = apply(xbvn.exp, 1, min), model = 'exp',
       k = 30, q1 = 0)
W.diag(xdat = xbvn, model = 'exp', k = 30, q1 = 0)
#inverse rate parametrization
W.diag(xdat = apply(xbvn.exp, 1, min), model = 'invexp',
       k = 30, q1 = 0)

## End(Not run)

mev documentation built on April 20, 2023, 5:10 p.m.