Package 'mev'

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, Swedish Meteorological and Hydrological Institute, distributed under a Creative Commons Attribution 4.0 SE license

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>

Description

Estimation of the bivariate angular dependence function

Usage

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

Arguments

Value

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

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

• coef point estimates of the angular dependence function

• lower level% confidence interval for lambda (lower bound)

• upper level% confidence interval for lambda (upper bound)

References

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

Examples

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

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) thresh, specified in terms of quantiles of the radial component R or marginal quantiles. Only complete tuples are kept.

Usage

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

Arguments

Details

The empirical likelihood weighted mean problem is implemented for all thresholds, while the Euclidean likelihood is only supported for diagonal thresholds 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(xdat = x, Rnorm = 'l1', Anorm = 'l1', marg = 'frechet', wgt = 'empirical')
wts2 <- angmeas(xdat = x, Rnorm = 'l2', Anorm = 'l2', marg = 'pareto')

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(
  xdat,
  thresh,
  Rnorm = c("l1", "l2", "linf"),
  Anorm = c("l1", "l2", "linf", "arctan"),
  marg = c("frechet", "pareto"),
  wgt = c("empirical", "euclidean"),
  region = c("sum", "min", "max"),
  is.angle = FALSE,
  ...
)

Arguments

Details

The cross-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 = 2L, param = 0.5, model = 'log')
out <- angmeasdir(x)

Description

Given a time series of observations in xdat, compute the maximum of blocks of size block ($b$), and then order them by further blocks of size m, increasing by row from left to right. If the length of xdat is not a multiple of block, the last observations are discarded without warning.

Usage

build.blocks(xdat, block = 1L, m = 2L)

Arguments

Value

a matrix with $\lfloor n/b \rfloor$ observations, ordered by row, with m columns.

Description

Daily measurements of wind speed during the month of January and including February 1st, from the Cheeseboro (California) weather station. between 1996 and 2026.

Usage

cheeseborowind

Format

A data frame with 992 rows and 3 variables:

date

date of observation

direction

angle (in degrees) of the wind

gust

maximum daily wind speed (in meters per second)

Source

Raw US Climate Archive, https://raws.dri.edu/cgi-bin/rawMAIN.pl?caCCHB, maintained by the Western Regional Climate Center, Desert Research Institute based in Reno, Nevada

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"),
  boundary = FALSE,
  ...
)

Arguments

Value

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

Description

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

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

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

Usage

dgeoaniso(loc, theta)

Arguments

Value

a d by d square matrix of pairwise distance

References

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

Description

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

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

Arguments

Details

For return levels, 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, while egp.retlev returns a plot of the return levels if plot=TRUE and a list with tail probabilities p, return levels retlev, thresholds thresh and model name model.

Usage

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

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

Author(s)

Leo Belzile

References

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

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

Examples

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

Description

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

Usage

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

Arguments

Value

an object of class eprof

Description

Density function, distribution function, quantile function and random number generation for various extended generalized Pareto distributions

Usage

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

degp(
  x,
  scale,
  shape,
  kappa,
  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",
    "logist"),
  log = FALSE
)

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

regp(
  n,
  scale,
  shape,
  kappa,
  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",
    "logist")
)

Arguments

References

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

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

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.

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.

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)

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.

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

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 = 0,
  model = c("pt-beta", "pt-gamma", "pt-power", "gj-tnorm", "gj-beta", "exptilt",
    "logist"),
  start = NULL,
  method = c("Nelder", "nlminb", "BFGS"),
  fpar = NULL,
  show = FALSE,
  ...
)

Arguments

Details

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

Value

fit.egp outputs the list returned 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 = "pt-gamma",
  show = TRUE)
thresh <- mev::qgp(seq(0.1, 0.5, by = 0.05), 0, 1, 0.5)
tstab.egp(
   xdat = xdat,
   thresh = thresh,
   model = "pt-gamma")
xdat <- regp(
  n = 100,
  scale = 1,
  shape = 0.1,
  kappa = 0.5,
  model = "pt-power"
)
fit.egp(
 xdat = xdat,
 model = "pt-power",
 show = TRUE,
 fpar = list(kappa = 1),
 method = "Nelder"
)

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, fit.gpd(rain)$est),
  rounded = 0.1,
  confint = TRUE,
  R = 20
)

