Type:	Package
Title:	Modelling of Extreme Values
Version:	1.17
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, Rcpp (≥ 0.12.16), Rsolnp, stats,
Suggests:	boot, cobs, evd, knitr, MASS, mvPot (≥ 0.1.4), mvtnorm, gmm, revdbayes, rmarkdown, ismev, tinytest, TruncatedNormal (≥ 1.1)
LinkingTo:	Rcpp, RcppArmadillo
LazyData:	true
RoxygenNote:	7.3.2
VignetteBuilder:	knitr
Encoding:	UTF-8
ByteCompile:	true
NeedsCompilation:	yes
Packaged:	2024-07-09 03:31:35 UTC; lbelzile
Author:	Leo Belzile [aut, cre], Jennifer L. Wadsworth [aut], Paul J. Northrop [aut], Scott D. Grimshaw [aut], Jin Zhang [ctb], Michael A. Stephens [ctb], Art B. Owen [ctb], Raphael Huser [aut]
Maintainer:	Leo Belzile <belzilel@gmail.com>
Repository:	CRAN
Date/Publication:	2024-07-09 04:20:02 UTC

Contrast matrix

Description

Produces a contrast matrix with (1,-1) elements running down the two diagonals

Usage

.C1(k)

Arguments

Value

a k-1 x k contrast matrix

Joint maximum likelihood estimation for exponential model

Description

Calculates the MLEs of the rate parameter, and joint asymptotic covariance matrix of these MLEs over 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

Joint maximum likelihood for the non-homogeneous Poisson Process

Description

Calculates the MLEs of the parameters (\mu, \sigma, \xi), and joint asymptotic covariance matrix of these MLEs over a range of thresholds as supplied by the user.

Usage

.Joint_MLE_NHPP(x, u = NULL, k, q1, q2 = 1, par, M)

Arguments

Value

a list with components

• mle matrix of MLEs above the supplied thresholds; columns are (\mu, \sigma, \xi)

• Cov.all joint asymptotic covariance matrix of all MLEs

• Cov.mu joint asymptotic covariance matrix of MLEs for \mu

• Cov.sig joint asymptotic covariance matrix of MLEs for \sigma

• Cov.xi joint asymptotic covariance matrix of MLEs for \xi

Author(s)

Jennifer L. Wadsworth

Algebraic calculation of the expected information

Description

Algebraic calculation of the expected information

Usage

.exp_info_algebraic(x, w, v, m)

Arguments

GP fitting function of Grimshaw (1993)

Description

Function for estimating parameters k and a for a random sample from a GPD.

Usage

.fit.gpd.grimshaw(x)

Arguments

Value

a list with the maximum likelihood estimates of components a and k

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 as fit.gpd, in addition to

• estimate: optimal bias-robust estimates of the scale and shape 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 x based 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 at x, depending on the value of type

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 the gpd.fit routine from ismev.

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 of ydat 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 components sigl and shl.

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 by optim. 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 of xdat).

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 returns the 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 possible 2^d-1 subsets of the power set

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

Description

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

Usage

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

Arguments

Value

an n sample from a multivariate Normal distribution

Generate from Brown-Resnick process Y \sim {P_x}, where P_{x} is probability of extremal function

Description

Generate from Brown-Resnick process Y \sim {P_x}, where P_{x} is probability of extremal function

Usage

.rPBrownResnick(index, Sigma_chol, Sigma)

Arguments

Value

a d-vector from P_x

Generate from extremal Husler-Reiss distribution Y \sim {P_x}, where P_{x} is probability of extremal function

Description

Generate from extremal Husler-Reiss distribution Y \sim {P_x}, where P_{x} is probability of extremal function

Usage

.rPHuslerReiss(index, cholesky, Sigma)

Arguments

Value

a d-vector from P_x

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

Description

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

Usage

.rPSmith(index, Sigma_chol, loc)

Arguments

Value

a d-vector from P_x

Generate from bilogistic Y \sim {P_x}, where P_{x} is the probability of extremal functions

Description

Generate from bilogistic Y \sim {P_x}, where P_{x} is the probability of extremal functions

Usage

.rPbilog(d, index, alpha)

Arguments

Value

a d-vector from P_x

Generate from extremal Dirichlet Y \sim {P_x}, where P_{x} is the probability of extremal functions from the Dirichlet model of Coles 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

a d-vector from P_x

Generate from extremal Dirichlet Y \sim {P_x}, where P_{x} is the probability of extremal functions from a Dirichlet mixture

Description

