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

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

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

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)

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

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.

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

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.

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

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)

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

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

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)

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

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)

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)

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

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)

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

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

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

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

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