Distributional regression models for Extended Generalized Pareto distributions

The EGPD (Naveau et al., 2016) accounts simultaneously and in a parsimonious way for both extreme and non extreme values of a random variable $Y$. Its cumulative distribution function (CDF) reads:

The function $G(\cdot)$ is a continuous CDF on the unit interval. It is constrained by several assumptions: 1) we need to ensure that the upper-tail behavior of $F$ follows a GPD; 2) the low values of $Y$, seen as the upper tail of $-Y$ when $y\rightarrow 0$, follows a GPD as well.

This leads to the following consequences for any candidate CDF $G(\cdot)$:

1. a.

   the upper-tail behavior of $1-F(y)$ is equivalent to the original GPD tail used to build $F(y)$ ;

2. b.

   when $y\rightarrow 0$, $\displaystyle F(y)\sim\frac{c}{\psi^{s}}y^{s}$ where $\displaystyle c=\lim_{y\rightarrow 0}\frac{G(y)}{y^{s}}$ is positive and finite for some real $s$.

Four parametric families, $G(\cdot,\kappa)$, $\kappa=(\kappa_{1},\ldots,\kappa_{k})^{T}$, that satisfies the above-mentioned constraints are proposed in Naveau et al. (2016).

1. 1.

   $G(u;\kappa)=u^{\kappa_{1}}$, $\kappa_{1}>0$

2. 2.

   $G(u;\kappa)=\kappa_{3}u^{\kappa_{1}}+(1-\kappa_{3})u^{\kappa_{2}}$, $\kappa_{2}\geq\kappa_{1}>0$ and $\kappa_{3}\in[0,1]$

3. 3.

   $G(u;\kappa)=1-Q_{\kappa_{1}}\{(1-u)^{\kappa_{1}}\}$, $\kappa_{1}>0$ where $Q_{\kappa_{1}}$ is the CDF of a Beta random variable with parameters $1/\kappa_{1}$ and $2$, that is:

   $$Q_{\kappa_{1}}(u)=\frac{1+\kappa_{1}}{\kappa_{1}}u^{1/\kappa_{1}}\Big{(}1-\frac{u}{1+\kappa_{1}}\Big{)}$$

4. 4.

   $G(u;\kappa)=[1-Q_{\kappa_{1}}\{(1-u)^{\kappa_{1}}\}]^{\kappa_{2}/2}$, $\kappa_{1},\kappa_{2}>0$

In model 1, $\xi$ controls the rate of upper tail decay, $\kappa_{1}$ the shape of the lower tail and $\psi$ is a scale parameter. $\kappa_{1}=1$ allows us to recover the GPD model, while varying it gives flexibility in the description of the lower tail.

Model 2, as a mixture of power laws, allows to increase the flexibility provided by model 1. Therein, $\kappa_{1}$ controls the lower tail behavior while $\kappa_{2}$ modifies the shape of the density in its central part.

Model 3 is connected to the work of Falk et al. (2010). $\xi$ keeps controlling the upper extreme tail, $\kappa_{1}$ the central part of the distribution in a threshold tuning way. However, this model imposes a behavior type $\chi^{2}$ for the lower tail (null density at $0$) instead of allowing the data to constrain it.

Hence Model 4, which adds an extra parameter $\kappa_{2}$ to circumvent this limitation. In this last model, $\kappa_{2},\kappa_{1},\xi$ control respectively the lower, moderate and upper parts of the distribution. By construction, the lower and upper tails are GPD of shape parameters $\kappa_{2}$ and $\xi$ respectively.

In the sequel we denote the vector of parameters $(\xi,\psi,\kappa^{T})^{T}$ with $\theta=(\theta_{1},...,\theta_{d})^{T}$ and the density of a parametrized EGPD with $EGPD(y;\theta)$. To simulate from Equation (2), the user randomly draws a uniform variable $U$ in $[0,1]$ and then applies the quantile function $Y=F^{-1}(U)$, given for $p\in[0,1]$ by:

In practice, the parameters of the distribution of $Y$ may depend on some covariates $\boldsymbol{x}=(x_{1},\ldots,x_{m})^{T}$, i.e.

where $\theta(\boldsymbol{x})=(\theta_{1}(\boldsymbol{x}),\ldots,\theta_{d}(\boldsymbol{x}))^{T}$. The previous specification is an instance of a distributional regression model (Stasinopoulos et al., 2018).

For relating the different distributional parameters $(\theta_{1}(\boldsymbol{x}),\ldots,\theta_{d}(\boldsymbol{x}))$ to the covariates, we rely on additive predictors of the form

where $f_{i1}(\cdot),\ldots,f_{iJ_{i}}(\cdot)$ are smooth functions of the covariates $\boldsymbol{x}$.

The representation (5) is particularly flexible and allows to capture complex dependence patterns between distribution parameters and covariates $\boldsymbol{x}$. For instance a semiparametric model with two covariates $(x_{1},x_{2})$ can be written as

