**1. Introduction**

In this paper, we continue the research we started in [1,2]. We develop the mathematical models for statistical regularities in precipitation proposed in the papers mentioned above. We consider the models for the statistical regularities in the duration of a wet period, maximum daily precipitation within a wet period and total precipitation volume per wet period. The base for the models is the generalized negative binomial (GNB) introduced in the recent paper [3]. The GNB distribution is a mixed Poisson distribution, the mixing distribution being generalized gamma (GG). The results of fitting the GNB distribution to real data are presented and demonstrate excellent concordance of the GNB model with the empirical distribution of the duration of wet periods measured in days. Based on this GNB model, asymptotic approximations are proposed for the distributions of the maximum daily precipitation volume within a wet period and of the total precipitation volume for a wet period. The asymptotic distribution of the maximum daily precipitation volume within a wet period turns out to be a tempered scale mixture of the gamma distribution in which the scale factor has the Weibull distribution, whereas the asymptotic approximation for the total precipitation volume for a wet period turns out to be the GG distribution. These asymptotic approximations are deduced using limit theorems for statistics constructed from samples with random sizes having the GNB distribution. The bounds for the accuracy of the proposed approximations are discussed theoretically and illustrated statistically. The proposed approximations appear to be very accurate. Based on these models, two approaches are proposed to the definition of abnormally extremal precipitation. These approaches can be regarded as a further development of those proposed in [2].

The importance of the problem of modeling statistical regularities in extreme precipitation is indisputable. Understanding climate variability and trends at relatively large time horizons is of crucial importance for long-range business, say, agricultural projects and forecasting of risks of water floods, dry spells and other natural disasters. Modeling regularities and trends in heavy and extreme daily precipitation is important for understanding climate variability and change at relatively small or medium time horizons. However, these models are much more uncertain as compared to those derived for mean precipitation or total precipitation during a wet period. In [4], a detailed review of this phenomenon is presented and it is noted that, at least for the European continent, most results hint at a growing intensity of heavy precipitation over the last decades.

In [2], we proposed a rather reasonable approach to the unambiguous (algorithmic) determination of extreme or abnormally heavy total precipitation for a wet period. This approach was based on the NB model for the duration of wet periods measured in days, and, as a consequence, on the distribution of the total precipitation volume during a wet period. This approach has some advantages. First, estimates of the parameters of the total precipitation are weakly affected by the accuracy of the daily records and are less sensitive to missing values. Second, the corresponding mathematical models are theoretically based on limit theorems of probability theorems that yield unambiguous asymptotic approximations, which are used as adequate mathematical models. Third, this approach gives an unambiguous algorithm for the determination of extreme or abnormally heavy total precipitation that does not involve statistical significance problems owing to the low occurrence of such (relatively rare) events.

The problem of the construction of a statistical test for the precipitation volume to be abnormally large can be mathematically formalized as follows. Let *m* 2 be a natural number and consider a sample of *m* positive observations *X*1, *X*2, ... , *Xm*. With finite *m*, among *Xi*'s there is always an extreme observation, say, *X*1, such that *X*1 *Xi*, *i* = 1, 2, ... , *m*. Two cases are possible: (i) *X*1 is a 'typical' observation and its extreme character is conditioned by purely stochastic circumstances (there must be an extreme observation within a finite homogeneous sample) and (ii) *X*1 is abnormally large so that it is an 'outlier' and its extreme character is due to some exogenous factors.

To construct a test for distinguishing between these two cases for abnormally extreme daily precipitation, we use the fact that the distribution of the maximum daily precipitation per wet period is a tempered scale mixture of the gamma distribution in which the scale factor has the Weibull distribution. According to this model, a daily precipitation volume is considered to be abnormally extremal if it exceeds a certain (pre-defined) quantile of this distribution.

As regards testing for anomalous extremeness of total precipitation volume during a wet period, we use the GG distribution as the model of statistical regularities of its behavior. The theoretical grounds for this model are provided by the law of large numbers for random sums in which the number of summands has the GNB distribution. It turns out that, as compared to the ordinary negative binomial (NB) model (see [2]), the additional exponent power parameter in the corresponding GG distribution matches slow global climate trends. Hence, the hypothesis that the total precipitation volume during a certain wet period is abnormally large can be re-formulated as the homogeneity hypothesis of a sample from the GG distribution. Two equivalent tests are proposed for testing this hypothesis. One of them is based on the beta distribution whereas the second is based on the

Snedecor–Fisher distribution. Both of these tests deal with the relative contribution of the total precipitation volume for a wet period to the considered set (sample) of successive wet periods. Within the second approach, it is possible to introduce the notions of relatively abnormal and absolutely abnormal precipitation volumes. These tests are scale-free and depend only on the easily estimated shape parameter of the GNB distribution and the time-scale parameter determining the denominator in the fractional contribution of a wet period under consideration. The tests appeared to be applicable not only to total precipitation volumes over wet periods but also to the precipitation intensities (the ratios of total precipitation volumes per wet periods to the durations of the corresponding wet periods measured in days).

#### **2. Generalized Negative Binomial Model for the Duration of Wet Periods**

The main results of this paper strongly rely on the GNB model, a wide and flexible family of discrete distributions that are mixed Poisson laws with the mixing GG distribution. Namely, we say that a random variable *Nr*,*γ*,*<sup>μ</sup>* (*r* > 0, *γ* ∈ R and *μ* > 0) has the generalized negative binomial distribution, if

$$\mathbb{P}(N\_{r,\gamma\_r\mu} = k) = \frac{1}{k!} \int\_0^\infty e^{-z} z^k g^\*(z; r, \gamma\_r \mu) dz, \quad k = 0, 1, 2... \tag{1}$$

where *g*<sup>∗</sup>(*<sup>z</sup>*;*r*, *γ*, *μ*) is the density of GG distribution:

$$\log^\*(\mathbf{x}; r, \gamma, \mu) = \frac{|\gamma| \mu^r}{\Gamma(r)} \mathbf{x}^{\gamma r - 1} \mathbf{e}^{-\mu \mathbf{x}^\gamma}, \quad \mathbf{x} \gg \mathbf{0}, \tag{2}$$

with *γ* ∈ R, *μ* > 0, *r* > 0. The GNB distributions seem to be very promising in the statistical description of many real phenomena, being very convenient and almost universal models.

