Movatterモバイル変換

Type:	Package
Title:	Modelling of Extreme Values
Version:	2.1
Description:	Various tools for the analysis of univariate, multivariate and functional extremes. Exact simulation from max-stable processes (Dombry, Engelke and Oesting, 2016, <doi:10.1093/biomet/asw008>, R-Pareto processes for various parametric models, including Brown-Resnick (Wadsworth and Tawn, 2014, <doi:10.1093/biomet/ast042>) and Extremal Student (Thibaud and Opitz, 2015, <doi:10.1093/biomet/asv045>). Threshold selection methods, including Wadsworth (2016) <doi:10.1080/00401706.2014.998345>, and Northrop and Coleman (2014) <doi:10.1007/s10687-014-0183-z>. Multivariate extreme diagnostics. Estimation and likelihoods for univariate extremes, e.g., Coles (2001) <doi:10.1007/978-1-4471-3675-0>.
License:	GPL-3
URL:	https://lbelzile.github.io/mev/,https://github.com/lbelzile/mev/
BugReports:	https://github.com/lbelzile/mev/issues/
Depends:	R (≥ 3.4)
Imports:	alabama, methods, nleqslv, numDeriv, Rcpp (≥ 0.12.16),Rsolnp, stats,
Suggests:	boot, cobs, evd, expint, knitr, MASS, mvPot (≥ 0.1.4),mvtnorm, gmm, revdbayes, rmarkdown, ismev, tinytest,TruncatedNormal (≥ 1.1)
LinkingTo:	Rcpp, RcppArmadillo
LazyData:	true
RoxygenNote:	7.3.3
VignetteBuilder:	knitr
Encoding:	UTF-8
ByteCompile:	true
NeedsCompilation:	yes
Packaged:	2025-11-10 18:08:19 UTC; lbelzile
Author:	Leo Belzile [aut, cre], Jennifer L. Wadsworth [aut], Paul J. Northrop [aut], Raphael Huser [aut], Scott D. Grimshaw [aut], Jin Zhang [ctb], Michael A. Stephens [ctb], Art B. Owen [ctb]
Maintainer:	Leo Belzile <belzilel@gmail.com>
Repository:	CRAN
Date/Publication:	2025-11-11 06:10:27 UTC

Joint maximum likelihood estimation for exponential model

Description

Calculates the MLEs of the rate parameter, and joint asymptotic covariance matrix of these MLEsover a range of thresholds as supplied by the user.

Usage

.Joint_MLE_Expl(x, u = NULL, k, q1, q2 = 1, param)

Arguments

Value

a list with

• mle vector of MLEs above the supplied thresholds

• cov joint asymptotic covariance matrix of these MLEs

Author(s)

Jennifer L. Wadsworth

GP fitting function of Grimshaw (1993)

Description

Function for estimating parametersk anda for a random sample from a GPD.

Usage

.fit.gpd.grimshaw(x)

Arguments

Value

a list with the maximum likelihood estimates of componentsa andk

Author(s)

Scott D. Grimshaw

Robust threshold selection of Dupuis

Description

The optimal bias-robust estimator (OBRE) for the generalized Pareto.This function returns robust estimates and the associated weights.

Usage

.fit.gpd.rob(dat, thresh, k = 4, tol = 1e-05, show = FALSE)

Arguments

Value

a list with the same components asfit.gpd,in addition to

• estimate: optimal bias-robust estimates of thescale andshape parameters.

• weights: vector of OBRE weights.

References

Dupuis, D.J. (1998). Exceedances over High Thresholds: A Guide to Threshold Selection,Extremes,1(3), 251–261.

See Also

fit.gpd

Examples

dat <- rexp(100).fit.gpd.rob(dat, 0.1)

Posterior predictive distribution and density for the GEV distribution

Description

This function calculates the posterior predictive density at points xbased on a matrix of posterior density parameters

Usage

.gev.postpred(x, posterior, Nyr = 100, type = c("density", "quantile"))

Arguments

Value

a vector of values for the posterior predictive density or quantile atx, depending on the value oftype

Maximum likelihood method for the generalized Pareto Model

Description

Maximum-likelihood estimation for the generalized Pareto model, including generalized linear modelling of each parameter. This function was adapted by Paul Northrop to include the gradient in thegpd.fit routine fromismev.

Usage

.gpd_2D_fit(  xdat,  threshold,  npy = 365,  ydat = NULL,  sigl = NULL,  shl = NULL,  siglink = identity,  shlink = identity,  siginit = NULL,  shinit = NULL,  show = TRUE,  method = "Nelder-Mead",  maxit = 10000,  ...)

Arguments

Details

For non-stationary fitting it is recommended that the covariates within the generalized linear models are (at least approximately) centered and scaled (i.e. the columns ofydat should be approximately centered and scaled).

The form of the GP model used follows Coles (2001) Eq (4.7). In particular, the shape parameter is defined so that positive values imply a heavy tail and negative values imply a bounded upper value.

Value

a list with components

nexc

scalar giving the number of threshold exceedances.

nllh

scalar giving the negative log-likelihood value.

mle

numeric vector giving the MLE's for the scale and shape parameters, resp.

rate

scalar giving the estimated probability of exceeding the threshold.

se

numeric vector giving the standard error estimates for the scale and shape parameter estimates, resp.

trans

logical indicator for a non-stationary fit.

model

list with componentssigl andshl.

link

character vector giving inverse link functions.

threshold

threshold, or vector of thresholds.

nexc

number of data points above the threshold.

data

data that lie above the threshold. For non-stationary models, the data are standardized.

conv

convergence code, taken from the list returned byoptim. A zero indicates successful convergence.

nllh

negative log likelihood evaluated at the maximum likelihood estimates.

vals

matrix with three columns containing the maximum likelihood estimates of the scale and shape parameters, and the threshold, at each data point.

mle

vector containing the maximum likelihood estimates.

rate

proportion of data points that lie above the threshold.

cov

covariance matrix.

se

numeric vector containing the standard errors.

n

number of data points (i.e., the length ofxdat).

npy

number of observations per year/block.

xdata

data that has been fitted.

References

Coles, S., 2001. An Introduction to Statistical Modeling of Extreme Values. Springer-Verlag, London.

Is the matrix conditionally negative semi-definite?Function adapted from 'is.CNSD' in the CEGO package, v 2.1.0

Description

Is the matrix conditionally negative semi-definite?Function adapted from 'is.CNSD' in the CEGO package, v 2.1.0

Usage

.is.CNSD(X, tol = 1e-08)

Arguments

Author(s)

Martin Zaefferer

Internal function

Description

Takes a list of asymmetry parameters with an associated dependence vector and returnsthe corresponding asymmetry matrix for the asymmetric logistic and asymmetric negative logistic models

Usage

.mvasym.check(asy, dep, d, model = c("alog", "aneglog", "maxlin"))

Arguments

Details

This function is extracted from the evd package and modified(C) Alec Stephenson

Value

a matrix of asymmetry components, enumerating all possible2^d-1 subsets ofthe power set

Generate from Brown-Resnick processY \sim {P_x}, whereP_{x} is probability of extremal function

Description

Generate from Brown-Resnick processY \sim {P_x}, whereP_{x} is probability of extremal function

Usage

.rPBrownResnick(index, Sigma_chol, Sigma)

Arguments

Value

ad-vector fromP_x

Generate from extremal Husler-Reiss distributionY \sim {P_x}, whereP_{x} is probability of extremal function

Description

Generate from extremal Husler-Reiss distributionY \sim {P_x}, whereP_{x} is probability of extremal function

Usage

.rPHuslerReiss(index, cholesky, Sigma)

Arguments

Value

ad-vector fromP_x

Generate from Smith model (moving maxima)Y \sim {P_x}, whereP_{x} is probability of extremal function

Description

Generate from Smith model (moving maxima)Y \sim {P_x}, whereP_{x} is probability of extremal function

Usage

.rPSmith(index, Sigma_chol, loc)

Arguments

Value

ad-vector fromP_x

Generate from bilogisticY \sim {P_x}, whereP_{x} is the probability of extremal functions

Description

Generate from bilogisticY \sim {P_x}, whereP_{x} is the probability of extremal functions

Usage

.rPbilog(d, index, alpha)

Arguments

Value

ad-vector fromP_x

Generate from extremal DirichletY \sim {P_x}, whereP_{x} is the probability of extremal functions from the Dirichlet model ofColes and Tawn.

Description

Note: we generate from the Dirichlet rather than the Gamma distribution, since the former is parallelized

Usage

.rPdir(d, index, alpha, irv = FALSE)

Arguments

Value

ad-vector fromP_x

Generate from extremal DirichletY \sim {P_x}, whereP_{x} is the probability of extremal functions from a Dirichlet mixture

Description

Generate from extremal DirichletY \sim {P_x}, whereP_{x} is the probability of extremal functions from a Dirichlet mixture

Usage

.rPdirmix(d, index, alpha, weight)

Arguments

Value

ad-vector fromP_x

Generate from extremal Student-tY \sim {P_x}, whereP_{x} is probability of extremal function

Description

Generate from extremal Student-tY \sim {P_x}, whereP_{x} is probability of extremal function

Usage

.rPexstud(index, cholesky, sigma, al)

Arguments

Value

ad-vector fromP_x

Generate from logisticY \sim {P_x}, whereP_{x} is probability of extremal function scaled by a Frechet variate

Description

Generate from logisticY \sim {P_x}, whereP_{x} is probability of extremal function scaled by a Frechet variate

Usage

.rPlog(d, index, theta)

Arguments

Value

ad-vector fromP_x

Generate from negative logisticY \sim {P_x}, whereP_{x} is probability of extremal function scaled by a Frechet variate

Description

Generate from negative logisticY \sim {P_x}, whereP_{x} is probability of extremal function scaled by a Frechet variate

Usage

.rPneglog(d, index, theta)

Arguments

Value

ad-vector fromP_x

Samples from exceedances at site (scaled extremal function definition)

Description

Samples from exceedances at site (scaled extremal function definition)

Usage

.rPsite(n, j, d, par, model, Sigma, loc)

Arguments

Value

an byd matrix containing the sample

Generates fromQ_i, the spectral measure of the bilogistic model

Description

Simulation algorithm of Boldi (2009) for the bilogistic model

Usage

.rbilogspec(n, alpha)

Arguments

Value

ann byd sample from the spectral distribution

References

Boldi (2009). A note on the representation of parametric modelsfor multivariate extremes.Extremes12, 211–218.

Generates fromQ_i, the spectral measure of the Brown-Resnick model

Description

Simulation algorithm of Dombry et al. (2015)

Usage

.rbrspec(n, Sigma_chol, Sigma)

Arguments

Value

ann byd sample from the spectral distribution

References

Dombry, Engelke and Oesting (2016). Exact simulation of max-stable processes,Biometrika,103(2), 303–317.

Generates fromQ_i, the spectral measure of the Dirichlet mixture model

Description

Simulation algorithm of Dombry et al. (2015)

Usage

.rdirmixspec(n, d, alpha, weight)

Arguments

Value

ann byd sample from the spectral distribution

Generates fromQ_i, the spectral measure of the extremal Dirichletmodel

Description

This model was introduced in Coles and Tawn (1991); thepresent method uses the simulation algorithm of Boldi (2009) for the extremal Dirichlet model

Usage

.rdirspec(n, d, alpha, irv = FALSE)

Arguments

Value

ann byd sample from the spectral distribution

References

Boldi (2009). A note on the representation of parametric modelsfor multivariate extremes.Extremes12, 211–218.

Generates fromQ_i, the spectral measure of the extremal Student model

Description

Generates fromQ_i, the spectral measure of the extremal Student model

Usage

.rexstudspec(n, sigma, al)

Arguments

Value

ann byd sample from the spectral distribution

Generates fromQ_i, the spectral measure of the Husler-Reiss model

Description

Generates fromQ_i, the spectral measure of the Husler-Reiss model

Usage

.rhrspec(n, Lambda)

Arguments

Value

ann byd sample from the spectral distribution

Generates fromQ_i, the spectral measure of the logistic model

Description

Simulation algorithm of Dombry et al. (2015)

Usage

.rlogspec(n, d, theta)

Arguments

Value

ann byd sample from the spectral distribution

References

Dombry, Engelke and Oesting (2016). Exact simulation of max-stable processes,Biometrika,103(2), 303–317.

Multivariate extreme value distribution sampling algorithm via angular measure

Description

This algorithm corresponds to Algorithm 1 in Dombry, Engelke and Oesting (2016),using the formulation of the Dirichlet mixture of Coles and Tawn (1991)as described and derived in Boldi (2009) for the bilogistic and extremalDirichlet model. Models currently implemented include logistic, negativelogistic, extremal Dirichlet and bilogistic MEV.

Usage

.rmevA1(n, d, par, model, Sigma, loc)

Arguments

Value

an byd matrix containing the sample

Multivariate extreme value distribution sampling algorithm via extremal functions

Description

Code implementing Algorithm 2 in Dombry, Engelke and Oesting (2016)

Usage

.rmevA2(n, d, par, model, Sigma, loc)

Arguments

Value

an byd matrix containing the sample

Random samples from asymmetric logistic distribution

Description

Simulation algorithm of Stephenson (2003), using exact-samples from the logistic

Usage

.rmevasy(n, d, par, asym, ncompo, Sigma, model)

Arguments

Value

an byd matrix containing the sample

References

Stephenson, A. G. (2003) Simulating multivariate extreme value distributions of logistic type.Extremes,6(1), 49–60.

Joe, H. (1990). Families of min-stable multivariate exponential and multivariateextreme value distributions,9, 75–81.

Random sampling from spectral distribution on l1 sphere

Description

Generate fromQ_i, the spectral measure of a given multivariate extreme value model

Usage

.rmevspec_cpp(n, d, par, model, Sigma, loc)

Arguments

Value

an byd matrix containing the sample

References

Dombry, Engelke and Oesting (2016). Exact simulation of max-stable processes,Biometrika,103(2), 303–317.

Boldi (2009). A note on the representation of parametric models for multivariate extremes.Extremes12, 211–218.

Multivariate Normal distribution sampler (Rcpp version), derived using the eigendecompositionof the covariance matrix Sigma. The function utilizes the arma random normal generator

Description

Multivariate Normal distribution sampler (Rcpp version), derived using the eigendecompositionof the covariance matrix Sigma. The function utilizes the arma random normal generator

Usage

.rmnorm_arma(n, Mu, Xmat, eigen = TRUE)

Arguments

Value

ann sample from a multivariate Normal distribution

Generates fromQ_i, the spectral measure of the negative logistic model

Description

Simulation algorithm of Dombry et al. (2016)

Usage

.rneglogspec(n, d, theta)

Arguments

Value

ann byd sample from the spectral distribution

References

Dombry, Engelke and Oesting (2016). Exact simulation of max-stable processes,Biometrika,103(2), 303–317.

Generates fromQ_i, the spectral measure of the pairwise Beta model

Description

This model was introduced in Cooley, Davis and Naveau (2010).The sample is drawn from a mixture and the algorithm follows from the proof of Theorem 1 in Ballani and Schlather (2011)and is written in full in Algorithm 1 of Sabourin et al. (2013).

Usage

.rpairbetaspec(n, d, alpha, beta)

Arguments

Value

a matrix of samples from the angular distribution

ann byd sample from the spectral distribution

Author(s)

Leo Belzile

References

Cooley, D., R.A. Davis and P. Naveau (2010). The pairwise beta distribution: A flexible parametric multivariate model for extremes,Journal of Multivariate Analysis,101(9), 2103–2117.

Ballani, D. and M. Schlather (2011). A construction principle for multivariate extreme value distributions,Biometrika,98(3), 633–645.

Sabourin, A., P. Naveau and A. Fougeres (2013). Bayesian model averaging for extremes,Extremes,16, 325–350.

Generates fromQ_i, the spectral measure of the pairwise exponential model

Description

The sample is drawn from a mixture and the algorithm follows from the proof of Theorem 1 in Ballani and Schlather (2011).

Usage

.rpairexpspec(n, d, alpha, beta)

Arguments

Value

a matrix of samples from the angular distribution

ann byd sample from the spectral distribution

Author(s)

Leo Belzile

References

Ballani, D. and M. Schlather (2011). A construction principle for multivariate extreme value distributions,Biometrika,98(3), 633–645.

Generates fromQ_i, the spectral measure of the Smith model (moving maxima)

Description

Simulation algorithm of Dombry et al. (2015)

Usage

.rsmithspec(n, Sigma_chol, loc)

Arguments

Value

ann byd sample from the spectral distribution

References

Dombry, Engelke and Oesting (2016). Exact simulation of max-stable processes,Biometrika,103(2), 303–317.

Generates fromQ_i, the spectral measure of the weighted Dirichlet model

Description

This model was introduced in Ballani and Schlather (2011).

Usage

.rwdirbsspec(n, d, alpha, beta)

Arguments

Value

a matrix of samples from the angular distribution

ann byd sample from the spectral distribution

Author(s)

Leo Belzile

References

Ballani, D. and M. Schlather (2011). A construction principle for multivariate extreme value distributions,Biometrika,98(3), 633–645.

Generates fromQ_i, the spectral measure of the weighted exponential model

Description

This model was introduced in Ballani and Schlather (2011).

Usage

.rwexpbsspec(n, d, alpha, beta)

Arguments

Value

a matrix of samples from the angular distribution

ann byd sample from the spectral distribution

Author(s)

Leo Belzile

References

Ballani, D. and M. Schlather (2011). A construction principle for multivariate extreme value distributions,Biometrika,98(3), 633–645.

Weighted empirical distribution function

Description

Compute an empirical distribution function with weightsw at each ofx

Usage

.wecdf(x, w)

Arguments

Value

a function of classecdf

Transform variogram matrix to covariance of conditional random field

Description

The matrixLambda is half the semivariogram matrix. The functionreturns the conditional covariance with respect to entries inco,restricted to thesubA rows and thesubB columns.

Usage

Lambda2cov(Lambda, co, subA, subB)

Arguments

Value

a covariance matrix

Score and likelihood ratio tests fit of equality of shape over multiple thresholds

Description

The function returns a P-value path for the score testand/or likelihood ratiotest for equality of the shape parameters overmultiple thresholds under the generalized Pareto model.

Usage

NC.diag(  xdat,  u,  GP.fit = c("Grimshaw", "nlm", "optim", "ismev"),  do.LRT = FALSE,  size = NULL,  plot = TRUE,  ...,  xi.tol = 0.001)

Arguments

Details

The default method is'Grimshaw' using the reduction of the parameters to a one-dimensionalmaximization. Other options are one-dimensional maximization of the profile thenlm function oroptim.Two-dimensional optimisation using 2D-optimizationismev using the routinefromgpd.fit from theismev library, with the addition of the algebraic gradient.The choice ofGP.fit should make no difference but the options were kept.Warning: the function will not recover from failure of the maximization routine, returning various error messages.

Value

a plot of P-values for the test at the different thresholdsu

Author(s)

Paul J. Northrop and Claire L. Coleman

References

Grimshaw (1993). Computing Maximum Likelihood Estimates for the GeneralizedPareto Distribution,Technometrics,35(2), 185–191.

Northrop & Coleman (2014). Improved threshold diagnostic plots for extreme valueanalyses,Extremes,17(2), 289–303.

Wadsworth & Tawn (2012). Likelihood-based procedures for thresholddiagnostics and uncertainty in extreme value modelling,J. R. Statist. Soc. B,74(3), 543–567.

Examples

## Not run: data(nidd)u <- seq(65,90, by = 1L)NC.diag(nidd, u, size = 0.05)## End(Not run)

Extreme U-statistic Pickands estimator

Description

[Deprecated function]Given a random sample ofn exceedances, the estimatorreturns an estimator of the shape parameter or extremevalue index using a kernel of order 3, based onk largest exceedances ofxdat. Oorschot et al. (2023) parametrize the model in terms ofm=n-k, so thatm=n corresponds to using only the three largest observations.

Usage

PickandsXU(xdat, m)

Arguments

Details

The calculations are based on the recursions provided in Lemma 4.3 of Oorschot et al.

References

Oorschot, J, J. Segers and C. Zhou (2023), Tail inference using extreme U-statistics,Electronic Journal of Statistics, 17(1): 1113-1159.doi:10.1214/23-EJS2129

Stein's vector of weights

Description

Computes a vector of decreasing weights for exceedances.

Usage

Stein_weights(n, gamma = 1)

Arguments

Value

a vector of positive weights of lengthn

References

Stein, M.L. (2023). A weighted composite log-likelihood approach to parametric estimation of the extreme quantiles of a distribution.Extremes26, 469-507 <doi:10.1007/s10687-023-00466-w>

Wadsworth's univariate and bivariate exponential threshold diagnostics

Description

Function to produce diagnostic plots and test statistics for thethreshold diagnostics exploiting structure of maximum likelihood estimatorsbased 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,  transform = TRUE,  ...)