Generate from extremal Dirichlet Y \sim {P_x}, where P_{x} is the probability of extremal functions from a Dirichlet mixture

Usage

.rPdirmix(d, index, alpha, weight)

Arguments

Value

a d-vector from P_x

Generate from extremal Student-t Y \sim {P_x}, where P_{x} is probability of extremal function

Description

Generate from extremal Student-t Y \sim {P_x}, where P_{x} is probability of extremal function

Usage

.rPexstud(index, cholesky, sigma, al)

Arguments

Value

a d-vector from P_x

Generate from logistic Y \sim {P_x}, where P_{x} is probability of extremal function scaled by a Frechet variate

Description

Generate from logistic Y \sim {P_x}, where P_{x} is probability of extremal function scaled by a Frechet variate

Usage

.rPlog(d, index, theta)

Arguments

Value

a d-vector from P_x

Generate from negative logistic Y \sim {P_x}, where P_{x} is probability of extremal function scaled by a Frechet variate

Description

Generate from negative logistic Y \sim {P_x}, where P_{x} is probability of extremal function scaled by a Frechet variate

Usage

.rPneglog(d, index, theta)

Arguments

Value

a d-vector from P_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

a n by d matrix containing the sample

Generates from Q_i, the spectral measure of the bilogistic model

Description

Simulation algorithm of Boldi (2009) for the bilogistic model

Usage

.rbilogspec(n, alpha)

Arguments

Value

an n by d sample from the spectral distribution

References

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

Generates from Q_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

an n by d sample from the spectral distribution

References

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

Generates from Q_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

an n by d sample from the spectral distribution

Generates from Q_i, the spectral measure of the extremal Dirichlet model

Description

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

Usage

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

Arguments

Value

an n by d sample from the spectral distribution

References

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

Generates from Q_i, the spectral measure of the extremal Student model

Description

Generates from Q_i, the spectral measure of the extremal Student model

Usage

.rexstudspec(n, sigma, al)

Arguments

Value

an n by d sample from the spectral distribution

Generates from Q_i, the spectral measure of the Husler-Reiss model

Description

Generates from Q_i, the spectral measure of the Husler-Reiss model

Usage

.rhrspec(n, Lambda)

Arguments

Value

an n by d sample from the spectral distribution

Generates from Q_i, the spectral measure of the logistic model

Description

Simulation algorithm of Dombry et al. (2015)

Usage

.rlogspec(n, d, theta)

Arguments

Value

an n by d 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 extremal Dirichlet model. Models currently implemented include logistic, negative logistic, extremal Dirichlet and bilogistic MEV.

Usage

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

Arguments

Value

a n by d 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

a n by d 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

a n by d 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 multivariate extreme value distributions, 9, 75–81.

Random sampling from spectral distribution on l1 sphere

Description

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

Usage

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

Arguments

Value

a n by d 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. Extremes 12, 211–218.

Generates from Q_i, the spectral measure of the negative logistic model

Description

Simulation algorithm of Dombry et al. (2016)

Usage

.rneglogspec(n, d, theta)

Arguments

Value

an n by d sample from the spectral distribution

References

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

Generates from Q_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

an n by d 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 from Q_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

an n by d 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 from Q_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

an n by d sample from the spectral distribution

References

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

Generates from Q_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

an n by d 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 from Q_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

an n by d 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.

Algebraic score

Description

Algebraic score

Usage

.score_algebraic(y, x, w, v, m)

Arguments

Weighted empirical distribution function

Description

Compute an empirical distribution function with weights w at each of x

Usage

.wecdf(x, w)

Arguments

Value

a function of class ecdf

Transform variogram matrix to covariance of conditional random field

Description

The matrix Lambda is half the semivariogram matrix. The function returns the conditional covariance with respect to entries in co, restricted to the subA rows and the subB columns.

Usage

Lambda2cov(Lambda, co, subA, subB)

Arguments

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 ratio test for equality of the shape parameters over multiple 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-dimensional maximization. Other options are one-dimensional maximization of the profile the nlm function or optim. Two-dimensional optimisation using 2D-optimization ismev using the routine from gpd.fit from the ismev library, with the addition of the algebraic gradient. The choice of GP.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 thresholds u

Author(s)

Paul J. Northrop and Claire L. Coleman

References

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

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

Wadsworth & Tawn (2012). Likelihood-based procedures for threshold diagnostics 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

Given a random sample of exceedances, the estimator returns an estimator of the shape parameter or extreme value index using a kernel of order 3, based on m largest exceedances of xdat.

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, Electron. J. Statist. 17(1): 1113-1159. doi:10.1214/23-EJS2129

