Block Maxima Model

The stationary BlockMaxima model is illustrated using the annual maximum sea-levels recorded at Port Pirie in South Australia from 1923 to 1987, studied by Coles (2001) in Chapter 3. The annual maxima are assumed independent and identically distributed.

The Extremes.jl package supports maximum likelihood inference, Bayesian inference and inference based on the probability weigthed moments. For the GEV parameter estimation, the following functions can be used:

gevfitpwm: estimation with the probability weighted moments;
gevfit: maximum likelihood estimation;
gevfitbayes: Bayesian estimation.

Log-scale paremeter

These functions shows the estimate of the log-scale parameter $\phi = \log \sigma$ instead of the scale parameter.

Load the data

Loading the annual maximum sea-levels at Port Pirie:

data = Extremes.dataset("portpirie")
first(data,5)

5×2 DataFrame

Row	Year	SeaLevel
	Int64	Float64
1	1923	4.03
2	1924	3.83
3	1925	3.65
4	1926	3.88
5	1927	4.01

The annual maxima can be shown as a function of the year using the Gadfly package:

set_default_plot_size(12cm, 8cm)
plot(data, x=:Year, y=:SeaLevel, Geom.line)

Maximum likelihood inference

GEV parameter estimation

The GEV parameter estimation with maximum likelihood is performed with the gevfit function:

julia> fm = gevfit(data, :SeaLevel)MaximumLikelihoodAbstractExtremeValueModel
model :
	BlockMaxima{GeneralizedExtremeValue}
	data :		Vector{Float64}[65]
	location :	μ ~ 1
	logscale :	ϕ ~ 1
	shape :		ξ ~ 1

θ̂  :	[3.874750223091266, -1.6192723640210762, -0.05010719929448139]

Note

In this example, the gevfit function is called using the data DataFrame structure as the first argument. The function can also be called using the vector of maxima as the first argument, e.g. gevfit(data[:,:SeaLevel]).

The vector of the parameter estimates (location scale and shape) can be extracted with the function params:

julia> params(fm)1-element Vector{Tuple{Float64, Float64, Float64}}:
 (3.874750223091266, 0.19804274969909974, -0.05010719929448139)

The location parameter with the function location:

julia> location(fm)1-element Vector{Float64}:
 3.874750223091266

The scale parameter with the function Extremes.scale:

julia> scale(fm)1-element Vector{Float64}:
 0.19804274969909974

The shape parameter with the function shape:

julia> shape(fm)1-element Vector{Float64}:
 -0.05010719929448139

Type-stable function

These functions return a unit dimension vector for the return level and a vector containing only one vector for the confidence interval. The reason is that the functions always return the same type in the stationary and non-stationary case. The functions are therefore type-stable allowing better performance of code execution.

The approximate covariance matrix of the parameter estimates can be obtained with the parametervar function:

julia> parametervar(fm)3×3 Matrix{Float64}:
  0.000780204   0.000995016  -0.0010741
  0.000995016   0.0104541    -0.00392576
 -0.0010741    -0.00392576    0.00965404

Confidence intervals for the parameters are obtained with the cint function:

julia> cint(fm)3-element Vector{Vector{Float64}}:
 [3.820004234825991, 3.929496211356541]
 [-1.819669858589598, -1.4188748694525544]
 [-0.24268345866324303, 0.14246906007428023]

For instance, the shape parameter 95% confidence interval is as follows:

julia> cint(fm)[3]2-element Vector{Float64}:
 -0.24268345866324303
  0.14246906007428023

Diagnostic plots

Several diagnostic plots for assessing the accuracy of the GEV model fitted to the Port Pirie data are can be shown with the diagnosticplots function:

set_default_plot_size(21cm ,16cm)
diagnosticplots(fm)

The diagnostic plots consist in the probability plot (upper left panel), the quantile plot (upper right panel), the density plot (lower left panel) and the return level plot (lower right panel). These plots can be displayed separately using respectively the probplot, qqplot, histplot and returnlevelplot functions.

Return level estimation

T-year return level estimate can be obtained using the returnlevel function. For example, the 100-year return level for the Port Pirie block maxima model is computed as follows:

julia> r = returnlevel(fm, 100)ReturnLevel
returnperiod :	100
value :		Vector{Float64}[1]

The return level can be accessed as follows:

julia> r.value1-element Vector{Float64}:
 4.688403360432851

The corresponding confidence interval can be computed with the cint function:

julia> c = cint(r)1-element Vector{Vector{Real}}:
 [4.377121171613512, 4.99968554925219]

To get the scalar return level in the stationary case, the following command can be used:

julia> r.value[]4.688403360432851

To get the scalar confidence interval in the stationary case, the following command can be used:

julia> c[]2-element Vector{Real}:
 4.377121171613512
 4.99968554925219

Bayesian Inference

Most functions described in the previous sections also work in the Bayesian context. To reproduce exactly the results, the seed should be fixed as follows:

import Random
Random.seed!(4786)

GEV parameter estimation

The Bayesian GEV parameter estimation is performed with the gevfitbayes function:

julia> fm = gevfitbayes(data, :SeaLevel)
Progress:   0%|                                         |  ETA: 0:06:38
Progress: 100%|█████████████████████████████████████████| Time: 0:00:02
BayesianAbstractExtremeValueModel
model :
	BlockMaxima{GeneralizedExtremeValue}
	data :		Vector{Float64}[65]
	location :	μ ~ 1
	logscale :	ϕ ~ 1
	shape :		ξ ~ 1