## End(Not run)

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))

Description

Given a matrix of n ordered samples of m order statistics from a postulated GEV, fit the parameters of the latter based on the marginal likelihood of the first m-1 order statistics using maximum likelihood.

Usage

fit.gevblock(
  xdat,
  marginal = FALSE,
  constraint = TRUE,
  rounding = 0,
  lb = NULL,
  start = NULL,
  vcov = FALSE
)

Arguments

Details

One can set constraint to TRUE to add a support constraint to the optimization to ensure that all values of xdat are in the support of the resulting distribution (only for the marginal likelihood).

Value

(constrained) maximum likelihood estimator of location, scale and shape parameters

Examples

set.seed(2026)
xdat <- build.blocks(mev::rgev(n = 200, shape = 0.1), m = 4)
fit.gevblock(xdat, marginal = TRUE)
fit.gevblock(round(xdat, 1), marginal = TRUE, lb = NULL, rounding = 0.1)
fit.gevblock(round(xdat, 1), marginal = TRUE, lb = -2, rounding = 0.1)
fit.gevblock(xdat, marginal = TRUE, lb = -2)
fit.gevblock(xdat)
fit.gevblock(round(xdat, 1), lb = NULL, rounding = 0.1)
fit.gevblock(round(xdat, 1), lb = -2, rounding = 0.1)
fit.gevblock(xdat, lb = -2)

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 = c("Grimshaw", "auglag", "nlm", "optim", "ismev", "zs", "zhang", "obre", "pwm"),
  show = FALSE,
  MCMC = NULL,
  k = 4,
  tol = 1e-08,
  fpar = NULL,
  warnSE = FALSE,
  returnsamp = TRUE,
  ...
)

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)

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)

Description

Estimator of the second order tail index parameter

Usage

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

Arguments

Examples

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

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)

Description

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

Usage

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

Arguments

Value

a data frame with the number of order statistics k and the shape parameter estimate shape, or a single numeric value if k is a scalar.

Description

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

Usage

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

Arguments

Value

a list with components

• estimate a vector containing the scale and shape 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 logical indicator of convergence.

• weights vector of weights for exceedances.

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

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")

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>.

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.

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

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')

Description

Given an object of class mev_gev, returns a matrix of parameter values to mimic the estimation uncertainty.

Usage

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

Arguments

Details

Two options are available: a normal approximation to the location, scale and shape based on the maximum likelihood estimates and the observed information matrix. This method uses forward sampling to simulate from a trivariate 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 3 whose columns contain scale and shape parameters

Examples

set.seed(2025)
xdat <- rgev(100, loc = 0, scale = 2, shape = -0.1)
fgev <- fit.gev(xdat)
pairs(gev.boot(fgev, method = "post"))
pairs(gev.boot(fgev, method = "norm"))

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)

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

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)

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)

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

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

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

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.

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

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')

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

Examples

set.seed(2025)
xdat <- rgev(100, loc = 0, scale = 2, shape = -0.1)
fgp <- fit.gpd(xdat)
plot(
 gpd.boot(fgp, method = "post")
)
points(
 gpd.boot(fgp, method = "norm"),
 col = 2,
 pch = 20
)

Description

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

Usage

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

Arguments

Value

a vector of length two with the scale and shape estimates

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)

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", "retlev", "VaR", "ES", "Nmean", "Nquant"),
  mod = "profile",
  mle = NULL,
  dat,
  m = NULL,
  N = NULL,
  p = NULL,
  q = NULL,
  correction = TRUE,
  thresh = NULL,
  plot = TRUE,
  ...
)

Arguments

Details

The 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"

• thresh: 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)

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", "retlev", "ES", "Nmean", "Nquant"),
  psi = NULL,
  m = NULL,
  thresh = 0,
  n.psi = 50,
  N = NULL,
  p = NULL,
  q = NULL,
  plot = FALSE,
  correction = TRUE,
  ...
)

Arguments

Details

As of version 1.11, this function is a wrapper 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)

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