Examples

samp <- rgp(n = 1000, shape = 0.2)
PickandsXU(samp, m = 3)

Wadsworth's univariate and bivariate exponential threshold diagnostics

Description

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

Usage

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

Arguments

Details

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

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

Value

plots of the requested diagnostics and an invisible list with components

• MLE: maximum likelihood estimates from all thresholds

• Cov: joint asymptotic covariance matrix for \xi, \eta or 1/\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

• cthresh: vector of candidate thresholds

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

## End(Not run)

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>

Bivariate angular dependence function for extrapolation based on rays

Description

The scale parameter g(w) in the Ledford and Tawn approach is estimated empirically for x large as

\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 a common unit Pareto scale. The coefficient \eta is estimated using maximum likelihood as the shape parameter of a generalized Pareto distribution on \min(X_P, Y_P).

Usage

angextrapo(dat, qu = 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 of w

• 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, transforming the data to pseudo-observations, than marginally to unit Frechet or unit Pareto. Empirical or Euclidean weights are computed and returned alongside with the angular and radial sample for values above threshold(s) th, specified in terms of quantiles of the radial component R or marginal quantiles. Only complete tuples are kept.

Usage

angmeas(
  x,
  th,
  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 specified via region=sum.

Value

a list with arguments ang for the d-1 pseudo-angular sample, rad with the radial component and possibly wts if Rnorm='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(x=x, th=0, Rnorm='l1', Anorm='l1', marg='Frechet', wgt='Empirical')
wts2 <- angmeas(x=x, Rnorm='l2', Anorm='l2', marg='Pareto', th=0)

Dirichlet mixture smoothing of the angular measure

Description

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

Usage

angmeasdir(
  x,
  th,
  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-validation bandwidth 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\},

where f is the density of the Dirichlet distribution, p_{k, -i} is the Euclidean weight obtained from estimating the Euclidean likelihood problem without observation i.

Value

an invisible list with components

• nu bandwidth parameter obtained by cross-validation;

• dirparmat n by d matrix of Dirichlet parameters for the mixtures;

• wts mixture weights.

Examples

set.seed(123)
x <- rmev(n=100, d=2, param=0.5, model='log')
out <- angmeasdir(x=x, th=0, Rnorm='l1', Anorm='l1', marg='Frechet', wgt='Empirical')

Automated mean residual life plots

Description

This function implements the automated proposal from Section 2.2 of Langousis et al. (2016) for mean residual life plots. It returns the threshold that minimize the weighted mean square error and moment estimators for the scale and shape parameter based 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 the NOAA NCDC daily rainfall database, Water Resources Research, 52, 2659–2681.

Parametric estimates of \bar{\chi}

Description

The function fits a generalized Pareto distribution to minima of Pareto variates, using the representation

\Pr(\min(X) > x) = \frac{L(x)}{x^{1/\eta}},

where \bar{\chi}=2\eta-1. The data are transformed to the unit Pareto scale and a generalized Pareto variable is fitted to the minimum. The parameter \eta corresponds to the shape of the latter. The confidence intervals can be based either on the delta-method, a profile likelihood or a tangent exponential model approximation.

Usage

chibar(dat, confint = c("delta", "profile", "tem"), qu = 0, level = 0.95)

Arguments

Value

a named vector of length 3 containing the point estimate, the lower and the upper confidence intervals

See Also

chiplot for empirical estimates of \chi and \bar{\chi}.

Examples

## Not run: 
set.seed(765)
# Max-stable model, chibar = 1
dat <- rmev(n = 1000, model = "log", d = 2, param = 0.5)
chibar(dat, 'profile', qu = 0.5)
s <- seq(0.05,1, length = 30)
chibar_est <- t(sapply(s, function(keep){chibar(dat, 'delta', qu = keep)}))
matplot(s, chibar_est, type = 'l', col = c(1, 2, 2),  lty = c(1, 2, 2),
 ylab = expression(bar(chi)), xlab = 'p')
abline(h = 1, lty = 3, col = 'grey')
# Multivariate normal sample, chibar = 0 - strong asymptotic independence at penultimate level
dat <- mvrnorm(n = 1000, mu = c(0, 0), Sigma = cbind(c(1, 0.75), c(0.75, 1)))
chibar(dat, 'tem', q = 0.1)
chibar_est <- t(sapply(s, function(keep){chibar(dat, 'profile', qu = keep)}))
matplot(s, chibar_est, type = 'l', col = c(1, 2, 2),  lty = c(1, 2, 2),
 ylab = expression(bar(chi)), xlab = 'p')
abline(h = 1, lty = 3, col = 'grey')

## End(Not run)

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 is loc, \sigma is scale and \xi is shape.

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 and NA indicating whether observations dat fall below the mthreshold mthresh

• cl cluster instance created by makeCluster (default to NULL)

• 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 with B1

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

Value

the value of the log-likelihood with attributes expme, 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 derive level 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"),
  ...
)

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 and computes a weighted statistic comparing the squared distance with the theoretical coefficient variation corresponding to a specific shape parameter (estimated from the data using a moment estimator as the value minimizing the test statistic, or using maximum likelihood). The procedure stops if there are no more than 10 exceedances above the highest 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

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

• cthresh: sorted vector of candidate thresholds

• cindex: index of selected threshold among cthresh or NA if none returned

• pval: bootstrap p-values, with NA if lazy and the p-value exceeds level at lower thresholds

• shape: shape parameter estimates

• nexc: number of exceedances of each threshold cthresh

• 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. The parametrization assumes there is a scale parameter, say \sigma, so that scale is the distortion for the second component only. The angle rho 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

a d by d square matrix of pairwise distance

Extended generalised Pareto families

Description

This function provides the log-likelihood and quantiles for the three different families presented in Papastathopoulos and Tawn (2013). The latter include an additional parameter, \kappa. All three families share the same tail index as the generalized Pareto distribution, while allowing for lower thresholds. In the case \kappa=1, the models reduce to the generalised Pareto.

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

Arguments

Details

For return levels, the p argument can be related to T year exceedances as follows: if there are n_y observations per year, than take p to equal 1/(Tn_y) to obtain the T-years return level.

Value

egp.ll returns the log-likelihood value.

egp.retlev returns a plot of the return levels if plot=TRUE and a matrix of return levels.

Usage

egp.ll(xdat, thresh, par, model=c('egp1','egp2','egp3'))

egp.retlev(xdat, thresh, par, model=c('egp1','egp2','egp3'), 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 Inference 143(3), 131–143.

Examples

set.seed(123)
xdat <- mev::rgp(1000, loc = 0, scale = 2, shape = 0.5)
par <- fit.egp(xdat, thresh = 0, model = 'egp3')$par
p <- c(1/1000, 1/1500, 1/2000)
#With multiple thresholds
th <- c(0, 0.1, 0.2, 1)
opt <- tstab.egp(xdat, th, model = 'egp1')
egp.retlev(xdat, opt$thresh, opt$par, 'egp1', p = p)
opt <- tstab.egp(xdat, th, model = 'egp2', plots = NA)
egp.retlev(xdat, opt$thresh, opt$par, 'egp2', p = p)
opt <- tstab.egp(xdat, th, model = 'egp3', plots = NA)
egp.retlev(xdat, opt$thresh, opt$par, 'egp3', p = p)

Extended generalised Pareto families of Papastathopoulos and Tawn (functions)

Description

This function provides the log-likelihood and quantiles for the three different families presented in Papastathopoulos and Tawn (2013). The latter include an additional parameter, \kappa. All three families share the same tail index than the GP model, while allowing for lower thresholds.

Usage

egp.ll(xdat, thresh, par, model = c("egp1", "egp2", "egp3"))

egp.retlev(
  xdat,
  thresh,
  par,
  model = c("egp1", "egp2", "egp3"),
  p,
  plot = TRUE
)

Arguments

Fit of extended GP models and parameter stability plots

Description

This function is an alias of fit.egp.

Usage

egp.fit(xdat, thresh, model = c("egp1", "egp2", "egp3"), 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
)

Fit an extended generalized Pareto distribution of Naveau et al.

Description

Deprecated function name to fit an extended generalized Pareto family. The user should call fit.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

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 length d).

• wts n 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 daily cumulated 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 list par must contain different arguments depending on the model. For the Brown–Resnick model, the user must supply the conditionally negative definite matrix Lambda following the parametrization in Engelke et al. (2015) or the covariance matrix Sigma, following Wadsworth and Tawn (2014). For the Husler–Reiss model, the user provides the mean and covariance matrix, m and Sigma. For the extremal student, the covariance matrix Sigma and the degrees of freedom df. For the logistic model, the strictly positive dependence parameter alpha.

Examples

## Not run: 
# Extremal Student
Sigma <- 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 model
D <- 5L
loc <- 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] / 2
expme(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 based on interexceedance time and gap of exceedances

Description

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

Usage

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

Arguments

Details

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

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

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

Value

a vector or matrix of estimated extremal index of dimension length(method) by length(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 climatology for extreme hourly precipitation in Switzerland. Theoretical and Applied Climatology, 120(3), 403-416.

Examples

set.seed(234)
#Moving maxima model with theta=0.5
a <- 1; theta <-  1/(1+a)
sim <- rgev(10001, loc=1/(1+a),scale=1/(1+a),shape=1)
x <- pmax(sim[-length(sim)]*a,sim[-1])
q <- seq(0.9,0.99,by=0.01)
ext.index(xdat=x,q=q,method=c('wls','mle'))

Estimators of the extremal coefficient

Description

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

Usage

extcoef(
  dat,
  coord = NULL,
  thresh = NULL,
  estimator = c("schlather", "smith", "fmado"),
  standardize = TRUE,
  method = c("nonparametric", "parametric"),
  prob = 0,
  plot = TRUE,
  ...
)

Arguments

Details

The Smith estimator: suppose Z(x) is simple max-stable vector (i.e., with unit Frechet marginals). Then 1/Z is unit exponential and 1/\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 using the reciprocal mean.

The Schlather and Tawn estimator: the likelihood of the naive estimator for a pair of two sites A 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 and z is a threshold on the unit Frechet scale. The search for the maximum likelihood estimate for every pair A is 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 case the Smith estimator if we do not censor the data (p = 0) and do not standardize observations by their harmonic mean.

The F-madogram estimator is a non-parametric estimate based on a stationary process Z; 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 significant for large data sets. If there are no missing observations, the routine fmadogram from the SpatialExtremes package should be prefered as it is noticeably faster.

The data will typically consist of max-stable vectors or block maxima. Both of the Smith and the Schlather–Tawn estimators require unit Frechet margins; the margins will be standardized to the unit Frechet scale, either parametrically or nonparametrically unless standardize = FALSE. If method = "parametric", a parametric GEV model is fitted to each column of dat using maximum likelihood estimation and transformed back using the probability integral transform. If method = "nonparametric", using the empirical distribution function. The latter is the default, as it is appreciably faster.

Value

an invisible list with vectors dist if coord is non-null or else a matrix of pairwise indices ind, extcoef and the supplied estimator, fmado and binned. If estimator == "schlather", an additional matrix with 2 columns containing the binned distance binned with the h 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(dat = dat, coord = coord)
# Extremal Student extremal coefficient function

XT.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 coef
plot(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 function
BR.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(dat = 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 Naveau et al. (2016) retain the tail index of the distribution while being compliant with the theoretical behavior of extreme low 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 carrier G(u)=u^\kappa.

• type 2 corresponds to a three parameters family, with carrier G(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 of type 2, with G(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 number generation 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 pairwise chi estimates and plots them as a function of the distance between sites.

Usage

extremo(dat, 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 grid
grid <- 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 grid
eta <- 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 variogram
Lambda <- ((dianiso/30)^0.7)/4
# Regular Euclidean distance between sites
di <- distg(coord, 1, 0)
# Simulate generalized max-Pareto field
set.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 variogram
mindistpvario <- 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,
                        fn = function(par){
                          mindistpvario(par, emp = extdat[,'prob'], coord = coord)})
stopifnot(opt$kkt1, opt$kkt2)
# Plotting the extremogram in the deformed space
distfa <- distg(loc = coord, opt$par[3], opt$par[4])
plot(
 x = c(distfa[lower.tri(distfa)]), 
 y = extdat[,2], 
 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)), 
  col = 2, lwd = 2,
  y = 2*(1-pnorm(sqrt(power.vario(distvec, alpha = opt$par[1],
                               scale = opt$par[2])/2))))

## End(Not run)

Parameter stability plot and maximum likelihood routine for extended GP models

Description

The function tstab.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, model = c("egp1", "egp2", "egp3"), init, show = FALSE)

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

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 by optim, 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 Inference 143(3), 131–143.

Examples

xdat <- mev::rgp(
  n = 100,
  loc = 0,
  scale = 1,
  shape = 0.5)
fitted <- fit.egp(
  xdat = xdat,
  thresh = 1,
  model = "egp2",
  show = TRUE)
thresh <- mev::qgp(seq(0.1, 0.5, by = 0.05), 0, 1, 0.5)
tstab.egp(
   xdat = xdat,
   thresh = thresh,
   model = "egp2",
   plots = 1:3)

Fit an extended generalized Pareto distribution of Naveau et al.

Description

This is a wrapper function to obtain PWM or MLE estimates for the extended GP models of Naveau et al. (2016) for rainfall intensities. The function calculates confidence intervals by means of nonparametric percentile bootstrap and returns histograms and QQ plots of the 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 carrier G(u)=u^\kappa.

• model 2 corresponds to a three parameters family, with carrier G(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 of type 2, with G(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, gp.fit(rain, 0)$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 class mev_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 parameter estimate.

• convergence components taken from the list returned by auglag. Values other than 0 indicate that the algorithm likely did not converge.

• counts components taken from the list returned by auglag.

• xdat vector of data

Examples

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

Maximum likelihood estimation for the generalized Pareto distribution

Description

Numerical optimization of the generalized Pareto distribution for data exceeding threshold. This function returns an object of class mev_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 optimal B-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 using nlm and optim. Method 'ismev' performs the two-dimensional optimization routine gpd.fit from the ismev 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 importance sampling assuming a GPD prior is placed on the parameter of the profile likelihood.

Value

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

• estimate a vector containing the scale and shape parameters (optimized and fixed).

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

• vcov the variance covariance matrix, obtained as the numerical inverse of the observed information matrix. For method = "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 parameter estimate.

• nat number of points lying above the threshold.

• pat proportion of points lying above the threshold.

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

• counts components taken from the list returned by optim.

• exceedances excess over the threshold.

Additionally, if method = "obre", a vector of OBRE weights.

Otherwise, a list containing

• threshold the threshold.

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

• 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 neither FALSE, nor NULL

• 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 within gp.fit, but users can call internal functions should they wish to use these features.

Author(s)

Scott D. Grimshaw for the Grimshaw option. Paul J. Northrop and Claire L. Coleman for the methods optim, nlm and ismev. J. Zhang and Michael A. Stephens (2009) and Zhang (2010) for the zs and zhang approximate methods and L. Belzile for methods auglag and obre, the wrapper and MCMC samplers.

If show = TRUE, the optimal B robust estimated weights for the largest observations are printed alongside with the p-value of the latter, obtained from the empirical distribution of the weights. This diagnostic can be used to guide threshold selection: small weights for the r-largest order statistics indicate that the robust fit is driven by the lower tail and that the threshold should perhaps be increased.

References

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

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

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

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

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

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

See Also

fpot and gpd.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 above threshold is modelled using the limiting point process of 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 parameter npp controls the frequency of observations. If data are recorded on a daily basis, using a value of npp = 365.25 yields location and scale parameters that correspond to those of the generalized 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 parameter estimate.

• nat number of points lying above the threshold.

• pat proportion of points lying above the threshold.

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

• counts components taken from the list returned by optim.

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)

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

Description

This uses a constrained optimization routine to return the maximum likelihood estimate based on an n by r matrix of observations. Observations should be ordered, i.e., the r-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 parameter estimate.

• convergence components taken from the list returned by auglag. Values other than 0 indicate that the algorithm likely did not converge.

• counts components taken from the list returned by auglag.

• xdat an n by r matrix of data

Examples

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

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

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

Source

European Climate Assessment and Dataset project https://www.ecad.eu/

References

Klein Tank, A.M.G. and Coauthors, 2002. Daily dataset of 20th-century surface air temperature and precipitation series for the European 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, see http://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 derivative for 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 shape in 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 are n_y observations per year, the L-year return level is obtained by taking N equal to n_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 using Firth's modified score function or implicit bias subtraction.

Usage

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

Arguments

Details

Method subtractsolves

\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 by gev.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}}),

where U is the score vector, b is the first order bias and i is either the observed or Fisher information.

The routine uses the MLE (bias-corrected) as starting values and proceeds to find the solution using a root finding algorithm. Since the bias-correction is not valid for \xi < -1/3, any solution that is unbounded will return a vector of NA 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)$estimate
gev.bcor(par, dat, 'subtract')
gev.bcor(par, dat, 'firth') #observed information
gev.bcor(par, dat, 'firth','exp')

