Non-Stationary Block Maxima Model

In the non-stationary block maxima model, the GEV parameters are allowed to be functions of multiple covariates:

\[ μ = X₁ × β₁ \\ ϕ = X₂ × β₂ \\ ξ = X₃ × β₃ \]

where $(X₁,β₁)$, $(X₂,β₂)$ and $(X₃,β₃)$ are respectively the design matrix and the corresponding coefficient parameter vector of $μ$, $ϕ$ and $ξ$.

Intercept

An intercept is included in all parameter functions by default.

The non-stationary BlockMaxima model is illustrated using the annual maximum sea-levels recorded at Fremantle in West Australia from 1897 to 1989, studied by Coles (2001) in Chapter 6.

Load the data

Loading the annual maximum sea-levels at Fremantle:

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

5×3 DataFrame

Row	Year	SeaLevel	SOI
	Int64	Float64	Float64
1	1897	1.58	-0.67
2	1898	1.71	0.57
3	1899	1.4	0.16
4	1900	1.34	-0.65
5	1901	1.43	0.06

The annual maxima can be plotted as function of the year:

plot(data, x=:Year, y=:SeaLevel, Geom.line,
    Coord.cartesian(xmin=1895, xmax=1990), Guide.xticks(ticks=1895:10:1990))

and as function of the Southern Oscillation Index:

plot(data, x=:SOI, y=:SeaLevel, Geom.point)

Both variables can be included in the block maxima model. Parameter estimation can be performed either by maximum likelihood or by the Bayesian approach. Probability weighted moment estimation cannot be used in the non-stationary case.

Maximum likelihood inference

GEV parameter estimation

The GEV parameter estimation with maximum likelihood is performed with the gevfit function. The parameter estimate vector $\mathbf{θ̂} = (\mathbf{β̂₁},\, \mathbf{β̂₂},\, \mathbf{β̂₃})^\top$ is contained in the field θ̂ of the returned structure.

Several non-stationary model can be fitted.

The stationary model

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

θ̂  :	[1.4823411063972953, -1.957072887171364, -0.21742418585608772]

The location parameter varying as a linear function of the year

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

θ̂  :	[-2.4727880606816433, 0.00203216246147196, -2.084849395454593, -0.12530315168965786]

The location parameter varying as a linear function of the year and the SOI

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

θ̂  :	[-2.62599181296525, 0.002114032686702418, 0.054517690075681226, -2.114182020887254, -0.14998552392836614]

Both the location and logscale parameters varying as a linear function of the year and the SOI

julia> fm₃ = gevfit(data, :SeaLevel, locationcovid = [:Year, :SOI], logscalecovid = [:Year, :SOI])MaximumLikelihoodAbstractExtremeValueModel
model :
	BlockMaxima{GeneralizedExtremeValue}
	data :		Vector{Float64}[86]
	location :	μ ~ 1 + Year + SOI
	logscale :	ϕ ~ 1 + Year + SOI
	shape :		ξ ~ 1

θ̂  :	[-2.215213912601999, 0.0019067329692207917, 0.06592556257024981, 6.760912614047801, -0.004564874869909601, 0.26512520665238254, -0.22292401254869737]

As show by Coles (2001), the best model is the one where the location parameter varies as a linear function of the year and the SOI, the fm₂ in the present section.

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

julia> params(fm₂)86-element Vector{Tuple{Float64, Float64, Float64}}:
 (1.3478013413585301, 0.12073200543257417, -0.14998552392836614)
 (1.4175173097390779, 0.12073200543257417, -0.14998552392836614)
 (1.3972790894947504, 0.12073200543257417, -0.14998552392836614)
 (1.3552337932201515, 0.12073200543257417, -0.14998552392836614)
 (1.396055385860587, 0.12073200543257417, -0.14998552392836614)
 (1.4226357041650213, 0.12073200543257417, -0.14998552392836614)
 (1.4203883216456696, 0.12073200543257417, -0.14998552392836614)
 (1.3041989668681433, 0.12073200543257417, -0.14998552392836614)
 (1.414258025904695, 0.12073200543257417, -0.14998552392836614)
 (1.4228475064841544, 0.12073200543257417, -0.14998552392836614)
 ⋮
 (1.5651780007967813, 0.12073200543257417, -0.14998552392836614)
 (1.4855154983699614, 0.12073200543257417, -0.14998552392836614)
 (1.5148883760945049, 0.12073200543257417, -0.14998552392836614)
 (1.560616560841752, 0.12073200543257417, -0.14998552392836614)
 (1.5665468318337523, 0.12073200543257417, -0.14998552392836614)
 (1.555031442001534, 0.12073200543257417, -0.14998552392836614)
 (1.4944501311012035, 0.12073200543257417, -0.14998552392836614)
 (1.6165030819544042, 0.12073200543257417, -0.14998552392836614)
 (1.6120749918320252, 0.12073200543257417, -0.14998552392836614)

The location parameter with the function location:

