Skip to contents

Smith's penultimate approximations

Source: R/penultimate.R
penultimate.Rd

The function takes as arguments the distribution and density functions. There are two options: method='bm' yields block maxima and method='pot' threshold exceedances. For method='bm', the user should provide in such case the block sizes via the argument m, whereas if method='pot', a vector of threshold values should be provided. The other argument (thresh or m depending on the method) is ignored.

Usage

penultimate(family, method = c("bm", "pot"), thresh, qlev, m, ...)

Arguments

family

the name of the parametric family. Will be used to obtain dfamily, pfamily, qfamily

method

either block maxima ('bm') or peaks-over-threshold ('pot') are supported

thresh

vector of thresholds for method 'pot'

qlev

vector of quantile levels for method 'pot', e.g., 0.9, 0.95, ... Ignored if argument thresh is provided.

m

vector of block sizes for method 'bm'

...

additional arguments passed to densF and distF

Value

a data frame containing

  • loc: location parameters (method='bm')

  • scale: scale parameters

  • shape: shape parameters

  • thresh: thresholds (if method='pot'), percentile corresponding to threshold (if method='pot')

  • m: block sizes (if method='bm')

Details

Alternatively, the user can provide functions densF, quantF and distF for the density, quantile function and distribution functions, respectively. The user can also supply the derivative of the density function, ddensF. If the latter is missing, it will be approximated using finite-differences.

For method = "pot", the function computes the reciprocal hazard and its derivative on the log scale to avoid numerical overflow. Thus, the density function should have argument log and the distribution function arguments log.p and lower.tail, respectively.

References

Smith, R.L. (1987). Approximations in extreme value theory. Technical report 205, Center for Stochastic Process, University of North Carolina, 1–34.

Author

Leo Belzile

Examples

# Threshold exceedance for Normal variables
quants <- seq(1, 5, by = 0.02)
penult <- penultimate(
   family = "norm",
   method = 'pot',
   thresh = quants,
   ddensF = function(x){-x*dnorm(x)}, # optional argument
   )
plot(x = quants,
     y = penult$shape,
     type = 'l',
     xlab = 'quantile',
    ylab = 'Penultimate shape',
    ylim = c(-0.5, 0))

# Block maxima for Gamma variables
# User must provide arguments for shape (or rate), for which there is no default
m <- seq(30, 3650, by = 30)
penult <- penultimate(family = 'gamma', method = 'bm', m = m, shape = 0.1)
plot(x = m,
     y = penult$shape,
     type = 'l',
     xlab = 'quantile',
     ylab = 'penultimate shape')


# Comparing density of GEV approximation with true density of maxima
m <- 100 # block of size 100
p <- penultimate(
  family = 'norm',
  ddensF = function(x){-x*dnorm(x)},
  method = 'bm',
  m = m)
x <- seq(1, 5, by = 0.01)
plot(
  x = x,
  y = m * dnorm(x) * exp((m-1) * pnorm(x, log.p = TRUE)),
  type = 'l',
  ylab = 'density',
  main = 'Distribution of the maxima of\n 100 standard normal variates')
lines(x, mev::dgev(x, loc = p$loc, scale = p$scale, shape = 0), col = 2)
lines(x, mev::dgev(x, loc = p$loc, scale = p$scale, shape = p$shape), col = 4)
legend(
 x = 'topright',
 lty = c(1, 1, 1),
 col = c(1, 2, 4),
 legend = c('exact', 'ultimate', 'penultimate'),
 bty = 'n')