Cox-Snell first order bias for the GEV distribution

Description

Bias vector for the GEV distribution based on an n sample. Due to numerical instability, values of the information matrix and the bias are 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

The function returns the expected or observed information matrix.

Usage

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

Arguments

log likelihood for the generalized extreme value distribution

Description

Function returning the density of an n 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 the fit.gev routine on the sample of block maxima and returns maximum likelihood estimates for all quantities of interest, including location, scale and shape parameters, quantiles and mean and quantiles of maxima of N 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 arguments args

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 corrections based on Severini's (1999) empirical covariance and the Fraser and Reid tangent exponential model 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 additional mod available are tem, the tangent exponential model (TEM) approximation and modif for the penalized profile likelihood based on p^* 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 of psi.max

• psi: vector of parameter \psi given in psi

• pll: values of the profile log likelihood at psi

• maxpll: value of maximum profile log likelihood

In addition, if mod includes tem

• 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, if mod includes modif

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

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

• 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 at psi

• 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 the 1-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

Tangent exponential model approximation for the GEV distribution

Description

The function gev.tem provides a tangent exponential model (TEM) approximation for 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 for psi that will cover 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 class fr (see tem in package hoa) 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 the log 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 derivative for the generalized extreme value distribution parametrized in terms of the quantiles/mean of N-block maxima parametrization z, 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 and random number generation for the generalized extreme value distribution.

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 parameters loc = \mu, scale = \sigma and shape = \xi is

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