Another example is when we want to model space-time data observed at site $(s_{1},s_{2})$ and at time $t$. In this case $\boldsymbol{x}=(s_{1},s_{2},t)^{T}$

In analogy to Generalized Linear Model the predictors are linked to the distributional parameters via known monotonic and twice differentiable link functions $h_{i}(\cdot)$.

In the case of model 1, i.e. $G(u,\kappa)=u^{\kappa}_{1}$, common link functions are

The functions $f_{ij}$in (5) are approximated in terms of basis function expansions

where $B_{k}(\boldsymbol{x})$ are the basis functions and $\beta_{ij,k}$ denote the corresponding basis coefficients. These basis can be of different types, thoroughly presented e.g. in Stasinopoulos et al. (2017) or Wood (2017).

The basis function expansions can be written as

where $\boldsymbol{z}_{ij}(\boldsymbol{x})$ is still a vector of transformed covariates that depends on the basis functions. and $\boldsymbol{\beta}_{ij}=(\beta_{ij,1},\ldots,\beta_{ij,K_{ij}})^{T}$ is a parameter vector to be estimated.

To ensure regularization of the functions $f_{ij}(\boldsymbol{x})$ so-called penalty terms are added to the objective log-likelihood function. Usually, the penalty for each function $f_{ij}(\boldsymbol{x})$ are quadratic penalty $\lambda\boldsymbol{\beta}_{ij}^{T}\boldsymbol{G}_{ij}(\boldsymbol{\lambda}_{ij})\boldsymbol{\beta}_{ij}$ where $\boldsymbol{G}_{ij}(\boldsymbol{\lambda}_{ij})$ is a known semi-definite matrix and the vector $\boldsymbol{\lambda}_{ij}$ regulates the amount of smoothing needed for the fit. A special case is when $\boldsymbol{G}_{ij}(\boldsymbol{\lambda}_{ij})=\lambda_{ij}\boldsymbol{G}_{ij}$. Therefore the type and properties of the smoothing functions are controlled by the vectors $\boldsymbol{z}_{ij}(\boldsymbol{x})$ and the matrices $\boldsymbol{G}_{ij}(\boldsymbol{\lambda}_{ij})$.

The penalized log-likelihood function for the latter models reads:

where $l$ represents the log-likelihood function for the model (4).

We illustrate our add-on to package gamlss with the analysis of the hourly precipitations recorded over the north-west area of France from 2016 to 2018 both included. This database, containing precipitations as well as other meteorological variables whose values are registered every 6 minutes over the 3 years of interest for each of the 274 stations of the network (see Figure 1), is freely accessible (see Larvor et al., 2020, fro more details).

Let us note $Y(\boldsymbol{x})$ the random variable for strictly positive hourly precipitations (dry and wet events are traditionally treated separately when it comes to rainfall simulation (Ailliot et al., 2015)), with covariates $\boldsymbol{x}=(x_{l},x_{L},x_{t})$, where $x_{l}$ and $x_{L}$ represents the position in space (longitude, latitude) and $x_{t}$ a cyclic time index, namely the day or the month of the year. The work of Naveau et al. (2016) and our tests show that EGPD fitting is improved by censoring the data, so in practice we remove values $Y(\boldsymbol{x})<0.5$ (mm). Other thresholds ($\leq 0.2$ and $\leq 0.1$mm) have been tested but $Y(\boldsymbol{x})<0.5$ (mm) gives the best results, in particular for the upper tails. Besides, as noted in Evin et al. (2018); Naveau et al. (2016), negative shape parameters $\xi$ are not expected for rainfall distributions at these time scales, so we use a link function constraining the parameter $\xi$ to positive values. Finally, as performed in Naveau et al. (2016), we only keep one every three hourly records of precipitation so as to ensure independency between samples.

We fit different specifications of EGPD-based models, where EGPD takes form 1, 3 or 4, and compare them with equivalently specified Gamma-based models. The Gamma (GA) distribution is a classical choice that works well for modeling the bulk of precipitations at a given site (Katz, 1977; Vrac and Naveau, 2007; Vlček and Huth, 2009). However non-negligible deviations are observed when one is interested in capturing extreme rainfalls (Katz et al., 2002), typically underestimated due to the lightness of the GA tail. For a marginal model M=EGPD1,EGPD3,EGPD4,GA, the different specifications read:

M.0 where

$$\eta_{a}(\boldsymbol{x})=\beta_{0},\qquad a=\mu,\sigma,\nu,\tau,$$

i.e. the parameters $\theta(\boldsymbol{x})$ do not depend on the covariates ;

M.t where

$$\eta_{a}(\boldsymbol{x})=\beta_{0}+s_{t}(x_{t}),\qquad a=\mu,\sigma,\nu,\tau,$$

uses cyclic cubic splines $s_{t}$ to model the nonlinear dependency in time while considering no variation in space ;

M.tnomu where

