# Space-Time Trend Detection and Dependence Modeling in Extreme Event Approaches by Functional Peaks-Over-Thresholds: Application to Precipitation in Burkina Faso

Authors: Sawadogo Béwentaoré; Diakarya Barro
Journal: International Journal of Mathematics and Mathematical Sciences
 (2022)

---

## Abstract

In this paper, we propose a new method for estimating trends in extreme spatiotemporal processes using both information from marginal distributions and dependence structure. We combine two statistical approaches of an extreme value theory: the temporal and spatial nonstationarities are handled via a tail trend function in the marginal distributions. The spatial dependence structure is modeled by a latent spatial process using generalizedℓ-Pareto processes. This methodology for trend analysis of extreme events is applied to precipitation data from Burkina Faso. We show that a significant increasing trend for the 50 and 100 year return levels in some parts of the country. We also show that extreme precipitation is spatially correlated with distance for a radius of approximately 200 km.


---

## Body

## 1. Introduction

In the framework of climate change, the modeling and accurate prediction of the magnitude and extent of extreme events that occur in space and time of climate variables is a particularly difficult task, since it implies taking into account spatial and temporal nonstationarities. Nowadays, there is a general consensus in the scientific community that climate change has accelerated in recent decades and that the climate will continue to change in the coming decades, mainly due to natural and anthropogenic changes (IPCC 2007, 2018, 2019). This change is manifested in most regions of the world by a resurgence of heavy rainfall, heat waves, and pollution peaks with very significant economic and social consequences.In sub-Saharan Africa, works on this topic have shown an increasing trend in the occurrence of extremes in meteorological parameters [1–4]. In Burkina Faso, a study conducted in the framework of the National Action Program for Adaptation to climate change (PANA, 2007) showed that precipitation is highly variable in space and time. There is a decreasing trend in cumulative and daily rainfall. According to recent publications from the National Meteorology Agency, precipitation has returned to humid periods since the late 1980s and decades of 1990 and 2000. The return of rains is more related to a high frequency of high intensity rainfall events than to an increase in the number of rainy days. For example, we can cite the floods of 1st September 2009 that affected Ouagadougou and its outskirts with a record of 261.3 mm of rainfall not registered since 1919. The floods have affected more than 150000 peoples, damaged several bridges, and flooded more than 9300 hectares of crops throughout the country. From a statistical point of view, this raises the question of trend detection in extreme events.The classical extreme value theory extended both to nonstationary and non-independent observations provides a rigorous mathematical framework to deal with this question of trend detection in extremes [5–8]. This issue was first generally studied in extreme value literature by parametric pointwise approaches in which an extreme value model (GEV or GPD) is fitted to the data at each site in turn, leaving the parameters of the marginal distributions to evolve with time or other significant covariates [9–13]. Although it is relatively simple to construct nonstationary models in the univariate framework, it is more difficult to account for spatial and temporal trends in these univariate models. The spatial and temporal dimension of extreme events was first developed for block maxima in a stationary framework and modeled by max-stable processes [14–18]. These spatial models have been readapted to handle nonstationarities induced by global warming [19, 20]. Although attractive, these models are expensive in computing time when tending toward large dimensions and are not adapted for modeling threshold exceedances.The threshold exceedance approach introduced by [21, 22] was first extended to multivariate environments [23] before being generalized to functional data [24–27] to give birth to the family of generalized ℓ-Pareto processes. This family of processes are good candidates for modeling the spatial dependence structure of exceedance data. However, these approaches do not sufficiently take into account marginal nonstationarities and spatial dependence between margins. In this paper, we propose a new method to capture nonstationarities in marginal distributions, while taking into account the spatial dependence structure. To reach this goal, we combine two statistical approaches to the extreme value theory. First temporal and spatial nonstationarities are controlled via a tail trend function in marginal distributions [5–7]. Subsequently, the spatial dependence structure is taken into account by a hidden auxiliary spatial process using generalized ℓ-Pareto processes [25, 26]. Finally, we use the model to simulate and predict future extreme events.The article is organized as follows: Section2 details the methodology implemented and the methods used to estimate the parameters. In Section 3, we present our main results from the data analysis, and in particular, we calculate the nonstationary return levels from the developed approaches. Section 4 concludes the study.


## 2. Methodology

### 2.1. Space-Time Trends Modeling

LetX=Xts,s∈S,t∈T be a continuous nonstationary space-time stochastic process with sample paths in the family of continuous functions CS×T equipped with the uniform norm ∥.∥∞, where S×T⊂ℝd×ℝ+ and C+S×T its restriction to nonnegative functions deprived of the null function. In practice X is observed at each stations s1,…,sm and at given dates t=1,…,n. Let Ft,s be the continuous univariate marginal distribution with a common right endpoint xF and Z=Zs,s∈S an unobserved latent spatial stochastic process with sample paths in CS×T satisfying the proportional tail condition such that(1)limx⟶xFPXtsj>xPZsj>x=cθtn,sj,with,1m∑j=1m∫01cθu,sjdu=1,u∈0,1,where cθ:0,1×S⟶0,∞ is assumed to be a continuous and positive function depending on a parameter vector θ∈Θ⊂ℝ called a tail trend function or skedasis function [5, 7]. The skedasis function describes the evolution of extreme events jointly in space and time. In the framework of the model (1), Mefleh et al. [8] shows that the empirical point measure converges in distribution in the space of point measure ℳp=ℳp0,1×0,∞ to a Poisson point process with intensity measure cuduγz−γ+1dz on 0,1×0,∞. In this case, the times of exceedances for high threshold x and the value of exceedances are asymptotically independent with distributions respectively equal to the trend density function cθu,s,s∈S and the Pareto distribution of tail index γ.Moreover, we assume that the continuous marginal distributionsFZ of the latent process is in the maximum domain of attraction condition for some constant γ∈ℝ and appropriate real normalization constants a˜n>0, b˜n∈ℝ. Thanks to the equation (1) and the convergence of Z exceedances to a GPD distribution we deduce a sample of Zts,s∈S as mentioned in [5, 7] in the following manner:(2)Ztsj=cθtn,sj−γXtsj−cθt/n,sjγ−1γa˜n−γb˜n,j=1,…,m;t=1,…,n,where γ,a˜n,b˜n, and cθ are parameters to be estimated. Parameter inference techniques will be discussed in Section 2.4.1.