for 1 + \xi (x - \mu) / \sigma > 0. If \xi = 0 the distribution 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 the annual 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 Statistical Modeling 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 derivative for the generalized extreme value distribution parametrized in terms of the return level z, 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 function fit.gpd is a wrapper around gp.fit

Usage

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

Arguments

Generalized Pareto distribution

Description

Likelihood, score function and information matrix, bias, approximate ancillary statistics and sample space derivative for 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 of log(scale) and shape in 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 local exponential 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 using Firth's modified score function or implicit bias subtraction.

Usage

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

Arguments

Details

Method subtract 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 by gpd.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}}),

where U is the score vector, b is the first order bias and i is either the observed or Fisher information.

The routine uses the MLE as starting value and proceeds to find the solution using a root finding algorithm. Since the bias-correction is not valid for \xi < -1/3, any solution that is unbounded will return a vector of NA 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)$estimate
gpd.bcor(par,dat, 'subtract')
gpd.bcor(par,dat, 'firth') #observed information
gpd.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 an n 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 class mev_gpd, returns a matrix of parameter values to mimic the estimation uncertainty.

Usage

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

Arguments

Details

Two options are available: a normal approximation to the scale and shape based on the maximum likelihood estimates and the observed information matrix. This method uses forward sampling to simulate from a bivariate normal distribution that satisfies the support and positivity constraints