It is necessary to explain why this combination of the mixed and mixing distributions is considered. First of all, the Poisson kernel is used as mixed for the following reasons. Pure Poisson processes can be regarded as the best models of stationary (time-homogeneous) chaotic flows of events [5]. Recall that the attractiveness of a Poisson process as a model of homogeneous discrete stochastic chaos is due to at least two circumstances. First, Poisson processes are point processes characterized by the time intervals between successive points that are independent random variables (r.v.'s) with one and the same exponential distribution, and, as is well known, the exponential distribution possesses the maximum differential entropy among all absolutely continuous distributions concentrated on the nonnegative half-line with finite expectations, whereas the entropy is a natural and convenient measure of uncertainty. Second, the points forming the Poisson process are uniformly distributed along the time axis in the sense that for any finite time interval [*<sup>t</sup>*1, *<sup>t</sup>*2], *t*1 < *t*2, the conditional joint distribution of the points of the Poisson process that fall into the interval [*<sup>t</sup>*1, *t*2] under the condition that the number of such points is fixed and equals, say, *n*, coincides with the joint distribution of the order statistics constructed from an independent sample of size *n* from the uniform distribution on [*<sup>t</sup>*1, *<sup>t</sup>*2], whereas the uniform distribution possesses the maximum differential entropy among all absolutely continuous distributions concentrated on finite intervals and very well corresponds to the conventional impression of an absolutely unpredictable random variable (see, e g., [5,6]). However, in actual practice, as a rule, the parameters of the chaotic stochastic processes are influenced by poorly predictable «extrinsic» factors, which can be regarded as stochastic so that most reasonable probabilistic models of non-stationary (time-non-homogeneous) chaotic point processes are doubly stochastic Poisson processes, also called Cox processes (see, e.g., [5,7,8]). These processes are defined as Poisson processes with stochastic intensities. Such processes proved to be adequate models in insurance [5,7,8], financial mathematics [9], physics [10] and many other fields. Their one-dimensional distributions are mixed Poisson.

In order to have a flexible model of a mixing distribution that is "responsible" for the description of statistical regularities of the manifestation of external stochastic factors, we sugges<sup>t</sup> to use the GG distributions defined by the density (2). The class of GG distributions was first described as a unitary

family in 1962 by E. Stacy [11] as the class of probability distributions simultaneously containing both Weibull and gamma distributions. The family of GG distributions contains practically all the most popular absolutely continuous distributions concentrated on the non-negative half-line. In particular, the family of GG distributions contains:

	- ◦ The exponential distribution ( *γ* = 1, *r* = 1),
	- ◦ The Erlang distribution ( *γ* = 1, *r* ∈ N),
	- ◦ The chi-square distribution ( *γ* = 1, *μ* = 1 2 );
	- ◦ The Lévy distribution (the one-sided stable distribution with the characteristic exponent 1 2 – the distribution of the first hit time of the unit level by the Brownian motion) ( *γ* = −1, *r* = 1 2 );

and other laws. The limit point of the class of GG-distributions is

• The log-normal distribution (*r* → ∞).

GG distributions are widely applied in many practical problems. There are dozens of papers dealing with the application of GG-distributions as models of regularities observed in practice. Apparently, the popularity of GG-distributions is due to the fact that most of them can serve as adequate asymptotic approximations, since all the representatives of the class of GG-distributions listed above appear as limit laws in various limit theorems of probability theory in rather simple limit schemes. Below we will formulate a general limit theorem (an analog of the law of large numbers) for random sums of independent r.v.'s in which the GG-distributions are limit laws. It is worth noting that the GG distribution and its limit cases give a general form of the exponential distribution of rank 1 for the scale parameter.

In [1], the data registered in so climatically different points as Potsdam (Brandenburg, Germany) and Elista (Kalmykia, Russia) was analyzed, and it was demonstrated that the fluctuations of the numbers of successive wet days with very high confidence fit the NB distribution with shape parameters *r* = 0.847 and *r* = 0.876, respectively. In the same paper, a schematic attempt was undertaken to explain this phenomenon by the fact that NB distributions can be represented as mixed Poisson laws with mixing gamma distributions whereas, as it already has been mentioned, the Poisson distribution is the best model for the discrete stochastic chaos and the mixing distribution accumulates the stochastic influence of factors that can be assumed exogenous with respect to the local system under consideration.

The NB distributions are special cases of the GNB distributions. This family of discrete distributions is very wide and embraces Poisson distributions (as limit points corresponding to a degenerate mixing distribution), NB (Polya) distributions including geometric distributions

(corresponding to the gamma mixing distribution, see [12]), Sichel distributions (corresponding to the inverse gamma mixing distributions, see [13,14]), Weibull–Poisson distributions (corresponding to the Weibull mixing distributions, see [15]) and many other types supplying descriptive statistics with many flexible models. More examples of mixed Poisson laws can be found in [8,16].

It is quite natural to expect that, having introduced one more free parameter into the pure negative binomial model, namely, the power parameter in the exponent of the original gamma mixing distribution, instead of the negative binomial model one might obtain a more flexible GNB model that provides an even better fit with the statistical data of the durations of wet days. The analysis of the real data shows that this is indeed so.

In Figures 1 and 2 there are the histograms constructed from real data of 3323 wet periods in Potsdam and 2937 wet periods in Elista. On the same pictures, there are the graphs of the fitted NB distribution (that is, the GNB distribution with *γ* = 1) and the fitted GNB distribution with additionally adjusted scale and power parameters. For vividness, in the GNB model, the value of the shape parameter *r* was taken the same as that obtained for the NB model and equal to 0.876 for Elista and 0.847 for Potsdam. For the "fine tuning" of the GNB models with these fixed values of *r*, the minimization of the 1-norm of the difference between the histogram and the fitted GNB model was used. In Appendix A, the Algorithm A1 of for the computation of GNB probabilities by the minimization of the 1, 2 and ∞-norms of the difference between the histogram and the fitted GNB model is presented.

**Figure 1.** The histograms constructed from real data of 3320 wet periods in Potsdam and the fitted negative binomial (NB) and generalized negative binomial (GNB) models, 1-distance minimization.

The analytic and asymptotic properties of the GNB distributions were studied in [3]. In particular, it was shown in that paper that the GNB distribution with shape parameter and exponent power parameter less than one is actually mixed geometric. The mixed geometric distributions were introduced and studied in [17] (also see [15,18]). A mixed geometric distribution can be interpreted in terms of the Bernoulli trials as follows. First, as a result of some "preliminary" experiment the value of some r.v. taking values in [0, 1] is determined, which is then used as the probability of success in the sequence of Bernoulli trials in which the original "unconditional" mixed Poisson r.v. is nothing else than the "conditionally" geometrically distributed r.v. having the sense of the number of trials up to the first failure. This makes it possible to assume that the sequence of wet/dry days is not independent but is conditionally independent and the random probability of success is determined by some outer stochastic factors. As such, we can consider the seasonality or the type of the cause of a wet period. So, since the GG-distribution is a more general and, hence, a more flexible model than the "pure" gamma

distribution, there arises a hope that the GNB distribution could provide an even better goodness of fit to the statistical regularities in the duration of wet periods than the "pure" NB binomial distribution.

**Figure 2.** The histograms constructed from real data of 2937 wet periods in Elista and the fitted NB and GNB models, 1-distance minimization.

#### **3. Notation, Definitions and Mathematical Preliminaries**

In the paper, conventional notation is used. The symbols *d*= and =⇒ denote the coincidence of distributions and convergence in distribution, respectively.

In what follows, for brevity and convenience, the results will be presented in terms of r.v.'s with the corresponding distributions. It will be assumed that all the r.v.'s are defined on the same probability space (<sup>Ω</sup>, F, <sup>P</sup>).

An r.v. having the gamma distribution with shape parameter *r* > 0 and scale parameter *μ* > 0 will be denoted *Gr*,*μ*,

$$\mathbb{P}(G\_{r,\mu} < \mathbf{x}) = \int\_0^\mathbf{x} \mathbf{g}(z; r, \mu) dz,\text{ with }\mathbf{g}(\mathbf{x}; r, \mu) = \frac{\mu^r}{\Gamma(r)} \mathbf{x}^{r-1} e^{-\mu \mathbf{x}},\text{ }\mathbf{x} \gg 0,$$

where <sup>Γ</sup>(*r*) is Euler's gamma-function, <sup>Γ</sup>(*r*) = ∞0 *<sup>x</sup>r*−1*e*<sup>−</sup>*xdx*, *r* > 0.

In this notation, obviously, *G*1,1 is an r.v. with the standard exponential distribution: <sup>P</sup>(*<sup>G</sup>*1,1 < *x*) = )1 − *e*<sup>−</sup>*<sup>x</sup>*\***1**(*x* 0) (here and in what follows **1**(*A*) is the indicator function of a set *A*).

A GG-distribution is the absolutely continuous distribution defined by the density (Equation (2)). The distribution function (d.f.) corresponding to the density *g*<sup>∗</sup>(*<sup>x</sup>*;*r*, *γ*, *μ*) will be denoted *<sup>F</sup>*<sup>∗</sup>(*<sup>x</sup>*;*r*, *γ*, *μ*).

The properties of GG-distributions are described in [11,19]. An r.v. with the density *g*<sup>∗</sup>(*<sup>x</sup>*;*r*, *γ*, *μ*) will be denoted *Gr*,*γ*,*μ*. It can be easily made sure that

$$\overline{\mathbf{G}}\_{r,\gamma,\mu} \stackrel{d}{=} \mathbf{G}\_{r,\mu}^{1/\gamma} \,. \tag{3}$$

and hence,

$$(\overline{G}\_{r,\gamma,\mu})^{\gamma} \stackrel{d}{=} G\_{r,\mu} \tag{4}$$

For convenience, for an r.v. with the Weibull distribution, a particular case of GG-distributions corresponding to the density *g*<sup>∗</sup>(*x*; 1, *γ*, 1) and the d.f. )1 − *e*<sup>−</sup>*x<sup>γ</sup>* \***1**(*x* 0) with *γ* > 0, we will use a special notation *W<sup>α</sup>*, that is, *<sup>W</sup>γ d*= *<sup>G</sup>*1,*γ*,1. Thus, *G*1,1 *d*= *W*1. The density *g*<sup>∗</sup>(*x*; 1, *α*, 1) with *α* < 0 defines the Fréchet or inverse Weibull distribution. It is easy to see that

$$\mathcal{W}\_1^{1/\gamma} \stackrel{d}{=} \mathcal{W}\_{\gamma}.\tag{5}$$

An r.v. *Nr*,*<sup>p</sup>* is said to have the negative binomial (NB) distribution with parameters *r* > 0 (shape) and *p* ∈ (0, 1) (success probability), if

$$\mathbb{P}(N\_{r,p} = k) = \frac{\Gamma(r+k)}{k!\Gamma(r)} \cdot p^r (1-p)^k, \quad k = 0, 1, 2, \dots$$

A particular case of the NB distribution corresponding to the value *r* = 1 is the geometric distribution. Let *p* ∈ (0, 1) and let *<sup>N</sup>*1,*p* be the r.v. having the geometric distribution with parameter *p* :

$$\mathbb{P}(N\_{1,p} = k) = p(1-p)^k, \quad k = 0, 1, 2, \dots$$

This means that for any *m* ∈ N

$$P(N\_{1,p} \gtrsim m) = \sum\_{k=m}^{\infty} p(1-p)^k = (1-p)^m.$$

Let *Y* be an r.v. taking values in the interval (0, 1). Moreover, let for all *p* ∈ (0, 1) the r.v. *Y* and the geometrically distributed r.v. *<sup>N</sup>*1,*p* be independent. Let *M* = *N*1,*Y*, that is, *<sup>M</sup>*(*ω*) = *<sup>N</sup>*1,*<sup>Y</sup>*(*ω*)(*ω*) for any *ω* ∈ Ω. The distribution

$$\mathbb{P}(M \ni m) = \int\_0^1 (1 - y)^m d\mathbb{P}(Y < y), \quad m \in \mathbb{N}\_{\prime}$$

of the r.v. *M* will be called *Y*-mixed geometric [17].

It is well known that the negative binomial distribution is a mixed Poisson distribution with the gamma mixing distribution [12] (also see [20]): for any *r* > 0, *p* ∈ (0, 1) and *k* ∈ {0} 5 N we have

$$\frac{\Gamma(r+k)}{k!\Gamma(r)} \cdot p^r (1-p)^k = \frac{1}{k!} \int\_0^\infty e^{-z} z^k g(z; r, \mu) dz,\tag{6}$$

where *μ* = *p*/(<sup>1</sup> − *p*).

The d.f. and the density of a strictly stable distribution with the characteristic exponent *α* and shape parameter *θ* defined by the characteristic function (ch.f.)

$$f(t; \mathfrak{a}, \theta) = \exp\left\{-|t|^{\mathfrak{a}} \exp\{-\frac{1}{2}i\pi\theta \mathrm{asigm} t\}\right\}, \quad t \in \mathbb{R},$$

where 0 < *α* 2, |*θ*| min{1, 2*α* − <sup>1</sup>}, will be respectively denoted *<sup>F</sup>*(*x*; *α*, *θ*) and *f*(*x*; *α*, *θ*) (see, e.g., [21]). An r.v. with the d.f. *<sup>F</sup>*(*x*; *α*, *θ*) will be denoted *S<sup>α</sup>*,*θ*.

To symmetric, strictly stable distributions, there corresponds the value *θ* = 0. To one-sided strictly stable distributions concentrated on the nonnegative halfline, there correspond the values *θ* = 1 and 0 < *α* 1. The pairs *α* = 1, *θ* = ±1 correspond to the distributions degenerate in ±1, respectively. All the other strictly stable distributions are absolutely continuous. Stable densities cannot explicitly be represented via elementary functions with four exceptions: the normal distribution (*α* = 2, *θ* = 0), the Cauchy distribution (*α* = 1, *θ* = 0), the Lévy distribution (*α* = 12 , *θ* = 1) and the distribution symmetric to the Lévy law (*α* = 12 , *θ* = −1). Expressions of stable densities in terms of the Fox functions (generalized Meijer G-functions) can be found in [22,23].

The standard normal d.f. will be denoted <sup>Φ</sup>(*x*),

$$\Phi(x) = \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{x} e^{-y^2/2} dy, \quad y \in \mathbb{R}.$$

An r.v. with the d.f. <sup>Φ</sup>(*x*) will be denoted *X*. The folded or half- normal distribution is the distribution of the r.v. |*X*|. It can be easily verified that *S*2,0 *d*= √2*X*.

In [24,25], it was demonstrated that if *γ* ∈ (0, 1], then