### 2.2. Nonstationary Return Period and Return Level

The concept of return period and level becomes very ambiguous when we leave the stationary context to the nonstationary framework. For example, in the stationary framework, an extreme event with a 100-year return period is likely to occur on average once every 100 years, but with an annual exceedance probability of1% in a given year. In this case, it is assumed that the underlying probability distribution remains unchanged over time. More formally, for a stationary random variable Z with distribution function FZ, the m-year return period is expressed as m=1/FZZ>zm, where zm is the return level associated with the return period. However, in the nonstationary context, the underlying probability distribution is no longer constant but is supposed to evolve with time. In this paper, we have chosen to follow return period based approaches, i.e., expected number of exceedances (ENE) [12] and expected waiting time (EWT) [9, 11]. In the ENE approach, the number of times the variable Xts,s∈S exceeds the return level value xm in m years is defined by Nm=∑t=1nxm1Xts>xm,s∈S under nonstationary context. The return level xm can be defined as the value for which the expected number of events exceeding xm in m years equals to one, i.e., the return level xm is the solution of the following equation:(3)1=∑t=1nxm1−Ft,sxm|θts,s∈S,where nx is the number of days in the year and θts the vector of time-dependent marginal parameters or other covariates. Parey et al. [12] uses the ENE method in a pointwise POT model where the parameters of the distribution of exceedances and the intensity of extreme event occurrences are described as polynomial functions of time.The EWT method was first proposed by Olsen et al. [11], and then derived by Salas and Obeysekera [13] using a geometric distribution with time-varying parameters. Under nonstationary conditions, the distribution describing waiting time Y before the first occurrence of an event exceeding the return level xm is as follows:(4)fsy=PYs=y=py,s∏i=1y−11−pi,s,y=1,2,…,ymax,s∈S,where variable Y is the day of the first occurrence of an event exceeding the quantile and pt,s=1−Ft,sxm|θts is daily exceedance probability varying with time step t. ymax is the time when the daily exceedance probability pt,s is equal to 1 for an increasing-trend series or is equal to 0 for a decreasing-trend series. Reused and simplified by [9], the EWT approach defines the m-year return level xm, as the value for which the expected waiting time until an exceedance of this level is m years, i.e, xm is the solution of the equation:(5)nxm=EY=1+∑y=1∞∏t=1yFt,sxm|θts,s∈S.Using the relationship (1),PXts>xm can be rewritten ascθt/n,sPZs>uPZs>xm|Zs>u for u<xm and we obtain the following results.Proposition 1.
LetXts,s∈S,t∈Tbe a nonstationary stochastic process defined on a regionS⊂ℝdandZs,s∈Sa latent spatial process satisfying equation (1). Given a return periodmand thresholdu<xm, the return levelxmfor alls∈Sis a solution of the following two equations:(i)
Return period as expected number of events(6)1=∑t=t0t0+nxmcθtn,sϕusF¯Zs|ℓZ>uxm−u.(ii)
Return period as expected waiting time(7)nxm=1+∑i=1∞∏t=1i1−cθtn,sϕusF¯Zs|ℓZ>uxm−u, where F¯Zs|ℓZ>u is a survival of generalized Pareto distribution of ℓ-exceedances at position s∈S; ϕus=PZs>u is the probability of exceedances, cθ is the tail trend function and nx is the number of days in the year.The nonstationary return levels obtained from (6) and (7) are evaluated by numerical algorithms taking into account information from the extrapolation of the trend function cθ. To derive the nonstationary return level at points where we have not observations, we use spatial marginal model (15) described in Section 2.4.1 to extrapolate the parameters at these points. An alternative approach is to calculate the nonstationary return levels from the following result.Proposition 2.
LetXts,s∈S,t∈Tbe a nonstationary spatio-temporal stochastic process defined on a regionS⊂ℝdandZs,s∈Sa latent spatial process. The nonstationary return levelxms,s∈Sof the nonstationary processXis deduced from the return levelzmof the latent spatial processZsuch that(8)xms=zmsc^θtm,sγ^+c^θtm,sγ^−1γ^a˜^n−γ^b˜^n,s∈S,wherea˜^n,b˜^n,c^θ, andγ^are the respective estimators ofa˜n,b˜n,cθ, andγ.tm=1+nxm/nwithnxthe number of days in the year and n the size of the sample observed.This result is a consequence of the (2). In practice, the return level zm of the latent spatial process is first computed at grid points where we have not observations using (15). Then, the nonstationary return level xm is deduced using the Proposition 2.


### 2.3. Spatial POT Modeling

After removing the nonstationarity, thanks to the (2), the modeling is then focused on the evaluation of the extreme spatial dependence structure in Z using functional POT [24–26]. In the multivariate and spatial framework a threshold exceedances for a random function Z=Zs,s∈S is defined by [26] to be an event of the form ℓZ>u for some u≥0, where ℓ:CS⟶ℝ+ is a continuous and homogeneous nonnegative risk function, i.e., there exists α>0 such that ℓλy=λαℓy when y∈C+S and λ>0. The risk function ℓ determines the type of extreme events of interest. For example, such a function can be the maximum, minimum, average, or value at a specific point s0∈S. Under minimal assumptions on the risk function, the conditional distribution of ℓ-exceedance for some threshold u≥0 of the process Z−bn/an can be approximated by a generalized ℓ-Pareto process, for n large enough [25].(9)PZ−bnℓan∈A|ℓZ−bnℓan≥u⟶PWℓ∈A,n⟶∞,where Wℓ is a nondegenerate stochastic process over S and belongs to the family of generalized ℓ-Pareto processes with tail index γ (see appendix B for more details). Specifically, Wℓ is a stochastic process taking values in z∈C+S:ℓz−b/ℓa≥0 and defined by the following:(10)Wℓ=aYℓγ−1γ+b,γ≠0,alogYℓ+b,γ=0,where a>0 and b are continuous functions on S, respectively, scale and location functions and Yℓ is a stochastic process whose probability measure is completely determined by the limit measure Λ. More details on generalized ℓ-Pareto processes can be found in [24, 25].For the modeling of spatial dependence structure of the latent process, we useYℓ whose margins are in the Frechet domain of attraction with tail index γ=1, as the process of reference. A pseudo-polar decomposition of Yℓ in (11) leads to the following formulation [24, 26].(11)Yℓ=RQ,where R is a unit Pareto random variable of index γR representing the intensity of process and Q is stochastic process denoted the angular component with state space S and taking values in S=y∈C+S:y1=1 whose probability measure is characterized by limit measure Λ.