The second approximation uses the ratio-of-uniforms method to obtain samples from the posterior distribution with uninformative priors, thus mimicking the joint distribution of maximum likelihood. The benefit of the latter is that it is more reliable in 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 an n sample from the GP distribution. gpd.ll.optim returns the negative log likelihood parametrized in terms of log(scale) and shape in 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

Generalized Pareto maximum likelihood estimates for various quantities of interest

Description

This function calls the fit.gpd routine on the sample of excesses and returns maximum likelihood estimates for all quantities of interest, including scale and shape parameters, quantiles and value-at-risk, expected shortfall and mean and quantiles of maxima of N 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 arguments args

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 the p^* formula. There are two small-sample corrections that use a proxy for \ell_{\lambda; \hat{\lambda}}, which are based on Severini's (1999) empirical covariance and the Fraser and Reid tangent exponential model approximation.

Usage

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

Arguments

Details

The three mod available are profile (the default), tem, the tangent exponential model (TEM) approximation and modif for the penalized profile likelihood based on p^* 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 of psi.max

• psi: vector of parameter \psi given in psi

• pll: values of the profile log likelihood at psi

• maxpll: value of maximum profile log likelihood

• family: a string indicating "gpd"

• threshold: value of the threshold, by default zero

In addition, if mod includes tem

