Non-Stationary Threshold Exceedance Model

In the non-stationary threshold exceedance model, the GP parameters are allowed to be functions of multiple covariates:

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

where $(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 ThresholdExceedance model is illustrated using the daily rainfall accumulations at a location in south-west England from 1914 to 1962, studied by Coles (2001) in Chapter 6.

Load the data

Loading the daily rainfall accumulations:

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

5×2 DataFrame

Row	Date	Rainfall
	Date	Float64
1	1914-01-01	0.0
2	1914-01-02	2.3
3	1914-01-03	1.3
4	1914-01-04	6.9
5	1914-01-05	4.6

Extract the exceedances over the threshold of 30 mm:

threshold = 30.0
df = filter(row -> row.Rainfall > threshold, data)
df[!,:Exceedance] = df[:,:Rainfall] .- threshold
df[!,:Year] = year.(df[:,:Date])
plot(df, x=:Date, y=:Exceedance, Geom.point)

Non-stationary 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.

The GP parameter estimation with maximum likelihood is performed with the gpfit function. The parameter estimate vector $\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₀ = gpfit(df, :Exceedance)MaximumLikelihoodAbstractExtremeValueModel
model :
	ThresholdExceedance
	data :		Vector{Float64}[152]
	logscale :	ϕ ~ 1
	shape :		ξ ~ 1

θ̂  :	[2.006896498380506, 0.1844926991237574]

The logscale parameter varying as a linear function of the year

julia> fm₁ = gpfit(df, :Exceedance, logscalecovid = [:Year])MaximumLikelihoodAbstractExtremeValueModel
model :
	ThresholdExceedance
	data :		Vector{Float64}[152]
	logscale :	ϕ ~ 1 + Year
	shape :		ξ ~ 1

θ̂  :	[-11.205671606002209, 0.006804339309339082, 0.1976859276224044]

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

julia> cint(fm₁)3-element Vector{Vector{Float64}}:
 [-38.83131900584391, 16.419975793839495]
 [-0.007421606369680473, 0.021030284988358634]
 [-0.0018240741601000532, 0.39719592940490883]

In particular, the 95% confidence interval for the rise in the log-scale parameter per year is as follows:

julia> cint(fm₁)[2]2-element Vector{Float64}:
 -0.007421606369680473
  0.021030284988358634

Diagnostic plots

Several diagnostic plots for assessing the accuracy of the GP model fitted to the rainfall 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)^[1].

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

julia> nobs = size(data,1)17531
julia> nobsperblock = 365365
julia> r = returnlevel(fm₁, threshold, nobs, nobsperblock, 100)ReturnLevel
returnperiod :	100
value :		Vector{Float64}[152]

The effective return levels can be accessed as follows:

julia> r.value152-element Vector{Float64}:
  96.07236738849316
  96.07236738849316
  96.07236738849316
  96.07236738849316
  96.523479213636
  96.523479213636
  96.523479213636
  96.97767102345072
  96.97767102345072
  96.97767102345072
   ⋮
 119.74262810056645
 119.74262810056645
 119.74262810056645
 119.74262810056645
 119.74262810056645
 120.35534961117199
 120.35534961117199
 120.97225450327086
 120.97225450327086

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

julia> c = cint(r)152-element Vector{Vector{Real}}:
 [55.894149909256036, 136.25058486773028]
 [55.894149909256036, 136.25058486773028]
 [55.894149909256036, 136.25058486773028]
 [55.894149909256036, 136.25058486773028]
 [56.49568315762775, 136.55127526964424]
 [56.49568315762775, 136.55127526964424]
 [56.49568315762775, 136.55127526964424]
 [57.086045158813825, 136.86929688808763]
 [57.086045158813825, 136.86929688808763]
 [57.086045158813825, 136.86929688808763]
 ⋮
 [61.3830404173303, 178.1022157838026]
 [61.3830404173303, 178.1022157838026]
 [61.3830404173303, 178.1022157838026]
 [61.3830404173303, 178.1022157838026]
 [61.3830404173303, 178.1022157838026]
 [60.873037890298875, 179.83766133204512]
 [60.873037890298875, 179.83766133204512]
 [60.335753531546054, 181.60875547499566]
 [60.335753531546054, 181.60875547499566]

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_plot = DataFrame(Year = year.(df[:,:Date]), r = r.value, rmin = rmin, rmax = rmax)

