{'fraction_non_alphanumeric': 0.044865813146205816, 'fraction_numerical': 0.010749878690605053, 'mean_word_length': 4.2938154514918, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling space-time trends and dependence in extreme precipitations of Burkina Faso by the approach of the Peaks-Over-Threshold

Modeling extremes of climate variables in the framework of climate change is a particularly difficult task, since it implies taking into account spatio-temporal nonstationarities. In this paper, we propose a new method for estimating extreme precipitation at the points where we have not observations using information from marginal distributions and dependence structure. To reach this goal we combine two statistical approaches of extreme values theory allowing on the one hand to control temporal and spatial non-stationarities via a tail trend function with a spatio-temporal structure in the marginal distributions and by modeling on the other hand the dependence structure by a latent spatial process using generalized `-Pareto processes. This new methodology for trend analysis of extreme events is applied to rainfall data from Burkina Faso. We show that extreme precipitation is spatially and temporally correlated for distances of approximately 200 km. Locally, extreme rainfall has more of an upward than downward trend.


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 non-stationarities. 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(IPCC , 2018(IPCC ,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. From a statistical point of view, the hypothesis of stationarity becomes untenable and debatable, and thus raises the question of detecting trends in extreme events, which has been the subject of much attention over the last 30 years. The classical Theory of Extreme Values extended both to non-stationary and non-independent observations provides a rigorous mathematical framework to deal with this question of trend detection in extremes ( [13], [16], [19], [5]).
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( [8], [20], [21], [7], [24]). Although it is relatively simple to construct non-stationary 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 ( [11], [17], [4], [3,2]). These spatial models have been readapted to handle non-stationarities induced by global warming ( [18], [6]). Although attractive, these models are expensive in computing time when tending towards large dimensions and are not adapted for modeling threshold exceedances.
The threshold exceedance approach introduced by [1] and [22] was first extended to multivariate environments ( [23]) before being generalized to functional data ( [15], [12], [9,10]) to give birth to the family of generalized -Pareto processes. This family of processes are better adapted to model the spatial dependence structure of exceedances 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 non-stationarities in marginal distributions, while taking into account the spatial dependence structure. To reach this goal we combine two statistical approaches of extreme value theory allowing on the one hand to control temporal and spatial non-stationarities via a tail trend function with a spatio-temporal structure in marginal distributions( [13], [16], [5]) and by modeling on the other hand the dependence structure by a hidden stationary auxiliary process using generalized -Pareto processes [10].
The article is organized as follows: section 2 provides a brief overview on generalized -Pareto processes. Section 3 details the methodology implemented and the methods used to estimate the parameters. In section 4 we present our main results, and in particular we compare the stationary and non-stationary return levels computed from the developed approaches.

Overview on generalized -Pareto process
Let S be a compact subset of R d and T be a compact subset of R + denoting the spatial and temporal domain respectively. We note by C(S) the set of continuous real functions on S with the uniform norm . ∞ ; C + (S) its restriction to non-negative functions deprived of the null function and M {C + (S)} the class of measures associated with C + (S).
Let {Z(s), s ∈ S} be a latent spatial stochastic process indexed by s ∈ S ⊂ R d with sample paths in C(S) of continuous marginal distribution F Z and common right endpoint z F . As in ( [15], [10]), we assume that the stationary process Z is a general functional regular variation, noted Z ∈ GRV (γ, a n , b n , Λ), i.e that there exists suitable sequences of continuous functions a n : S −→ R + , b n : S −→ R and γ ∈ R such that where Λ is a non-zero measure in M {C + (S)} and homogeneous of order −1, with Λ (tA) = t −1 Λ (A) for any positive real t > 0 and Borel set A ⊂ C + (S).
In the multivariate and spatial framework a threshold exceedances for a random function Z = {Z(s), s ∈ S} is defined by Dombry and Ribatet [12] to be an event of the form { (Z) > u} for some u ≥ 0, where : C(S) −→ R + 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 s 0 ∈ S. Moreover, as in ( [10], [14]) we assume that there exists a continuous and positive real function A Z such that the sequence a n and the risk function satisfy the following asymptotic decomposition lim n→∞ sup s∈S a n (s) (a n ) − A Z (s) = 0, i.e a n (s) ≈ (a n )A Z (s), n → ∞ The marginal distributions of Z are supposed to belong to a location-scale family, thus ensuring a constant γ ∈ R shape parameter over any S, i.e. that there exist a H distribution function such that where A Z : S → (0; ∞[ verifies asymptotic decomposition(2) and B Z : S → R are continuous functions. In particular H must belong to the domain of attraction of a generalized extreme value (GEV) distribution of index γ ∈ R, i.e., F Z ∈ D(G γ ). In other words, there are appropriate normalization real sequencesã n > 0,b n ∈ R such that lim n→∞ H n ã n z +b n = G γ (z). The condition of the max-domain of attraction is equivalent to Under these conditions and the assumption of the general functional regular variation the suitable sequence of continuous functions a n (s) and b n (s), satisfy: Functions A Z (s) and B Z (s) as well as the spatial dependence structure of Z are assumed to belong to a family of parametric functions Under minimal assumptions on the risk function and assumptions (1, 2) of Z, the conditional distribution of -exceedance for some threshold u ≥ 0 of the process (Z − b n )/ (a n ) can be approximated by a generalized -Pareto process, for n large enough( [10]) where z = max (z, W is a non-degenerate stochastic process over S and belongs to the family of generalized -Pareto processes with tail index γ, zero location, scaling function A Z and limit measure Λ. Specifically, a generalized -Pareto process W associated to the limit measure Λ and tail index γ ∈ R is a stochastic process taking values in {z ∈ C + (S) : {(z − b)/ (a)} ≥ 0} and defined by: where a = (a)A Z > 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 Λ. For the modelling of spatial dependence structure of latent process, we use Y whose margins are in the Frechet max-domain of attraction with tail index γ = 1, as the process of reference. A pseudo-polar decomposition of Y in (7) leads to the following formulation ( [12], [9]) 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) : y 1 = 1} whose probability measure is characterized by limit measure Λ. More details on generalized -Pareto processes can be found in ( [15], [12], [9,10]).

Space-time trends detection
Let X = {X t (s), s ∈ S, t ∈ T } be a continuous non-stationary space-time stochastic process with sample paths in the family of continuous functions C(S × T ) equipped with the uniform norm . ∞ , where S × T ⊂ R d × R + and C + (S × T ) its restriction to non-negative functions deprived of the null function. In practice X is observed at each stations s 1 , · · · , s m and at given dates t = 1, · · · , n. Let F t,s be the continuous univariate marginal distribution with a common right endpoint x F and Z = {Z(s), s ∈ S} an unobserved latent spatial stochastic process with sample paths in C(S × T ) satisfying the properties described in the section (2) and the proportional tail condition such that where c θ : [0, 1] × S −→ (0, ∞[ is an assumed continuous and positive function depending on a parameter vector θ ∈ Θ ⊂ R, called tail scale function or skedasis function ( [13], [15,16], [5]). In added we assume that the continuous marginal distributions F Z have a some common right endpoint as F t,s . The skedasis function describes the evolution of extreme events jointly in space and time. The tail trend function c θ is designed to ensure uniqueness of c θ such that: In the framework of the model (9), Mefleh et al. [19] shows that the empirical point measure converges in distribution in the space of point measure In such a 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 γ. Several models are eligible to model the function c θ , but 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. We are interested in monotonic log-linear and simple linear trend models of c θ . We introduce a spatial structure in the tail function by letting the parameter θ evolve as a function of geographic coordinates (longitude, latitude), that is, The parameter vector θ = (θ 0 , θ 1 , θ 2 ) is then estimated using the maximum likelihood method and multiple linear regression.

Marginal model of latent process
The marginal parameters γ, a n , A Z , b n and B Z of the latent spatial process Z are estimated under the constraint that u n = (b n ) + u (a n ), following the modeling assumptions described in section (2). For simplicity purpose we choose u = 0 so that -exceedances of Z are defined as events for which (Z) ≥ (b n ). In this case a natural choice of u n = (b n ) is a high quantile of the variable (Z), for example u n = q 0.95 { (Z)}. In general, a parametric model may be necessary for a n and b n , as in Engelke et al. [14], but we consider in this work that a n (s j ) = a j > 0 and b n (s j ) = b j ∈ R for any j = 1, · · · , m. The threshold stability of the generalized Pareto distributions does not allow us to identify the b n function without additional assumptions, so under the assumption that whereb n = (b n ) = u n (see 2 and 5) and u q {Z(s)} is an empirical quantile of the order q of the -exceedances at each location s, where q is chosen such that (b n ) =b n in order to impose the identifiability of the parameters. Thus, the tail index γ ∈ R and the scale parameters a n (s) > 0 are estimated by maximizing the independent log-likelihood; that is The estimate of parameters A Z and B Z are deducted from (5) under the assumptions of parameter identifiability, i.e, (A Z ) = 1 and (B Z ) = 0, which implies that (a n ) =ã n . Thanks to the equation (9) and the convergence of Z exceedances to a GPD distribution we deduce a sample of {Z t (s)} from observations of {X t (s)} ( [16] and [5]) in the following manner: where γ, ã n , b n and c θ are respectively consistent estimators of γ,ã n ,b n and c θ described in section (3.1).

Spatial extremal dependence
After removing the non-stationarity the modeling is focused on the evaluation of the remaining extreme spatial dependence structure in Z. Thus we approach the limit distribution of -exceedances of Z by a generalized -Pareto process [10]. In the framework of the Brown-Resnick model, the angular component of Q is a log-Gaussian random function whose underlying Gaussian process has stationary increments, which allows in particular to use classical dependence structures from the geostatistical literature to characterize the spatial dependence. In order to better capture this dependence structure, we use a flexible parametric semi-variogram belonging to the class of power semi-variograms. We estimate the parameter vector θ Λ of the dependence structure using the gradient score method or the censored likelihood method [9].

Non-stationary return period and return level
The concept of return period and level becomes very ambiguous when we leave the stationary context to the non-stationary framework. In this paper, we have chosen to follow return period based approaches,i.e., Expected Number of Exceedances(ENE) [21] and Expected Waiting Time (EWT) ( [20], [7]). In the ENE approach, the number of times the variable X t (s), s ∈ S exceeds the return level value x m in m years is defined by N m = nxm t=1 1 {Xt(s)>xm,s∈S} under non stationary context. The return level x m can be defined as the value for which the expected number of events exceeding x m in m years equals to one, i.e., the return level x m is the solution of the following equation: where n x is the number of days in the year and θ t (s) the vector of time-dependent marginal parameters or other covariates. Parey et al. [21] 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. [20], and then derived by Salas and Obeysekera [24] using a geometric distribution with time-varying parameters. Under non stationary conditions, the distribution describing waiting time Y before the first occurrence of an event exceeding the return level x m is (1 − p i,s ) , y = 1, 2, · · · , y max and s ∈ S, where variable Y is the day of the first occurrence of an event exceeding the quantile x m , p t,s = 1 − F t,s (x m | θ t (s)) is daily exceedance probability varying with time step t. y max is the time when the daily exceedance probability p t,s is equal to 1 for an increasing-trend series or is equal to 0 for a decreasing-trend series. Reused and simplified by [7], the EWT approach defines the m-year return level x m , as the value for which the expected waiting time until an exceedance of this level is m years, i.e, x m is the solution of the equation: Using the relationship (9), P (X t (s) > x m ) can be rewritten c θ ( t n , s)P (Z(s) > u)P (Z(s) > x m | Z(s) > u) for u < x m and we obtain the following results Proposition 3.1 Let {X t (s), s ∈ S, t ∈ T } be a non-stationary stochastic process defined on a region S ⊂ R d and {Z(s), s ∈ S} a latent spatial process satisfying the equation 9.
Given a return period m and threshold u < x m , the return level x m for all s ∈ S is a solution of the following two equations: i) Return period as expected number of events ii) Return period as expected waiting time whereF Z,s is a survival of generalized Pareto distribution of -exceedances at position s; φ u (s) = P (Z(s) > u) is the probability of exceedances, c θ is the tail trend function and n x is the number of days in the year.
The return levels derived from the equations (18 & 19) are evaluated by numerical algorithms taking into account the information from the extrapolation of the trend function c θ on the one hand and the spatial dependence structure of the latent stationary process Z on the other hand. To derive the return level at points where we have not observations, we assume, the spatial process Z is locally stationary in space. Under these new assumptions and using the knowledge of the spatial dependence structure of the latent process Z and the equation (14), we propose the following result: Proposition 3.2 Let {X t (s), s ∈ S, t ∈ T } be a non-stationary spatio-temporal stochastic process defined on a region S ⊂ R d and {Z(s), s ∈ S} a latent spatial process. The nonstationary return level x m (s), s ∈ S of the non-stationary process X is deduced from the return level z m of the latent spatial process Z such that: where ã n , b n c θ and γ are the respective estimators ofã n ,b n , c θ and γ. t m = 1 + n x m n with n x the number of days in the year and n the size of the sample observed.
This result is a consequence of the equation (14) and will be used to derive the nonstationary return levels of process X, using the z m computed from the estimated parameters. Thus, to calculate the z m at grid points where we have not observations, we use a spatial model from marginal parameters to build a spatial map on the domain of study based on covariates such as longitude and latitude.

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. The figure1 gives the spatial distribution of the synoptic 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 during this period that the most rainfall is recorded in Burkina Faso. Thus, a serie of 60 by 184 days 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.

Marginal characteristics
The scale a n (s), location b n (s) and parameters of log-linear and linear trend function c θ are estimated for any s ∈ S and are shown respectively on maps (a), (b), (c) and (d) of figure (2). In practice we deduce from the generalized -Pareto model fitted to the data,  the marginal tail behavior of the survival distributionF Z, γ(s), σ(s) for any s ∈ S from the equation (21) given a sufficiently large u q > 0 threshold.
with σ(s) = a n (s) + γ u q (s) − b n (s) , s ∈ S, where a n , b n , and γ are the marginal parameters estimators described in section (3.1). We can use it to check the quality of the marginal adjustment of the stationary process Z. Figure 3 shows the QQplots 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 zones. Globally, the fit of the marginal models seems convincing, as most of the observations remain within the confidence intervals. Furthermore, the sign of the θ parameter of the trend function tells us that the frequencies of extreme precipitation are locally quite variable throughout the study region. They have an increasing trend (θ > 0) in areas such as Ouahigouya, Bogande, Boromo, Gaoua, and Po. On the other hand, at the Ouagadougou and Fada stations, extreme rainfall frequencies tend to decrease (θ < 0). Figure 4 gives us the details of the adjustment of the tail trend function by a log-linear model.

Estimated spatial dependence model
We use a spatial model of the marginal parameters dependent on significant covariates to obtain the marginal parameters at the different grid points using a generalized linear regression model, in order to prepare the ground for calculations of non-stationary return levels at locations where we have no observations.
where N s 0 is a small neighborhood around s 0 . In what follows, these regional neighborhoods will be determined by a small number, D 0 of nearest stations from the site s 0 , thus we will write N s 0 = N s 0 ,D 0 . Obviously, the choice of neighbourhood is important; the assumed stationary marginal parameters could be a poor approximation for large neighbourhoods (i.e., for large D 0 ), while the simulation of the process could be cumbersome for small neighbourhoods characterized by a small D 0 number of stations. In principle, the choice of D 0 should be such that the spatial dependence structure and marginal parameters is approximately stationary within each selected neighborhood N s 0 ,D 0 around s 0 ∈ G. We obtain the relatively homogeneous, non-overlapping sub-regions using the k-means clustering method centered on the reference stations. 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 ∈ R 2 , κ ∈ ]0, 2]. We check the adequacy of the fit of the model to the data using the extremogram and the variogram. Figures (5a) and (5b) show respectively the good quality of the fit of the dependence measures such as the semi-variogram and extremogram. In Figure(5a), 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 5b is the variogram model fitted to the data as a function of distance. It is noted that extreme precipitation is spatially correlated for distances of the order of 200km. This spatial dependence decreases gradually as the distance increases before stabilizing for example for a value of π(h) = 0.0625. This reflects the very localized nature of extreme rainfall and it comforts us in our analysis because this property is generally present in rainfall data.

Non-stationary return level results
The return levels x 50 and x 100 are first estimated punctually (see Figure7) using the equations (18 & 19) before being spatially interpolated to the points where we have no observations using the dependence structure estimated from the observed data and parameters of trend function spatially extrapolated to the section ( period m, we first compute the return levels z m on each point of grid from the latent spatial process Z using the estimated spatial model, before deriving the non-stationary return levels x m using the equation 20 of proposition 3.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 6). 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 Figures (7).

Conclusion
In this study we proposed a new flexible methodology for trend detection in the extremes capable of capturing marginal non-stationarities and the dependence structure between margins using generalized -Pareto processes. We computed the non-stationary return  Figure 7: Maps of the non-stationary 50-years return levels(left) and 100-years return level (right) obtained by extrapolating the information from the log-linear tail trend function and the dependence structure.
levels of rainfall in Burkina Faso at points where we have no observations and we showed that these extreme rainfall events are spatially correlated around a radius of 200km. This spatial dependence decreases progressively as the distance increases. In sum, we set up a non-stationary stochastic generator of extreme rainfall in Burkina Faso. with σ(s) = a n (s) + γ(s)(u(s) − b n (s)) ii) The result (19) is shown in a similar way

B Asymptotic distribution of -Pareto Process
For a threshold vector u = (u 1 , · · · , u d ) ∈ R d + and for any regularly varying stochastic process Z ∈ GRV {γ, a n , b n , Λ}, the density function of the -excess of a -Pareto process is obtained by renormalizing suitably the intensity function λ found by taking the partial derivatives of the previously defined measure Λ (1) by Λ (A (u)): where while λ is the intensity function and A (u) = z ∈ R d + : (z/u) ≥ 1 the region of exceedances. In order to make inferences and to model the dependence structure of -Pareto processes, we focus on the Brown-Resnick model for which the formulas of Λ and λ are available. The d-dimensional intensity function of the Brown-Resnick model is given: wherez i = log z i z 1 + γ(s i − s 1 ) and Σ the covariance matrix. To better capture the possible dependence structure we use anisotropic semi-variograms whose parameters change with time and other covariates. The different parameters are estimated using the gradient scoring rule method [9]. For any z ∈ A (u), the log-density function is given by: where w is a differentiable weighting function that vanishing on the boundaries of A (u) and θ is an element in the Θ ⊂ R p parameter space.