### 2.4. Statistical Inference

#### 2.4.1. Marginal Model of the Spatial Latent Process

Marginal parametersan, bn, and γ of latent process can be estimated into fixing b˜n at a high quantile of ℓZ, i.e., b˜^n=q1−αℓZ. In general, a parametric model may be necessary for an and bn, as in Engelke et al. [28]. However, some forms of parameterization can lead to problems of parameter identifiability and inference. For these reasons and simplicity purpose, we consider in this work that ansj=aj>0 and bnsj=bj∈ℝ for any j=1,…,m. The threshold stability of the generalized Pareto distributions does not allow us to identify the bn function without additional assumptions, so we set(12)bns=uq′Zs−b˜n,s∈S,where uq′Zs is an empirical quantile of the order q′ of the ℓ-exceedances at each location s∈S, where q′ is chosen such that ℓbn=b˜n in order to impose the identifiability of the parameters. Thus, the tail index γ∈ℝ and the scale parameters ans>0 are estimated by maximizing the independent log-likelihood; that is,(13)ℓindepγ,ansj=∑t=1n∑j=1m1xtsj≥bnsjlog1ansj1+γxtsj−bnsjansj+−1/γ−1.The assumptions of parameter identifiability (see appendix B) implies thatℓan=a˜n.The trend function parametercθ is estimated using the maximum likelihood method under the assumption that the dates of exceeding the marginal thresholds b˜^n are asymptotically independent and identically distributed, from the density function cθt/n,s [8]. Several models are eligible to model the function cθ. However, in this study, we opt for parametric models because of their flexibility in trend detection and in order to make extrapolations of the trend beyond the observed data. To facilitate inference by the maximum likelihood method, we are interested in monotonic log-linear and simple linear trend models of cθ such that(14)cθ1tn,s=θsexpθs−1expθstn,θs∈ℝ;tn∈0,1cθ2tn,s=θs2tn−1+1,θs∈−1;1tn∈0,1,where θs,s∈S is a characteristic parameter of the trend function. θs>0 suggests that extreme are getting more frequent, while θs<0 indicates that extremes are getting less frequent.


#### 2.4.2. Spatial Extremal Dependence

The spatial dependence structure is taken into account in the marginal parameters of spatial processZ by letting these parameters evolve as a function of space covariates. More precisely, we use a spatial model of the marginal parameters as functions of the significant space covariates using a generalized linear regression model [9, 12, 13].(15)ans=ganvs>0,bns=gbnvs∈ℝθs=gθvs∈ℝγs=γ0∈ℝ,,s∈S,where v is a vector of covariates which can contain for example the space coordinates, or any other spatialized data explaining significantly the weather variable studied.Furthermore, after removing the nonstationarity, modeling is focused on the evaluation of the remaining spatial dependence structure inZ. Thus, we approach the limit distribution of ℓ-exceedances of spatial latent process by a generalized ℓ-Pareto process (9). To estimate the parameters of spatial dependence structure, we have chosen to model its angular component Q (11) by a log-Gaussian process with stationary increments and within the framework of the Brown–Resnick model. In this case, we have access to the limit measure Λ. In order to better capture this dependence structure, we use a flexible parametric semi-variogram belonging to the class of power semi-variograms. The parameter of the dependence structure is estimated using the gradient score method or the censored likelihood method [24, 29].


## 2.1. Space-Time Trends Modeling

LetX=Xts,s∈S,t∈T be a continuous nonstationary space-time stochastic process with sample paths in the family of continuous functions CS×T equipped with the uniform norm ∥.∥∞, where S×T⊂ℝd×ℝ+ and C+S×T its restriction to nonnegative functions deprived of the null function. In practice X is observed at each stations s1,…,sm and at given dates t=1,…,n. Let Ft,s be the continuous univariate marginal distribution with a common right endpoint xF and Z=Zs,s∈S an unobserved latent spatial stochastic process with sample paths in CS×T satisfying the proportional tail condition such that(1)limx⟶xFPXtsj>xPZsj>x=cθtn,sj,with,1m∑j=1m∫01cθu,sjdu=1,u∈0,1,where cθ:0,1×S⟶0,∞ is assumed to be a continuous and positive function depending on a parameter vector θ∈Θ⊂ℝ called a tail trend function or skedasis function [5, 7]. The skedasis function describes the evolution of extreme events jointly in space and time. In the framework of the model (1), Mefleh et al. [8] shows that the empirical point measure converges in distribution in the space of point measure ℳp=ℳp0,1×0,∞ to a Poisson point process with intensity measure cuduγz−γ+1dz on 0,1×0,∞. In this case, the times of exceedances for high threshold x and the value of exceedances are asymptotically independent with distributions respectively equal to the trend density function cθu,s,s∈S and the Pareto distribution of tail index γ.Moreover, we assume that the continuous marginal distributionsFZ of the latent process is in the maximum domain of attraction condition for some constant γ∈ℝ and appropriate real normalization constants a˜n>0, b˜n∈ℝ. Thanks to the equation (1) and the convergence of Z exceedances to a GPD distribution we deduce a sample of Zts,s∈S as mentioned in [5, 7] in the following manner:(2)Ztsj=cθtn,sj−γXtsj−cθt/n,sjγ−1γa˜n−γb˜n,j=1,…,m;t=1,…,n,where γ,a˜n,b˜n, and cθ are parameters to be estimated. Parameter inference techniques will be discussed in Section 2.4.1.