$$\mathcal{W}\_{\gamma} \stackrel{d}{=} \mathcal{W}\_1 \cdot \mathbb{S}\_{\gamma, 1}^{-1} \tag{7}$$

with the r.v.'s on the right-hand side being independent.

For *r* ∈ (0, 1) let *Gr*, 1 and *G*1−*r*, 1 be independent gamma-distributed r.v.'s. Let *μ* > 0. Introduce the r.v.

$$Z\_{r, \mu} = \frac{\mu(G\_{r, 1} + G\_{1-r, 1})}{G\_{r, 1}} \stackrel{d}{=} \mu Z\_{r, 1} \stackrel{d}{=} \mu \left(1 + \frac{1-r}{r} Q\_{1-r, r}\right),\tag{8}$$

where *Q*<sup>1</sup>−*r*,*<sup>r</sup>* is the r.v. with the Snedecor–Fisher distribution defined by the probability density

$$q(\mathbf{x}; 1-r, r) = \frac{(1-r)^{1-r}r^r}{\Gamma(1-r)\Gamma(r)} \cdot \frac{1}{\mathbf{x}^r[r + (1-r)\mathbf{x}]'} \quad \mathbf{x} \gtrless 0. \tag{9}$$

In the paper [26], it was shown that any gamma distribution with shape parameter no greater than one is mixed exponential. For convenience, we formulate this result as the following lemma.

**Lemma 1** ([26])**.** *The density of a gamma distribution g*(*<sup>x</sup>*;*r*, *μ*) *with* 0 < *r* < 1 *can be represented as*

$$g(x;r,\mu) = \int\_0^\infty z e^{-zx} p(z;r,\mu) dz.$$

*where*

$$p(z;r,\mu) = \frac{\mu^r}{\Gamma(1-r)\Gamma(r)} \cdot \frac{\mathbf{1}(z \rhd \mu)}{(z-\mu)^r z}$$

*is the density of the r.v. Zr*,*<sup>μ</sup> introduced above. In other words, if* 0 < *r* < 1*, then*

$$G\_{r, \mu} \stackrel{d}{=} \frac{\mathcal{W}\_1}{Z\_{r, \mu}} \, \tag{10}$$

*where the random variables W*1 *and Zr*,*<sup>μ</sup> are independent. Moreover, a gamma distribution with shape parameter r* > 1 *cannot be represented as a mixed exponential distribution.*

Let *r* > 0, *γ* ∈ R and *μ* > 0. Let the r.v. *Nr*,*γ*,*<sup>μ</sup>* have the GNB distribution. Its d.f. will be denoted *FGNB*(*<sup>x</sup>*;*r*, *γ*, *μ*).

Along with the arguments given above in favor of the adequacy of the GNB models for the duration of wet periods based on their definition as mixed Poisson distributions, this effect can also be explained (at least in part) by their one more important property of being mixed geometric formulated as the following theorem.

**Theorem 1** ([3])**.** *If r* ∈ (0, 1]*, γ* ∈ (0, 1] *and μ* > 0*, then a GNB distribution is a Yr*,*γ*,*μ-mixed geometric distribution:*

$$\mathbb{P}(N\_{r,\gamma,\mu} = k) = \int\_0^1 y(1-y)^k d\mathbb{P}(Y\_{r,\gamma,\mu} < y), \quad k = 0, 1, 2... \tag{11}$$

*where*

$$\mathcal{Y}\_{r,\gamma,\mu} \stackrel{d}{=} \frac{\mathcal{S}\_{\gamma,1} \mathcal{Z}\_{r,\mu}^{1/\gamma}}{1 + \mathcal{S}\_{\gamma,1} \mathcal{Z}\_{r,\mu}^{1/\gamma}} \stackrel{d}{=} \frac{\mu^{1/\gamma} \mathcal{S}\_{\gamma,1} (\mathcal{G}\_{r,1} + \mathcal{G}\_{1-r,1})^{1/\gamma}}{\mathcal{G}\_{r,1}^{1/\gamma} + \mu^{1/\gamma} \mathcal{S}\_{\gamma,1} (\mathcal{G}\_{r,1} + \mathcal{G}\_{1-r,1})^{1/\gamma}},\tag{12}$$

*where the r.v.'s <sup>S</sup>γ*,<sup>1</sup> *and <sup>Z</sup>μ*,*<sup>r</sup> or <sup>S</sup>γ*,1*, Gr*,<sup>1</sup> *and G*1−*r*,<sup>1</sup> *are independent.*

#### **4. The Asymptotic Approximation to the Probability Distribution of Extremal Daily Precipitation within a Wet Period**

In this section, the probability distribution of extremal daily precipitation within a wet period will be deduced as an asymptotic approximation. We will require some auxiliary statements formulated as lemmas.

The following asymptotic property of the GNB distribution will play the fundamental role in the construction of asymptotic approximations to the distributions of extreme daily precipitation within a wet period and the total precipitation volume per wet period and the corresponding statistical tests for precipitation to be abnormally heavy.

**Lemma 2** ([3])**.** *For r* > 0*, γ* ∈ R*, μ* > 0 *let Nr*,*γ*,*<sup>μ</sup> be an r.v. with the GNB distribution. We have*

$$
\mu^{1/\gamma} N\_{r,\gamma\_r \mu} \implies \overline{G}\_{r,\gamma\_r 1} \stackrel{d}{=} G\_{r,1}^{1/\gamma} \tag{13}
$$

*as μ* → 0*. If, moreover, r* ∈ (0, 1] *and γ* ∈ (0, 1]*, then the limit law can be represented as*

$$\overline{\mathbf{G}}\_{r,\gamma,1} \stackrel{d}{=} \frac{\mathbf{W}\_1}{\mathbf{S}\_{\gamma,1} \mathbf{Z}\_{r,1}^{1/\gamma}} \stackrel{d}{=} \frac{\mathbf{W}\_1^{1/\gamma}}{\mathbf{Z}\_{r,1}^{1/\gamma}} \stackrel{d}{=} \left(\frac{\mathbf{W}\_1 \mathbf{G}\_{r,1}}{\mathbf{G}\_{r,1} + \mathbf{G}\_{1-r,1}}\right)^{1/\gamma} \stackrel{d}{=} \mathbf{W}\_1^{1/\gamma} \cdot \left(1 + \frac{\mathbf{1} - \mathbf{r}}{r} \mathbf{Q}\_{1-r,r}\right)^{-1/\gamma},\tag{14}$$

*where the r.v.'s W*1*, <sup>S</sup>γ*,<sup>1</sup> *and Zr*,<sup>1</sup> *are independent as well as the r.v.'s W*1 *and Zr*,1*, or the r.v.'s W*1*, Gr*,<sup>1</sup> *and G*1−*r*,1*, and the r.v. Q*<sup>1</sup>−*r*,*<sup>r</sup> has the Snedecor–Fisher distribution with parameters* 1 − *r and r, see Equation* (9)*.*

Let *μ* > 0, *γ* > 0. Instead of an infinitesimal parameter *μ*, in order to construct asymptotic approximations with "large" sample size, introduce an auxiliary "infinitely large" parameter *n* ∈ N and assume that *μ* = *μn* = *μn*<sup>−</sup>*γ*. It can be easily made sure that in this case

$$
\overline{G}\_{r,\gamma,\mu/n^{\gamma}} \stackrel{d}{=} n \overline{G}\_{r,\gamma,\mu}.\tag{15}
$$

Then for *r* > 0, *μ* > 0 for any *n* ∈ N, we have

$$n^{-1}\overline{\mathcal{G}}\_{r,\gamma\_{\gamma}\lambda/n^{\gamma}} \stackrel{d}{=} \overline{\mathcal{G}}\_{r,\gamma\_{\gamma}\lambda} \stackrel{d}{=} \lambda^{-1/\gamma}\overline{\mathcal{G}}\_{r,\gamma\_{\gamma}1} \stackrel{d}{=} \lambda^{-1/\gamma}\mathcal{G}\_{r,1}^{1/\gamma}.\tag{16}$$

The standard Poisson process (the Poisson process with unit intensity) will be denoted *<sup>P</sup>*(*t*), *t* 0.

**Lemma 3** ([27])**.** *Let* Λ1, Λ2, ... *be a sequence of positive r.v.'s such that for any n* ∈ N *the r.v.* Λ*n is independent of the standard Poisson process <sup>P</sup>*(*t*)*, t* 0*. The convergence*

$$n^{-1}P(\Lambda\_\mathbb{N}) \implies \Lambda$$

*as n* → ∞ *to some nonnegative r.v.* Λ *takes place if and only if*

$$n^{-1} \Lambda\_{\mathfrak{n}} \implies \Lambda\_{\mathfrak{r}} \quad n \to \infty. \tag{17}$$

Lemma 3 can be regarded as a special case of the following result. Consider a sequence of r.v.'s *W*1, *W*2, ... Let *N*1, *N*2, ... be natural-valued r.v.'s such that for every *n* ∈ N the r.v. *Nn* is independent of the sequence *W*1, *W*2, ... In the following statement, the convergence is meant as *n* → ∞.

**Lemma 4** ([28,29])**.** *Assume that there exists an infinitely increasing* (*convergen<sup>t</sup> to zero*) *sequence of positive numbers* {*bn*}*n*<sup>1</sup> *and an r.v. W such that*

$$b\_n^{-1} \mathcal{W}\_n \implies \mathcal{W}\_\*$$

*If there exist an infinitely increasing* (*convergen<sup>t</sup> to zero*) *sequence of positive numbers* {*dn*}*n*<sup>1</sup> *and an r.v. N such that*

$$d\_n^{-1} b\_{N\_n} \Longrightarrow N\_\prime \tag{18}$$

*then*

$$d\_n^{-1} \mathcal{W}\_{\mathcal{N}\_n} \Longrightarrow \mathcal{W} \cdot \mathcal{N}\_\* \tag{19}$$

*where the r.v.'s on the right-hand side of Equation* (19) *are independent. If, in addition, Nn* −→ ∞ *in probability and the family of scale mixtures of the d.f. of the r.v. W is identifiable, then Condition* (18) *is not only sufficient for Equation* (19)*, but is necessary as well.*

Consider a sequence of independent identically distributed (i.i.d.) r.v.'s *X*1, *X*2, .... Let *N*1, *N*2, ... be a sequence of natural-valued r.v.'s such that for each *n* ∈ N the r.v. *Nn* is independent of the sequence *X*1, *X*2, . . .. Denote *Mn* = max{*<sup>X</sup>*1,..., *XNn* }.

**Lemma 5** ([30])**.** *Let* Λ1, Λ2, ... *be a sequence of positive r.v.'s such that for each n* ∈ N *the r.v.* Λ*n is independent of the Poisson process <sup>P</sup>*(*t*)*, t* 0*. Let Nn* = *<sup>P</sup>*(<sup>Λ</sup>*n*)*. Assume that there exists a nonnegative r.v.* Λ *such that Convergence* (17) *takes place. Let X*1, *X*2, ... *be i.i.d. r.v.'s a common d.f. <sup>F</sup>*(*x*)*. Assume also that* sup{*x* : *<sup>F</sup>*(*x*) < 1} = ∞ *and there exists a number α* > 0 *such that for each x* > 0

$$\lim\_{y \to \infty} \frac{1 - F(xy)}{1 - F(y)} = x^{-a}. \tag{20}$$

*Then*

$$\lim\_{N \to \infty} \sup\_{x \geqslant 0} \left| \mathbb{P} \left( \frac{M\_{\text{II}}}{F^{-1} \left( 1 - \frac{1}{n} \right)} < x \right) - \int\_0^\infty e^{-zx^{-a}} d\mathbb{P} \left( \Lambda < z \right) \right| = 0.1$$

Now we turn to the main results of this section. The principal role in our reasoning will be played by Lemma 5. In order to justify its applicability, we need to make sure that the daily precipitation volumes satisfy Condition (20). A thorough statistical analysis shows that, although being rather adequate and, in general, acceptable model, the traditional gamma distribution (used, e.g., in [4]) is not the best model for statistical regularities in daily precipitation. The analysis of meteorological data (daily precipitation volumes) registered over 60 years at two geographic points with a very different climate: Potsdam (Brandenburg, Germany) with a mild climate influenced by the closeness to the ocean with warm Gulfstream flow and Elista (Kalmykia, Russia) with a radically continental climate convincingly suggests the Pareto-type model for the distribution of daily precipitation volumes, see Figures 3 and 4. For comparison, on these figures, the graphs of the best gamma-densities there are also presented. It can be seen that the gamma model fits the histograms in a noticeably worse way than the Pareto distribution.

**Theorem 2.** *Let n* ∈ N*, γ* > 0*, μ* > 0 *and let Nr*,*γ*,*μ<sup>n</sup> be an r.v. with the GNB distribution with parameters r* > 0*, γ* > 0 *and μn* = *μ*/*nγ. Let X*1, *X*2, ... *be i.i.d. r.v.'s with a common d.f. <sup>F</sup>*(*x*)*. Assume that* rext(*F*) = ∞ *and there exists a number α* > 0 *such that Relation* (20) *holds for any x* > 0*. Then*

$$\lim\_{n \to \infty} \sup\_{x \geqslant 0} \left| \mathbb{P} \left( \frac{\max \{ X\_1, \dots, X\_{N\_{r, \gamma, \mu\_0}} \}}{F^{-1} \left( 1 - \frac{1}{n} \right)} < x \right) - F(x; r, \kappa, \gamma, \mu) \right| = 0, \nu$$

*where*

$$F(\mathbf{x}; r, \mathbf{a}, \gamma, \mu) = \int\_0^\infty e^{-\lambda x^{-\mathfrak{a}}} g^\*(\lambda; r, \gamma, \mu) d\lambda \equiv \mathbb{P}(M\_{r, \mathbf{a}, \gamma, \mu} < \mathbf{x}), \quad \mathbf{x} \in \mathbb{R}.$$

*The limit r.v. Mr*,*α*,*γ*,*<sup>μ</sup> admits the following product representations:*

$$M\_{r,a,\gamma,\mu} \stackrel{d}{=} \frac{\overline{G}\_{r,a\gamma,\mu}}{W\_a} \stackrel{d}{=} \left(\frac{\overline{G}\_{r,\gamma,\mu}}{W\_1}\right)^{1/a} \stackrel{d}{=} \mu^{-1/a\gamma} \left(\frac{G\_{r,1}}{W\_\gamma}\right)^{1/a\gamma} \tag{21}$$

*and in each term, the involved random variables are independent.*

**Proof.** By definition, the GNB distribution is a mixed Poisson distribution with the GG mixing distribution. So, *Nr*,*γ*,*μ<sup>n</sup> d*= *P*-*Gr*,*γ*,*μ<sup>n</sup>* . Therefore, from Equation (16), Lemma 3 with Λ*n* = *Gr*,*γ*,*μ<sup>n</sup>* and Lemma 5 with the account of the absolute continuity of the limit distribution it immediately follows that

$$\lim\_{n \to \infty} \sup\_{x \ge 0} \left| \mathbb{P} \left( \frac{\max \{ X\_1, \dots, X\_{N\_{\tau/q, n}} \}}{F^{-1} (1 - \frac{1}{n})} < x \right) - \int\_0^\infty e^{-zx^{-a}} g^\*(z; r, \gamma, \mu) dz \right| = 0.1$$

Since the Fréchet (inverse Weibull) d.f. *e*<sup>−</sup>*<sup>x</sup>* −*α* with *α* > 0 corresponds to the r.v. *W*−<sup>1</sup> *α* , it is easy to make sure

$$F(\mathbf{x}; r, \mathfrak{a}, \gamma, \mu) \equiv \int\_0^\infty e^{-z\mathbf{x}^{-\mathfrak{a}}} g^\*(z; r, \gamma, \mu) dz = \mathbb{P}\left(\frac{\overline{G}\_{r, \gamma, \mu}^{1/\mathfrak{a}}}{W\_\mathfrak{a}} < x\right).$$

Moreover, using relation *Gr*,*γ*,*<sup>μ</sup> d*= *G*1/*<sup>γ</sup> <sup>r</sup>*,*μ* , it is easy to see that

$$\frac{\overline{G}\_{r,\gamma,\mu}^{1/a}}{\mathcal{W}\_{\boldsymbol{\alpha}}} \stackrel{d}{=} \frac{\overline{G}\_{r,\boldsymbol{a}\gamma,\mu}}{\mathcal{W}\_{\boldsymbol{a}}} \stackrel{d}{=} \left(\frac{\overline{G}\_{r,\boldsymbol{\gamma},\mu}}{\mathcal{W}\_{1}}\right)^{1/a} \stackrel{d}{=} \mu^{-1/a\gamma} \left(\frac{G\_{r,1}}{\mathcal{W}\_{\boldsymbol{\gamma}}}\right)^{1/a\gamma}$$

where in each term the involved random variables are independent. The theorem is proved.

If *γ* = 1, then the limit distribution *<sup>F</sup>*(*<sup>x</sup>*;*r*, *α*, 1, *μ*) corresponds to the results of [31,32].

**Theorem 3.** *The distribution of the r.v. Mr*,*α*,*γ*,*<sup>μ</sup> admits the following representations.*

*(i) If r* ∈ (0, 1]*, it is the scale mixture of the distribution of the ratio of two independent Weibull-distributed r.v.'s:*

$$M\_{r, \alpha, \gamma, \mu} \overset{d}{=} \left(\mu Z\_{r, 1}\right)^{-1/\alpha\gamma} \cdot \frac{\mathcal{W}\_{\alpha\gamma}}{\mathcal{W}\_{\gamma}}.$$

*where all the involved random variables are independent and the r.v. Zr*,<sup>1</sup> *is defined in Equation* (8)*.*

*(ii) If γ* ∈ (0, 1]*, it is the scale mixture of the tempered Snedecor–Fisher distribution with parameters r and* 1*:*

$$\mathcal{M}\_{r,\alpha,\gamma,\mu} \stackrel{d}{=} \left(\frac{S\_{\gamma,1}}{\mu r} \cdot Q\_{r,1}\right)^{1/\alpha\gamma}\prime$$

*where <sup>S</sup>γ*,<sup>1</sup> *is a positive strictly stable r.v. with characteristic exponent γ independent of the r.v. Qr*,<sup>1</sup> *with the Snedecor–Fisher distribution in Equation* (9) *with parameters r and* 1*.*

*(iii) If γ* ∈ (0, 1] *and r* ∈ (0, 1]*, it is the scale mixture of the Pareto laws:*

$$M\_{r, \mathfrak{a}, \gamma, \mu} \overset{d}{=} \Pi\_{\mathfrak{a}} \left( S\_{\gamma, 1} Z\_{r, 1}^{1/\gamma} \right)^{-1/a} \prime$$

*where* <sup>P</sup>(<sup>Π</sup>*α* > *x*)=(*x<sup>α</sup>* + <sup>1</sup>)−1*, x* 0*.*

*(iv) If r* ∈ (0, 1] *and αγ* ∈ (0, 1]*, it is the scale mixture of the folded normal laws:*

$$\mathcal{M}\_{\mathcal{T},\boldsymbol{\alpha},\boldsymbol{\gamma},\mu} \overset{d}{=} |X| \cdot \frac{\sqrt{2W\_1}}{\mu^{1/\alpha\gamma} \mathcal{W}\_{\boldsymbol{\alpha}} \mathcal{S}\_{\boldsymbol{\alpha}\gamma,1} Z\_{r,1}^{1/\alpha\gamma}},$$

*where all the involved r.v.'s are independent.*

**Proof.** To prove (i) it suffices to consider the rightmost term in Equation (21), apply relations *W*1/*<sup>γ</sup>* 1 *d*= *<sup>W</sup>γ* and *Gr*,*<sup>μ</sup> d*= *W*1 *Zr*,*<sup>μ</sup>* (here 0 < *r* < 1 and the r.v.'s *W*1 and *Zr*,*<sup>μ</sup>* are independent (for details, see Lemma 1).

To prove (ii) it suffices to transform the rightmost term in (21) with the account of representation in Equation (7) and use the definition of the Snedecor–Fisher distribution as the distribution of the ratio of two independent gamma-distributed r.v.'s (see, e.g., Section 27 in [33]).

To prove (iii) it suffices to transform the second term in Equation (21) with the account of Equation (14) and notice that the distribution of the ratio of two independent exponentially distributed r.v.'s coincides with that of the random variable Π1.

To prove (iv) it suffices to transform the second term in (21) with the account of (14) and notice that *W*1 *d*= |*X*|√2*W*1 with the r.v.'s on the right-hand side being independent (see, e.g., [25]). The theorem is proved.

**Figure 3.** The histogram of daily precipitation volumes in Potsdam and the fitted Pareto and gamma distributions.

**Figure 4.** The histogram of daily precipitation volumes in Elista and the fitted Pareto and gamma distributions.

The product representations for the random value *Mr*,*α*,*γ*,*<sup>μ</sup>* established in Theorem 3 can be used for computer simulation.

**Theorem 4.** *If r* ∈ (0, 1]*, μ* > 0 *and αγ* ∈ (0, 1]*, then the d.f. <sup>F</sup>*(*<sup>x</sup>*;*r*, *α*, *γ*, *μ*) *is mixed exponential:*

$$1 - F(\mathfrak{x}; r, \mathfrak{a}, \gamma, \mathfrak{y}, \mathfrak{y}) = \int\_0^\infty e^{-\mathfrak{u}\mathfrak{x}} dA(\mathfrak{u}), \quad \mathfrak{x} \not\ni 0,$$

*where <sup>A</sup>*(*u*) = P*<sup>μ</sup>*1/*αγWαSαγ*,1*Z*1/*αγ <sup>r</sup>*,1 < *<sup>u</sup>, u* 0*, and all the involved r.v.'s are independent.*

**Proof.** To prove this statement, it suffices to transform the second term in Equation (21) with the account of Equation (14) and obtain

$$\mathcal{M}\_{r,\alpha,\gamma,\mu} \stackrel{d}{=} \frac{\mathcal{W}\_1}{\mu^{1/a\gamma}\mathcal{W}\_a\mathcal{S}\_{\alpha\gamma,1}Z\_{r,1}^{1/a\gamma}}.$$

**Corollary 1.** *Let r* ∈ (0, 1]*, αγ* ∈ (0, 1]*, μ* > 0*. Then the distribution function <sup>F</sup>*(*<sup>x</sup>*;*r*, *α*, *γ*, *μ*) *is infinitely divisible.*

**Proof.** This statement immediately follows from Theorem 3 and the result of [34] stating that the product of two independent non-negative r.v.'s is infinitely divisible, if one of the two is exponentially distributed.

It is possible to deduce explicit expressions for the moments of the r.v. *Mr*,*α*,*γ*,*μ*.

**Theorem 5.** *Let* 0 < *δ* < *α. Then*

$$\mathbb{E}M\_{r,\alpha,\gamma,\mu}^{\delta} = \frac{\Gamma\left(r + \frac{\delta}{\alpha\gamma}\right)\Gamma\left(1 - \frac{\delta}{a}\right)}{\mu^{\delta/\alpha}\gamma\Gamma(r)}.$$

**Proof.** From Equation (14) it follows that <sup>E</sup>*Mδ<sup>r</sup>*,*α*,*γ*,*<sup>μ</sup>* = *μ*<sup>−</sup>*δ*/*αγ*E*Gδ*/*αγ <sup>r</sup>*,1 · E*W*−*δ*/*<sup>α</sup>* 1 . It is easy to verify that E*Gδ*/*αγ <sup>r</sup>*,1 = Γ*r* + *δαγ* /Γ(*r*), E*W*−*δ*/*<sup>α</sup>* 1 = Γ-1 − *δα* . Hence follows the desired result.

Consider the bounds for the rate of convergence in Theorem 2. For this purpose, we will use one more auxiliary statement.

**Lemma 6.** *Let λ* > 0*, X*1, *X*2, ... *be i.i.d. r.v.'s with a common d.f. <sup>F</sup>*(*x*)*, P*(*t*) *be the standard Poisson process independent of X*1, *X*2,...*. Assume that there exists a d.f. <sup>H</sup>*(*x*) *such that for any x* ∈ R

$$\lim\_{n \to \infty} \mathbb{P}\left(\frac{1}{F^{-1}\left(1 - \frac{1}{n}\right)} \max\_{1 \le k \le n} X\_k < \infty\right) = H(\mathbf{x}).\tag{22}$$

.

*Then for any n* ∈ N

$$\left| \mathbb{P} \left( \frac{1}{F^{-1} \left( 1 - \frac{1}{n} \right)} \max\_{1 \le k \le P(n\lambda)} X\_k < \mathbf{x} \right) - H^{\lambda} (\mathbf{x}) \right| \leqslant \left| n \left[ 1 - F \left( \mathbf{x} F^{-1} (1 - \frac{1}{n}) \right) \right] - \log H(\mathbf{x}) \right| \lambda H^{\lambda} (\mathbf{x}).$$

**Proof.** This statement is a special case of Corollary 2 in [35].

**Theorem 6.** *Let n* ∈ N*, γ* > 0*, μ* > 0 *and let Nr*,*γ*,*μ<sup>n</sup> be an r.v. with the GNB distribution with parameters r* > 0*, γ* > 0 *and μn* = *μ*/*nγ. Let X*1, *X*2,... *be i.i.d. r.v.'s with the common Pareto d.f.*

$$F(x) = 1 - \frac{c}{a x^a + c'}, \quad x \gg 0,$$

*with a*, *c*, *α* > 0*,*

$$F(x;r,\varkappa,\gamma,\mu) = \int\_0^\infty e^{-\lambda x^{-x}} g^\*(\lambda; r, \gamma, \mu) d\lambda, \quad x \in \mathbb{R}.$$

*Then for any x* ∈ R

$$\left| \mathbb{P} \left( \left[ \frac{a}{c(n-1)} \right]^{1/\gamma} \max\_{1 \le k \le N\_{r,\gamma,\mu\_n}} X\_k < \infty \right) - F(x; r, \alpha, \gamma, \mu) \right| \lesssim \varepsilon$$

$$\leq \left| \frac{\mathbf{x}^{a} - \mathbf{1}}{\mathbf{x}^{a}(n - 1) + 1} \right| \cdot \int\_{0}^{\infty} \lambda e^{-\lambda \mathbf{x}^{-a}} \mathbf{g}^{\*} (\lambda ; r\_{\prime} \gamma\_{\prime} \mu) d\lambda \\ \lesssim \left| \frac{\mathbf{x}^{a} - \mathbf{1}}{\mathbf{x}^{a}(n - 1) + 1} \right| \cdot \frac{\Gamma(r + \frac{1}{\gamma})}{\mu^{1/\gamma} \Gamma(r)}.$$

**Proof.** First of all, check Condition (20). We have

$$\frac{1 - F(xy)}{1 - F(y)} = \frac{ay^{\alpha} + c}{a\mathbf{x}^{\alpha}y^{\alpha} + c} \longrightarrow \mathbf{x}^{-\alpha}$$

as *y* → <sup>∞</sup>, that is, Condition (20) holds implying Equation (22) with *<sup>H</sup>*(*x*) = *e*<sup>−</sup>*<sup>x</sup>* −*α* in accordance with the classical theory of extremes (see, e.g., [36]). Second, note that in the case under consideration *<sup>F</sup>*−<sup>1</sup>(<sup>1</sup> − 1*n* )=[ *<sup>c</sup>*(*<sup>n</sup>*−<sup>1</sup>) *a* ]1/*α* so that *FxF*−<sup>1</sup>(<sup>1</sup> − 1*n* ) = 1 − [*x<sup>α</sup>*(*n* − 1) + <sup>1</sup>]−<sup>1</sup> and

$$n\left[1 - F\left(\mathbf{x}F^{-1}\left(1 - \frac{1}{n}\right)\right)\right] - \log H(\mathbf{x}) = \frac{n}{\mathbf{x}^a(n-1) + 1} - \frac{1}{\mathbf{x}^a} = \frac{\mathbf{x}^a - 1}{\mathbf{x}^a(n-1) + 1}.$$

Third, from Equation (15) it follows that *Nr*,*γ*,*μ*/*n<sup>γ</sup> d*= *<sup>P</sup>*(*nGr*,*γ*,*<sup>μ</sup>*) with independent *P*(*t*) and *Gr*,*γ*,*μ*. Therefore, by Lemma 6 we have

$$\begin{split} & \left| \mathsf{P} \left( \left[ \frac{a}{c(n-1)} \right]^{1/\gamma} \max\_{1 \le k \le k\_{r\gamma,\mu\_{\mathbb{P}}}} X\_k < \infty \right) - \mathsf{F} (\mathsf{x}; r, a, \gamma, \mu) \right| \leqslant \\ & \leqslant \int\_0^\infty \left| \mathsf{P} \left( \left[ \frac{a}{c(n-1)} \right]^{1/\gamma} \max\_{1 \le k \le \mathbb{P}(n\lambda)} X\_k < \infty \right) - e^{-\lambda x^{-\alpha}} \right| \mathsf{g}^\*(\lambda; r, \gamma, \mu) d\lambda \leqslant \\ & \leqslant \left| \frac{\mathsf{x}^a - 1}{\mathsf{x}^a (n-1) + 1} \right| \cdot \int\_0^\infty \lambda e^{-\lambda \mathsf{x}^a - \mathsf{g}^\*} \mathsf{g}^\*(\lambda; r, \gamma, \mu) d\lambda \leqslant \left| \frac{\mathsf{x}^a - 1}{\mathsf{x}^a (n-1) + 1} \right| \cdot \mathsf{E} \overline{\mathsf{G}}\_{r, \gamma, \mu} = \\ & = \left| \frac{\mathsf{x}^a - 1}{\mathsf{x}^a (n-1) + 1} \right| \cdot \frac{\Gamma(r + \frac{1}{\gamma})}{\mu^{1/\gamma} \Gamma(r)}. \end{split}$$

The theorem is proved.

Actually, Theorem 6 states that the rate of convergence in Theorem 2 is *<sup>O</sup>*(*μ*1/*<sup>γ</sup> n* ) as *μn* → 0.

The results of this section serve as a theoretical base for the construction of a test for abnormally extreme daily precipitation. The distribution of the maximum daily precipitation per wet period can be assumed to be a tempered scale mixture of the gamma distribution in which the scale factor has the Weibull distribution. According to the typical construction of a test, a daily precipitation volume is considered to be abnormally extremal, if it exceeds a certain (pre-defined) quantile of this distribution. A detailed description of this test and algorithm of estimation of the parameters of the distribution mentioned above deserve a separate study as well as its application to real data.

#### **5. The Asymptotic Approximation to the Probability Distribution of Total Precipitation over a Wet Period. Generalized R ényi Theorem for Gnb Random Sums**

As far ago as in the 1950s, being interested in modeling rare events, A. Rényi studied rarefaction of renewal point processes and proved his famous theorem on convergence of rarefied renewal processes to the Poisson process [37,38]. The Rényi theorem states that the distribution of a geometric sum

(i.e., a sum of a random number of i.i.d. r.v.'s in which the number of summands is a r.v. with the geometric distribution independent of the summands) normalized by its expectation converges to the exponential law as the expectation of the sum infinitely increases. The normalization of a sum by its expectation is typical for laws of large numbers. Therefore, the Rényi theorem can be regarded as the law of large numbers for geometric sums. A general law of large numbers for random sums of independent identically distributed (i.i.d.) random variables (r.v.'s) was proved in [28]. It was demonstrated there that the distribution of a random sum normalized by its expectation converges to some distribution, if and only if the distribution of the random index (the number of summands) converges to the same distribution (up to a scale parameter) under the same normalization. In [3] the law of large numbers for GNB random sums was proved. However, a direct application of this result to modeling the probability distribution of total precipitation over a wet period is hampered by the following very interesting practical observation.

One might have expected that successive daily precipitation volumes *X*1, *X*2, ... satisfy the classical law of large numbers, that is, the arithmetic mean 1 *n* (*<sup>X</sup>*1 + ... + *Xn*) converges to some number *a* almost surely as *n* infinitely grows, as it was done in [2]. However, a thorough analysis of real data shows that this not quite so. In Figure 5, there are the graphs of the averaged daily precipitation volumes in Potsdam and Elista demonstrating the slowly decreasing trend for Potsdam and slowly increasing trend for Elista.

**Figure 5.** Stabilization of the cumulative averages of daily precipitation volumes as *n* grows in Potsdam (continuous line) and Elista (dash line).

This means that, in order to match the stabilization of the averages at some level *a*, it is required to normalize the sum *X*1 + ... + *Xn* not by *n*, but by a somewhat more complicated function of *n* that can match the influence of slow global trends. As such, a function of *n*, consider a power function *nβ* with *β* > 0 and assume that not necessarily i.i.d. r.v.'s *X*1, *X*2, . . . satisfy the condition

$$\frac{1}{m^{\beta}}\sum\_{j=1}^{n}X\_{j} \implies a \in (0,\infty) \tag{23}$$

as *n* → ∞. The parameters *a* and *β* can be rather reliably estimated by the least squares technique.

Let *X*1, *X*2, ... , *Xn* be the observed values of successive nonzero daily precipitation volumes, *n* ∈ N be the total number of available observations. For a natural *k* = 1, ... , *n* denote *sk* = *X*1 + ... + *Xk*. If Condition (23) holds, then for *k* large enough (1 *m k n*), the following estimates of the parameters *a* and *β* in Relation (23) can be used:

$$\widetilde{\alpha} = \exp\left\{ \frac{\sum\_{k=m}^{n} \log s\_k \cdot \sum\_{k=m}^{n} (\log k)^2 - \sum\_{k=m}^{n} \log k \cdot \sum\_{k=m}^{n} \left( \log k \cdot \log s\_k \right)}{(n-m+1) \sum\_{k=m}^{n} (\log k)^2 - \left( \sum\_{k=m}^{n} \log k \right)^2} \right\},\tag{24}$$

$$\widetilde{\beta} = \frac{\sum\_{k=m}^{n} \log s\_k - (n-m+1)\log \widetilde{\alpha}}{\sum\_{k=m}^{n} \log k}.\tag{25}$$

Indeed, if Condition (23) holds, the following approximate equality can be written:

$$\frac{T\_k}{k^{\beta}} \approx a \Leftrightarrow -\beta \log k + \log T\_k \approx \log a.$$

Therefore, the estimates of the parameters *a* and *β* can be found as the solution of the least squares problem

$$\sum\_{k=m}^n (\log T\_k - \beta \log k - \log a)^2 \longrightarrow \min\_{\beta, \log a} \dots$$

This solution can be found explicitly and has the form

$$\begin{split} \widehat{\log a} &= \frac{\sum\_{k=m}^{n} \log T\_k \cdot \sum\_{k=m}^{n} (\log k)^2 - \sum\_{k=m}^{n} \log k \cdot \sum\_{k=m}^{n} \left( \log k \cdot \log T\_k \right)}{\left( n - m + 1 \right) \sum\_{k=m}^{n} (\log k)^2 - \left( \sum\_{k=m}^{n} \log k \right)^2}, \\ &\tilde{\beta} = \frac{\sum\_{k=m}^{n} \log T\_k - (n - m + 1) \widehat{\log a}}{\sum\_{k=m}^{n} \log k}, \end{split}$$

that leads to Formulas (24) and (25). This least squares method for estimation of *a* and *β* is realized by Algorithm A2 (see Appendix A).

The application of Equation (23) to real data from Potsdam and Elista with *a* and *β* estimated by Equations (24) and (25) is illustrated in Figure 6. It can be seen that the cumulative averages stabilize at the level *a* = 4.087 with *β* = 0.981 for Potsdam and at the level *a* = 0.96 with *β* = 1.146 for Elista.

**Figure 6.** Stabilization of the cumulative averages of daily precipitation volumes as *n* grows with *β* = 1.139 for Potsdam (solid line) and with *β* = 0.981 for Elista (dashed line).

So, to construct the asymptotic approximation to the probability distribution of total precipitation over a wet period, we should prove a generalized Rényi theorem for GNB random sums improving an analogous statement proved in [3]. It must be especially noted that in the following theorem, the r.v.'s *X*1, *X*2, ... are not assumed to be i.i.d.

**Theorem 7.** *Assume that the nonzero daily precipitation volumes X*1, *X*2, ... *satisfy Condition* (23) *with some β* > 0 *and a* > 0*. Let the numbers r* > 0*, γ and μ* > 0 *be arbitrary. For each n* ∈ N*, let the r.v. Nr*,*γ*,*μ<sup>n</sup> have the GNB distribution with parameters r, γ and μn* = *μ*/*nγ. Assume that for each n* ∈ N *the r.v. Nr*,*γ*,*μ<sup>n</sup> is independent of the sequence X*1, *X*2, ... *Then*

$$\frac{a\mu^{\otimes/\gamma}}{n^{\beta}}\sum\_{j=1}^{N\_{r,\gamma,\mu/n^{\gamma}}}X\_j \implies \overline{G}\_{r,\gamma/\beta,1} \overset{d}{=} G\_{r,1}^{\beta/\gamma}$$

*as n* → ∞*.*

**Proof.** The proof is based on Lemma 4 and Equation (16). From Equation (16) it follows that

$$\frac{\mu^{1/\gamma}}{n} \cdot \mathcal{N}\_{r,\gamma,\mu'/n^{\gamma}} \implies \overline{\mathcal{G}}\_{r,\gamma,1} \tag{26}$$

as *n* → ∞. By virtue of Condition (23), in Lemma 4 let *bn* = *<sup>n</sup>β*/*<sup>a</sup>*. As *Nn* in Lemma 4 take *Nr*,*γ*,*μ*/*n<sup>γ</sup>* . Then *bNn* = 1*a <sup>N</sup>β<sup>r</sup>*,*γ*,*μ*/*n<sup>γ</sup>* . From Equation (26) it follows that, as *n* → <sup>∞</sup>,

$$\frac{1}{4} \mathcal{N}\_{r,\gamma,\mu/n^{\gamma}}^{\mathcal{S}} \cdot \frac{\mu^{\mathcal{S}/\gamma}}{n^{\mathcal{S}}} \implies \frac{1}{4} \overline{\mathcal{G}}\_{r,\gamma,1}^{\mathcal{S}} \stackrel{d}{=} \frac{1}{4} \overline{\mathcal{G}}\_{r,\gamma/\beta,1} \stackrel{d}{=} \frac{1}{a} \mathcal{G}\_{r,1}^{\mathcal{S}/\gamma}. \tag{27}$$

Therefore, as *dn* we can take *dn* = *<sup>n</sup>β*/*μβ*/*γ*. So, using Equation (27) in the role of Equation (18) in Lemma 4, we obtain Equation (19) in the form

$$\frac{\mu^{\otimes/\gamma}}{n^{\otimes}}\sum\_{j=1}^{N\_{r,\gamma,\mu/n^{\gamma}}}X\_j \implies \frac{1}{a}\mathsf{G}\_{r,\gamma/\beta,1} \triangleq \frac{1}{a}\mathsf{G}\_{r,1}^{\beta/\gamma},\tag{28}$$

whence follows the desired result. The theorem is proved.

Theorem 7 presents a good tool for the account of the parameters *β* and *γ* characterizing the deviation from traditional NB and arithmetic mean models due to the influence of possible (slow) global trends. If in Theorem 7 *r* = *γ* = *β* = 1, then we obtain a version of the Rényi theorem [39] generalized to non-identically distributed and not necessarily independent summands. If in Theorem 7 *β* = 1, then we obtain the law of large numbers for GNB random sums (see [3]). If in Theorem 7 *γ* = 1, then we obtain the law of large numbers for NB random sums modified for the case *β* -= 1.

Therefore, if daily precipitation volumes *X*1, *X*2, ... (of course, being non-identically distributed and not independent), with the account of the excellent fit of the GNB model for the duration of a wet period (see Figure 1), with rather small *μ*, the GG distribution can be regarded as an adequate and theoretically well-based model for the total precipitation volume over a (long enough) wet period.

As regards the bounds for the rate of convergence in Theorem 7, consider a special case of *β* = 1 and i.i.d. *X*1, *X*2,...

As a measure of the distance between probability distributions, consider the *ζ*-metric proposed by V. M. Zolotarev in [40,41] (also see [42], p. 44). Let *s* > 0. There exists a unique representation of the number *s* as *s* = *m* + *α* where *m* is an integer and 0 < *α* 1. By F*s* we denote the set of all real-valued bounded functions *f* on R that are *m* times differentiable and | *f* (*m*)(*x*) − *f* (*m*)(*y*)| |*x* − *y*|*<sup>α</sup>*. Let *X* and *Y* be two r.v.'s in which the distribution functions will be denoted *FX*(*x*) and *FY*(*x*), respectively. The *ζ*-metric *ζs*(*<sup>X</sup>*, *Y*) ≡ *ζs*(*FX*, *FY*) in the space of probability distributions is defined by the equality *ζs*(*<sup>X</sup>*, *Y*) = sup E- *f*(*X*) − *f*(*Y*) : *f* ∈ <sup>F</sup>*s*. In particular,

$$\mathcal{J}\_1(X,Y) = \int\_{\mathbb{R}} |F\_X(\mathbf{x}) - F\_Y(\mathbf{x})| d\mathbf{x}.$$

In [43], it was shown that in the case *β* = 1 and i.i.d. *X*1, *X*2, . . . for 1 *s* 2 we have

$$\mathbb{Z}\_s \left( \frac{\mu^{1/\gamma}}{na} \sum\_{j=1}^{N\_{r,\gamma,\mu/n^\gamma}} X\_{j\prime} \, \big| \, \mathbb{G}\_{r,\gamma,1} \right) \leqslant \frac{(\mathbb{E}X\_1^2)^{s/2}}{n^{s/2} |\mathbb{E}X\_1|^s} \cdot \frac{\Gamma(1+\gamma)\Gamma(r+\frac{s}{2\gamma})}{\Gamma(1+s)\Gamma(r)}.$$

In particular,

$$\zeta\_2 \left( \frac{\mu^{1/\gamma}}{na} \sum\_{j=1}^{N\_{r,\gamma,\mu/n^{\gamma}}} \mathbf{X}\_{j\prime} \, \overline{\mathbf{G}}\_{r,\gamma,\mu} \right) \lessdot \frac{\mathbb{E}X\_1^2}{2n(\mathbb{E}X\_1)^2} \cdot \frac{\Gamma(r + \frac{1}{\gamma})}{\Gamma(r)}.$$

The results presented above justify the GG models for the probability distribution of total precipitation volume over a wet period improving the models considered in [2]. Statistical tests for the detection of anomalously extreme total volumes will be considered below.

#### **6. Statistical Tests for Anomalously Extreme Total Precipitation Volumes**

Now we turn to the construction of the tests for the total precipitation volume during a wet period to be abnormally large.

In what follows, based on the results of the preceding section, we will assume that the total precipitation volume during a wet period has the GG distribution with some parameters *r* > 0, *γ* > 0 and *μ* > 0.

Let *m* ∈ N and *G*(1) *<sup>r</sup>*,*γ*,*μ*, *G*(2) *<sup>r</sup>*,*γ*,*μ*, ... , *G*(*m*) *<sup>r</sup>*,*γ*,*μ* be independent r.v.'s having the same GG distribution with parameters *r* > 0, *γ* and *μ* > 0. Also, let *G*(1) *<sup>r</sup>*,*μ* , *G*(2) *<sup>r</sup>*,*μ* , ... , *G*(*m*) *<sup>r</sup>*,*μ* be i.i.d. r.v.'s having the same gamma distribution with parameters *r* > 0 and *μ* > 0.

The base for the first step in the construction of the desired test is the following obvious conclusion: if the r.v.'s *G*(1) *<sup>r</sup>*,*γ*,*μ*, *G*(2) *<sup>r</sup>*,*γ*,*μ*, ... , *G*(*m*) *<sup>r</sup>*,*γ*,*μ* are identically distributed (that is, the sample *G*(1) *<sup>r</sup>*,*γ*,*μ*, *G*(2) *<sup>r</sup>*,*γ*,*μ*, ... , *G*(*m*) *<sup>r</sup>*,*γ*,*μ* is homogeneous), then the r.v.'s -*G*(1) *<sup>r</sup>*,*γ*,*<sup>μ</sup><sup>γ</sup>*, -*G*(2) *<sup>r</sup>*,*γ*,*<sup>μ</sup><sup>γ</sup>*, ... , -*G*(*m*) *<sup>r</sup>*,*γ*,*<sup>μ</sup><sup>γ</sup>* are also identically distributed (that is, the sample -*G*(1) *<sup>r</sup>*,*γ*,*<sup>μ</sup><sup>γ</sup>*, -*G*(2) *<sup>r</sup>*,*γ*,*<sup>μ</sup><sup>γ</sup>*,..., -*G*(*m*) *<sup>r</sup>*,*γ*,*<sup>μ</sup><sup>γ</sup>* is homogeneous.

Consider the relative contribution of the r.v. -*G*(1) *<sup>r</sup>*,*γ*,*<sup>μ</sup><sup>γ</sup>* to the sum -*G*(1) *<sup>r</sup>*,*γ*,*<sup>μ</sup><sup>γ</sup>* + -*G*(2) *<sup>r</sup>*,*γ*,*<sup>μ</sup><sup>γ</sup>* + ... + -*G*(*m*) *<sup>r</sup>*,*γ*,*<sup>μ</sup><sup>γ</sup>*:

$$R = \frac{(\overline{\mathbf{G}}\_{r,\gamma,\mu}^{(1)})^\gamma}{(\overline{\mathbf{G}}\_{r,\gamma,\mu}^{(1)})^\gamma + (\overline{\mathbf{G}}\_{r,\gamma,\mu}^{(2)})^\gamma + \dots + (\overline{\mathbf{G}}\_{r,\gamma,\mu}^{(m)})^\gamma} . \tag{29}$$

From Equation (4), it obviously follows that

$$R \stackrel{d}{=} \frac{G\_{r,\mu}^{(1)}}{G\_{r,\mu}^{(1)} + G\_{r,\mu}^{(2)} + \dots + G\_{r,\mu}^{(m)}} \stackrel{d}{=} \frac{G\_{r,1}^{(1)}}{G\_{r,1}^{(1)} + G\_{r,1}^{(2)} + \dots + G\_{r,1}^{(m)}} \stackrel{d}{=} R^\*$$

(see Equation(29)).

So, the r.v. *R* characterizes the relative precipitation volume for one (long enough) wet period with respect to the total precipitation volume registered for *m* wet periods.

Note that

$$R = \left(1 + \frac{1}{G\_{r,\mu}^{(1)}} (G\_{r,\mu}^{(2)} + \dots + G\_{r,\mu}^{(m)})\right)^{-1} \stackrel{d}{=} \left(1 + \frac{G\_{(m-1)r,\mu}}{G\_{r,\mu}}\right)^{-1} \stackrel{d}{\to}$$

where the gamma-distributed r.v.'s on the right hand side are independent. The distribution of the r.v. *R* was described in [2] where it was demonstrated that *R d*= -1 + *kr Qk*,*<sup>r</sup>*−<sup>1</sup> where *Qk*,*<sup>r</sup>* is the r.v. having the Snedecor–Fisher distribution determined for *k* > 0, *r* > 0 by the Lebesgue density

$$f\_{k,r}(x) = \frac{\Gamma(k+r)}{\Gamma(k)\Gamma(r)} \left(\frac{k}{r}\right)^k \frac{x^{k-1}}{(1+\frac{k}{r}x)^{k+r'}}, \quad x \gg 0. \tag{30}$$

It should be noted that the particular value of the scale parameter is insignificant. For convenience, it is assumed equal to one. It can be easily made sure by standard calculation using Equation (30), the distribution of the r.v. *R* is determined by the density

$$p(x;k,r) = \frac{\Gamma(k+r)}{\Gamma(r)\Gamma(k)}(1-x)^{k-1}x^{r-1}, \quad 0 \le x \ll 1, \ldots$$

that is, it is the beta distribution with parameters *k* = (*m* − <sup>1</sup>)*r* and *r*.

Then the test for the homogeneity of an independent sample of size *m* consisting of the GG distributed observations of total precipitation volumes during *m* wet periods with known *γ* based on the r.v. *R* looks as follows. Let *V*1, ... , *Vm* be the total precipitation volumes during *m* wet periods and, moreover, *V*1 *Vj* for all *j* 2. Calculate the quantity

$$SR = \frac{V\_1^{\gamma}}{V\_1^{\gamma} + \dots + V\_m^{\gamma}}$$

(*SR* means «*S*ample *R*»). From what was said above, it follows that under the hypothesis *H*0: «the precipitation volume *V*1 under consideration is not abnormally large» the r.v. *SR* has the beta distribution with parameters *k* = (*m* − <sup>1</sup>)*r* and *r*. Let *ε* ∈ (0, 1) be a small number, *β<sup>k</sup>*,*<sup>r</sup>*(<sup>1</sup> − *ε*) be the (1 − *ε*)-quantile of the beta distribution with parameters *k* = (*m* − <sup>1</sup>)*r* and *r*. If *SR* > *β<sup>k</sup>*,*<sup>r</sup>*(<sup>1</sup> − *ε*), then the hypothesis *H*0 must be rejected, that is, the volume *V*1 of precipitation during one wet period must be regarded as abnormally large. Moreover, the probability of erroneous rejection of *H*0 is equal to *ε*.

Instead of *R*, the quantity

$$R\_0 = \frac{(m-1)\left(\overline{G}\_{r,\gamma,\mu}^{(1)}\right)^\gamma}{\left(\overline{G}\_{r,\gamma,\mu}^{(2)}\right)^\gamma + \dots + \left(\overline{G}\_{r,\gamma,\mu}^{(m)}\right)^\gamma} \stackrel{d}{=} \frac{(m-1)G\_{r,\mu}^{(1)}}{G\_{r,\mu}^{(2)} + \dots + G\_{r,\mu}^{(m)}} \stackrel{d}{=} \frac{k}{r} \frac{G\_{r,\mu}}{G\_{k,\mu}} \stackrel{d}{=} \frac{k}{r} \frac{G\_{r,1}}{G\_{k,1}} \stackrel{d}{=} Q\_{r,k}$$

can be considered. Then, as is easily seen, the r.v.'s *R* and *R*0 are related by the one-to-one correspondence

$$R = \frac{R\_0}{m - 1 + R\_0} \text{ or } R\_0 = \frac{(m - 1)R}{1 - R}.$$

so that the homogeneity test for a sample from the GG distribution equivalent to the one described above and, correspondingly, the test for a precipitation volume during a wet period to be abnormally large, can be based on the r.v. *R*0, which has the Snedecor–Fisher distribution with parameters *r* and *k* = (*m* − <sup>1</sup>)*<sup>r</sup>*.

Namely, again let *V*1, ... , *Vm* be the total precipitation volumes during *m* wet periods and, moreover, *V*1 *Vj* for all *j* 2. Calculate the quantity

$$SR\_{GG} = \frac{(m-1)V\_1^{\gamma}}{V\_2^{\gamma} + \dots + V\_m^{\gamma}}.\tag{31}$$

(*SR*0 means «*S*ample *R*0»). From what was said above, it follows that under the hypothesis *H*0: «the precipitation volume *V*1 under consideration is not abnormally large» the r.v. *SR* has the Snedecor–Fisher distribution with parameters *r* and *k* = (*m* − <sup>1</sup>)*<sup>r</sup>*. Let *ε* ∈ (0, 1) be a small number, *qr*,*<sup>k</sup>*(<sup>1</sup> − *ε*) be the (1 − *ε*)-quantile of the Snedecor–Fisher distribution with parameters *r*

and *k* = (*m* − <sup>1</sup>)*<sup>r</sup>*. If *SR*0 > *qr*,*<sup>k</sup>*(<sup>1</sup> − *ε*), then the hypothesis *H*0 must be rejected, that is, the volume *V*1 of precipitation during one wet period must be regarded as abnormally large. Moreover, the probability of erroneous rejection of *H*0 is equal to *ε*.

Let *l* be a natural number, 1 *l* < *m*. It is worth noting that, unlike the test based on the statistic *R*, the test based on *R*0 can be modified for testing the hypothesis *<sup>H</sup>*0: «the precipitation volumes *Vi*1 , *Vi*2 , ... , *Vil* do not make an abnormally large cumulative contribution to the total precipitation volume *V*1 + ... + *Vm*». For this purpose denote

$$T\_l^\gamma = V\_{i\_1}^\gamma + V\_{i\_2}^\gamma + \dots + V\_{i\_l}^\gamma, \quad T^\gamma = V\_1^\gamma + V\_2^\gamma + \dots + V\_m^\gamma$$

and consider the quantity

$$SR\_0' = \frac{(m-l)T\_l^\gamma}{l(T^\gamma - T\_l^\gamma)}.$$

In the same way as it was done above, it is easy to make sure that

$$SR\_0' \stackrel{d}{=} \frac{(m-l)G\_{lr,l}}{lG\_{(m-l)r,1}} \stackrel{d}{=} Q\_{lr,(m-l)r} \cdot 1$$

Let *ε* ∈ (0, 1) be a small number, *qlr*,(*<sup>m</sup>*−<sup>1</sup>)*r*(<sup>1</sup> − *ε*) be the (1 − *ε*)-quantile of the Snedecor–Fisher distribution with parameters *lr* and *k* = (*m* − *l*)*<sup>r</sup>*. If *SR*0 > *qlr*,(*<sup>m</sup>*−*l*)*r*(<sup>1</sup> − *ε*), then the hypothesis *H*0 must be rejected, that is, the cumulative contribution of the precipitation volumes *Vi*1 , *Vi*2 , ... , *Vil* into the total precipitation volume *V*1 + ... + *Vm* must be regarded as abnormally large. Moreover, the probability of erroneous rejection of *H*0is equal to *ε*.

#### **7. Comparison of Tests for Anomalously Extreme Precipitation Volumes Based on Gamma and Gg Distributions**

In this section, the results of the application of the test based on the statistic *R* in Equation (29) to the analysis of the time series of daily precipitation observed in Potsdam and Elista from 1950 to 2007 are considered and compared with similar results for the case of gamma distributed total precipitation volumes during wet periods [2].

The results of the application of the tests for a total precipitation volume during one wet period to be abnormally large based on GG and gamma models in the moving mode are shown in Figures 7 and 8 (Potsdam) and Figures 9 and 10 (Elista).

If *m* is the window width (the number of observations in a moving window). A fixed sample point falls in exactly *m* windows. One of the following cases can occur for a fixed observation:


Algorithm A3 (see Appendix A) realizes the method based on the statistic *R* described above. For the sake of vividness on these figures the time horizon equals 90 and 360 days and the significance level *α* of the tests is 0.01. The absolutely, intermediate and relatively abnormal precipitation volumes are marked by downward-pointing triangles, circles and squares, respectively, for the test based on the gamma model, whereas the corresponding test based on the statistic *R* based on the GG distribution are marked by upward-pointing triangles, diamonds and right-pointing triangles, respectively. It is worth noting that MATLAB's notations are used here for these markers.

**Figure 7.** Abnormal precipitation volumes (Potsdam, 90 days).

**Figure 8.** Abnormal precipitation volumes (Potsdam, 360 days).

Figures 7–10 demonstrate non-trivial values of the parameter *γ*, that is, *γ* -= 1. For Potsdam *γ* = 1.286, whereas for Elista *γ* equals 1.279. At the same time, the results of the two methods are quite close, although the approach based on the GG distribution demonstrates a higher quality of determining potentially extreme observations. The same conclusions are valid for smaller window sizes.

**Figure 9.** Abnormal precipitation volumes (Elista, 90 days).

**Figure 10.** Abnormal precipitation volumes (Elista, 360 days).

#### **8. Comparison of GG-Based Statistical Test and Peaks over Threshold Methodology for Extreme Precipitation Intensities**

One important precipitation indicator is the precipitation intensity that is defined as the ratio of the total precipitation volume over a wet period to the duration of this wet period measured in days. The extreme precipitation volumes and intensities are relevant to various problems of climatology and hydrometeorology (see, for example, [44–47]). Traditionally, these phenomena are investigated for different geographical regions or countries [48,49]. In particular, the issue of determining threshold values, the excess of which leads to the extreme events, for example, in daily rainfalls or their intensities, is the key point of the study. Precipitation intensities are important not only for forecasting floods but also for solving problems such as runoff and soil erosion [50,51]. It can be explained by the contemporary climate change scenarios that predict a significant increase in the frequency of high intensity rainfall events, primarily in the dry areas. Moreover, precipitation can induce shallow landslides [52,53] and debris flows [54].

Statistical analysis of real data shows that the probability distribution of the precipitation intensity can be approximated by the gamma distribution with very high accuracy. In [55], some theoretical arguments were presented to justify the gamma model for the distribution of precipitation intensities. So, the statistical approach described in Section 6 and in [2] can be also used for identification of abnormally large intensities. For the analysis, the precipitation intensities in Potsdam and Elista (verified samples without missing values) are used as the initial data. This section presents a comparison of a non-parametric approach based on the extreme value theory as well as modified Peaks over Threshold (PoT) methodology [56] with the parametric approach that significantly involves testing parametric statistical hypotheses to determine extreme intensities of wet periods (see Section 6). The classical version of PoT [57] is quite popular for solving a wide range of climatic problems. In particular, the following results can be mentioned: the inverse Weibull distribution as an extreme wind speed model [58], a time-dependent versions of the PoT model for severe storm waves [59] and daily temperatures [60], probability model for rainfalls of high magnitude [61], analysis of precipitation extremes in a changing climate [62]. Most applications of the extreme value theory assume stationarity, but it is well-known that real events are not stationary. So, the generalized results analogous to Theorem 7 are required. All the numerical methods are implemented as a MATLAB program. Algorithm A4 demonstrates the method based on the PoT and GG-test (see Appendix A).

Figures 11 and 12 present the results obtained by the modified PoT algorithm in which the Weibull distribution is considered as the distribution of time between extreme events. Starting from the maximum threshold value that coincides with the maximum of the analyzed data, the hypothesis that the time intervals between the moments of excess of a certain threshold have the Weibull distribution is tested. The corresponding *P*-value is saved, and the threshold is shifted down by a certain (small) step (in this case, 0.01). It is worth noting that a similar procedure was suggested in [56] for precipitation volumes under the assumption that the time intervals between excesses have the exponential distribution. For a given significance level (in both cases, *α* is chosen as 0.01), the corresponding hypothesis is not rejected for all thresholds for which the *P*-values are located to the right of the red vertical line in the upper graphs in Figures 11 and 12. The lower graphs show the parameters of the fitted Weibull distribution.

On Figures 13 and 14 the results of the test (see Section 6) for both the GG and usual gamma distribution are compared with those of the PoT method based on the exponential and Weibull distributions for the intensities in Potsdam and Elista. The following notation is used:


**Figure 11.** P-values that correspond to the various thresholds (Potsdam).

**Figure 12.** P-values that correspond to the various thresholds (Elista).

The green-filled downward-pointing triangles mark the intensities, which are classified as absolutely abnormal based on the GG test (see Section 6). The black upward-pointing triangles correspond to the decision based on the classical gamma distribution test (that is, *γ* = 1, see (31)). The circles denote intermediate extreme observations, and the squares mark relatively extreme ones. This classification is described in Section 7. It is worth noting that for the GG test in Potsdam the value of *γ* is 1.0775 and for Elista *γ* = 1.1257.

For Potsdam, the results of gamma and GG tests are good and close. In addition, the PoT method is also effective in the case where the threshold is chosen with the maximum *P*-value. However, for Elista, with less rainfalls with lower intensities, the results are quite different. Indeed, the decisions of the PoT method are close for the exponential and Weibull cases (the thresholds differ by only 0.29). However, a statistical test based on the gamma distribution identifies only four intensities as absolutely extreme, while the GG test identifies more absolute extremes, including those below the thresholds mentioned above.

**Figure 13.** Comparison of statistical tests andpPeaks over threshold methodology for extreme precipitation intensities (Potsdam).

**Figure 14.** Comparison of statistical tests and peaks over threshold methodology for extreme precipitation intensities (Elista).

#### **9. Conclusions and Discussion**

In the paper, asymptotic models for some precipitation characteristics based on GNB distributions were considered. Also, a statistical test based on the GG distribution was proposed for the determination of the type of precipitation extremes. The GG and GNB distributions are not quite widespread, so the methods for the estimation of their parameters are, as a rule, not implemented in standard statistical packages. Therefore, the implementation of appropriate procedures requires the creation of specialized software solutions, for example, based on the functional approach, as it was done in the study described in this paper using the MATLAB programming language. However, as was demonstrated in the paper, the results of fitting such distributions to real data turned out to be better as compared to conventional models. Therefore, for processing spatial meteorological data from a large number of stations, the proposed methods and models can be effectively implemented as high-performance computing services.

**Author Contributions:** Conceptualization, V.K., A.G.; formal analysis, V.K., A.G.; funding acquisition, A.G.; investigation, A.G., V.K.; methodology, V.K., A.G.; project administration, V.K., A.G.; resources, A.G.; software, A.G.; supervision, V.K.; validation, A.G.; visualization, A.G.; writing—original draft, V.K., A.G.; writing—review and editing, A.G., V.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** The results presented in Section 8 were obtained by Andrey Gorshenin to meet the research goals of gran<sup>t</sup> No 18-71-00156 of the Russian Science Foundation.

**Acknowledgments:** Authors thank the reviewers for their valuable comments that helped to improve the presentation of the material.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Appendix A. Algorithms**

This section presents the above-described algorithms in the form pseudo-code. All of these algorithms have been implemented using Matlab programming language without any tools specific to this development environment. They can be successfully implemented, for example, with the Python programming language taking into account some minor changes.

#### **Algorithm A1.** GNB approximations


8: HISTOGRAMS(*γl*1,*μl*1,*errl*1,*γl*2,*μl*2,*errl*2,*γl*∞,*μl*∞,*errl*∞);

#### **Algorithm A2.** Stabilization of averages



 17:level=level-step;