• 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, if mod includes modif

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

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

• 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 at psi

• 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 function gpd.tem provides a tangent exponential model (TEM) approximation for higher order likelihood inference for a scalar parameter for the generalized Pareto distribution. Options include scale 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 for psi that will cover the range of options, but the fit may fail and return an error. In such cases, the user can try to find good grid of starting values and provide them to the routine.

Usage

gpd.tem(
  dat,
  param = c("scale", "shape", "quant", "VaR", "ES", "Nmean", "Nquant"),
  psi = NULL,
  m = NULL,
  threshold = 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 around gpd.pll.

The interpretation for m is as follows: if there are on average m_y observations per year above the threshold, then m = Tm_y corresponds to T-year return level.

Value

an invisible object of class fr (see tem in package hoa) 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 plots
m1 <- 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 psi
psi <- 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 the log 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 derivative for the generalized Pareto distribution parametrized in terms of average maximum of N exceedances.

The parameter N 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 of N 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 derivative for the generalized Pareto distribution parametrized in terms of expected shortfall.

The parameter m corresponds to \zeta_u/(1-\alpha), where \zeta_u is the rate of exceedance over the threshold u 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 expected shortfall 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) - optimization The 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 and random number generation for the generalized Pareto distribution.

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 Statistical Modeling 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 derivative for 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 for m is as follows: if there are on average m_y observations per year above the threshold, then m=Tm_y corresponds to T-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 of log(scale) and shape in order to perform unconstrained optimization