## 2.2. Nonstationary Return Period and Return Level

The concept of return period and level becomes very ambiguous when we leave the stationary context to the nonstationary framework. For example, in the stationary framework, an extreme event with a 100-year return period is likely to occur on average once every 100 years, but with an annual exceedance probability of1% in a given year. In this case, it is assumed that the underlying probability distribution remains unchanged over time. More formally, for a stationary random variable Z with distribution function FZ, the m-year return period is expressed as m=1/FZZ>zm, where zm is the return level associated with the return period. However, in the nonstationary context, the underlying probability distribution is no longer constant but is supposed to evolve with time. In this paper, we have chosen to follow return period based approaches, i.e., expected number of exceedances (ENE) [12] and expected waiting time (EWT) [9, 11]. In the ENE approach, the number of times the variable Xts,s∈S exceeds the return level value xm in m years is defined by Nm=∑t=1nxm1Xts>xm,s∈S under nonstationary context. The return level xm can be defined as the value for which the expected number of events exceeding xm in m years equals to one, i.e., the return level xm is the solution of the following equation:(3)1=∑t=1nxm1−Ft,sxm|θts,s∈S,where nx is the number of days in the year and θts the vector of time-dependent marginal parameters or other covariates. Parey et al. [12] uses the ENE method in a pointwise POT model where the parameters of the distribution of exceedances and the intensity of extreme event occurrences are described as polynomial functions of time.The EWT method was first proposed by Olsen et al. [11], and then derived by Salas and Obeysekera [13] using a geometric distribution with time-varying parameters. Under nonstationary conditions, the distribution describing waiting time Y before the first occurrence of an event exceeding the return level xm is as follows:(4)fsy=PYs=y=py,s∏i=1y−11−pi,s,y=1,2,…,ymax,s∈S,where variable Y is the day of the first occurrence of an event exceeding the quantile and pt,s=1−Ft,sxm|θts is daily exceedance probability varying with time step t. ymax is the time when the daily exceedance probability pt,s is equal to 1 for an increasing-trend series or is equal to 0 for a decreasing-trend series. Reused and simplified by [9], the EWT approach defines the m-year return level xm, as the value for which the expected waiting time until an exceedance of this level is m years, i.e, xm is the solution of the equation:(5)nxm=EY=1+∑y=1∞∏t=1yFt,sxm|θts,s∈S.Using the relationship (1),PXts>xm can be rewritten ascθt/n,sPZs>uPZs>xm|Zs>u for u<xm and we obtain the following results.Proposition 1.
LetXts,s∈S,t∈Tbe a nonstationary stochastic process defined on a regionS⊂ℝdandZs,s∈Sa latent spatial process satisfying equation (1). Given a return periodmand thresholdu<xm, the return levelxmfor alls∈Sis a solution of the following two equations:(i)
Return period as expected number of events(6)1=∑t=t0t0+nxmcθtn,sϕusF¯Zs|ℓZ>uxm−u.(ii)
Return period as expected waiting time(7)nxm=1+∑i=1∞∏t=1i1−cθtn,sϕusF¯Zs|ℓZ>uxm−u, where F¯Zs|ℓZ>u is a survival of generalized Pareto distribution of ℓ-exceedances at position s∈S; ϕus=PZs>u is the probability of exceedances, cθ is the tail trend function and nx is the number of days in the year.The nonstationary return levels obtained from (6) and (7) are evaluated by numerical algorithms taking into account information from the extrapolation of the trend function cθ. To derive the nonstationary return level at points where we have not observations, we use spatial marginal model (15) described in Section 2.4.1 to extrapolate the parameters at these points. An alternative approach is to calculate the nonstationary return levels from the following result.Proposition 2.
LetXts,s∈S,t∈Tbe a nonstationary spatio-temporal stochastic process defined on a regionS⊂ℝdandZs,s∈Sa latent spatial process. The nonstationary return levelxms,s∈Sof the nonstationary processXis deduced from the return levelzmof the latent spatial processZsuch that(8)xms=zmsc^θtm,sγ^+c^θtm,sγ^−1γ^a˜^n−γ^b˜^n,s∈S,wherea˜^n,b˜^n,c^θ, andγ^are the respective estimators ofa˜n,b˜n,cθ, andγ.tm=1+nxm/nwithnxthe number of days in the year and n the size of the sample observed.This result is a consequence of the (2). In practice, the return level zm of the latent spatial process is first computed at grid points where we have not observations using (15). Then, the nonstationary return level xm is deduced using the Proposition 2.


## 2.3. Spatial POT Modeling

After removing the nonstationarity, thanks to the (2), the modeling is then focused on the evaluation of the extreme spatial dependence structure in Z using functional POT [24–26]. In the multivariate and spatial framework a threshold exceedances for a random function Z=Zs,s∈S is defined by [26] to be an event of the form ℓZ>u for some u≥0, where ℓ:CS⟶ℝ+ is a continuous and homogeneous nonnegative risk function, i.e., there exists α>0 such that ℓλy=λαℓy when y∈C+S and λ>0. The risk function ℓ determines the type of extreme events of interest. For example, such a function can be the maximum, minimum, average, or value at a specific point s0∈S. Under minimal assumptions on the risk function, the conditional distribution of ℓ-exceedance for some threshold u≥0 of the process Z−bn/an can be approximated by a generalized ℓ-Pareto process, for n large enough [25].(9)PZ−bnℓan∈A|ℓZ−bnℓan≥u⟶PWℓ∈A,n⟶∞,where Wℓ is a nondegenerate stochastic process over S and belongs to the family of generalized ℓ-Pareto processes with tail index γ (see appendix B for more details). Specifically, Wℓ is a stochastic process taking values in z∈C+S:ℓz−b/ℓa≥0 and defined by the following:(10)Wℓ=aYℓγ−1γ+b,γ≠0,alogYℓ+b,γ=0,where a>0 and b are continuous functions on S, respectively, scale and location functions and Yℓ is a stochastic process whose probability measure is completely determined by the limit measure Λ. More details on generalized ℓ-Pareto processes can be found in [24, 25].For the modeling of spatial dependence structure of the latent process, we useYℓ whose margins are in the Frechet domain of attraction with tail index γ=1, as the process of reference. A pseudo-polar decomposition of Yℓ in (11) leads to the following formulation [24, 26].(11)Yℓ=RQ,where R is a unit Pareto random variable of index γR representing the intensity of process and Q is stochastic process denoted the angular component with state space S and taking values in S=y∈C+S:y1=1 whose probability measure is characterized by limit measure Λ.