$$\eta_{a}(\boldsymbol{x})=\beta_{0}+s_{t}(x_{t}),\qquad a=\sigma,\nu,\tau,$$

does the same as before but considers as constant the parameter $\mu$ ;

M.st where

$$\eta_{a}(\boldsymbol{x})=\beta_{0}+TP(x_{l},x_{L})+s_{t}(x_{t}),\qquad a=\mu,\sigma,\nu,\tau,$$

uses smooth surface fitting (thin-plate splines) to model the dependency in space as well as cyclic cubic splines for the dependency in time ;

M.st2mu where

$$\eta_{a}(\boldsymbol{x})=\beta_{0}+TP(x_{l},x_{L})+s_{t}(x_{t}),\qquad a=\mu,\sigma,$$

does the same as before but considers constant in space the parameters that are not $\mu$ or $\sigma$ ;

M.st2nomu where

$$\eta_{a}(\boldsymbol{x})=\beta_{0}+TP(x_{l},x_{L})+s_{t}(x_{t}),\qquad a=\sigma,\nu,$$

does the same as before but considers constant in space the parameters $\mu$ and $\tau$ ;

M.st3nomu where

$$\eta_{a}(\boldsymbol{x})=\beta_{0}+TP(x_{l},x_{L})+s_{t}(x_{t}),\qquad a=\sigma,\nu,\tau,$$

does the same as before but applies to M=EGPD4 only and considers parameter $\mu$ as a constant over space.

These variations within each model class will allow us to assess whether the smooth regression of distribution parameters over covariates is justified or not. Below we reproduce the code to fit megpd1.0, megpd1.tnomu and mga.st on the training set (training), that contains 60 % of the stations, as reported on Figure 1. The dataset consisting of the remaining stations form the validation set (validation). The gamlss() function reads:

In order to evaluate the goodness of fit on the training set, we use summary statistics: Global Deviance (GD), Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). The GAIC(modelName,k=a) function allows to compute GD (minus twice the fitted log-likelihood, $a=0$), the AIC ($a=2$) and the BIC ($a=\log(\mbox{number of observations})$) for each fitted model modelName and to compare them as performed on Figure 2. Besides, Table 1 indicates the effective degree of freedom (df) computed for each model. The function call is reproduced below for the models of class EGPD1.

From Figure 2, we can see that, as could be expected, for each classes of models (where EGPD1-based models represent one class, EGPD3-based ones another, etc) the GD ranks first models with the full space-time description M.st. All other variations are rather close in fitting performance, at the exception of the constant-parameters one M.0, that is largely behind. This shows the clear interest of allowing a smooth variation of distribution parameters as a function of the chosen covariates. However, when it comes to modeling efficiency, where the trade-off between number of parameters and performances is taken into account (through the AIC and BIC), the full space-time model M.st becomes much less desirable, falling just above and occasionally behind the constant parameter models M.0. This is all the more through that we consider the BIC indice, that penalizes more the number of parameters than the AIC, and a model with a large number of parameters (EGPD4 versus GA). The best trade-off, according to the BIC, are the time-only full model M.t for the GA class and the time-only without the shape parameter $\mu$ otherwise. This may also be an indication that the seasonal variation of distribution parameters is more important than their spatial variation over the area of interest.

The ranking of model variations within a same model class allows us to draw conclusions w.r.t. the advantages or challenges of each class. Thus, for EGPD4 the GD ranks model variations following decreasing complexity: megpd4.st followed by megpd4.st3nomu, etc. We can deduce that the estimation of $\mu$ is not that sensitive when it comes to the fitting set (at least with enough data). On the contrary, megpd1.st2nomu is preferred to megpd1.st and megpd1.st2mu which can be interpreted as $\mu$ being a rather sensitive parameter to estimate for EGPD1. Considering it constant in space consequently improves the global fitting.

Performances on the validation set are obtained by means of the getTGD() and TGD() functions, providing the validation global deviance (GD-v), that is the GD evaluated using predictive values for the parameters at the validation sample.

Thus, when it comes to extrapolating parameters for new covariate values, the GD-v ranking displayed in Figure 2 reveals that models considering only the dependence in time (M.t, M.tnomu) are systematically preferred, followed by space-time models not considering $\mu$ (M.st2nomu), before the full space-time models (M.st). This supports the fact that the rainfall distributions are relatively close over the studied geographical area, compared to their time evolution that shows stronger variations. Consequently, given the limited size of the dataset, keeping low the number of parameters to regress on covariates (by considering the seasonality only) and in particular avoiding in this category the shape parameter $\mu$, makes estimations more robust. This agrees with observations from Evin et al. (2018), who noted that fitting the shape parameter is the most challenging part of (extended)-GPD modeling. They concluded that mono-site fitting is not robust, leading e.g. to negative shape parameters that do not match observations on rainfall intensities. To circumvent this issue, the authors use constant parameters within regions of influence defined by homogeneity tests.