library(mev)
# library(xts)
# library(lubridate)
data(frwind, package = "mev")
<- with(frwind,
lyon ::xts(x = S2, order.by = date))
xts# Create series of yearly maximum
<- xts::apply.yearly(lyon, max) ymax
Likelihood-based inference
The mev
package provides gradient-based optimization routines for fitting univariate extreme value models, either block maxima or threshold exceedances, using one of four likelihoods: that of the generalized extreme value distribution, the generalized Pareto distribution, and the inhomogeneous Poisson point process and the
Relative to other packages such as evd
or ismev
, the package functions include analytic expressions for the score and observed informations, with careful interpolation when mev
does not handle generalized linear or generalized additive models for the parameters, to avoid having as many inequality constraints in the optimization as there are observations times the number of covariates.
Basic theory
Let
By definition, the maximum likelihood estimator solves the score equation, i.e.
Statistical inference
This section presents some test statistics that can easily be computed using some of the functionalities of mev
, as well as confidence intervals for parameters and common functionals, based on the profile likelihood.
The three main type of test statistics for likelihood-based inference are the Wald, score and likelihood ratio tests. The three main classes of statistics for testing a simple null hypothesis
Oftentimes, we are interested in a functional of the parameter vector
We can define score and information in the usual fashion: for example, the observed profile information function is [j_() =- = {j^{}(, _{})}^{-1}. ]
We can turn tests and their asymptotic distribution into confidence intervals. For the hypothesis
Likelihoods
There are four basic likelihoods for univariate extremes: the likelihood of the generalized extreme value (GEV) distribution for block maxima, the likelihood for the generalized Pareto distribution and that of the non-homogeneous Poisson process (NHPP) for exceedances above a threshold
Generalized extreme value distribution
The generalized extreme value (GEV) distribution with location parameter
The max-stability property allows one to extrapolate the distribution beyond observed levels: one can show that the distribution of the maximum of a larger block (or
The GEV distribution is suitable for maximum of a large number of observations: the larger the block size, the closer the approximation will be, but the smaller the sample size
The fit.gev
function includes two optimization routines: either use the PORT methods from nlminb
, or Broyden-Fletcher-Goldfarb-Shanno algorithm (BFGS
) inside a constrained optimization algorithm (augmented Lagrangian). The default option is nlminb
, which sometimes returns diagnostics indicating false convergence when the model is near the maximum likelihood estimate.
As for other model, parameters can be fixed and nested models can be compared using the anova
S3 method. For these, we distinguish between estimated coefficients (estimate
) or with the coef
method, and the full vector of parameters, param
.
Numerical example
We consider in the tutorial daily mean wind speed data, measured in km/h, at four weather stations located in the south of France. We first take the data for Lyon’s airport (station
We can then fit a GEV distribution via maximum likelihood to the yearly maximum, extract the coefficients
<- mev::fit.gev(xdat = ymax, show = TRUE) opt_gev
Log-likelihood: -141.6626
Estimates
loc scale shape
36.18449 3.94287 -0.01124
Standard Errors
loc scale shape
0.6589 0.4881 0.1318
Optimization Information
Convergence: successful
Function Evaluations: 27
Gradient Evaluations: 11
<- coef(opt_gev)
mle isTRUE(all.equal(rep(0,3),
::gev.score(par = mle, dat = ymax),
mevcheck.attributes = FALSE,
tolerance = 1e-5))
[1] TRUE
Having found the MLE, we can compute the covariance matrix of the parameters from the observed information matrix. The standard errors are the square root of the elements on the diagonal. While mev
uses exact formulae, these can be approximated by computing the hessian via finite differences.
# Compute observed information matrix
<- mev::gev.infomat(par = mle, dat = ymax)
jmat # Compute standard errors
<- sqrt(diag(solve(jmat)))
se_mle # Compare with 'mev' output
isTRUE(all.equal(se_mle, opt_gev$std.err))
[1] TRUE
Even if we have parameter estimates, there is no guarantee that the model is adequate. Standard visual goodness-of-fit diagnostics can be obtained with the plot
method. To see other methods, query methods(class = "mev_gev")
.
# PP and QQ plots
par(mfrow = c(1,2))
plot(opt_gev)
graphics.off()
The Gumbel distribution, which corresponds to a GEV with shape
<- mev::fit.gev(xdat = ymax,
opt_gumb fpar = list(shape = 0))
anova(opt_gev, opt_gumb)
Analysis of Deviance Table
npar Deviance Df Chisq Pr(>Chisq)
opt_gev 3 283.32
opt_gumb 2 283.33 1 0.0073 0.9321
None of the parameters are of interest in themselves. We may be interested rather by a risk summary, which is a function of parameters. For example, we could get the parameters of the GEV for 50 years via max-stability and return the average if
All of these are invariant to reparametrization, so we can use the formula and plug-in the parameter values. For inference, we reparametrize the model in terms of this quantity, then vary over a grid of values of the 50-year average maximum and compute profile-likelihood-based confidence intervals at level 95%.
gev.mle(xdat = ymax, args = "Nmean", N = 50)
Nmean
53.4114
# Compute profile log-likelihood
<- mev::gev.pll(param = "Nmean", dat = ymax, N = 50) prof
# Extract confidence intervals
confint(prof)) (
Estimate Lower CI Upper CI
53.41140 47.79346 73.63881
Generalized Pareto distribution
The generalized Pareto (GP) distribution with scale
Except for this boundary case, the maximum likelihood estimator solves the score equation
If mev
does not report them and prints an optional warning.
The figure shows the profile likelihood for
# Only keep data from September to April
<- with(
windlyon
frwind, ::month(date) <= 4 |
S2[lubridate::month(date) >= 9])
lubridate# Keep only 100 largest points (fewer because of ties)
<- quantile(windlyon,
u probs = 1-100/length(windlyon))
# Fit generalized Pareto via ML
<- fit.gpd(
fitted_gp xdat = windlyon,
threshold = u,
show = TRUE)
Method: Grimshaw
Log-likelihood: -207.5276
Threshold: 33.84
Number Above: 90
Proportion Above: 0.0079
Estimates
scale shape
3.57863 0.03088
Standard Errors
scale shape
0.6091 0.1337
Optimization Information
Convergence: successful
# P-P and Q-Q diagnostic plots
par(mfrow = c(1, 2))
plot(fitted_gp)
graphics.off()
# Fit exponential by passing a list with a fixed parameter
<- fit.gpd(windlyon,
reduced_gp threshold = u,
fpar = list(shape = 0))
# The MLE is sample mean of exceedances - check this
isTRUE(coef(reduced_gp) == mean(windlyon))
[1] TRUE
# Compare nested models using likelihood ratio test
anova(fitted_gp, reduced_gp)
Analysis of Deviance Table
npar Deviance Df Chisq Pr(>Chisq)
fitted_gp 2 415.06
reduced_gp 1 502.50 1 87.448 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The mev
package includes alternative routines for estimation, including the optimal bias-robust estimator of Dupuis (1999) and the approximate Bayesian estimators of Zhang & Stephens (2009) and Zhang (2010). The latter two are obtained by running a Markov chain Monte Carlo algorithm, but only the posterior mean and standard deviation are returned to reduce the memory footprint of the returned object, and these are calculated on the fly using running mean and variance estimators.
# Bayesian point estimates (based on MAP)
fit.gpd(windlyon,
threshold = u,
show = TRUE,
MCMC = TRUE,
method = "zhang")
Method: Zhang
Threshold: 33.84
Number Above: 90
Proportion Above: 0.0079
Approximate posterior mean estimates
scale shape
3.456 0.066
Posterior mean estimates
scale shape
3.5594 0.0495
Monte Carlo standard errors
scale shape
0.453 0.130
Estimates based on an adaptive MCMC
Runs: 10000
Burnin: 2000
Acceptance rate: 0.43
Thinning: 1
If the sample is small, maximum likelihood estimators are biased for the generalized Pareto distribution (the shape parameter is negatively biased, regardless of the true value for subtract
), which solves the implicit equation
# First-order bias corrected estimates
<- gpd.bcor(par = coef(fitted_gp),
corr_coef dat = windlyon,
corr = "firth")
Error in bcor.st$value : $ operator is invalid for atomic vectors
Inhomogeneous Poisson process
Let
Consider a sample of
The model includes additional arguments, np
and npp
(number of observations per period). 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. Alternatively, one can specify instead the number of periods np
, akin to npp*np
theoretically equal to the number of exceedances.
The tuning parameters impact the convergence of the estimation since the dependence between parameters becomes very strong: Sharkey & Tawn (2017) suggest to pick a value of np
that near-orthogonalize the parameters. Wadsworth:2011 recommended picking this to be the number of observations (so npp=1
), but Moins et al. (2023) show that a better choice leads to orthogonalization.
Another option is to fit the generalized Pareto model: if the probability of exceeding threshold np
. With the point estimates of the generalized Pareto model, say
The log likelihood of the
We can simulate from the
Risk measures
Two typical questions in extreme values are: given the intensity of an extreme event, what is its recurrence period? and what is a typical worst-case scenario over a given period of time? For the latter, suppose for simplicity that the daily observations are blocked into years, so that inference is based on
Quantiles, mean and return levels of