julia> location(fm₂)86-element Vector{Float64}:
 1.3478013413585301
 1.4175173097390779
 1.3972790894947504
 1.3552337932201515
 1.396055385860587
 1.4226357041650213
 1.4203883216456696
 1.3041989668681433
 1.414258025904695
 1.4228475064841544
 ⋮
 1.5651780007967813
 1.4855154983699614
 1.5148883760945049
 1.560616560841752
 1.5665468318337523
 1.555031442001534
 1.4944501311012035
 1.6165030819544042
 1.6120749918320252

The scale parameter with the function Extremes.scale:

julia> scale(fm₂)86-element Vector{Float64}:
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 ⋮
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417
 0.12073200543257417

The shape parameter with the function shape:

julia> shape(fm₂)86-element Vector{Float64}:
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
  ⋮
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614
 -0.14998552392836614

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

julia> parametervar(fm₂)5×5 Matrix{Float64}:
  1.02181      -0.000524501   6.71532e-5    0.0127779    -0.0209508
 -0.000524501   2.69285e-7   -1.61708e-8   -6.4457e-6     1.05896e-5
  6.71532e-5   -1.61708e-8    0.00038545    0.000239274  -0.000355337
  0.0127779    -6.4457e-6     0.000239274   0.00703707   -0.00256783
 -0.0209508     1.05896e-5   -0.000355337  -0.00256783    0.00444377

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

julia> cint(fm₂)5-element Vector{Vector{Float64}}:
 [-4.607209165945568, -0.6447744599849323]
 [0.0010969546332637297, 0.0031311107401411067]
 [0.01603794386329232, 0.09299743628807014]
 [-2.2785980179345966, -1.9497660238399117]
 [-0.2806398353854793, -0.019331212471252968]

In particular, the 95% confidence interval for the rise in annual maximum sea-levels per year is as follows:

julia> cint(fm₂)[2]2-element Vector{Float64}:
 0.0010969546332637297
 0.0031311107401411067

Diagnostic plots

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

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

The diagnostic plots consist in the residual probability plot (upper left panel), the residual quantile plot (upper right panel) and the residual density plot (lower left panel) of the standardized data (see Chapter 6 of Coles, 2001). These plots can be displayed separately using respectively the probplot, qqplot, histplot and returnlevelplot functions.

Return level estimation

Since the model parameters vary in time, the quantiles also vary in time. Therefore, a T-year return level can be estimated for each year. This set of return levels are referred to as effective return levels as proposed by Katz et al. (2002).

The 100-year effective return levels for the fm₂ model can be computed using the returnlevel function:

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

The effective return levels can be accessed as follows:

julia> r.value86-element Vector{Float64}:
 1.7489936832950526
 1.8187096516756003
 1.7984714314312729
 1.756426135156674
 1.7972477277971095
 1.8238280461015437
 1.821580663582192
 1.7053913088046657
 1.8154503678412175
 1.8240398484206768
 ⋮
 1.9663703427333037
 1.8867078403064839
 1.9160807180310273
 1.9618089027782744
 1.9677391737702747
 1.9562237839380565
 1.895642473037726
 2.0176954238909266
 2.013267333768548

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

julia> c = cint(r)86-element Vector{Vector{Real}}:
 [1.6464829066065283, 1.8515044599835768]
 [1.7251036134260933, 1.9123156899251073]
 [1.7048688798346274, 1.8920739830279183]
 [1.6550379488049871, 1.8578143215083607]
 [1.7037278450906306, 1.8907676105035884]
 [1.7317757216930512, 1.9158803705100362]
 [1.7297984243804045, 1.9133629027839796]
 [1.5809572007674173, 1.8298254168419141]
 [1.723770734795733, 1.907130000886702]
 [1.7329937388835963, 1.9150859579577573]
 ⋮
 [1.859790953836958, 2.0729497316296497]
 [1.7587332274998884, 2.0146824531130796]
 [1.797245351179525, 2.03491608488253]
 [1.8523402194326652, 2.0712775861238835]
 [1.858201200934169, 2.0772771466063804]
 [1.8440901683740472, 2.068357399502066]
 [1.7646548289729513, 2.0266301171025005]
 [1.9073298034620776, 2.1280610443197756]
 [1.9027214984142238, 2.123813169122872]

The effective return levels along with their confidence intervals can be plotted as follows:

rmin = [c[i][1] for i in eachindex(c)]
rmax = [c[i][2] for i in eachindex(c)]
df = DataFrame(Year = data[:,:Year], r = r.value, rmin = rmin, rmax = rmax)

set_default_plot_size(12cm, 9cm)
plot(df, x=:Year, y=:r, ymin=:rmin, ymax=rmax, Geom.line, Geom.ribbon,
    Coord.cartesian(xmin=1895, xmax=1990), Guide.xticks(ticks=1895:10:1990),
    Guide.ylabel("100-year Effective Return Level"))

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)

Most functions described in the previous sections also work in the Bayesian context.

GEV parameter estimation

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