## 2.4. Statistical Inference

### 2.4.1. Marginal Model of the Spatial Latent Process

Marginal parametersan, bn, and γ of latent process can be estimated into fixing b˜n at a high quantile of ℓZ, i.e., b˜^n=q1−αℓZ. In general, a parametric model may be necessary for an and bn, as in Engelke et al. [28]. However, some forms of parameterization can lead to problems of parameter identifiability and inference. For these reasons and simplicity purpose, we consider in this work that ansj=aj>0 and bnsj=bj∈ℝ for any j=1,…,m. The threshold stability of the generalized Pareto distributions does not allow us to identify the bn function without additional assumptions, so we set(12)bns=uq′Zs−b˜n,s∈S,where uq′Zs is an empirical quantile of the order q′ of the ℓ-exceedances at each location s∈S, where q′ is chosen such that ℓbn=b˜n in order to impose the identifiability of the parameters. Thus, the tail index γ∈ℝ and the scale parameters ans>0 are estimated by maximizing the independent log-likelihood; that is,(13)ℓindepγ,ansj=∑t=1n∑j=1m1xtsj≥bnsjlog1ansj1+γxtsj−bnsjansj+−1/γ−1.The assumptions of parameter identifiability (see appendix B) implies thatℓan=a˜n.The trend function parametercθ is estimated using the maximum likelihood method under the assumption that the dates of exceeding the marginal thresholds b˜^n are asymptotically independent and identically distributed, from the density function cθt/n,s [8]. Several models are eligible to model the function cθ. However, in this study, we opt for parametric models because of their flexibility in trend detection and in order to make extrapolations of the trend beyond the observed data. To facilitate inference by the maximum likelihood method, we are interested in monotonic log-linear and simple linear trend models of cθ such that(14)cθ1tn,s=θsexpθs−1expθstn,θs∈ℝ;tn∈0,1cθ2tn,s=θs2tn−1+1,θs∈−1;1tn∈0,1,where θs,s∈S is a characteristic parameter of the trend function. θs>0 suggests that extreme are getting more frequent, while θs<0 indicates that extremes are getting less frequent.


### 2.4.2. Spatial Extremal Dependence

The spatial dependence structure is taken into account in the marginal parameters of spatial processZ by letting these parameters evolve as a function of space covariates. More precisely, we use a spatial model of the marginal parameters as functions of the significant space covariates using a generalized linear regression model [9, 12, 13].(15)ans=ganvs>0,bns=gbnvs∈ℝθs=gθvs∈ℝγs=γ0∈ℝ,,s∈S,where v is a vector of covariates which can contain for example the space coordinates, or any other spatialized data explaining significantly the weather variable studied.Furthermore, after removing the nonstationarity, modeling is focused on the evaluation of the remaining spatial dependence structure inZ. Thus, we approach the limit distribution of ℓ-exceedances of spatial latent process by a generalized ℓ-Pareto process (9). To estimate the parameters of spatial dependence structure, we have chosen to model its angular component Q (11) by a log-Gaussian process with stationary increments and within the framework of the Brown–Resnick model. In this case, we have access to the limit measure Λ. In order to better capture this dependence structure, we use a flexible parametric semi-variogram belonging to the class of power semi-variograms. The parameter of the dependence structure is estimated using the gradient score method or the censored likelihood method [24, 29].


## 2.4.1. Marginal Model of the Spatial Latent Process

Marginal parametersan, bn, and γ of latent process can be estimated into fixing b˜n at a high quantile of ℓZ, i.e., b˜^n=q1−αℓZ. In general, a parametric model may be necessary for an and bn, as in Engelke et al. [28]. However, some forms of parameterization can lead to problems of parameter identifiability and inference. For these reasons and simplicity purpose, we consider in this work that ansj=aj>0 and bnsj=bj∈ℝ for any j=1,…,m. The threshold stability of the generalized Pareto distributions does not allow us to identify the bn function without additional assumptions, so we set(12)bns=uq′Zs−b˜n,s∈S,where uq′Zs is an empirical quantile of the order q′ of the ℓ-exceedances at each location s∈S, where q′ is chosen such that ℓbn=b˜n in order to impose the identifiability of the parameters. Thus, the tail index γ∈ℝ and the scale parameters ans>0 are estimated by maximizing the independent log-likelihood; that is,(13)ℓindepγ,ansj=∑t=1n∑j=1m1xtsj≥bnsjlog1ansj1+γxtsj−bnsjansj+−1/γ−1.The assumptions of parameter identifiability (see appendix B) implies thatℓan=a˜n.The trend function parametercθ is estimated using the maximum likelihood method under the assumption that the dates of exceeding the marginal thresholds b˜^n are asymptotically independent and identically distributed, from the density function cθt/n,s [8]. Several models are eligible to model the function cθ. However, in this study, we opt for parametric models because of their flexibility in trend detection and in order to make extrapolations of the trend beyond the observed data. To facilitate inference by the maximum likelihood method, we are interested in monotonic log-linear and simple linear trend models of cθ such that(14)cθ1tn,s=θsexpθs−1expθstn,θs∈ℝ;tn∈0,1cθ2tn,s=θs2tn−1+1,θs∈−1;1tn∈0,1,where θs,s∈S is a characteristic parameter of the trend function. θs>0 suggests that extreme are getting more frequent, while θs<0 indicates that extremes are getting less frequent.


## 2.4.2. Spatial Extremal Dependence

