Generalized Pareto Distribution

This is a location-scale family distribution.

Parameters¶

Interpretation¶

We can interpret the shape parameters as follows:

$\kappa=0$. This corresponds to a type 1, short tail distribution with exponential decay.

$\kappa>0$. This corresponds to a type 2, heavy tail distribution with a slow power-law decay.

$\kappa<0$. This corresponds to a type 3, thin-tailed, polynomial decay with a finite upper bound.

Probability Density Function¶

This is denoted as the probability that our rv $Y$ will be equivalent to some specific value, $y$, conditioned on the fact that our values are greater than some threshold $y_0$.

We can define the probability density function as

where $a_+=\text{max}(a,0)$.

Cumulative Distribution Function¶

This is denoted as the probability that our rv $Y$ will be less than or equal to some specific value $y$ conditioned on the fact that our values are greater than some threshold $y_0$.

We can define the cumulative density function is defined as:

Survival Function¶

This is the exceedence probability of $Y$ above some value $y$. However, in this case, this probability is conditioned on the probability of a threshold value, $y_0$.

We denote as the survival function of the GPD. This is simply 1 minus the CDF of the GPD given as:

We can plug in the CDF function into this equation defined as:

Quantile Function¶

This is also known as the Point-Percentile-Function or the inverse CDF. This function maps an input threshold, $y_0$, to a value $y$ st the probability of $Y$ being less than or equal to $y$ is $y_p$.

We can take the inverse of this function to see that it is the inverse CDF which we denote as the quantile function.

where $y_p\in[0,1]$ is the data within the probability transform domain. These can be computed in closed form

# function for kappa neq 0
def gpd_quantile(yp, loc, scale, shape):
    # loc - scale / shape * (1 - (- log(1 - p)) ** (- shape))
    return loc + scale / shape * ((1 - yp) ** (- shape) - 1)
# function for kappa = 0
def quantile(yp, loc, scale):
    return loc - scale * log(1 - yp)

# location scale inverse
def location_scale_inverse(z, location, scale):
    return location + scale * z

def location_scale_forward(x, location, scale):
    return (x - location) / scale

def gpd_quantile_shape(y_p, shape):
    # get boolian array
    is_shape_zero = jnp.equal(shape, 0.0)
    # calculate the true shape
    safe_shape = jnp.where(is_shape_zero, 0.0, shape)
    # calculate negative log term, -log(1-yp)
    shape_zero_term = - jnp.log1p(- y_p)
    # calculate k neq 0, exp(-k log(1-yp)) / k
    shape_nonzero_term = jnp.expm1(shape * neglog1mp) / safe_shape
    # zero 
    term = jnp.where(is_shape_zero, shape_zero_term,  shape_nonzero_term)
    return term

# function for kappa > 0
def gpd_quantile(y_p, loc, scale, shape):
    # calculate the shape parameter
    z = gpd_quantile_shape(y_p, shape)
    # calculate the inverse of location-scale
    return location_scale_inverse(z, loc, scale)

Inverse Survival Function¶

The inverse survival function maps the input value $y_s$ to some probability value $y_p$ which represents the probability that there is some value greater than said input value.

So recall the survival function is the negative of the CDF. This implies that the inverse survival function is the negative of the quantile function.

We can take the inverse of this function to see that it is the inverse CDF which we denote as the quantile function.

where $y_s\in[0,1]$ is the survival probability within the probability transform domain. These can be computed in closed form

where $y_s=1-y_p$ is the survival probability.

# function for kappa > 0
def gpd_isf(y_s, loc, scale, shape):
    # calculate the shape parameter
    z = gpd_quantile_shape(1 - y_s, shape)
    # calculate the inverse of location-scale
    return location_scale_inverse(z, loc, scale)

Joint Distribution¶

We can write the likelihood that the observations, $y$, follow the GEVD distribution. So, given some observations, $\mathcal{D}=\{y_n\}_{n=1}^{N}$, which we believe follow the GEVD distribution, we can write the joint distribution decomposition as

This implies that the global prior parameters come from some distribution

and that these parameters get passed through our data likelihood term

Log Probability¶

Let’s say we are given some samples.

where $N$ are the number of exceedances above our threshold, $y_0$. Recall the GPD PDF for our iid samples is

We can add the log term to get

which reduces to

where $z_n=(y_n - y_0)/\sigma$ and $[1 + \kappa z_n]_+ = \text{max}(1 + \kappa z_n,0)$.