julia> fm = gevfitbayes(data, :SeaLevel, locationcovid = [:Year, :SOI])
Progress:   6%|██▍                                      |  ETA: 0:00:02
Progress:  11%|████▌                                    |  ETA: 0:00:02
Progress:  19%|███████▉                                 |  ETA: 0:00:01
Progress:  27%|███████████                              |  ETA: 0:00:01
Progress:  35%|██████████████▎                          |  ETA: 0:00:01
Progress:  42%|█████████████████▍                       |  ETA: 0:00:01
Progress:  50%|████████████████████▍                    |  ETA: 0:00:01
Progress:  58%|███████████████████████▋                 |  ETA: 0:00:01
Progress:  65%|██████████████████████████▉              |  ETA: 0:00:00
Progress:  73%|██████████████████████████████           |  ETA: 0:00:00
Progress:  80%|█████████████████████████████████        |  ETA: 0:00:00
Progress:  88%|████████████████████████████████████▎    |  ETA: 0:00:00
Progress:  96%|███████████████████████████████████████▌ |  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:01
BayesianAbstractExtremeValueModel
model :
	BlockMaxima{GeneralizedExtremeValue}
	data :		Vector{Float64}[86]
	location :	μ ~ 1 + Year + SOI
	logscale :	ϕ ~ 1
	shape :		ξ ~ 1

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

Prior

Currently, only the improper uniform prior is implemented, i.e. \[ f_{(β₁,β₂,β₃)}(β₁,β₂,β₃) ∝ 1. \]

Inference for the non-stationary Gumbel distribution

The package aslo provides functions for the inference of the non-stationary Gumbel model. Documentation on the Gumbel model can be found here in the Block Maxima section.

Example on the annual maximum sea-levels recorded at Fremantle

The location parameter varying as a linear function of the year and the SOI

julia> fm₂ = gumbelfit(data, :SeaLevel, locationcovid = [:Year, :SOI])MaximumLikelihoodAbstractExtremeValueModel
model :
	BlockMaxima{Gumbel}
	data :		Vector{Float64}[86]
	location :	μ ~ 1 + Year + SOI
	logscale :	ϕ ~ 1

θ̂  :	[-3.112447859214716, 0.002358975977102567, 0.045892020491916104, -2.147213954852745]

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

julia> cint(fm₂)4-element Vector{Vector{Float64}}:
 [-4.825555304377321, -1.3993404140521108]
 [0.0014786231231302162, 0.0032393288310749177]
 [0.011612282486563964, 0.08017175849726824]
 [-2.306675317663155, -1.987752592042335]

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

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

The 100-year effective return levels for the fm₂ model can be computed using the returnlevel function:

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

The effective return levels can be accessed as follows:

julia> r.value86-element Vector{Float64}:
 1.8691213924827486
 1.9283864738698273
 1.9119297214452446
 1.8771161608238944
 1.9120584713502575
 1.9355921517061487
 1.9342797660438973
 1.8370530575535422
 1.930278234104639
 1.9386675476981967
 ⋮
 2.1007765495184634
 2.0342974947576913
 2.059602480980752
 2.0986750733513877
 2.1042464907629244
 2.0951324616170472
 2.044715614028447
 2.148037035087765
 2.1448889686058377

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

julia> c = cint(r)86-element Vector{Vector{Real}}:
 [1.7583457823872852, 1.9798970025782119]
 [1.8180216271016414, 2.038751320638013]
 [1.8037675690928139, 2.0200918737976754]
 [1.7676627712661332, 1.9865695503816556]
 [1.8048195914347755, 2.0192973512657395]
 [1.8275998081282645, 2.043584495284033]
 [1.827056655321212, 2.0415028767665824]
 [1.7161554213121801, 1.9579506937949043]
 [1.8244859886440883, 2.03607047956519]
 [1.833241669248867, 2.0440934261475263]
 ⋮
 [2.0000054363435025, 2.2015476626934243]
 [1.9265267632441145, 2.142068226271268]
 [1.9566630009071906, 2.1625419610543135]
 [1.9977907834878426, 2.1995593632149326]
 [2.002961830717341, 2.2055311508085076]
 [1.9939002689389977, 2.1963646542950968]
 [1.9356352812226434, 2.1537959468342502]
 [2.0402732588868613, 2.255800811288669]
 [2.0381263835124295, 2.251651553699246]

The effective return levels along with their confidence intervals can be plotted as follows:

rmin = [c[i][1] for i in eachindex(c)]
rmax = [c[i][2] for i in eachindex(c)]
df = DataFrame(Year = data[:,:Year], r = r.value, rmin = rmin, rmax = rmax)

set_default_plot_size(12cm, 9cm)
plot(df, x=:Year, y=:r, ymin=:rmin, ymax=rmax, Geom.line, Geom.ribbon,
    Coord.cartesian(xmin=1895, xmax=1990), Guide.xticks(ticks=1895:10:1990),
    Guide.ylabel("100-year Effective Return Level"))