• gpdr.score: score vector

• gpdr.infomat: observed information matrix for GPD parametrized in terms of rate of m-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 of m-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 as dat with unit Pareto observations

Interpret bivariate threshold exceedance models

Description

This is an adaptation of the evir package interpret.gpdbiv function. interpret.fbvpot deals with the output of a call to fbvpot from the evd 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 list fitted must contain

• model a string; see bvevd from package evd for options

• param a named vector containing the parameters of the model, as well as parameters scale1, shape1,scale2 and shape2, corresponding to marginal GPD parameters.

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

• pat the proportion of observations above the corresponding threshold

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 package evir

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 based on the difference between the expected quadratic score and the second derivative of the log-likelihood. The asymptotic distribution for each threshold u and gap K is asymptotically \chi^2 with one degree of freedom. The approximation is good for N>80 and conservative for smaller sample sizes. The test assumes independence between gaps.

Usage

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

Arguments

Details

The procedure proposed in Suveges & Davison (2010) was corrected for erratas. The maximum likelihood is based on the limiting mixture distribution of the intervals between exceedances (an exponential with a point mass at zero). The condition D^{(K)}(u_n) should be checked by the user.

Fukutome et al. (2015) propose an ad hoc automated procedure

1. Calculate the interexceedance times for each K-gap and each threshold, along with the number of clusters

2. Select the (u, K) pairs for which IMT < 0.05 (corresponding to a P-value of 0.82)

3. Among those, select the pair (u, K) for which the number of clusters is the largest

Value

an invisible list of matrices containing

• IMT a matrix of test statistics

• pvals a matrix of approximate p-values (corresponding to probabilities under a \chi^2_1 distribution)

• mle a matrix of maximum likelihood estimates for each given pair (u, K)

• loglik a matrix of log-likelihood values at MLE for each given pair (u, K)

• threshold a vector of thresholds based on empirical quantiles at supplied levels.

• q the vector q supplied by the user

• K the largest gap number, supplied by the user

Author(s)

Leo Belzile

References

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

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

White (1982), Maximum Likelihood Estimation of Misspecified Models. Econometrica, 50(1), 1-25.

Examples

infomat.test(xdat = rgp(n = 10000),
             q = seq(0.1, 0.9, length = 10),
             K = 3)

Intensity function for the Brown-Resnick model

Description

The intensity function includes the normalizing constants

Usage

intensBR(tdat, Lambda, cholPrecis = NULL)

Arguments

Value

log intensity contribution

Intensity function for the extremal Student model

Description

The intensity function includes the normalizing constants

Usage

intensXstud(tdat, df, Sigma, cholPrecis = NULL)

Arguments

Value

log intensity contribution

Jacobian of the transformation from generalized Pareto to unit Pareto distribution

Description

If dat is a vector, the arguments loc, scale and shape should be numericals of length 1.

Usage

jac(dat, loc = 0, scale, shape, lambdau = 1, censored)

Arguments

Value

log-likelihood contribution for the Jacobian

Estimation of the bivariate angular dependence function of Wadsworth and Tawn (2013)

Description

Estimation of the bivariate angular dependence function of Wadsworth and Tawn (2013)

Usage

lambdadep(dat, qu = 0.95, method = c("hill", "mle", "bayes"), plot = TRUE)

Arguments

Value

a plot of the lambda function if plot=TRUE, plus an invisible list with components

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

• lambda point estimates of lambda

• lower.confint 95

• upper.confint 95

Examples

set.seed(12)
dat <- mev::rmev(n = 1000, d = 2, model = "log", param = 0.1)
lambdadep(dat, method = 'hill')
## Not run: 
lambdadep(dat, method = 'bayes')
lambdadep(dat, method = 'mle')
# With independent observations
dat <- matrix(runif(n = 2000), ncol = 2)
lambdadep(dat, method = 'hill')

## End(Not run)

Likelihood for multivariate peaks over threshold models

Description

Likelihood for the 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 is loc, \sigma is scale and \xi is shape.

Usage

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

Arguments