def gpd_logpdf(x, location, scale, shape):
    # calculate location scale: z=(y−μ​)/σ
    z = (x - location) / scale
    # calculate t(z) = max(1+κz)
    t = max(1.0 + shape * z, 0)
    # term 1: −log σ
    t1 = - np.log(scale)
    # term 2: - (1+1/κ)log(1+κz)
    t2 = - (1.0 / shape + 1.0) * np.log(t) 
    return  t1 + t2

y: Array["T"] = ...
params: PyTree = ...
# apply vectorized operation
nll: Array["T"] = vectorize(gpd_logpdf, y, *params)
# take the sume
nll: Scalar = sum(nll)

Return Period¶

We can calculate the RP using equation (8).

We assume that the occurrence of events over the threshold, $y_0$, is given by the Poisson process.

We can calculate the RP using equation (8). Practically, we set this to the survival function of the GEVD (equation ).

If we rearrange this equation, we get

To make things simpler, we can simply use the quantile function in equation (11) and set the probability to

However, if we expand this out, we get

Average Recurrence Interval¶

We can calculate the RP using equation (14).

Removing the exponential components, we are left with

We assume that the occurrence of events over the threshold, $y_0$, is given by the Poisson process.

We can calculate the RP using equation (14). Practically, we set this to the survival function of the GEVD (equation ).

If we rearrange this equation, we get

To make things simpler, we can simply use the quantile function in equation (11) and set the probability to

However, if we expand this out, we get

Reparameterization¶

In this instance, we are assuming that there is a threshold parameter, $y_0$. We can write the reparameterization of the GPD in terms of the GEV distribution

We can also re-parameterize the $\lambda_{y_0}$ in terms of the GPD and GEV parameters.

Lastly, we can do something similar like so

There has also been a similar reparameterization found in

Similarly, another reparameterization is

Background¶

where $y>0$. If we let $y_0\rightarrow\infty$, then this leads to an approximate family of distributions given by the GPD

Marginal Survival Function¶

We are interested in the marginal probability of occurrence above an arbitrary maximum value, $y$. We can write the joint distribution of both quantities to be factored as follows.

The first term is the rate of exceedences above some quantity $y$ given some threshold, $y_0$. The second term is the probability that an event is above some threshold, $y_0$. We could also describe it as the rate of exceedences above some threshold, $y_0$.

Let’s define an arrival rate λ to be the average number of events per year larger than a threshold, $y_0$. This is analogous to a Poisson distribution.

where $k$ is the number of occurrences within some period $T$. We can write down this distribution as

We know that the expected value is simply the parameter λ.

We can calculate this approximately by summing the number of events, over the threshold, $y_0$, and then we divide them by the total number of events, $N_{y_0}$

To relate this back to our function, we would need the rate parameter λ in units as events per year

where $t$ is some conversion factor from whatever time unit to years.

Literature Review¶

There are many cases where the GPD is used within the literature. We split this section into theory and applications. However, most of the theoretic work that presents this method is based in applied settings, in particular hydrology.

Theory. From a more theoretic perspective, Davison & Smith, 1990Martins & Stedinger, 2001Chavez-Demoulin * et al., 2005 introduce the Poisson-GPD method where they also relate the AEP and ARI. In Coles, 2001 (Chapter 4.3.3 - Appendix A.1), there is a staple chapter where they introduce the Poisson-GPD method. In Wang & Holmes, 2020, the authors discuss some of the key differences between the annual exceedence probability (AEP) and the average recurrence interval (ARI). Birkhäuser Basel (2007) is a another good book that has a few chapters dedicated to BM, POTs and PPs.

Applications. In Nemukula & Sigauke, 2020, they apply this to study maximum daily temperature in South Africa where they apply a non-stationary Poisson-GPD. Thiombiano et al. (2016) study how the Arctic Oscillation and Pacific North American covariates are related to the extreme daily precipitation months in Southeastern Canada where they apply a Poisson-GPD method with spline functions to map the covariates to the GPD parameters. Silva et al. (2015) study how El Niño-Southern Oscillation can effect the flood regime in the Itajaí river basin in Southern Brasil where they apply a non-stationary Poisson-GPD. Cid et al., 2015 look at how the North Atlantic Oscillation effect storm surges in the Atlantic and North Atlantic regions where they apply a non-stationary Poisson-GPD. Katz et al. (2005) investigate ecological disturbances in extremes for paleoecology where they apply a non-stationary Poisson-GPD.

Algorithms. Silva et al. (2015) showcase how we can use the delta-method (aka Laplace approximation) to calculate uncertainty intervals on the parameters. MacDonald et al. (2011) consider a Kernel Density estimator for the event occurrences parameterization.