Arguments

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 uniformityof the p-value distribution under theLR test.

There are two options for the bivariate diagnostic: either provide pairwise minimum of marginallyexponentially distributed margins or provide an times 2 matrix with the original data, whichis transformed to exponential margins using the empirical distribution function.

Value

plots of the requested diagnostics and an invisible list with components

• MLE: maximum likelihood estimates from all thresholds

• Cov: joint asymptotic covariance matrix for\xi,\eta or1/\eta.

• WN: values of the white noise process

• LRT: values of the likelihood ratio test statistic vs threshold

• pval:P-value of the likelihood ratio test

• k: final number of thresholds used

• thresh: threshold selected by the likelihood ratio procedure

• qthresh: quantile level of threshold selected by the likelihood ratio procedure

• thresh0: vector of candidate thresholds

• qthresh0: quantile level of candidate thresholds

• mle.u: maximum likelihood estimates for the selected threshold

• model: model 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 onlyW.diag(xdat = abs(rnorm(5000)), model = 'nhpp',       k = 30, q1 = 0, plots = "PS")W.diag(rexp(1000), model = 'nhpp', k = 20, q1 = 0)xbvn <- rmnorm(n = 6000,                mu = rep(0, 2),                Sigma = cbind(c(1, 0.7), c(0.7, 1)))# Transform margins to exponential manuallyxbvn.exp <- -log(1 - pnorm(xbvn))#rate parametrizationW.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 parametrizationW.diag(xdat = apply(xbvn.exp, 1, min), model = 'invexp',       k = 30, q1 = 0)## End(Not run)

Abisko rainfall

Description

Daily non-zero rainfall measurements in Abisko (Sweden) from January 1913 until December 2014.

Arguments

Format

a data frame with 15132 rows and two variables

Source

Abisko Scientific Research Station

References

A. Kiriliouk, H. Rootzen, J. Segers and J.L. Wadsworth (2019),Peaks over thresholds modeling with multivariate generalized Pareto distributions, Technometrics,61(1), 123–135, <doi:10.1080/00401706.2018.1462738>

Estimation of the bivariate angular dependence function

Description

Estimation of the bivariate angular dependence function

Usage

adf(  xdat,  qlev = 0.95,  estimator = c("hill", "mle", "bayes"),  level = 0.95,  ties.method = "random",  angles = seq(0, 1, by = 0.02),  plot = TRUE)

Arguments

Value

a plot of the angular dependence function ifplot=TRUE, plus an invisible list with components

• angle the sequence of angles in (0,1) at which thelambda values are evaluated

• coef point estimates of the angular dependence function

• lowerlevel% confidence interval for lambda (lower bound)

• upperlevel% confidence interval for lambda (upper bound)

References

J.L. Wadsworth and J.A. Tawn (2013). A new representation for multivariate tail probabilities,Bernoulli, 19(5B), 2689-2714.

Examples

set.seed(12)dat <- mev::rmev(n = 1000, d = 2, model = "log", param = 0.1)adf(xdat = dat, estimator = 'hill')

Bivariate angular dependence function for extrapolation based on rays

Description

The scale parameterg(w) in the Ledford and Tawn approach is estimated empirically forx large as

\hat{g}(w) = \frac{\Pr(X_P>xw, Y_P>x(1-w))}{\Pr(X_P>x, Y_P>x)}

where the sample (X_P, Y_P) are observations on the unit Pareto scale.The coefficient\eta is estimated using maximum likelihood as theshape parameter of a generalized Pareto distribution on\min(X_P, Y_P).

Usage

angextrapo(xdat, qlev = 0.95, w = seq(0.05, 0.95, length = 20), ...)

Arguments

Value

a list with elements

• w: angles between zero and one

• g: scale function at a given value ofw

• eta: Ledford and Tawn tail dependence coefficient

References

Ledford, A.W. and J. A. Tawn (1996), Statistics for near independence in multivariate extreme values.Biometrika,83(1), 169–187.

Examples

angextrapo(rmev(n = 1000, model = 'log', d = 2, param = 0.5))

Rank-based transformation to angular measure

Description

The method uses the pseudo-polar transformation for suitable norms, transformingthe data to pseudo-observations, than marginally to unit Frechet or unit Pareto.Empirical or Euclidean weights are computed and returned alongside with the angular andradial sample for values above threshold(s)thresh, specified in terms of quantilesof the radial componentR or marginal quantiles. Only complete tuples are kept.

Usage

angmeas(  xdat,  thresh,  Rnorm = c("l1", "l2", "linf"),  Anorm = c("l1", "l2", "linf", "arctan"),  marg = c("frechet", "pareto"),  wgt = c("empirical", "euclidean"),  region = c("sum", "min", "max"),  is.angle = FALSE,  ...)

Arguments

Details

The empirical likelihood weighted mean problem is implemented for all thresholds,while the Euclidean likelihood is only supported for diagonal thresholds specifiedviaregion=sum.

Value

a list with argumentsang for thed-1 pseudo-angular sample,rad with the radial componentand possiblywts ifRnorm='l1' and the empirical likelihood algorithm converged. The Euclidean algorithm always returns weights even if some of these are negative.

a list with components

• ang matrix of pseudo-angular observations

• rad vector of radial contributions

• wts empirical or Euclidean likelihood weights for angular observations

Author(s)

Leo Belzile

References

Einmahl, J.H.J. and J. Segers (2009). Maximum empirical likelihood estimation of the spectral measure of an extreme-value distribution,Annals of Statistics,37(5B), 2953–2989.

de Carvalho, M. and B. Oumow and J. Segers and M. Warchol (2013). A Euclidean likelihood estimator for bivariate tail dependence,Comm. Statist. Theory Methods,42(7), 1176–1192.

Owen, A.B. (2001).Empirical Likelihood, CRC Press, 304p.

Examples

x <- rmev(n = 25, d = 3, param = 0.5, model = 'log')wts <- angmeas(xdat = x, Rnorm = 'l1', Anorm = 'l1', marg = 'frechet', wgt = 'empirical')wts2 <- angmeas(xdat = x, Rnorm = 'l2', Anorm = 'l2', marg = 'pareto')

Dirichlet mixture smoothing of the angular measure

Description

This function computes the empirical or Euclidean likelihoodestimates of the spectral measure and uses the points returned from a call toangmeas to compute the Dirichletmixture smoothing of de Carvalho, Warchol and Segers (2012), placing a Dirichlet kernel at each observation.

Usage

angmeasdir(  xdat,  thresh,  Rnorm = c("l1", "l2", "linf"),  Anorm = c("l1", "l2", "linf", "arctan"),  marg = c("frechet", "pareto"),  wgt = c("empirical", "euclidean"),  region = c("sum", "min", "max"),  is.angle = FALSE,  ...)

Arguments

Details

The cross-validationbandwidth is the solution of

\max_{\nu} \sum_{i=1}^n \log \left\{ \sum_{k=1,k \neq i}^n p_{k, -i} f(\mathbf{w}_i; \nu \mathbf{w}_k)\right\},

wheref is the density of the Dirichlet distribution,p_{k, -i} is the Euclidean weightobtained from estimating the Euclidean likelihood problem without observationi.