set_default_plot_size(12cm, 8cm)
plot(df_plot, x=:Year, y=:r, ymin=:rmin, ymax=rmax, Geom.line, Geom.ribbon,
    Coord.cartesian(xmin=1910, xmax=1965), Guide.xticks(ticks=1910:5:1965),
    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)

GP parameter estimation

The Bayesian GP parameter estimation is performed with the gpfitbayes function:

julia> fm = gpfitbayes(df, :Exceedance, logscalecovid = [:Year])
Progress:  12%|█████                                    |  ETA: 0:00:01
Progress:  24%|█████████▉                               |  ETA: 0:00:01
Progress:  38%|███████████████▌                         |  ETA: 0:00:01
Progress:  52%|█████████████████████▌                   |  ETA: 0:00:00
Progress:  67%|███████████████████████████▌             |  ETA: 0:00:00
Progress:  82%|█████████████████████████████████▋       |  ETA: 0:00:00
Progress:  97%|███████████████████████████████████████▋ |  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:00
BayesianAbstractExtremeValueModel
model :
	ThresholdExceedance
	data :		Vector{Float64}[152]
	logscale :	ϕ ~ 1 + Year
	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. \]

The model fit can be assessed with the diagnostic plots using the function diagnosticplots:

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

The empirical covariance matrix of the parameters and the credible intervals can be obtained with the functions parametervar and cint, respectively:

julia> parametervar(fm)3×3 Matrix{Float64}:
 210.232     -0.108249     -0.20519
  -0.108249   5.57417e-5    0.000101175
  -0.20519    0.000101175   0.0115637

julia> cint(fm, .95)3-element Vector{Vector{Float64}}:
 [-39.25043820184902, 16.596351149913975]
 [-0.007719583007889886, 0.02100563222999514]
 [0.035250436765367404, 0.44066637811066683]

Return level estimation

The 100-year effective return level estimates can be obtained using the function returnlevel:

julia> nobs = size(data,1)17531
julia> nobsperblock = 365365
julia> r = returnlevel(fm, threshold, nobs, nobsperblock, 100)ReturnLevel
returnperiod :	100
value :		Matrix{Float64}[456000]

The corresponding 95% credible interval can be computed using the function cint:

julia> c = cint(r, 0.95)152-element Vector{Vector{Real}}:
 [63.90230140290233, 168.0372574616761]
 [63.90230140290233, 168.0372574616761]
 [63.90230140290233, 168.0372574616761]
 [63.90230140290233, 168.0372574616761]
 [64.11418140628982, 168.3911173704221]
 [64.11418140628982, 168.3911173704221]
 [64.11418140628982, 168.3911173704221]
 [64.48520658159299, 168.09173468161583]
 [64.48520658159299, 168.09173468161583]
 [64.48520658159299, 168.09173468161583]
 ⋮
 [74.00108362294426, 222.6767296225598]
 [74.00108362294426, 222.6767296225598]
 [74.00108362294426, 222.6767296225598]
 [74.00108362294426, 222.6767296225598]
 [74.00108362294426, 222.6767296225598]
 [73.08731068757521, 224.4069539397245]
 [73.08731068757521, 224.4069539397245]
 [71.70757332656396, 225.85208724974973]
 [71.70757332656396, 225.85208724974973]

The 100-year effective return levels along with their 95% credible intervals can be illustrated as follows:

df_plot = DataFrame(Year=df.Year, ReturnLevel = vec(mean(r.value, dims=1)),
    LowerBound = first.(c), UpperBound=last.(c))

set_default_plot_size(12cm, 8cm)
plot(df_plot, x=:Year, y=:ReturnLevel, Geom.line,
    ymin=:LowerBound, ymax=:UpperBound, Geom.ribbon,
    Guide.ylabel("100-year effective return level"))

1Katz, R. W., M. B. Parlange, and P. Naveau (2002), Statistics of extremes in hydrology, Adv. Water Resour., 25, 1287–1304.