The spatial dependence structure is taken into account in the marginal parameters of spatial processZ by letting these parameters evolve as a function of space covariates. More precisely, we use a spatial model of the marginal parameters as functions of the significant space covariates using a generalized linear regression model [9, 12, 13].(15)ans=ganvs>0,bns=gbnvs∈ℝθs=gθvs∈ℝγs=γ0∈ℝ,,s∈S,where v is a vector of covariates which can contain for example the space coordinates, or any other spatialized data explaining significantly the weather variable studied.Furthermore, after removing the nonstationarity, modeling is focused on the evaluation of the remaining spatial dependence structure inZ. Thus, we approach the limit distribution of ℓ-exceedances of spatial latent process by a generalized ℓ-Pareto process (9). To estimate the parameters of spatial dependence structure, we have chosen to model its angular component Q (11) by a log-Gaussian process with stationary increments and within the framework of the Brown–Resnick model. In this case, we have access to the limit measure Λ. In order to better capture this dependence structure, we use a flexible parametric semi-variogram belonging to the class of power semi-variograms. The parameter of the dependence structure is estimated using the gradient score method or the censored likelihood method [24, 29].


## 3. Application to Extreme Precipitation in Burkina Faso

### 3.1. Data Set Analysis

This study uses time series of daily precipitation measurements from 1957 to 2016 provided by ten synoptic stations extracted from the Burkina Faso climatological database. These stations have been selected to ensure good spatial uniformity and representativeness of different climatic regimes and data quality. Figure1 gives the spatial distribution of selected stations in our study area. In order to limit the problems related to seasonal rainfall cycles on each station, we worked from the sub-series corresponding to rainy days. The period from May to October was chosen because it is during this period that the most rainfall is recorded in the country. Thus, a series of 60 years (1957-2016) by 184 days (May to October) is extracted to constitute the time series of daily rainfall. On these time series, we apply a run declustering procedure with a daily step (r = 1 day) to identify the groups of approximately independent extreme observations within the sample in order to avoid short-term dependencies in the time series.Figure 1
Spatial distribution of synoptic sampling stations associated with climatic regimes; source: National Meteorology Agency of Burkina Faso.


### 3.2. Marginal Characteristics

The first line of Figure2 shows the maps of the parameters an and bn estimated on our selected stations. Map (a) gives us the variability parameter of the precipitation distribution and shows a very high variability of extreme precipitation. While map (b) gives us the estimated level or threshold for each station using equation (12) and we see that the level is higher in the Sudanian climate and lower in the Sahelian climate.In practice, we deduce from the generalized ℓ-Pareto model fitted to the data, the marginal tail behavior of the survival distribution F¯Zs,γ^s,σ^s for any s∈S from the (16) given a sufficiently large uq>0 threshold.(16)PZs−uqs≥z|Zs≥uqs≈F¯Zs,γ^s,σ^sz,z≥0,with σ^s=a^ns+γ^uqs−b^ns, s∈S, where a^n, b^n, and γ^ are the marginal parameters estimators described in Section 2.4.1. We can use it to check the quality of the marginal adjustment of the stationary process Z. Figure 3 shows the qq-plots of the local tail distribution due to two stations by climatic zone. Columns 1, 2, and 3 of Figure 3 represent the respective adjustments for the Sudanian, Sudano-Sahelian, and Sahelian climate. Globally, the fit of the marginal models seems convincing, as most of the observations remain within the confidence intervals.Figure 2
Estimated functionsan (a) and bn (b) of the generalized ℓ-Pareto process for modelling extreme precipitation, and θ for Log-linear (c) and simple linear (d) trends. Estimates are obtained by shifting the local empirical quantiles u0.80ℓZs by b˜^n=33.4mm.
(a)(b)(c)(d)Figure 3
q-q plot of the distributions of the local tails of six synoptic positions including two stations for each climate region. (a) Sudanian climate; (b) Sudano-Sahelian climate; and (c) Sahelian climate. The threshold is taken equal to bns for each station.
(a)(b)(c)The second row of Figure2 shows the maps of parameter θ estimated for the log-linear and simple linear trend functions. We can observe a strong spatial variability in the occurrence of extreme precipitation. Thus, we observe a decreasing trend in the frequency of extreme rainfall in Ouaga, Dori, and Bobo. In contrast, in Boromo, Ouahigouya, and Pô, extreme precipitation events will become more and more frequent. Figure 4 gives us the details of the adjustment of the tail trend function by a log-linear model.Figure 4
Local fitting of the log-linear trend modelcθ1 for some of the selected sites with two stations for each climate region.


### 3.3. Estimated Spatial Dependence Model

We use a spatial marginal model (15) to estimate the marginal parameters at the different grid points. This is to prepare the ground for the calculation of nonstationary return levels at locations where we have no observations.(17)ans=a0+a1Long×Lats,bns=b0+b1Lat2s+b2Long×Lats,θs=θ0+θ1Long+θ2Lat,γs=γ0.,s∈S,An alternative to the spatial model of marginal parameters is to consider small subsetsNs0⊂S, where Ns0 is a small neighborhood around s0. In what follows, these regional neighborhoods will be determined by a small number D0 of nearest stations from the site s0, thus we will write Ns0=Ns0,D0. In principle, the choice of D0 should be such that the spatial dependence structure and marginal parameters is approximately stationary within each selected neighborhood Ns0,D0 around s0∈S. Obviously, the choice of neighborhood is important; the assumed stationary marginal parameters could be a poor approximation for large neighborhoods (i.e., for large D0), while the simulation of the process could be cumbersome for small neighborhoods characterized by a small D0 number of stations. We obtain the relatively homogeneous, nonoverlapping subregions using the k-means clustering method centered on the reference stations [30, 31]. This method is extensively used because it is computationally simple and produces accurate results, compared to other more complex clustering methods. The longitude and latitude of the grid points were used as input variables in the k-means clustering algorithm to form the ten clusters centered on the reference stations.In addition, to better take into account the dependence structure of precipitation data and capture the possible isotropic, we choose a flexible parametric semi-variogramν belonging to the power class models such that νh=h/τκ, with τ>0, h∈ℝ2, κ∈0,2. We check the adequacy of the fit of the model to the data using the extremogram and the variogram. Figures 5(a) and 5(b) show, respectively, the good quality of the fit of the dependence measures such as the semi-variogram and extremogram. In Figure 5(a), the points in gray represent the calculated empirical extremogram and the blue curve is the fitted empirical extremogram, while the red curve represents our fitted dependence model. The red curve in Figure 5(b) is the variogram model fitted to the data as a function of distance.Figure 5
Isotropic extremogramπh=limq⟶1PZs+h>uqs+h|Zs>uqs,ℓZ>u (a) and variogram (b) estimated for risk functional ℓZ=maxs∈SZs of distance between locations s and s+h. Blue line: empirical extremogram; red line: fitted dependence model.
(a)(b)It is noted that extreme precipitation is spatially correlated for distances of the order of 200 km. In other words, an extreme precipitation event is likely to affect an area within a 200 km radius. This spatial dependence decreases gradually as the distance increases before stabilizing for example for a value ofπh=0.0625. We find the very localized nature of extreme precipitation for long distances. This is generally a property present in rainfall data [32].