Value

an invisible list with components

• nu bandwidth parameter obtained by cross-validation;

• dirparmatn byd matrix of Dirichlet parameters for the mixtures;

• wts mixture weights.

Examples

set.seed(123)x <- rmev(n = 100, d = 2L, param = 0.5, model = 'log')out <- angmeasdir(x)

Automated mean residual life plots

Description

This function implements the automated proposal fromSection 2.2 of Langousis et al. (2016)for mean residual life plots. It returns the thresholdthat minimize the weighted mean square error andmoment estimators for the scale and shape parameterbased on weighted least squares.

Usage

automrl(xdat, kmax, thresh, plot = TRUE, ...)

Arguments

Details

The procedure consists in estimating the usual

Value

a list containing

• thresh: selected threshold

• scale: scale parameter estimate

• shape: shape parameter estimate

References

Langousis, A., A. Mamalakis, M. Puliga and R. Deidda (2016).Threshold detection for the generalized Pareto distribution:Review of representative methods and application to theNOAA NCDC daily rainfall database, Water Resources Research,52, 2659–2681.

Censored likelihood for multivariate peaks over threshold models

Description

Censored likelihoods for various parametric limiting models over region determined by

\{y \in F: \max_{j=1}^D \sigma_j \frac{y^\xi_j-1}{\xi_j}+\mu_j > u\};

where\mu isloc,\sigma isscale and\xi isshape.

Usage

clikmgp(  dat,  thresh,  mthresh = thresh,  loc,  scale,  shape,  par,  model = c("log", "neglog", "br", "xstud"),  likt = c("mgp", "pois", "binom"),  lambdau = 1,  ...)

Arguments

Details

Optional arguments can be passed to the function via...

• censored matrix of booleans andNA indicating whether observationsdat fall below the mthresholdmthresh

• cl cluster instance created bymakeCluster (default toNULL)

• ncors number of cores for parallel computing of the likelihood

• numAbovePerRow number of observations above mthreshold (non-missing) per row

• numAbovePerCol number of observations above mthreshold (non-missing) per column

• mmax maximum per column

• B1 number of replicates for quasi Monte Carlo integral for the exponent measure

• B2 number of replicates for quasi Monte Carlo integral for the censored intensity contribution

• genvec1 generating vector for the quasi Monte Carlo routine (exponent measure), associated withB1

• genvec2 generating vector for the quasi Monte Carlo routine (individual obs contrib), associated withB2

Value

the value of the log-likelihood withattributesexpme, giving the exponent measure

Note

The location and scale parameters are not identifiable unless one of them is fixed.

Confidence intervals for profile likelihood objects

Description

Computes confidence intervals for the parameter psi for profile likelihood objects.This function uses spline interpolation to derivelevel confidence intervals

Usage

## S3 method for class 'eprof'confint(  object,  parm,  level = 0.95,  prob = c((1 - level)/2, 1 - (1 - level)/2),  print = FALSE,  method = c("cobs", "smooth.spline"),  boundary = FALSE,  ...)

Arguments

Value

returns a 2 by 3 matrix containing point estimates, lower and upper confidence intervals based on the likelihood root and modified version thereof

Threshold selection via coefficient of variation

Description

This function computes the empirical coefficient of variation andcomputes a weighted statistic comparing the squared distance withthe theoretical coefficient variation corresponding to a specificshape parameter (estimated from the data using a moment estimatoras the value minimizing the test statistic, or using maximum likelihood).The procedure stops if there are no more than 10 exceedances above thehighest threshold

Usage

cvselect(  xdat,  thresh,  method = c("mle", "wcv", "cv"),  nsim = 999L,  nthresh = 10L,  level = 0.05,  lazy = FALSE)

Arguments

Value

a list with elements

• thresh0: value of threshold returned by the procedure,NA if the hypothesis is rejected at all thresholds

• thresh: sorted vector of candidate thresholds

• cindex: index of selected threshold amongthresh orNA if none returned

• pval: bootstrap p-values, withNA iflazy and the p-value exceeds level at lower thresholds

• shape: shape parameter estimates

• nexc: number of exceedances of each thresholdthresh

• method: estimation method for the shape parameter

References

del Castillo, J. and M. Padilla (2016).Modelling extreme values by the residual coefficient of variation, SORT, 40(2), pp. 303–320.

Distance matrix with geometric anisotropy

Description

The function computes the distance between locations, with geometric anisotropy.Consider real parameters\theta_1 and\theta_2, and the transformation\psi=\arctan(\theta_1/\theta_2)/2 andr=1 +\theta_1^2 + \theta_2^2.The dilation and rotation matrix is

\left(\begin{matrix} \sqrt{r}\cos(\rho) & -\sqrt{r}\sin(\rho) \\ \sin(\rho)/\sqrt{r} & \cos(\rho)/\sqrt{r} \end{matrix} \right).

The parametrization is convenient for optimization purposes, as the parameter vector is unconstrainedand the transformation has unit Jacobian.

Usage

dgeoaniso(loc, theta)

Arguments

Value

ad byd square matrix of pairwise distance

References

Rai, K. and Brown, P.E. (2025), A parameter transformation of the anisotropic Matérn covariance function. Canadian Journal of Statistics e11839.doi:10.1002/cjs.11839

Distance matrix with geometric anisotropy

Description

The function computes the distance between locations, with geometric anisotropy.The parametrization assumes there is a scale parameter, say\sigma, so thatscaleis the distortion for the second component only. The anglerho must lie in[-\pi/2, \pi/2]. The dilation and rotation matrix is

\left(\begin{matrix} \cos(\rho) & \sin(\rho) \\ - \sigma\sin(\rho) & \sigma\cos(\rho) \end{matrix} \right)

Usage

distg(loc, scale, rho)

Arguments

Value

ad byd square matrix of pairwise distance

Extended generalised Pareto families

Description

This function provides the log-likelihood and quantiles for the three different families presentedin Papastathopoulos and Tawn (2013) and the two proposals of Gamet and Jalbert (2022), plus exponential tilting. All of the models contain an additional parameter,\kappa \ge 0.All families share the same tail index as the generalized Pareto distribution, while allowing for lower thresholds.For most models, the distribution reduce to the generalised Pareto when\kappa=1 (for modelsgj-tnorm andlogist, on the boundary of the parameter space when\kappa \to 0).

egp.retlev gives the return levels for the extended generalised Pareto distributions

Arguments

Details

For return levels, thep argument can be related toT year exceedances as follows:if there aren_y observations per year, than takepto equal1/(Tn_y) to obtain theT-years return level.

Value

egp.ll returns the log-likelihood value, whileegp.retlev returns a plot of the return levels ifplot=TRUE and a list with tail probabilitiesp, return levelsretlev, thresholdsthresh and model namemodel.

Usage

egp.ll(xdat, thresh, model, par)

egp.retlev(xdat, thresh, par, model, p, plot=TRUE)

Author(s)

Leo Belzile

References

Papastathopoulos, I. and J. Tawn (2013). Extended generalised Pareto models for tail estimation,Journal of Statistical Planning and Inference143(3), 131–143, <doi:10.1016/j.jspi.2012.07.001>.

Gamet, P. and Jalbert, J. (2022). A flexible extended generalized Pareto distribution for tail estimation.Environmetrics,33(6), <doi:10.1002/env.2744>.

Examples

set.seed(123)xdat <- rgp(1000, loc = 0, scale = 2, shape = 0.5)par <- fit.egp(xdat, thresh = 0, model = 'gj-beta')$parp <- c(1/1000, 1/1500, 1/2000)# With multiple thresholdsth <- c(0, 0.1, 0.2, 1)opt <- tstab.egp(xdat, thresh = th, model = 'gj-beta')egp.retlev(xdat = xdat, thresh = th, model = 'gj-beta', p = p)opt <- tstab.egp(xdat, th, model = 'pt-power', plots = NA)egp.retlev(xdat = xdat, thresh = th, model = 'pt-power', p = p)

Extended generalised Pareto log likelihood

Description

This function provides the log-likelihood and quantiles for the different extended generalized Pareto distributions. All families share the same tail index as the generalized Pareto, and reduce to the latter whenkappa=1.

Usage

egp.ll(  xdat,  thresh,  par,  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",    "logist"))egp.retlev(  xdat,  thresh,  par,  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",    "logist"),  p,  confint = FALSE,  plot = TRUE)

Arguments

Fit of extended GP models and parameter stability plots

Description

This function is an alias offit.egp.

Usage

egp.fit(  xdat,  thresh,  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",    "logist"),  init,  show = FALSE)

Details

Supported for backward compatibility

Deprecated function for parameter stability plots

Description

Deprecated function for parameter stability plots

Usage

egp.fitrange(  xdat,  thresh,  model = c("egp1", "egp2", "egp3"),  plots = 1:3,  umin,  umax,  nint)

Profile log likelihood for extended generalized Pareto models

Description

Computes the profile log likelihood over a grid of values of\psi for various parameters, including return levels.

Usage

egp.pll(  psi,  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",    "logist"),  param = c("kappa", "scale", "shape", "retlev"),  mle = NULL,  xdat,  thresh = NULL,  plot = FALSE,  method = c("Nelder", "nlminb", "BFGS"),  p,  ...)

Arguments

Value

an object of classeprof

Fit an extended generalized Pareto distribution of Naveau et al.

Description

Deprecated function name to fit an extended generalized Pareto family. The user should callfit.extgp instead.

Usage

egp2.fit(  data,  model = 1,  method = c("mle", "pwm"),  init,  censoring = c(0, Inf),  rounded = 0,  CI = FALSE,  R = 1000,  ncpus = 1,  plots = TRUE)

Arguments

See Also

fit.extgp

Extended generalized Pareto distribution

Description

Density function, distribution function, quantile function andrandom number generation for various extended generalized Paretodistributions

Usage

pegp(  q,  scale,  shape,  kappa,  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",    "logist"),  lower.tail = TRUE,  log.p = FALSE)degp(  x,  scale,  shape,  kappa,  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",    "logist"),  log = FALSE)qegp(  p,  scale,  shape,  kappa,  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",    "logist"),  lower.tail = TRUE,  log.p = FALSE)regp(  n,  scale,  shape,  kappa,  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",    "logist"))

Arguments

References

Papastathopoulos, I. and J. Tawn (2013). Extended generalised Pareto models for tail estimation,Journal of Statistical Planning and Inference143(3), 131–143, <doi:10.1016/j.jspi.2012.07.001>.

Gamet, P. and Jalbert, J. (2022). A flexible extended generalized Pareto distribution for tail estimation.Environmetrics,33(6), <doi:10.1002/env.2744>.

Self-concordant empirical likelihood for a vector mean

Description

Self-concordant empirical likelihood for a vector mean

Usage

emplik(  dat,  mu = rep(0, ncol(dat)),  lam = rep(0, ncol(dat)),  eps = 1/nrow(dat),  M = 1e+30,  thresh = 1e-30,  itermax = 100)

Arguments

Value

a list with components

• logelr log empirical likelihood ratio.

• lam Lagrange multiplier (vector of lengthd).

• wtsn vector of observation weights (probabilities).

• conv boolean indicating convergence.

• niter number of iteration until convergence.

• ndec Newton decrement.

• gradnorm norm of gradient of log empirical likelihood.

Author(s)

Art Owen,C++ port by Leo Belzile

References

Owen, A.B. (2013). Self-concordance for empirical likelihood,Canadian Journal of Statistics,41(3), 387–397.

Eskdalemuir Observatory Daily Rainfall

Description

This dataset contains exceedances of 30mm for dailycumulated rainfall observations over the period 1970-1986.These data were aggregated from hourly series.

Format

a vector with 93 daily cumulated rainfall measurements exceeding 30mm.

Details

The station is one of the rainiest of the whole UK, with an average 1554m of cumulated rainfall per year.The data consisted of 6209 daily observations, of which 4409 were non-zero.Only the 93 largest observations are provided.

Source

Met Office.

Exponent measure for multivariate generalized Pareto distributions

Description

Integrated intensity over the region defined by[0, z]^c for logistic, Huesler-Reiss, Brown-Resnick and extremal Student processes.

Usage

expme(  z,  par,  model = c("log", "neglog", "hr", "br", "xstud"),  method = c("TruncatedNormal", "mvtnorm", "mvPot"))

Arguments

Value

numeric giving the measure of the complement of[0,z].

Note

The listpar must contain different arguments depending on the model. For the Brown–Resnick model, the user must supply the conditionally negative definite matrixLambda following the parametrization in Engelkeet al. (2015) or the covariance matrixSigma, following Wadsworth and Tawn (2014). For the Husler–Reiss model, the user provides the mean and covariance matrix,m andSigma. For the extremal student, the covariance matrixSigma and the degrees of freedomdf. For the logistic model, the strictly positive dependence parameteralpha.

Examples

## Not run: # Extremal StudentSigma <- stats::rWishart(n = 1, df = 20, Sigma = diag(10))[, , 1]expme(z = rep(1, ncol(Sigma)), par = list(Sigma = cov2cor(Sigma), df = 3), model = "xstud")# Brown-Resnick modelD <- 5Lloc <- cbind(runif(D), runif(D))di <- as.matrix(dist(rbind(c(0, ncol(loc)), loc)))semivario <- function(d, alpha = 1.5, lambda = 1) {  (d / lambda)^alpha}Vmat <- semivario(di)Lambda <- Vmat[-1, -1] / 2expme(z = rep(1, ncol(Lambda)), par = list(Lambda = Lambda), model = "br", method = "mvPot")Sigma <- outer(Vmat[-1, 1], Vmat[1, -1], "+") - Vmat[-1, -1]expme(z = rep(1, ncol(Lambda)), par = list(Lambda = Lambda), model = "br", method = "mvPot")## End(Not run)

Extremal index estimators

Description

The function implements estimators of the extremal index based on interexceedance time and gap of exceedances. The maximum likelihood estimator and iteratively reweighted leastsquare estimators of Suveges (2007) as well as the intervals estimator. The implementationdiffers from the presentation of the paper in that an iteration limit is enforced to make surethe iterative procedure terminates. Multiple thresholds can be supplied.