sim :
	MambaLite.Chains
	Iterations :		2001:5000
	Thinning interval :	1
	Chains :		1
	Samples per chain :	3000
	Value :			Array{Float64, 3}[3000,3,1]

Prior

Currently, only the improper uniform prior is implemented, i.e. \[ f_{(μ,ϕ,ξ)}(μ,ϕ,ξ) ∝ 1. \] It yields to a proper posterior as long as the sample size is larger than 3 (Northrop and Attalides, 2016).

Sampling scheme

Currently, the No-U-Turn Sampler extension (Hoffman and Gelman, 2014) to Hamiltonian Monte Carlo (Neel, 2011, Chapter 5) is implemented for simulating an autocorrelated sample from the posterior distribution.

The generated sample from the posterior distribution is contained in the field sim of the fitted structure. It is an object of type Chains from the Mamba.jl package.

Credible intervals for the parameters are obtained with the cint function:

julia> cint(fm)3-element Vector{Vector{Float64}}:
 [3.8182738930811233, 3.931583774977113]
 [-1.7897113975890682, -1.4005833922456459]
 [-0.2018902582086667, 0.1553392907141933]

Inference based on the probability weighted moments

Most functions described in the previous sections also work for the model fitted with the probability weighted moments.

GEV parameter estimation

The parameter estimation based on the probability weighted moments is performed with the gevfitpwm function:

julia> fm = gevfitpwm(data, :SeaLevel)pwmAbstractExtremeValueModel
model :
	BlockMaxima{GeneralizedExtremeValue}
	data :		Vector{Float64}[65]
	location :	μ ~ 1
	logscale :	ϕ ~ 1
	shape :		ξ ~ 1

θ̂  :	[3.8731723562720766, -1.5932320395836068, -0.051477125862911276]

The approximate covariance matrix of the parameter estimates using a bootstrap procedure can be obtained with the parametervar function:

julia> parametervar(fm)3×3 Matrix{Float64}:
  0.000871239   0.00116247  -0.00111277
  0.00116247    0.0106812   -0.00361624
 -0.00111277   -0.00361624   0.0072629

Confidence intervals on the parameter estimates using a bootstrap procedure can be obtained with the cint function:

julia> cint(fm)3-element Vector{Vector{Float64}}:
 [3.822276679146202, 3.9317746782451266]
 [-1.8146460804733326, -1.4222820715987998]
 [-0.224777189866786, 0.08648421062397504]

Inference for the Gumbel distribution

The inference for the Gumbel distribution is also provided for modeling the series of block maxima. The Gumbel distribution is a sub-family of the GEV distribution when the shape parameter $\xi = 0$. The library provides functions analogous to those for the GEV distribution for the adjustment of the parameters of the Gumbel distribution:

by maximum likelihood gumbelfit;
under the Bayesian paradigm gumbelfitbayes;
by the method of moments gumbelfitpwm.

The inference for the fitted model is performed with the same methods defined previously for the GEV distribution, notably with the cint and returnlevel methods. Diagnostic plot methods diagnosticplots are also implemented, such as probplot, qqplot, histplot and returnlevelplot.

Gumbel distribution for modelling block maxima

As extreme-value statisticians, we advocate avoiding the use of the Gumbel distribution for modelling the block maxima. This is because the choice of family is made with the data at hand, and when extrapolating to large quantiles, i.e. larger than the range of the data, the uncertainty associated with this choice is not taken into account. If the data suggest that the Gumbel family is the best one, this does not imply that the other families are not plausible. In applications, the confidence intervals on the shape parameter are often wide, representing the difficulty of discriminating the tail behavior using only the limited number of data. Therefore, we plead for the use of the GEV distribution for the block maxima model without choosing a subfamily such as the Gumbel. As Coles (2001) also argued in Page 64, this is "... the safest option is to accept there is uncertainty about the value of the shape parameter ... and to prefer the inference based on the GEV model. The larger measures of uncertainty generated by the GEV model then provide a more realistic quantification of the genuine uncertainties involved in model extrapolation".

Example on the annual maximum sea-levels recorded at Port Pirie

The Gumbel distribution parameter estimation with maximum likelihood can be performed with the gumbelfit function:

julia> fm = gumbelfit(data, :SeaLevel)MaximumLikelihoodAbstractExtremeValueModel
model :
	BlockMaxima{Gumbel}
	data :		Vector{Float64}[65]
	location :	μ ~ 1
	logscale :	ϕ ~ 1

θ̂  :	[3.8694436499474767, -1.6353194298573894]

The approximate covariance matrix of the parameter estimates can be obtained with the parametervar function:

julia> parametervar(fm)2×2 Matrix{Float64}:
 0.000649941  0.000783584
 0.000783584  0.00935991

Confidence intervals for the parameters are obtained with the cint function:

julia> cint(fm)2-element Vector{Vector{Float64}}:
 [3.819476453475028, 3.9194108464199253]
 [-1.824939298384268, -1.4456995613305108]

The diagnostic plots for assessing the accuracy of the Gumbel model can be shown with the diagnosticplots function:

set_default_plot_size(21cm ,16cm)
diagnosticplots(fm)

The 100-year return level for the Port Pirie block maxima model is computed with the returnlevel function as follows:

julia> r = returnlevel(fm, 100)ReturnLevel
returnperiod :	100
value :		Vector{Float64}[1]

The return level can be accessed as follows:

julia> r.value1-element Vector{Float64}:
 4.7659672278955645

The corresponding confidence interval can be computed with the cint function:

julia> c = cint(r)1-element Vector{Vector{Real}}:
 [4.574150814798266, 4.957783640992863]