### 3.4. Nonstationary Return Level Results

After estimating the marginal parameters and those of the dependence structure, we now evaluate the nonstationary return levels. The return levelsx50 and x100 are first estimated (see Figure 6) using the (6) and (7). Then, they are spatially interpolated at points where we have no observations using the spatially extrapolated trend function and dependence structure estimated from the observed data.Figure 6
Maps of the nonstationary 50-years return levels (a) and 100-years return level (b) obtained by extrapolating the information from the log-linear tail trend function and the dependence structure.
(a)(b)For a given return periodm, we first compute the return levels zm on each point of grid from the latent spatial process Z using the estimated spatial model. The nonstationary return levels xm is deducted using the (8) of Proposition 2. We repeat this operation for different values of the return period m and we deduce on each point of the grid, the associated return level (Figure 7). We can then produce maps of the return levels obtained from the dependence structure and the extrapolation of the trend function for a future period m. Thus, the extreme precipitation, likely to be observed on average at least once every 50 years (resp. 100 years), will be particularly intense in the Sudanian and Sudano-Sahelian zone and less intense in the Sahelian zone, with a potentially quite strong spatial dependence within a radius of 200 km. The southwest and eastern regions of the country will be most affected by extreme precipitations. The results of the return level for a return period m=50 and 100 years are displayed in Figure 6.Figure 7
Trends of stationary and nonstationary return levels at reference stations computed with a log-linear tail trend function.


## 3.1. Data Set Analysis

This study uses time series of daily precipitation measurements from 1957 to 2016 provided by ten synoptic stations extracted from the Burkina Faso climatological database. These stations have been selected to ensure good spatial uniformity and representativeness of different climatic regimes and data quality. Figure1 gives the spatial distribution of selected stations in our study area. In order to limit the problems related to seasonal rainfall cycles on each station, we worked from the sub-series corresponding to rainy days. The period from May to October was chosen because it is during this period that the most rainfall is recorded in the country. Thus, a series of 60 years (1957-2016) by 184 days (May to October) is extracted to constitute the time series of daily rainfall. On these time series, we apply a run declustering procedure with a daily step (r = 1 day) to identify the groups of approximately independent extreme observations within the sample in order to avoid short-term dependencies in the time series.Figure 1
Spatial distribution of synoptic sampling stations associated with climatic regimes; source: National Meteorology Agency of Burkina Faso.


## 3.2. Marginal Characteristics

The first line of Figure2 shows the maps of the parameters an and bn estimated on our selected stations. Map (a) gives us the variability parameter of the precipitation distribution and shows a very high variability of extreme precipitation. While map (b) gives us the estimated level or threshold for each station using equation (12) and we see that the level is higher in the Sudanian climate and lower in the Sahelian climate.In practice, we deduce from the generalized ℓ-Pareto model fitted to the data, the marginal tail behavior of the survival distribution F¯Zs,γ^s,σ^s for any s∈S from the (16) given a sufficiently large uq>0 threshold.(16)PZs−uqs≥z|Zs≥uqs≈F¯Zs,γ^s,σ^sz,z≥0,with σ^s=a^ns+γ^uqs−b^ns, s∈S, where a^n, b^n, and γ^ are the marginal parameters estimators described in Section 2.4.1. We can use it to check the quality of the marginal adjustment of the stationary process Z. Figure 3 shows the qq-plots of the local tail distribution due to two stations by climatic zone. Columns 1, 2, and 3 of Figure 3 represent the respective adjustments for the Sudanian, Sudano-Sahelian, and Sahelian climate. Globally, the fit of the marginal models seems convincing, as most of the observations remain within the confidence intervals.Figure 2
Estimated functionsan (a) and bn (b) of the generalized ℓ-Pareto process for modelling extreme precipitation, and θ for Log-linear (c) and simple linear (d) trends. Estimates are obtained by shifting the local empirical quantiles u0.80ℓZs by b˜^n=33.4mm.
(a)(b)(c)(d)Figure 3
q-q plot of the distributions of the local tails of six synoptic positions including two stations for each climate region. (a) Sudanian climate; (b) Sudano-Sahelian climate; and (c) Sahelian climate. The threshold is taken equal to bns for each station.
(a)(b)(c)The second row of Figure2 shows the maps of parameter θ estimated for the log-linear and simple linear trend functions. We can observe a strong spatial variability in the occurrence of extreme precipitation. Thus, we observe a decreasing trend in the frequency of extreme rainfall in Ouaga, Dori, and Bobo. In contrast, in Boromo, Ouahigouya, and Pô, extreme precipitation events will become more and more frequent. Figure 4 gives us the details of the adjustment of the tail trend function by a log-linear model.Figure 4
Local fitting of the log-linear trend modelcθ1 for some of the selected sites with two stations for each climate region.


## 3.3. Estimated Spatial Dependence Model