Usage

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

Arguments

Details

The iteratively reweighted least square is a procedure based on the gaps of exceedancesS_n=T_n-1The model is first fitted to non-zero gaps, which are rescaled to have unit exponential scale. The slopebetween the theoretical quantiles and the normalized gap of exceedances isb=1/\theta,with intercepta=\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 forthe 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 conditionD^{(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.

Value

a vector or matrix of estimated extremal index of dimensionlength(method) bylength(q).

Author(s)

Leo Belzile

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 climatologyfor extreme hourly precipitation in Switzerland.Theoretical and Applied Climatology,120(3), 403-416.

Examples

## Not run: set.seed(234)# Moving maxima model with theta=0.5a <- 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, quantile = q, method = c ('wls', 'mle'))## End(Not run)

Estimators of the extremal coefficient

Description

These functions estimate the extremal coefficient using an approximatesample from the Frechet distribution.

Usage

extcoef(  xdat,  coord = NULL,  thresh = NULL,  estimator = c("schlather", "smith", "fmado"),  margtrans = c("emp", "gev", "none"),  ties.method = "random",  prob = 0,  plot = TRUE,  ...)

Arguments

Details

TheSmith estimator: supposeZ(x) is simple max-stable vector (i.e., with unit Frechet marginals). Then1/Z is unit exponential and1/\max(Z(s_1), Z(s_2))is exponential with rate\theta = \max\{Z(s_1), Z(s_2)\}.The extremal index for the pair can therefore be calculated usingthe reciprocal mean.

TheSchlather and Tawn estimator: the likelihood of the naive estimator for a pairof two sitesA is

\mathrm{card}\left\{ j: \max_{i \in A} X_i^{(j)}\bar{X}_i)>z \right\}\log(\theta_A) - \theta_A \sum_{j=1}^n \left[ \max \left\{z, \max_{i \in A}(X_i^{(j)}\bar{X}_i)\right\}\right]^{-1},

where\bar{X}_i = n^{-1} \sum_{j=1}^n 1/X_i^{(j)} is the harmonic mean andzis a threshold on the unit Frechet scale.The search for the maximum likelihood estimate for every pairAis restricted to the interval[1,3].A binned version of the extremal coefficient cloud is also returned.The Schlather estimator is not self-consistent.The Schlather and Tawn estimator includes as special casethe Smith estimator if we do not censor the data (p = 0)and do not standardize observations by their harmonic mean.

TheF-madogram estimator is a non-parametric estimate based on a stationary processZ; the extremal coefficient satisfies

\theta(h)=\frac{1+2\nu(h)}{1-2\nu(h)},

where

\nu(h) = \frac{1}{2} \mathsf{E}[|F(Z(s+h)-F(Z(s))|]

The implementation only uses complete pairs to calculate the relative ranks.

All estimators are coded in plain R and computations are not optimized.The estimation time can therefore be large for large data sets.If there are no missing observations, the routinefmadogramfrom theSpatialExtremes package should be preferred as it isnoticeably faster.

The data will typically consist of max-stable vectors or block maxima.Both of the Smith and the Schlather–Tawn estimators require unitFrechet margins; the margins will be standardizedto the unit Frechet scale, either parametrically or nonparametrically unlessmargtrans="none".Ifmargtrans = "gev", a parametric GEV model is fitted to each column ofdat using maximum likelihoodestimation and transformed back using the probability integral transform. Ifmethod = "emp",using the empirical distribution function. The latter is the default, as it is appreciably faster.

Value

an invisible list with vectorsdist ifcoord is non-null or else a matrix of pairwise indicesind,extcoef and the suppliedestimator,fmado andbinned. Ifestimator == "schlather", an additional matrix with 2 columns containing the binned distancebinned with theh and the binned extremal coefficient.

References

Schlather, M. and J. Tawn (2003). A dependence measure for multivariate and spatial extremes,Biometrika,90(1), pp.139–156.

Cooley, D., P. Naveau and P. Poncet (2006). Variograms for spatial max-stable random fields, In: Bertail P., Soulier P., Doukhan P. (eds)Dependence in Probability and Statistics. Lecture Notes in Statistics, vol. 187. Springer, New York, NY

R. J. Erhardt, R. L. Smith (2012), Approximate Bayesian computing for spatial extremes,Computational Statistics and Data Analysis,56, pp.1468–1481.

Examples

## Not run: coord <- 10 * cbind(runif(50), runif(50))di <- as.matrix(dist(coord))dat <- rmev(  n = 1000,  d = 100,  param = 3,  sigma = exp(-di / 2),  model = 'xstud')res <- extcoef(xdat = dat, coord = coord)# Extremal Student extremal coefficient functionXT.extcoeffun <- function(h, nu, corrfun, ...) {  if (!is.function(corrfun)) {    stop('Invalid function \"corrfun\".')  }  h <- unique(as.vector(h))  rhoh <- sapply(h, corrfun, ...)  cbind(    h = h,    extcoef = 2 * pt(sqrt((nu + 1) * (1 - rhoh) / (1 + rhoh)), nu + 1)  )}#This time, only one graph with theoretical extremal coefplot(res$dist, res$extcoef, ylim = c(1, 2), pch = 20)abline(v = 2, col = 'gray')extcoefxt <- XT.extcoeffun(  seq(0, 10, by = 0.1),  nu = 3,  corrfun = function(x) {    exp(-x / 2)  })lines(  extcoefxt[, 'h'],  extcoefxt[, 'extcoef'],  type = 'l',  col = 'blue',  lwd = 2)# Brown--Resnick extremal coefficient functionBR.extcoeffun <- function(h, vario, ...) {  if (!is.function(vario)) {    stop('Invalid function \"vario\".')  }  h <- unique(as.vector(h))  gammah <- sapply(h, vario, ...)  cbind(h = h, extcoef = 2 * pnorm(sqrt(gammah / 4)))}extcoefbr <- BR.extcoeffun(  seq(0, 20, by = 0.25),  vario = function(x) {    2 * x^0.7  })lines(  extcoefbr[, 'h'],  extcoefbr[, 'extcoef'],  type = 'l',  col = 'orange',  lwd = 2)coord <- 10 * cbind(runif(20), runif(20))di <- as.matrix(dist(coord))dat <- rmev(  n = 1000,  d = 20,  param = 3,  sigma = exp(-di / 2),  model = 'xstud')res <- extcoef(  xdat = dat,  coord = coord,  estimator = "smith")## End(Not run)

Extended generalised Pareto families of Naveau et al. (2016)

Description

Density function, distribution function, quantile function and random generation for the extended generalized Pareto distribution (GPD) with scale and shape parameters.

Arguments

Details

The extended generalized Pareto families proposed in Naveauet al. (2016)retain the tail index of the distribution while being compliant with the theoretical behavior of extremelow rainfall. There are five proposals, the first one being equivalent to the GP distribution.

• type 0 corresponds to uniform carrier,G(u)=u.

• type 1 corresponds to a three parameters family, with carrierG(u)=u^\kappa.

• type 2 corresponds to a three parameters family, with carrierG(u)=1-V_\delta((1-u)^\delta).

• type 3 corresponds to a four parameters family, with carrier

  G(u)=1-V_\delta((1-u)^\delta))^{\kappa/2}

  .

• type 4 corresponds to a five parameter model (a mixture oftype 2, withG(u)=pu^\kappa + (1-p)*u^\delta

Usage

pextgp(q, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1)

dextgp(x, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1, log=FALSE)

qextgp(p, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1)

rextgp(n, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1, unifsamp=NULL, censoring=c(0,Inf))

Author(s)

Raphael Huser and Philippe Naveau

References

Naveau, P., R. Huser, P. Ribereau, and A. Hannart (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection,Water Resour. Res., 52, 2753-2769,doi:10.1002/2015WR018552.

Carrier distribution for the extended GP distributions of Naveau et al.

Description

Density, distribution function, quantile function and random numbergeneration for the carrier distributions of the extended Generalized Pareto distributions.

Arguments

Usage

pextgp.G(u, type=1, prob, kappa, delta)

dextgp.G(u, type=1, prob=NA, kappa=NA, delta=NA, log=FALSE)

qextgp.G(u, type=1, prob=NA, kappa=NA, delta=NA)

rextgp.G(n, prob=NA, kappa=NA, delta=NA,type=1, unifsamp=NULL, direct=FALSE, censoring=c(0,1))

Author(s)

Raphael Huser and Philippe Naveau

See Also

extgp

Pairwise extremogram for max-risk functional

Description

The function computes the pairwisechi estimates and plots them as a function of the distance between sites. The function also includes utilities for geometric anisotropy.

Usage

extremo(xdat, margp, coord, scale = 1, rho = 0, plot = FALSE, ...)

Arguments

Value

an invisible matrix with pairwise estimates of chi along with distance (unsorted)

Examples

## Not run: lon <- seq(650, 720, length = 10)lat <- seq(215, 290, length = 10)# Create a gridgrid <- expand.grid(lon,lat)coord <- as.matrix(grid)dianiso <- distg(coord, 1.5, 0.5)sgrid <- scale(grid, scale = FALSE)# Specify marginal parameters `loc` and `scale` over grideta <- 26 + 0.05*sgrid[,1] - 0.16*sgrid[,2]tau <- 9 + 0.05*sgrid[,1] - 0.04*sgrid[,2]# Parameter matrix of Huesler--Reiss# associated to power variogramLambda <- ((dianiso/30)^0.7)/4# Regular Euclidean distance between sitesdi <- distg(coord, 1, 0)# Simulate generalized max-Pareto fieldset.seed(345)simu1 <- rgparp(n = 1000, thresh = 50, shape = 0.1, riskf = "max",                scale = tau, loc = eta, sigma = Lambda, model = "hr")extdat <- extremo(dat = simu1, margp = 0.98, coord = coord,                  scale = 1.5, rho = 0.5, plot = TRUE)# Constrained optimization# Minimize distance between extremal coefficient from fitted variogrammindistpvario <- function(par, emp, coord){alpha <- par[1]; if(!isTRUE(all(alpha > 0, alpha < 2))){return(1e10)}scale <- par[2]; if(scale <= 0){return(1e10)}a <- par[3]; if(a<1){return(1e10)}rho <- par[4]; if(abs(rho) >= pi/2){return(1e10)}semivariomat <- power.vario(distg(coord, a, rho), alpha = alpha, scale = scale)  sum((2*(1-pnorm(sqrt(semivariomat[lower.tri(semivariomat)]/2))) - emp)^2)}hin <- function(par, ...){  c(1.99-par[1], -1e-5 + par[1],    -1e-5 + par[2],    par[3]-1,    pi/2 - par[4],    par[4]+pi/2)  }opt <- alabama::auglag(  par = c(0.7, 30, 1, 0),  hin = hin,  control.outer = list(trace = 0),  fn = function(par){   mindistpvario(par, emp = extdat$prob, coord = coord)})stopifnot(opt$kkt1, opt$kkt2)# Plotting the extremogram in the deformed spacedistfa <- distg(loc = coord, opt$par[3], opt$par[4])plot( x = c(distfa[lower.tri(distfa)]), y = extdat$prob, pch = 20, yaxs = "i", xaxs = "i", bty = 'l', xlab = "distance", ylab= "cond. prob. of exceedance", ylim = c(0,1))lines(  x = (distvec <- seq(0,200, length = 1000)),  y = 2*(1-pnorm(sqrt(power.vario(distvec, alpha = opt$par[1],                               scale = opt$par[2])/2))),  col = 2,  lwd = 2)## End(Not run)

Parameter stability plot and maximum likelihood routine for extended GP models

Description

The functiontstab.egp provides classical threshold stability plot for (\kappa,\sigma,\xi).The fitted parameter values are displayed with pointwise normal 95% confidence intervals.The function returns an invisible list with parameter estimates and standard errors, and p-values for the Wald test that\kappa=1.The plot is for the modified scale (as in the generalised Pareto model) and as such it is possible that the modified scale be negative.tstab.egp can also be used to fit the model to multiple thresholds.

Usage

fit.egp(  xdat,  thresh = 0,  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",    "logist"),  start = NULL,  method = c("Nelder", "nlminb", "BFGS"),  fpar = NULL,  show = FALSE,  ...)

Arguments

Details

fit.egp is a numerical optimization routine to fit the extended generalised Pareto models of Papastathopoulos and Tawn (2013),using maximum likelihood estimation.

Value

fit.egp outputs the list returned byoptim, which contains the parameter values, the hessian and in addition the standard errors

tstab.egp returns a plot(s) of the parameters fit over the range of provided thresholds, with pointwise normal confidence intervals; the function also returns an invisible list containing notably the matrix of point estimates (par) and standard errors (se).

Author(s)

Leo Belzile

References

Papastathopoulos, I. and J. Tawn (2013). Extended generalised Pareto models for tail estimation,Journal of Statistical Planning and Inference143(3), 131–143.

Examples

xdat <- mev::rgp(  n = 100,  loc = 0,  scale = 1,  shape = 0.5)fitted <- fit.egp(  xdat = xdat,  thresh = 1,  model = "pt-gamma",  show = TRUE)thresh <- mev::qgp(seq(0.1, 0.5, by = 0.05), 0, 1, 0.5)tstab.egp(   xdat = xdat,   thresh = thresh,   model = "pt-gamma")xdat <- regp(  n = 100,  scale = 1,  shape = 0.1,  kappa = 0.5,  model = "pt-power")fit.egp( xdat = xdat, model = "pt-power", show = TRUE, fpar = list(kappa = 1), method = "Nelder")

Fit an extended generalized Pareto distribution of Naveau et al.

Description

This is a wrapper function to obtain PWM or MLE estimates forthe extended GP models of Naveau et al. (2016) for rainfall intensities. The function calculates confidence intervalsby means of nonparametric percentile bootstrap and returns histograms and QQ plots ofthe fitted distributions. The function handles both censoring and rounding.

Usage

fit.extgp(  data,  model = 1,  method = c("mle", "pwm"),  init,  censoring = c(0, Inf),  rounded = 0,  confint = FALSE,  R = 1000,  ncpus = 1,  plots = TRUE)

Arguments

Details

The different models include the following transformations:

• model 0 corresponds to uniform carrier,G(u)=u.

• model 1 corresponds to a three parameters family, with carrierG(u)=u^\kappa.

• model 2 corresponds to a three parameters family, with carrierG(u)=1-V_\delta((1-u)^\delta).

• model 3 corresponds to a four parameters family, with carrier

  G(u)=1-V_\delta((1-u)^\delta))^{\kappa/2}

  .

• model 4 corresponds to a five parameter model (a mixture oftype 2, withG(u)=pu^\kappa + (1-p)*u^\delta

Author(s)

Raphael Huser and Philippe Naveau

References

Naveau, P., R. Huser, P. Ribereau, and A. Hannart (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection,Water Resour. Res., 52, 2753-2769,doi:10.1002/2015WR018552.

See Also

egp.fit,egp,extgp

Examples

## Not run: data(rain, package = "ismev")fit.extgp(  rain[rain > 0],  model = 1,  method = 'mle',  init = c(0.9, fit.gpd(rain)$est),  rounded = 0.1,  confint = TRUE,  R = 20)## End(Not run)

Maximum likelihood estimation for the generalized extreme value distribution

Description

This function returns an object of classmev_gev, with default methods for printing and quantile-quantile plots.The default starting values are the solution of the probability weighted moments.

Usage

fit.gev(  xdat,  start = NULL,  method = c("nlminb", "BFGS"),  show = FALSE,  fpar = NULL,  warnSE = FALSE)

Arguments

Value

a list containing the following components:

• estimate a vector containing the maximum likelihood estimates.

• std.err a vector containing the standard errors.

• vcov the variance covariance matrix, obtained as the numerical inverse of the observed information matrix.

• method the method used to fit the parameter.

• nllh the negative log-likelihood evaluated at the parameterestimate.

• convergence components taken from the list returned byauglag.Values other than0 indicate that the algorithm likely did not converge.

• counts components taken from the list returned byauglag.

• xdat vector of data

Examples

xdat <- mev::rgev(n = 100)fit.gev(xdat, show = TRUE)# Example with fixed parameterfit.gev(xdat, show = TRUE, fpar = list(shape = 0))

Maximum likelihood estimation for the generalized Pareto distribution

Description

Numerical optimization of the generalized Pareto distribution fordata exceedingthreshold.This function returns an object of classmev_gpd, with default methods for printing and quantile-quantile plots.

Usage

fit.gpd(  xdat,  threshold = 0,  method = "Grimshaw",  show = FALSE,  MCMC = NULL,  k = 4,  tol = 1e-08,  fpar = NULL,  warnSE = FALSE)

Arguments

Details

The default method is'Grimshaw', which maximizes the profile likelihood for the ratio scale/shape. Other options include'obre' for optimalB-robust estimator of the parameter of Dupuis (1998), vanilla maximization of the log-likelihood using constrained optimization routine'auglag', 1-dimensional optimization of the profile likelihood usingnlm andoptim. Method'ismev' performs the two-dimensional optimization routinegpd.fit from theismev library, with in addition the algebraic gradient.The approximate Bayesian methods ('zs' and'zhang') are extracted respectively from Zhang and Stephens (2009) and Zhang (2010) and consists of a approximate posterior mean calculated via importancesampling assuming a GPD prior is placed on the parameter of the profile likelihood.

Value

Ifmethod is neither'zs' nor'zhang', a list containing the following components:

• estimate a vector containing thescale andshape parameters (optimized and fixed).

• std.err a vector containing the standard errors. Formethod = "obre", these are Huber's robust standard errors.

• vcov the variance covariance matrix, obtained as the numerical inverse of the observed information matrix. Formethod = "obre",this is the sandwich Godambe matrix inverse.

• threshold the threshold.

• method the method used to fit the parameter. See details.

• nllh the negative log-likelihood evaluated at the parameterestimate.

• nat number of points lying above the threshold.

• pat proportion of points lying above the threshold.

• convergence components taken from the list returned byoptim.Values other than0 indicate that the algorithm likely did not converge (in particular 1 and 50).

• counts components taken from the list returned byoptim.

• exceedances excess over the threshold.

Additionally, ifmethod = "obre", a vector of OBREweights.

Otherwise, a list containing

• threshold the threshold.

• method the method used to fit the parameter. SeeDetails.

• nat number of points lying above the threshold.

• pat proportion of points lying above the threshold.

• approx.mean a vector containing containing the approximate posterior mean estimates.

and in addition if MCMC is neitherFALSE, norNULL

• post.mean a vector containing the posterior mean estimates.

• post.se a vector containing the posterior standard error estimates.

• accept.rate proportion of points lying above the threshold.

• niter length of resulting Markov Chain

• burnin amount of discarded iterations at start, capped at 10000.

• thin thinning integer parameter describing

Note

Some of the internal functions (which are hidden from the user) allow for modelling of the parameters using covariates. This is not currently implemented withingp.fit, but users can call internal functions should they wish to use these features.

Author(s)

Scott D. Grimshaw for theGrimshaw option. Paul J. Northrop and Claire L. Coleman for the methodsoptim,nlm andismev.J. Zhang and Michael A. Stephens (2009) and Zhang (2010) for thezs andzhang approximate methods and L. Belzile for methodsauglag andobre, the wrapper and MCMC samplers.

Ifshow = TRUE, the optimalB robust estimated weights for the largest observations are printed alongside with thep-value of the latter, obtained from the empirical distribution of the weights. This diagnostic can be used to guide threshold selection:small weights for ther-largest order statistics indicate that the robust fit is driven by the lower tailand that the threshold should perhaps be increased.

References

Davison, A.C. (1984). Modelling excesses over high thresholds, with an application, inStatistical extremes and applications, J. Tiago de Oliveira (editor), D. Reidel Publishing Co., 461–482.

Grimshaw, S.D. (1993). Computing Maximum Likelihood Estimates for the GeneralizedPareto Distribution,Technometrics,35(2), 185–191.

Northrop, P.J. and C. L. Coleman (2014). Improved threshold diagnostic plots for extreme valueanalyses,Extremes,17(2), 289–303.

Zhang, J. (2010). Improving on estimation for the generalized Pareto distribution,Technometrics52(3), 335–339.

Zhang, J. and M. A. Stephens (2009). A new and efficient estimation method for the generalized Pareto distribution.Technometrics51(3), 316–325.

Dupuis, D.J. (1998). Exceedances over High Thresholds: A Guide to Threshold Selection,Extremes,1(3), 251–261.

See Also

fpot andgpd.fit

Examples

data(eskrain)fit.gpd(eskrain, threshold = 35, method = 'Grimshaw', show = TRUE)fit.gpd(eskrain, threshold = 30, method = 'zs', show = TRUE)

Maximum likelihood estimation of the point process of extremes

Description

Data abovethreshold is modelled using the limiting point processof extremes.

Usage

fit.pp(  xdat,  threshold = 0,  npp = 1,  np = NULL,  method = c("nlminb", "BFGS"),  start = NULL,  show = FALSE,  fpar = NULL,  warnSE = FALSE)

Arguments

Details

The parameternpp controls the frequency of observations.If data are recorded on a daily basis, using a value ofnpp = 365.25yields location and scale parameters that correspond to those of thegeneralized extreme value distribution fitted to block maxima.

Value

a list containing the following components:

• estimate a vector containing all parameters (optimized and fixed).

• std.err a vector containing the standard errors.

• vcov the variance covariance matrix, obtained as the numerical inverse of the observed information matrix.

• threshold the threshold.

• method the method used to fit the parameter. See details.

• nllh the negative log-likelihood evaluated at the parameterestimate.

• nat number of points lying above the threshold.

• pat proportion of points lying above the threshold.

• convergence components taken from the list returned byoptim.Values other than0 indicate that the algorithm likely did not converge (in particular 1 and 50).

• counts components taken from the list returned byoptim.

References

Coles, S. (2001), An introduction to statistical modelling of extreme values. Springer : London, 208p.

Examples

data(eskrain)pp_mle <- fit.pp(eskrain, threshold = 30, np = 6201)plot(pp_mle)

Estimator of the second order tail index parameter

Description

Estimator of the second order tail index parameter

Usage

fit.rho(xdat, k, method = c("fagh", "dk", "ghp", "gbw"), ...)

Arguments

Examples

# Example with rho = -0.2n <- 1000xdat <- mev::rgp(n = n, shape = 0.2)kmin <- floor(n^0.995)kmax <- ceiling(n^0.999)rho_est <- fit.rho(   xdat = xdat,   k = n - kmin:kmax)rho_med <- mean(rho_est$rho)

Maximum likelihood estimates of point process for the r-largest observations

Description

This uses a constrained optimization routine to return the maximum likelihood estimatebased on ann byr matrix of observations. Observations should be ordered, i.e.,ther-largest should be in the last column.

Usage

fit.rlarg(  xdat,  start = NULL,  method = c("nlminb", "BFGS"),  show = FALSE,  fpar = NULL,  warnSE = FALSE)

Arguments

Value

a list containing the following components:

• estimate a vector containing all the maximum likelihood estimates.

• std.err a vector containing the standard errors.

• vcov the variance covariance matrix, obtained as the numerical inverse of the observed information matrix.

• method the method used to fit the parameter.

• nllh the negative log-likelihood evaluated at the parameterestimate.

• convergence components taken from the list returned byauglag.Values other than0 indicate that the algorithm likely did not converge.

• counts components taken from the list returned byauglag.

• xdat ann byr matrix of data

Examples

xdat <- rrlarg(n = 10, loc = 0, scale = 1, shape = 0.1, r = 4)fit.rlarg(xdat)

Shape parameter estimates

Description

Wrapper to estimate the tail index or shape parameter of an extreme value distribution. Each function has similar sets of arguments, a vector or scalar number of order statisticsk anda vector of positive observationsxdat. Themethod argument allows users to choose between different indicators, including the Hill estimator (hill, for positive observations and shape only), the moment estimator of Dekkers and de Haan (mom ordekkers), the de Vries estimator of de Haan and Peng (vries), the generalized jackknife estimator of Gomes et al. (genjack), the Beirlant, Vynckier and Teugels generalized quantile estimator (bvt orgenquant), the Pickands estimator (pickands), the extremeU-statistics estimator of Oorschot, Segers and Zhou (osz), or the exponential rgression model of Beirlant et al. (erm).

Usage

fit.shape(  xdat,  k,  method = c("hill", "rbm", "osz", "vries", "genjack", "mom", "dekkers", "genquant",    "pickands", "erm"),  ...)

Arguments

Value

a data frame with the number of order statisticsk and the shape parameter estimateshape, or a single numeric value ifk is a scalar.

Maximum likelihood estimation for weighted generalized Pareto distribution

Description

Weighted maximum likelihood estimation, with user-specified vector of weights.

Usage

fit.wgpd(xdat, threshold = 0, weightfun = Stein_weights, start = NULL, ...)

Arguments

Value

a list with components

• estimate a vector containing thescale andshape parameters (optimized and fixed).

• std.err a vector containing the standard errors.

• vcov the variance covariance matrix, obtained as the numerical inverse of the observed information matrix.

• threshold the threshold.

• method the method used to fit the parameter. See details.

• nllh the negative log-likelihood evaluated at the parameterestimate.

• nat number of points lying above the threshold.

• pat proportion of points lying above the threshold.

• convergence logical indicator of convergence.

• weights vector of weights for exceedances.

• exceedances excess over the threshold, sorted in decreasing order.

French wind data

Description

Daily mean wind speed (in km/h) at four stations in the south of France, namely Cap Cepet (S1), Lyon St-Exupery (S2), Marseille Marignane (S3) and Montelimar (S4).The data includes observations from January 1976 until April 2023; days containing missing values are omitted.

Format

A data frame with 17209 observations and 8 variables:

date

date of measurement

S1

wind speed (in km/h) at Cap Cepet

S2

wind speed (in km/h) at Lyon Saint-Exupery

S3

wind speed (in km/h) at Marseille Marignane

S4

wind speed (in km/h) at Montelimar

H2

humidity (in percentage) at Lyon Saint-Exupery

T2

mean temperature (in degree Celcius) at Lyon Saint-Exupery

Themetadata attribute includes latitude and longitude (in degrees, minutes, seconds), altitude (in m), station name and station id.

Source

European Climate Assessment and Dataset projecthttps://www.ecad.eu/

References

Klein Tank, A.M.G. and Coauthors, 2002. Daily dataset of 20th-century surface air temperature and precipitation series for theEuropean Climate Assessment. Int. J. of Climatol., 22, 1441-1453.

Examples

data(frwind, package = "mev")head(frwind)attr(frwind, which = "metadata")

Magnetic storms

Description

Absolute magnitude of 373 geomagnetic storms lasting more than 48h with absolute magnitude (dst) larger than 100 in 1957-2014.

Format

a vector of size 373

Note

For a detailed article presenting the derivation of the Dst index, seehttp://wdc.kugi.kyoto-u.ac.jp/dstdir/dst2/onDstindex.html

Source

Aki Vehtari

References

World Data Center for Geomagnetism, Kyoto, M. Nose, T. Iyemori, M. Sugiura, T. Kamei (2015),Geomagnetic Dst index, <doi:10.17593/14515-74000>.

Generalized extreme value distribution

Description

Likelihood, score function and information matrix, bias,approximate ancillary statistics and sample space derivativefor the generalized extreme value distribution

Arguments

Usage

gev.ll(par, dat)gev.ll.optim(par, dat)gev.score(par, dat)gev.infomat(par, dat, method = c('obs','exp'))gev.retlev(par, p)gev.bias(par, n)gev.Fscore(par, dat, method=c('obs','exp'))gev.Vfun(par, dat)gev.phi(par, dat, V)gev.dphi(par, dat, V)

Functions

• gev.ll: log likelihood

• gev.ll.optim: negative log likelihood parametrized in terms of location,log(scale) and shapein order to perform unconstrained optimization

• gev.score: score vector

• gev.infomat: observed or expected information matrix

• gev.retlev: return level, corresponding to the(1-p)th quantile

• gev.bias: Cox-Snell first order bias

• gev.Fscore: Firth's modified score equation

• gev.Vfun: vector implementing conditioning on approximate ancillary statistics for the TEM

• gev.phi: canonical parameter in the local exponential family approximation

• gev.dphi: derivative matrix of the canonical parameter in the local exponential family approximation

References

Firth, D. (1993). Bias reduction of maximum likelihood estimates,Biometrika,80(1), 27–38.

Coles, S. (2001).An Introduction to Statistical Modeling of Extreme Values, Springer, 209 p.

Cox, D. R. and E. J. Snell (1968). A general definition of residuals,Journal of the Royal Statistical Society: Series B (Methodological),30, 248–275.

Cordeiro, G. M. and R. Klein (1994). Bias correction in ARMA models,Statistics and Probability Letters,19(3), 169–176.

Firth's modified score equation for the generalized extreme value distribution

Description

Firth's modified score equation for the generalized extreme value distribution

Usage

gev.Fscore(par, dat, method = "obs")

Arguments

References

Firth, D. (1993). Bias reduction of maximum likelihood estimates,Biometrika,80(1), 27–38.

See Also

gev

N-year return levels, median and mean estimate

Description

N-year return levels, median and mean estimate

Usage

gev.Nyr(par, nobs, N, type = c("retlev", "median", "mean"), p = 1/N)

Arguments

Details

If there aren_y observations per year, theL-year return level is obtained by takingN equal ton_yL.

Value

a list with components

• est point estimate

• var variance estimate based on delta-method

• type statistic

Asymptotic bias of block maxima for fixed sample sizes

Description

Asymptotic bias of block maxima for fixed sample sizes

Usage

gev.abias(shape, rho)

Arguments

Value

a vector of length three containing the bias for location, scale and shape (in this order)

References

Dombry, C. and A. Ferreira (2017). Maximum likelihood estimators based on the block maxima method.https://arxiv.org/abs/1705.00465

Bias correction for GEV distribution

Description

Bias corrected estimates for the generalized extreme value distribution usingFirth's modified score function or implicit bias subtraction.

Usage

gev.bcor(par, dat, corr = c("subtract", "firth"), method = c("obs", "exp"))

Arguments

Details

Methodsubtractsolves

\tilde{\boldsymbol{\theta}} = \hat{\boldsymbol{\theta}} + b(\tilde{\boldsymbol{\theta}}

for\tilde{\boldsymbol{\theta}}, using the first order term in the bias expansion as given bygev.bias.

The alternative is to use Firth's modified score and find the root of

U(\tilde{\boldsymbol{\theta}})-i(\tilde{\boldsymbol{\theta}})b(\tilde{\boldsymbol{\theta}}),

whereU is the score vector,b is the first order bias andi is either the observed or Fisher information.

The routine uses the MLE (bias-corrected) as starting values and proceedsto find the solution using a root finding algorithm.Since the bias-correction is not valid for\xi < -1/3, any solution that is unboundedwill return a vector ofNA as the solution does not exist then.

Value

vector of bias-corrected parameters

Examples

set.seed(1)dat <- mev::rgev(n=40, loc = 1, scale=1, shape=-0.2)par <- mev::fit.gev(dat)$estimategev.bcor(par, dat, 'subtract')gev.bcor(par, dat, 'firth') #observed informationgev.bcor(par, dat, 'firth','exp')

Cox-Snell first order bias for the GEV distribution

Description

Bias vector for the GEV distribution based on ann sample.Due to numerical instability, values of the information matrix and the biasare linearly interpolated when the value of the shape parameter is close to zero.

Usage

gev.bias(par, n)

Arguments

See Also

gev

Information matrix for the generalized extreme value distribution

Description

Information matrix for the generalized extreme value distribution

Usage

gev.infomat(par, dat, method = c("obs", "exp"), nobs = length(dat))

Arguments

Value

a 3 by 3 matrix containing the expected or observed information

Log likelihood for the generalized extreme value distribution

Description

Function returning the density of ann sample from the GEV distribution.

Usage

gev.ll(par, dat)gev.ll.optim(par, dat)

Arguments

Details

gev.ll.optim returns the negative log likelihood parametrized in terms of location,log(scale) and shape in order to perform unconstrained optimization

See Also

gev

Generalized extreme value maximum likelihood estimates for various quantities of interest

Description

This function calls thefit.gev routine on the sample of block maxima and returns maximum likelihoodestimates for all quantities of interest, including location, scale and shape parameters, quantiles and mean andquantiles of maxima ofN blocks.

Usage

gev.mle(  xdat,  args = c("loc", "scale", "shape", "quant", "Nmean", "Nquant"),  N,  p,  q)

Arguments

Value

named vector with maximum likelihood estimated parameter values for argumentsargs

Examples

dat <- mev::rgev(n = 100, shape = 0.2)gev.mle(xdat = dat, N = 100, p = 0.01, q = 0.5)

Profile log-likelihood for the generalized extreme value distribution

Description

This function calculates the profile likelihood along with two small-sample correctionsbased on Severini's (1999) empirical covariance and the Fraser and Reid tangent exponentialmodel approximation.

Usage

gev.pll(  psi,  param = c("loc", "scale", "shape", "quant", "Nmean", "Nquant"),  mod = "profile",  dat,  N = NULL,  p = NULL,  q = NULL,  correction = TRUE,  plot = TRUE,  ...)

Arguments

Details

The two additionalmod available aretem, the tangent exponential model (TEM) approximation andmodif for the penalized profile likelihood based onp^* approximation proposed by Severini.For the latter, the penalization is based on the TEM or an empirical covariance adjustment term.

Value

a list with components

• mle: maximum likelihood estimate

• psi.max: maximum profile likelihood estimate

• param: string indicating the parameter to profile over

• std.error: standard error ofpsi.max

• psi: vector of parameter\psi given inpsi

• pll: values of the profile log likelihood atpsi

• maxpll: value of maximum profile log likelihood

In addition, ifmod includestem

• normal: maximum likelihood estimate and standard error of the interest parameter\psi

• r: values of likelihood root corresponding to\psi

• q: vector of likelihood modifications

• rstar: modified likelihood root vector

• rstar.old: uncorrected modified likelihood root vector

• tem.psimax: maximum of the tangent exponential model likelihood

In addition, ifmod includesmodif

• tem.mle: maximum of tangent exponential modified profile log likelihood

• tem.profll: values of the modified profile log likelihood atpsi

• tem.maxpll: value of maximum modified profile log likelihood

• empcov.mle: maximum of Severini's empirical covariance modified profile log likelihood

• empcov.profll: values of the modified profile log likelihood atpsi

• empcov.maxpll: value of maximum modified profile log likelihood

References

Fraser, D. A. S., Reid, N. and Wu, J. (1999), A simple general formula for tail probabilities for frequentist and Bayesian inference.Biometrika,86(2), 249–264.

Severini, T. (2000) Likelihood Methods in Statistics. Oxford University Press. ISBN 9780198506508.

Brazzale, A. R., Davison, A. C. and Reid, N. (2007) Applied asymptotics: case studies in small-sample statistics. Cambridge University Press, Cambridge. ISBN 978-0-521-84703-2

Examples

## Not run: set.seed(123)dat <- rgev(n = 100, loc = 0, scale = 2, shape = 0.3)gev.pll(psi = seq(0,0.5, length = 50), param = 'shape', dat = dat)gev.pll(psi = seq(-1.5, 1.5, length = 50), param = 'loc', dat = dat)gev.pll(psi = seq(10, 40, length = 50), param = 'quant', dat = dat, p = 0.01)gev.pll(psi = seq(12, 100, length = 50), param = 'Nmean', N = 100, dat = dat)gev.pll(psi = seq(12, 90, length = 50), param = 'Nquant', N = 100, dat = dat, q = 0.5)## End(Not run)

Return level for the generalized extreme value distribution

Description

This function returns the1-pth quantile of the GEV.

Usage

gev.retlev(par, p)

Arguments

See Also

gev

Score vector for the generalized extreme value distribution

Description

Score vector for the generalized extreme value distribution

Usage

gev.score(par, dat)

Arguments

Value

a vector of length containing the derivative of the

Tangent exponential model approximation for the GEV distribution

Description

The functiongev.tem provides a tangent exponential model (TEM) approximationfor higher order likelihood inference for a scalar parameter for the generalized extreme value distribution.Options include location scale and shape parameters as well as value-at-risk (or return levels).The function attempts to find good values forpsi that willcover the range of options, but the fail may fit and return an error.

Usage

gev.tem(  param = c("loc", "scale", "shape", "quant", "Nmean", "Nquant"),  dat,  psi = NULL,  p = NULL,  q = 0.5,  N = NULL,  n.psi = 50,  plot = TRUE,  correction = TRUE)

Arguments

Value

an invisible object of classfr (seetem in packagehoa) with elements

• normal: maximum likelihood estimate and standard error of the interest parameter\psi

• par.hat: maximum likelihood estimates

• par.hat.se: standard errors of maximum likelihood estimates

• th.rest: estimated maximum profile likelihood at (\psi,\hat{\lambda})

• r: values of likelihood root corresponding to\psi

• psi: vector of interest parameter

• q: vector of likelihood modifications

• rstar: modified likelihood root vector

• rstar.old: uncorrected modified likelihood root vector

• param: parameter

Author(s)

Leo Belzile

Examples

## Not run: set.seed(1234)dat <- rgev(n = 40, loc = 0, scale = 2, shape = -0.1)gev.tem('shape', dat = dat, plot = TRUE)gev.tem('quant', dat = dat, p = 0.01, plot = TRUE)gev.tem('scale', psi = seq(1, 4, by = 0.1), dat = dat, plot = TRUE)dat <- rgev(n = 40, loc = 0, scale = 2, shape = 0.2)gev.tem('loc', dat = dat, plot = TRUE)gev.tem('Nmean', dat = dat, p = 0.01, N=100, plot = TRUE)gev.tem('Nquant', dat = dat, q = 0.5, N=100, plot = TRUE)## End(Not run)

Tangent exponential model statistics for the generalized extreme value distribution

Description

Matrix of approximate ancillary statistics, sample space derivative of thelog likelihood and mixed derivative for the generalized extreme value distribution.

Usage

gev.Vfun(par, dat)gev.phi(par, dat, V)gev.dphi(par, dat, V)

Arguments

See Also

gev

Generalized extreme value distribution (quantile/mean of N-block maxima parametrization)

Description

Likelihood, score function and information matrix,approximate ancillary statistics and sample space derivativefor the generalized extreme value distribution parametrized in terms of thequantiles/mean of N-block maxima parametrizationz, scale and shape.

Arguments

Usage

gevN.ll(par, dat, N, q, qty = c('mean', 'quantile'))gevN.ll.optim(par, dat, N, q = 0.5, qty = c('mean', 'quantile'))gevN.score(par, dat, N, q = 0.5, qty = c('mean', 'quantile'))gevN.infomat(par, dat, qty = c('mean', 'quantile'), method = c('obs', 'exp'), N, q = 0.5, nobs = length(dat))gevN.Vfun(par, dat, N, q = 0.5, qty = c('mean', 'quantile'))gevN.phi(par, dat, N, q = 0.5, qty = c('mean', 'quantile'), V)gevN.dphi(par, dat, N, q = 0.5, qty = c('mean', 'quantile'), V)

Functions

• gevN.ll: log likelihood

• gevN.score: score vector

• gevN.infomat: expected and observed information matrix

• gevN.Vfun: vector implementing conditioning on approximate ancillary statistics for the TEM

• gevN.phi: canonical parameter in the local exponential family approximation

• gevN.dphi: derivative matrix of the canonical parameter in the local exponential family approximation

Author(s)

Leo Belzile

Information matrix of the generalized extreme value distribution (quantile/mean of N-block maxima parametrization)

Description

Information matrix of the generalized extreme value distribution (quantile/mean of N-block maxima parametrization)

Usage

gevN.infomat(  par,  dat,  method = c("obs", "exp"),  qty = c("mean", "quantile"),  N,  q = 0.5,  nobs = length(dat))

Arguments

See Also

gevN

Negative log likelihood of the generalized extreme value distribution (quantile/mean of N-block maxima parametrization)

Description

Negative log likelihood of the generalized extreme value distribution (quantile/mean of N-block maxima parametrization)

Usage

gevN.ll(par, dat, N, q = 0.5, qty = c("mean", "quantile"))

Arguments

See Also

gevN

This function returns the mean of N observations from the GEV.

Description

This function returns the mean of N observations from the GEV.

Usage

gevN.mean(par, N)

Arguments

See Also

gevN

Score of the generalized extreme value distribution (quantile/mean of N-block maxima parametrization)

Description

Score of the generalized extreme value distribution (quantile/mean of N-block maxima parametrization)

Usage

gevN.score(par, dat, N, q = 0.5, qty = c("mean", "quantile"))

Arguments

See Also

gevN

Tangent exponential model statistics for the generalized Pareto distribution (mean of maximum of N exceedances parametrization)

Description

Vector implementing conditioning on approximate ancillary statistics for the TEM

Usage

gevN.Vfun(par, dat, N, q = 0.5, qty = c("mean", "quantile"))gevN.phi(par, dat, N, q = 0.5, qty = c("mean", "quantile"), V)gevN.dphi(par, dat, N, q = 0.5, qty = c("mean", "quantile"), V)

Arguments

See Also

gevN

Generalized extreme value distribution

Description

Density function, distribution function, quantile function andrandom number generation for the generalized extreme valuedistribution.

Usage

qgev(p, loc = 0, scale = 1, shape = 0, lower.tail = TRUE, log.p = FALSE)rgev(n, loc = 0, scale = 1, shape = 0)dgev(x, loc = 0, scale = 1, shape = 0, log = FALSE)pgev(q, loc = 0, scale = 1, shape = 0, lower.tail = TRUE, log.p = FALSE)

Arguments

Details

The distribution function of a GEV distribution with parametersloc =\mu,scale =\sigma andshape =\xi is

F(x) = \exp\{-[1 + \xi (x - \mu) / \sigma] ^ {-1/\xi} \}

for1 + \xi (x - \mu) / \sigma > 0. If\xi = 0 thedistribution function is defined as the limit as\xi tends to zero.

The quantile function, when evaluated at zero or one,returns the lower and upper endpoint, whether the latter is finite or not.

Author(s)

Leo Belzile, with code adapted from Paul Northrop

References

Jenkinson, A. F. (1955) The frequency distribution of theannual maximum (or minimum) of meteorological elements.Quart. J. R. Met. Soc.,81, 158-171.Chapter 3:doi:10.1002/qj.49708134804

Coles, S. G. (2001)An Introduction to StatisticalModeling of Extreme Values, Springer-Verlag, London.doi:10.1007/978-1-4471-3675-0_3

Generalized extreme value distribution (return level parametrization)

Description

Likelihood, score function and information matrix,approximate ancillary statistics and sample space derivativefor the generalized extreme value distribution parametrized in terms of the return levelz, scale and shape.

Arguments

Usage

gevr.ll(par, dat, p)gevr.ll.optim(par, dat, p)gevr.score(par, dat, p)gevr.infomat(par, dat, p, method = c('obs', 'exp'), nobs = length(dat))gevr.Vfun(par, dat, p)gevr.phi(par, dat, p, V)gevr.dphi(par, dat, p, V)

Functions

• gevr.ll: log likelihood

• gevr.ll.optim: negative log likelihood parametrized in terms of return levels,log(scale) and shape in order to perform unconstrained optimization

• gevr.score: score vector

• gevr.infomat: observed information matrix

• gevr.Vfun: vector implementing conditioning on approximate ancillary statistics for the TEM

• gevr.phi: canonical parameter in the local exponential family approximation

• gevr.dphi: derivative matrix of the canonical parameter in the local exponential family approximation

Author(s)

Leo Belzile

Observed information matrix for GEV distribution (return levels)

Description

The information matrix is parametrized in terms of return level ((1-p)th quantile), scale and shape.

Usage

gevr.infomat(par, dat, method = c("obs", "exp"), p, nobs = length(dat))

Arguments

See Also

gevr

Negative log likelihood of the generalized extreme value distribution (return levels)

Description

Negative log likelihood of the generalized extreme value distribution (return levels)

Negative log likelihood parametrized in terms of location, log return level and shape in order to perform unconstrained optimization

Usage

gevr.ll(par, dat, p)gevr.ll.optim(par, dat, p)

Arguments

See Also

gevr

Score of the log likelihood for the GEV distribution (return levels)

Description

Score of the log likelihood for the GEV distribution (return levels)

Usage

gevr.score(par, dat, p)

Arguments

See Also

gevr

Tangent exponential model statistics for the GEV distribution (return level)

Description

Vector implementing conditioning on approximate ancillary statistics for the TEM

Usage

gevr.Vfun(par, dat, p)gevr.phi(par, dat, p, V)gevr.dphi(par, dat, p, V)

Arguments

See Also

gevr

Maximum likelihood estimate of generalized Pareto applied to threshold exceedances

Description

The functionfit.gpd is a wrapper aroundgp.fit

Usage

gp.fit(  xdat,  threshold,  method = c("Grimshaw", "auglag", "nlm", "optim", "ismev", "zs", "zhang"),  show = FALSE,  MCMC = NULL,  fpar = NULL,  warnSE = TRUE)

Arguments

Value

a list; seefit.gpd for details

Generalized Pareto distribution

Description

Likelihood, score function and information matrix, bias,approximate ancillary statistics and sample space derivativefor the generalized Pareto distribution

Arguments

Usage

gpd.ll(par, dat, tol=1e-5)gpd.ll.optim(par, dat, tol=1e-5)gpd.score(par, dat)gpd.infomat(par, dat, method = c('obs','exp'))gpd.bias(par, n)gpd.Fscore(par, dat, method = c('obs','exp'))gpd.Vfun(par, dat)gpd.phi(par, dat, V)gpd.dphi(par, dat, V)

Functions

• gpd.ll: log likelihood

• gpd.ll.optim: negative log likelihood parametrized in terms oflog(scale) and shapein order to perform unconstrained optimization

• gpd.score: score vector

• gpd.infomat: observed or expected information matrix

• gpd.bias: Cox-Snell first order bias

• gpd.Fscore: Firth's modified score equation

• gpd.Vfun: vector implementing conditioning on approximate ancillary statistics for the TEM

• gpd.phi: canonical parameter in the local exponential family approximation

• gpd.dphi: derivative matrix of the canonical parameter in the localexponential family approximation

Author(s)

Leo Belzile

References

Firth, D. (1993). Bias reduction of maximum likelihood estimates,Biometrika,80(1), 27–38.

Coles, S. (2001).An Introduction to Statistical Modeling of Extreme Values, Springer, 209 p.

Cox, D. R. and E. J. Snell (1968). A general definition of residuals,Journal of the Royal Statistical Society: Series B (Methodological),30, 248–275.

Cordeiro, G. M. and R. Klein (1994). Bias correction in ARMA models,Statistics and Probability Letters,19(3), 169–176.

Giles, D. E., Feng, H. and R. T. Godwin (2016). Bias-corrected maximum likelihood estimation of the parameters of the generalized Pareto distribution,Communications in Statistics - Theory and Methods,45(8), 2465–2483.

Firth's modified score equation for the generalized Pareto distribution

Description

Firth's modified score equation for the generalized Pareto distribution

Usage

gpd.Fscore(par, dat, method = c("obs", "exp"))

Arguments

References

Firth, D. (1993). Bias reduction of maximum likelihood estimates,Biometrika,80(1), 27–38.

See Also

gpd,gpd.bcor

Asymptotic bias of threshold exceedances for k order statistics

Description

The formula given in de Haan and Ferreira, 2007 (Springer). Note that the latter differs from that found in Drees, Ferreira and de Haan.

Usage

gpd.abias(shape, rho)

Arguments

Value

a vector of length containing the bias for scale and shape (in this order)

References

Dombry, C. and A. Ferreira (2017). Maximum likelihood estimators based on the block maxima method.https://arxiv.org/abs/1705.00465

Bias correction for GP distribution

Description

Bias corrected estimates for the generalized Pareto distribution usingFirth's modified score function or implicit bias subtraction.

Usage

gpd.bcor(par, dat, corr = c("subtract", "firth"), method = c("obs", "exp"))

Arguments

Details

Methodsubtract solves

\tilde{\boldsymbol{\theta}} = \hat{\boldsymbol{\theta}} + b(\tilde{\boldsymbol{\theta}}

for\tilde{\boldsymbol{\theta}}, using the first order term in the bias expansion as given bygpd.bias.

The alternative is to use Firth's modified score and find the root of

U(\tilde{\boldsymbol{\theta}})-i(\tilde{\boldsymbol{\theta}})b(\tilde{\boldsymbol{\theta}}),

whereU is the score vector,b is the first order bias andi is either the observed or Fisher information.

The routine uses the MLE as starting value and proceedsto find the solution using a root finding algorithm.Since the bias-correction is not valid for\xi < -1/3, any solution that is unboundedwill return a vector ofNA as the bias correction does not exist then.

Value

vector of bias-corrected parameters

Examples

set.seed(1)dat <- rgp(n=40, scale=1, shape=-0.2)par <- gp.fit(dat, threshold=0, show=FALSE)$estimategpd.bcor(par,dat, 'subtract')gpd.bcor(par,dat, 'firth') #observed informationgpd.bcor(par,dat, 'firth','exp')

Cox-Snell first order bias expression for the generalized Pareto distribution

Description

Bias vector for the GP distribution based on ann sample.

Usage

gpd.bias(par, n)

Arguments

References

Coles, S. (2001).An Introduction to Statistical Modeling of Extreme Values, Springer, 209 p.

Cox, D. R. and E. J. Snell (1968). A general definition of residuals,Journal of the Royal Statistical Society: Series B (Methodological),30, 248–275.

Cordeiro, G. M. and R. Klein (1994). Bias correction in ARMA models,Statistics and Probability Letters,19(3), 169–176.

Giles, D. E., Feng, H. and R. T. Godwin (2016). Bias-corrected maximum likelihood estimation of the parameters of the generalized Pareto distribution,Communications in Statistics - Theory and Methods,45(8), 2465–2483.

See Also

gpd,gpd.bcor

Bootstrap approximation for generalized Pareto parameters

Description

Given an object of classmev_gpd,returns a matrix of parameter values to mimicthe estimation uncertainty.

Usage

gpd.boot(object, B = 1000L, method = c("post", "norm"))

Arguments

Details

Two options are available: a normal approximation tothe scale and shape based on the maximum likelihoodestimates and the observed information matrix.This method uses forward sampling to simulatefrom a bivariate normal distribution that satisfiesthe support and positivity constraints

The second approximation uses the ratio-of-uniformsmethod to obtain samples from the posteriordistribution with uninformative priors, thusmimicking the joint distribution of maximum likelihood.The benefit of the latter is that it is more reliablein small samples and when the shape is negative.

Value

a matrix of size B by 2 whose columns contain scale and shape parameters

Information matrix for the generalized Pareto distribution

Description

The function returns the expected or observed information matrix.

Usage

gpd.infomat(par, dat, method = c("obs", "exp"), nobs = length(dat))

Arguments

See Also

gpd

Log likelihood for the generalized Pareto distribution

Description

Function returning the density of ann sample from the GP distribution.gpd.ll.optim returns the negative log likelihood parametrized in terms oflog(scale) and shapein order to perform unconstrained optimization. The function is coded in such a way that the log likelihood is infinite when\xi < -1.

Usage

gpd.ll(par, dat, tol = 1e-05)gpd.ll.optim(par, dat, tol = 1e-05)

Arguments

See Also

gpd

Estimation of generalized Pareto parameters via L-moments

Description

Given a sample of exceedances, compute the first four L-moments and use either the first two to obtain the scale and shape (default), or else use L-skewness and L-scale to compute the scale and shape of the generalized Pareto distribution

Usage

gpd.lmom(xdat, thresh, sorted = FALSE, Lskew = FALSE)

Arguments

Value

a vector of length two with the scale and shape estimates

Generalized Pareto maximum likelihood estimates for various quantities of interest

Description

This function calls thefit.gpd routine on the sample of excesses and returns maximum likelihoodestimates for all quantities of interest, including scale and shape parameters, quantiles and value-at-risk,expected shortfall and mean and quantiles of maxima ofN threshold exceedances

Usage

gpd.mle(  xdat,  args = c("scale", "shape", "quant", "VaR", "ES", "Nmean", "Nquant"),  m,  N,  p,  q)

Arguments

Value

named vector with maximum likelihood values for argumentsargs

Examples

xdat <- mev::rgp(n = 30, shape = 0.2)gpd.mle(xdat = xdat, N = 100, p = 0.01, q = 0.5, m = 100)

Profile log-likelihood for the generalized Pareto distribution

Description

This function calculates the (modified) profile likelihood based on thep^* formula.There are two small-sample corrections that use a proxy for\ell_{\lambda; \hat{\lambda}},which are based on Severini's (1999) empirical covarianceand the Fraser and Reid tangent exponential model approximation.

Usage

gpd.pll(  psi,  param = c("scale", "shape", "quant", "retlev", "VaR", "ES", "Nmean", "Nquant"),  mod = "profile",  mle = NULL,  dat,  m = NULL,  N = NULL,  p = NULL,  q = NULL,  correction = TRUE,  thresh = NULL,  plot = TRUE,  ...)

Arguments

Details

The threemod available areprofile (the default),tem, the tangent exponential model (TEM) approximation andmodif for the penalized profile likelihood based onp^* approximation proposed by Severini.For the latter, the penalization is based on the TEM or an empirical covariance adjustment term.

Value

a list with components

• mle: maximum likelihood estimate

• psi.max: maximum profile likelihood estimate

• param: string indicating the parameter to profile over

• std.error: standard error ofpsi.max

• psi: vector of parameter\psi given inpsi

• pll: values of the profile log likelihood atpsi

• maxpll: value of maximum profile log likelihood

• family: a string indicating "gpd"

• thresh: value of the threshold, by default zero

In addition, ifmod includestem

• normal: maximum likelihood estimate and standard error of the interest parameter\psi

• r: values of likelihood root corresponding to\psi

• q: vector of likelihood modifications

• rstar: modified likelihood root vector

• rstar.old: uncorrected modified likelihood root vector

• tem.psimax: maximum of the tangent exponential model likelihood

In addition, ifmod includesmodif

• tem.mle: maximum of tangent exponential modified profile log likelihood

• tem.profll: values of the modified profile log likelihood atpsi

• tem.maxpll: value of maximum modified profile log likelihood

• empcov.mle: maximum of Severini's empirical covariance modified profile log likelihood

• empcov.profll: values of the modified profile log likelihood atpsi

• empcov.maxpll: value of maximum modified profile log likelihood

Examples

## Not run: dat <- rgp(n = 100, scale = 2, shape = 0.3)gpd.pll(psi = seq(-0.5, 1, by=0.01), param = 'shape', dat = dat)gpd.pll(psi = seq(0.1, 5, by=0.1), param = 'scale', dat = dat)gpd.pll(psi = seq(20, 35, by=0.1), param = 'quant', dat = dat, p = 0.01)gpd.pll(psi = seq(20, 80, by=0.1), param = 'ES', dat = dat, m = 100)gpd.pll(psi = seq(15, 100, by=1), param = 'Nmean', N = 100, dat = dat)gpd.pll(psi = seq(15, 90, by=1), param = 'Nquant', N = 100, dat = dat, q = 0.5)## End(Not run)

Score vector for the generalized Pareto distribution

Description

Score vector for the generalized Pareto distribution

Usage

gpd.score(par, dat)

Arguments

See Also

gpd

Tangent exponential model approximation for the GP distribution

Description

The functiongpd.tem provides a tangent exponential model (TEM) approximationfor higher order likelihood inference for a scalar parameter for the generalized Pareto distribution. Options includescale and shape parameters as well as value-at-risk (also referred to as quantiles, or return levels)and expected shortfall. The function attempts to find good values forpsi that willcover the range of options, but the fit may fail and return an error. In such cases, the user can try to find goodgrid of starting values and provide them to the routine.

Usage

gpd.tem(  dat,  param = c("scale", "shape", "quant", "VaR", "retlev", "ES", "Nmean", "Nquant"),  psi = NULL,  m = NULL,  thresh = 0,  n.psi = 50,  N = NULL,  p = NULL,  q = NULL,  plot = FALSE,  correction = TRUE,  ...)

Arguments

Details

As of version 1.11, this function is a wrapper aroundgpd.pll.

The interpretation form is as follows: if there are on averagem_y observations per year above the threshold, thenm = Tm_y corresponds toT-year return level.

Value

an invisible object of classfr (seetem in packagehoa) with elements

• normal: maximum likelihood estimate and standard error of the interest parameter\psi

• par.hat: maximum likelihood estimates

• par.hat.se: standard errors of maximum likelihood estimates

• th.rest: estimated maximum profile likelihood at (\psi,\hat{\lambda})

• r: values of likelihood root corresponding to\psi

• psi: vector of interest parameter

• q: vector of likelihood modifications

• rstar: modified likelihood root vector

• rstar.old: uncorrected modified likelihood root vector

• param: parameter

Author(s)

Leo Belzile

Examples

set.seed(123)dat <- rgp(n = 40, scale = 1, shape = -0.1)#with plotsm1 <- gpd.tem(param = 'shape', n.psi = 50, dat = dat, plot = TRUE)## Not run: m2 <- gpd.tem(param = 'scale', n.psi = 50, dat = dat)m3 <- gpd.tem(param = 'VaR', n.psi = 50, dat = dat, m = 100)#Providing psipsi <- c(seq(2, 5, length = 15), seq(5, 35, length = 45))m4 <- gpd.tem(param = 'ES', dat = dat, m = 100, psi = psi, correction = FALSE)mev:::plot.fr(m4, which = c(2, 4))plot(fr4 <- spline.corr(m4))confint(m1)confint(m4, parm = 2, warn = FALSE)m5 <- gpd.tem(param = 'Nmean', dat = dat, N = 100, psi = psi, correction = FALSE)m6 <- gpd.tem(param = 'Nquant', dat = dat, N = 100, q = 0.7, correction = FALSE)## End(Not run)

Tangent exponential model statistics for the generalized Pareto distribution

Description

Matrix of approximate ancillary statistics, sample space derivative of thelog likelihood and mixed derivative for the generalized Pareto distribution.

Usage

gpd.Vfun(par, dat)gpd.phi(par, dat, V)gpd.dphi(par, dat, V)

Arguments

See Also

gpd

Generalized Pareto distribution (mean of maximum of N exceedances parametrization)

Description

Likelihood, score function and information matrix,approximate ancillary statistics and sample space derivativefor the generalized Pareto distribution parametrized in terms of average maximum ofN exceedances.

The parameterN corresponds to the number of threshold exceedances of interest over which the maxima is taken.z is the corresponding expected value of this block maxima.Note that the actual parametrization is in terms of excess expected mean, meaning expected mean minus threshold.

Arguments

Details

The observed information matrix was calculated from the Hessian using symbolic calculus in Sage.

Usage

gpdN.ll(par, dat, N, tol=1e-5)gpdN.score(par, dat, N)gpdN.infomat(par, dat, N, method = c('obs', 'exp'), nobs = length(dat))gpdN.Vfun(par, dat, N)gpdN.phi(par, dat, N, V)gpdN.dphi(par, dat, N, V)

Functions

• gpdN.ll: log likelihood

• gpdN.score: score vector

• gpdN.infomat: observed information matrix for GP parametrized in terms of mean of the maximum ofN exceedances and shape

• gpdN.Vfun: vector implementing conditioning on approximate ancillary statistics for the TEM

• gpdN.phi: canonical parameter in the local exponential family approximation

• gpdN.dphi: derivative matrix of the canonical parameter in the local exponential family approximation

Author(s)

Leo Belzile

Information matrix of the generalized Pareto distribution (mean of maximum of N exceedances parametrization)

Description

Information matrix of the generalized Pareto distribution (mean of maximum of N exceedances parametrization)

Usage

gpdN.infomat(par, dat, N, method = c("obs", "exp"), nobs = length(dat))

Arguments

See Also

gpdN

Negative log likelihood of the generalized Pareto distribution (mean of maximum of N exceedances parametrization)

Description

Negative log likelihood of the generalized Pareto distribution (mean of maximum of N exceedances parametrization)

Usage

gpdN.ll(par, dat, N)

Arguments

See Also

gpdN

This function returns the mean of N maxima from the GP.

Description

This function returns the mean of N maxima from the GP.

Usage

gpdN.mean(par, N)

Arguments

See Also

gpdN

This function returns the qth percentile of N maxima from the GP.

Description

This function returns the qth percentile of N maxima from the GP.

Usage

gpdN.quant(par, q, N)

Arguments

See Also

gpdN

Score vector of the generalized Pareto distribution (mean of maximum of N exceedances parametrization)

Description

Score vector of the generalized Pareto distribution (mean of maximum of N exceedances parametrization)

Usage

gpdN.score(par, dat, N)

Arguments

See Also

gpdN

Tangent exponential model statistics for the generalized Pareto distribution (mean of maximum of N exceedances parametrization)

Description

Vector implementing conditioning on approximate ancillary statistics for the TEM

Usage

gpdN.Vfun(par, dat, N)gpdN.phi(par, dat, N, V)gpdN.dphi(par, dat, N, V)

Arguments

See Also

gpdN

Generalized Pareto distribution (expected shortfall parametrization)

Description

Likelihood, score function and information matrix,approximate ancillary statistics and sample space derivativefor the generalized Pareto distribution parametrized in terms of expected shortfall.

The parameterm corresponds to\zeta_u/(1-\alpha), where\zeta_u is the rate of exceedance over the thresholdu and\alpha is the percentile of the expected shortfall.Note that the actual parametrization is in terms of excess expected shortfall, meaning expected shortfall minus threshold.

Arguments

Details

The observed information matrix was calculated from the Hessian using symbolic calculus in Sage.

Usage

gpde.ll(par, dat, m, tol=1e-5)gpde.ll.optim(par, dat, m, tol=1e-5)gpde.score(par, dat, m)gpde.infomat(par, dat, m, method = c('obs', 'exp'), nobs = length(dat))gpde.Vfun(par, dat, m)gpde.phi(par, dat, V, m)gpde.dphi(par, dat, V, m)

Functions

• gpde.ll: log likelihood

• gpde.ll.optim: negative log likelihood parametrized in terms of log expectedshortfall and shape in order to perform unconstrained optimization

• gpde.score: score vector

• gpde.infomat: observed information matrix for GPD parametrized in terms of rate of expected shortfall and shape

• gpde.Vfun: vector implementing conditioning on approximate ancillary statistics for the TEM

• gpde.phi: canonical parameter in the local exponential family approximation

• gpde.dphi: derivative matrix of the canonical parameter in the local exponential family approximation

Author(s)

Leo Belzile

Observed information matrix for the GP distribution (expected shortfall)

Description

The information matrix is parametrized in terms of excess expected shortfall and shape

Usage

gpde.infomat(par, dat, m, method = c("obs", "exp"), nobs = length(dat))

Arguments

See Also

gpde

Negative log likelihood of the generalized Pareto distribution (expected shortfall)

Description

Negative log likelihood of the generalized Pareto distribution (expected shortfall)

Negative log likelihood of the generalized Pareto distribution (expected shortfall) - optimizationThe negative log likelihood is parametrized in terms of log expected shortfall and shape in order to perform unconstrained optimization

Usage

gpde.ll(par, dat, m)gpde.ll.optim(par, dat, m)

Arguments

See Also

gpde

Score vector for the GP distribution (expected shortfall)

Description

Score vector for the GP distribution (expected shortfall)

Usage

gpde.score(par, dat, m)

Arguments

See Also

gpde

Tangent exponential model statistics for the generalized Pareto distribution (expected shortfall)

Description

Vector implementing conditioning on approximate ancillary statistics for the TEM

Usage

gpde.Vfun(par, dat, m)gpde.phi(par, dat, V, m)gpde.dphi(par, dat, V, m)

Arguments

See Also

gpde

Generalized Pareto distribution

Description

Density function, distribution function, quantile function andrandom number generation for the generalized Paretodistribution.

Usage

pgp(q, loc = 0, scale = 1, shape = 0, lower.tail = TRUE, log.p = FALSE)dgp(x, loc = 0, scale = 1, shape = 0, log = FALSE)qgp(p, loc = 0, scale = 1, shape = 0, lower.tail = TRUE)rgp(n, loc = 0, scale = 1, shape = 0)

Arguments

References

Coles, S. G. (2001)An Introduction to StatisticalModeling of Extreme Values, Springer-Verlag, London.doi:10.1007/978-1-4471-3675-0_3

Generalized Pareto distribution (return level parametrization)

Description

Likelihood, score function and information matrix,approximate ancillary statistics and sample space derivativefor the generalized Pareto distribution parametrized in terms of return levels.

Arguments

Details

The observed information matrix was calculated from the Hessian using symbolic calculus in Sage.

The interpretation form is as follows: if there are on averagem_y observations per year above the threshold, thenm=Tm_y corresponds toT-year return level.

Usage

gpdr.ll(par, dat, m, tol=1e-5)gpdr.ll.optim(par, dat, m, tol=1e-5)gpdr.score(par, dat, m)gpdr.infomat(par, dat, m, method = c('obs', 'exp'), nobs = length(dat))gpdr.Vfun(par, dat, m)gpdr.phi(par, V, dat, m)gpdr.dphi(par, V, dat, m)

Functions

• gpdr.ll: log likelihood

• gpdr.ll.optim: negative log likelihood parametrized in terms oflog(scale) and shapein order to perform unconstrained optimization

• gpdr.score: score vector

• gpdr.infomat: observed information matrix for GPD parametrized in terms of rate ofm-year return level and shape

• gpdr.Vfun: vector implementing conditioning on approximate ancillary statistics for the TEM

• gpdr.phi: canonical parameter in the local exponential family approximation

• gpdr.dphi: derivative matrix of the canonical parameter in the local exponential family approximation

Author(s)

Leo Belzile

Observed information matrix for GP distribution (return levels)

Description

The information matrix is parametrized in terms of rate ofm-year return level and shape

Usage

gpdr.infomat(par, dat, m, method = c("obs", "exp"), nobs = length(dat))

Arguments

See Also

gpdr

Negative log likelihood of the generalized Pareto distribution (return levels)

Description

Negative log likelihood of the generalized Pareto distribution (return levels)

Negative log likelihood parametrized in terms of log return level and shape in order to perform unconstrained optimization

Usage

gpdr.ll(par, dat, m)gpdr.ll.optim(par, dat, m)

Arguments

See Also

gpdr

Score of the profile log likelihood for the GP distribution (return levels parametrization)

Description

Score of the profile log likelihood for the GP distribution (return levels parametrization)

Usage

gpdr.score(par, dat, m)

Arguments

See Also

gpdr

Tangent exponential model statistics for the generalized Pareto distribution (return level)

Description

Vector implementing conditioning on approximate ancillary statistics for the TEM

Usage

gpdr.Vfun(par, dat, m)gpdr.phi(par, dat, V, m)gpdr.dphi(par, dat, V, m)

Arguments

See Also

gpdr

gpdr

Transformation from the generalized Pareto to unit Pareto

Description

Transformation from the generalized Pareto to unit Pareto

Usage

gpdtopar(dat, loc = 0, scale, shape, lambdau = 1)

Arguments

Value

a vector or matrix of the same dimension asdat with unit Pareto observations

Interpret bivariate threshold exceedance models

Description

This is an adaptation of theevir packageinterpret.gpdbiv function.interpret.fbvpot deals with the output of a call tofbvpot from theevd and to handle families other than the logistic distribution.The likelihood derivation comes from expression 2.10 in Smith et al. (1997).

Usage

ibvpot(fitted, q, silent = FALSE)

Arguments

Details

The listfitted must contain

• model a string; seebvevd from packageevd for options

• param a named vector containing the parameters of themodel, as well as parametersscale1,shape1,scale2 andshape2, corresponding to marginal GPD parameters.

• threshold a vector of length 2 containing the two thresholds.

• pat the proportion of observations above the correspondingthreshold

Value

an invisible numeric vector containing marginal, joint and conditional exceedance probabilities.

Author(s)

Leo Belzile, adapting original S code by Alexander McNeil

References

Smith, Tawn and Coles (1997), Markov chain models for threshold exceedances.Biometrika,84(2), 249–268.

See Also

interpret.gpdbiv in packageevir

Examples

if (requireNamespace("evd", quietly = TRUE)) {y <- rgp(1000,1,1,1)x <- y*rmevspec(n=1000,d=2,sigma=cbind(c(0,0.5),c(0.5,0)), model='hr')mod <- evd::fbvpot(x, threshold = c(1,1), model = 'hr', likelihood ='censored')ibvpot(mod, c(20,20))}

Information matrix test statistic and MLE for the extremal index

Description

The information matrix test (IMT), proposed by Suveges and Davison (2010), is basedon the difference between the expected quadratic score and the second derivative ofthe log-likelihood. The asymptotic distribution for each thresholdu and gapKis asymptotically\chi^2 with one degree of freedom. The approximation is good forN>80 and conservative for smaller sample sizes. The test assumes independence between gaps.

Usage

infomat.test(xdat, thresh, q, K, plot = TRUE, ...)

Arguments

Details