We use a spatial marginal model (15) to estimate the marginal parameters at the different grid points. This is to prepare the ground for the calculation of nonstationary return levels at locations where we have no observations.(17)ans=a0+a1Long×Lats,bns=b0+b1Lat2s+b2Long×Lats,θs=θ0+θ1Long+θ2Lat,γs=γ0.,s∈S,An alternative to the spatial model of marginal parameters is to consider small subsetsNs0⊂S, where Ns0 is a small neighborhood around s0. In what follows, these regional neighborhoods will be determined by a small number D0 of nearest stations from the site s0, thus we will write Ns0=Ns0,D0. In principle, the choice of D0 should be such that the spatial dependence structure and marginal parameters is approximately stationary within each selected neighborhood Ns0,D0 around s0∈S. Obviously, the choice of neighborhood is important; the assumed stationary marginal parameters could be a poor approximation for large neighborhoods (i.e., for large D0), while the simulation of the process could be cumbersome for small neighborhoods characterized by a small D0 number of stations. We obtain the relatively homogeneous, nonoverlapping subregions using the k-means clustering method centered on the reference stations [30, 31]. This method is extensively used because it is computationally simple and produces accurate results, compared to other more complex clustering methods. The longitude and latitude of the grid points were used as input variables in the k-means clustering algorithm to form the ten clusters centered on the reference stations.In addition, to better take into account the dependence structure of precipitation data and capture the possible isotropic, we choose a flexible parametric semi-variogramν belonging to the power class models such that νh=h/τκ, with τ>0, h∈ℝ2, κ∈0,2. We check the adequacy of the fit of the model to the data using the extremogram and the variogram. Figures 5(a) and 5(b) show, respectively, the good quality of the fit of the dependence measures such as the semi-variogram and extremogram. In Figure 5(a), the points in gray represent the calculated empirical extremogram and the blue curve is the fitted empirical extremogram, while the red curve represents our fitted dependence model. The red curve in Figure 5(b) is the variogram model fitted to the data as a function of distance.Figure 5
Isotropic extremogramπh=limq⟶1PZs+h>uqs+h|Zs>uqs,ℓZ>u (a) and variogram (b) estimated for risk functional ℓZ=maxs∈SZs of distance between locations s and s+h. Blue line: empirical extremogram; red line: fitted dependence model.
(a)(b)It is noted that extreme precipitation is spatially correlated for distances of the order of 200 km. In other words, an extreme precipitation event is likely to affect an area within a 200 km radius. This spatial dependence decreases gradually as the distance increases before stabilizing for example for a value ofπh=0.0625. We find the very localized nature of extreme precipitation for long distances. This is generally a property present in rainfall data [32].


## 3.4. Nonstationary Return Level Results

After estimating the marginal parameters and those of the dependence structure, we now evaluate the nonstationary return levels. The return levelsx50 and x100 are first estimated (see Figure 6) using the (6) and (7). Then, they are spatially interpolated at points where we have no observations using the spatially extrapolated trend function and dependence structure estimated from the observed data.Figure 6
Maps of the nonstationary 50-years return levels (a) and 100-years return level (b) obtained by extrapolating the information from the log-linear tail trend function and the dependence structure.
(a)(b)For a given return periodm, we first compute the return levels zm on each point of grid from the latent spatial process Z using the estimated spatial model. The nonstationary return levels xm is deducted using the (8) of Proposition 2. We repeat this operation for different values of the return period m and we deduce on each point of the grid, the associated return level (Figure 7). We can then produce maps of the return levels obtained from the dependence structure and the extrapolation of the trend function for a future period m. Thus, the extreme precipitation, likely to be observed on average at least once every 50 years (resp. 100 years), will be particularly intense in the Sudanian and Sudano-Sahelian zone and less intense in the Sahelian zone, with a potentially quite strong spatial dependence within a radius of 200 km. The southwest and eastern regions of the country will be most affected by extreme precipitations. The results of the return level for a return period m=50 and 100 years are displayed in Figure 6.Figure 7
Trends of stationary and nonstationary return levels at reference stations computed with a log-linear tail trend function.


## 4. Discussion and Conclusion

In summary, modeling the spatiotemporal trends of extreme precipitation in the sub-Saharan Africa is a challenge and a necessity. Climate models do not often agree on this issue. Advanced statistical methods remain a fairly powerful alternative to obtain pertinent information on possible spatiotemporal changes in the underlying precipitation distribution. A frequently used approach to dealing with nonstationarity in trend analysis of extreme precipitation is to include covariates, such as space and time, in the parameters of the underlying distribution [9, 10, 12, 13]. These approaches are built under certain important assumptions that the underlying process that generates the data is locally stationary and that the observations generated can be considered as independent approximations and generally do not take into account spatial dependence.In this study, we proposed a new flexible methodology for trend detection in extremes of spatiotemporal processes. This is capable of capturing marginal nonstationarities and the dependence structure between margins. Marginal nonstationarity is taken into account by a trend function called the skedasis function which models the frequency of extreme events jointly in space and time. While the spatial dependence structure is governed by a latent spatial process and modeled using generalizedℓ-Pareto processes. We calculated the nonstationary return levels of precipitation in Burkina Faso and exhibit some regions in which there may be a significant increase or decrease of extreme precipitation. We showed that these extreme rainfall events are spatially correlated around a radius of 200 km. This spatial dependence decreases progressively as the distance increases. In sum, we set up a nonstationary stochastic generator of extreme rainfall in Burkina Faso.Previous studies conducted in the sub-Saharan African region have shown trends in extreme precipitation similar to those we obtained using Burkina Faso data set. According to the latest publications of the National Meteorological Agency of Burkina Faso, rainfall has started to move back toward a wet period since the end of the 1980s and the decade 1990-2000. This return of rain is more related to an increase in the frequency of heavy rainfall events than to an increase in the number of rainy days. Our results are also in agreement with the results found by [2–4] in the region which indicate a similar intensification of rainfall events. Based on a point GEV model with time-dependent parameters, these authors show that intense precipitation events will become more frequently occurring in the region. Nkrumah et al. [33] and Mouhamed et al. [34] observed that the average annual precipitation has a decreasing trend in the same region due to less frequent but more intense precipitation during the last decade.Our method generates extreme events consistent enough with the data already observed. The return levels computed can be improved by introducing also a nonstationarity in the structure of spatial dependence.



---

*Source: 2608270-2022-06-08.xml*