VMSTA Modern Stochastics: Theory and Applications 2351-60542351-60462351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA216 10.15559/22-VMSTA216 Research Article A class of fractional Ornstein–Uhlenbeck processes mixed with a Gamma distribution https://orcid.org/0000-0001-7040-0366 BianchiLuigi Amedeoluigiamedeo.bianchi@unitn.it∗ https://orcid.org/0000-0002-5209-2678 BonaccorsiStefanostefano.bonaccorsi@unitn.it https://orcid.org/0000-0001-8007-3119 TubaroLucianoluciano.tubaro@unitn.it Università degli Studi di Trento, Via Sommarive 14, 38123 Povo (Trento), Italy Corresponding author. 2023 2111202210137571410202171120227112022 © 2023 The Author(s). Published by VTeX2023 Open access article under the CC BY license.

We consider a sequence of fractional Ornstein–Uhlenbeck processes, that are defined as solutions of a family of stochastic Volterra equations with a kernel given by the Riesz derivative kernel, and leading coefficients given by a sequence of independent Gamma random variables. We construct a new process by taking the empirical mean of this sequence. In our framework, the processes involved are not Markovian, hence the analysis of their asymptotic behaviour requires some ad hoc construction. In our main result, we prove the almost sure convergence in the space of trajectories of the empirical means to a given Gaussian process, which we characterize completely.

 Fractional Ornstein–Uhlenbeck processes empirical means Gamma mixing stochastic Volterra equations generalized Wright function 60G22 60G17 Introduction

Recently, there has been a certain interest in the literature about the convergence of the empirical means of a sequence of realizations of the solution of a stochastic differential equation (SDE) with random coefficients as a mean to construct new classes of stochastic processes. For instance, K. El-Sebaiy and co-authors, in a sequence of papers [10, 11, 8] have constructed several classes of Ornstein–Uhlenbeck type processes led by a fractional Brownian motion B H ( t ) or a Hermite process Z q H ( t ) of order q ≥ 1.

The construction of this kind kind of processes was introduced in [13] in the analysis of a stationary Ornstein–Uhlenbeck process d X k ( t ) = − α k X k ( t ) d t + d W ( t ) , with initial condition X k ( 0 ) = 0, where W = { W ( t ) , t ≥ 0 } is a real standard Brownian motion defined on a stochastic basis ( Ω , F , { F t , t ∈ R } , P ) and { α k , k ≥ 1 } is a sequence of independent, identically distributed random variables, independent of W. With no loss of generality, we may assume that the sequence of α k is defined on a different probability space ( Ω ¯ , F ¯ , P ¯ ).

The empirical means are defined by Y n ( t ) = 1 n ∑ k = 1 n X k ( t ) . Under the assumption that α k are independent, identically distributed, with a Gamma distribution Γ ( μ , λ ), with shape parameter μ > 1 2 \frac{1}{2}$]]> and scale parameter λ > 00$]]>, in [13] the authors prove that P ¯-a.s. the sequence ( Y n ( t ) ) n converges to a process Y ( t ) for any t ∈ R and in C ( [ a , b ] ) for arbitrary a < b, where the process Y ( t ) is defined in law by Y ( t ) = ∫ 0 t λ t − s + λ μ d W ( s ) .

In order to extend the above construction, we recall the notion of stochastic Volterra processes defined as convolution processes between a given kernel and a Brownian motion W X ( t ) = ∫ 0 t g ( t , s ) d W ( s ) , t ≥ 0 . Such processes are related to (and somehow coincide with a subclass of) ambit processes, as introduced by Barndorff-Nielsen [3, 4], and in particular to the so-called Brownian semistationary processes, see [2].

These processes cannot be expressed in the Itô differential form (1), except in the case of an exponential kernel g ( t , s ) = e λ ( t − s ) , due to the presence of a memory term. It is however possible to describe their behaviour through an integral stochastic (Volterra) equation. Although this process is not Markovian, several properties can be studied, see for instance [6] for a survey of known results and the relation with the fractional Brownian motion.

In analogy with the existing literature, we aim to construct a class of stochastic Volterra processes via an aggregation procedure (i.e. taking the limit of the empirical means). Thanks to this approach, we are able to prove some of its properties. We shall refer to our processes as fractional Ornstein–Uhlenbeck processes, since they are given as solutions of the integral equation X k ( t ) = − α k g ρ ∗ X k ( t ) + W ( t ) , k ≥ 1 , t > 0 , 0,\]]]> with initial condition X k ( 0 ) = 0, where ∗ denotes the convolution product and g ρ ( t ) = 1 Γ ( ρ ) t ρ − 1 , t > 0 , 0,\]]]> denotes the fractional integration kernel for any ρ > 00$]]>, cf. [18]. It is well known [6, 7] that the solution of (3) is a Gaussian process defined in terms of the scalar resolvent equation s α ( t ) + α ∫ 0 t g ρ ( t − τ ) s α ( τ ) d τ = 1 , t ≥ 0 , α ≥ 0 , by means of the formula X k ( t ) = ∫ 0 t s α k ( t − τ ) d W ( τ ) . The scalar resolvent function s α ( t ) can be defined in terms of the well known Mittag-Leffler function E ρ ( t ) [15]; this link will be exploited in order to prove finer properties of the function s α ( t ), see Section 2. Notice further that the choice ρ = 1 reduces the problem to the evolution equation studied in [13]. In this case, s α ( t ) = e − α t satisfies the semigroup property, X k is a Markov process and its behaviour at infinity is much easier to study than in the case 1 < ρ ≤ 2.

Presentation of the results and outline of the paper

In this paper, we assume that α k are i.i.d., with a Gamma distribution Γ ( μ , λ ) , μ , λ > 0 ; 0\text{;}\hspace{2.5pt}\]]]> and we assume that the fractional integration parameter ρ satisfies ρ ∈ ( 1 , 2 ] .

Notice that condition (6) is quite standard in the literature, compare, e.g., [13]. In our setting, we prove that the sequence of the empirical mean processes Y n (given in (2)) converges in L 2 ( P ), for fixed time, P ¯-a.s., to a stochastic process Y ( t ) given by the convolution Wiener integral Y ( t ) = ∫ 0 t E ¯ [ s α ( t − s ) ] d W ( s ) , where E ¯, the expectation with respect to P ¯, of the resolvent scalar kernel with respect to the random variable α can be explicitly computed and related to a generalized Wright–Fox function G ρ μ ( z ), see Section 2.2 below. In Section 3, we prove that Y n ( t ) converges pathwise on [ 0 , T ], P ¯-a.s., for arbitrary T > 00$]]>. First, in Theorem 2, we prove weak convergence of the laws in the space of continuous trajectories on [ 0 , T ], by proving a suitable tightness result. Then, in Theorem 3, we prove that the convergence holds almost surely, by using a particular application of the integration by parts formula.

Finally, in Section 4 we consider the asymptotic behaviour of the limit process Y ( t ) as t → ∞. Notice that in the previous section no condition was imposed on the values of the parameters μ and λ defining the α-distribution. In this section, we shall impose the bound μ > 1 2 ρ ; \frac{1}{2\rho };\]]]> notice that this condition is coherent with the condition μ > 1 2 \frac{1}{2}$]]> imposed in [13] for the evolutionary problem ( ρ = 1).

We show that the process Y ( t ) converges in distribution, for t → ∞, to some Gaussian random variable η and we give a characterization of the limit as “stationary” solution of the problem. The emphasis here comes from the fact that we cannot use the notion of invariant measures for SDEs, since our system is not Markovian. However, if we consider the solution of a modified problem with initial time equal to − ∞, we have that this process has the law equal to that of η for all times.

 Some special functions Mittag-Leffler function

We introduce the Mittag-Leffler function E ρ ( z ) defined by the series E ρ ( z ) = ∑ k = 0 ∞ z k Γ ( ρ k + 1 ) , z ∈ C , for any ρ > 00$]]>. Notice that this is a generalization of the exponential function, obtained as a special case when ρ = 1. This function was introduced by Mittag-Leffler in [15], and since then studied by many authors. Basic references for the next results are [9, 12].

The series in (9) converges in the whole convex plane for ρ > 00$]]>. For x ∈ R, x → ∞, we get the following asymptotic expansion [12, (3.4.15)]: for any m ≥ 1 we have E ρ ( − x ) = − ∑ k = 1 m ( − x ) − k Γ [ 1 − k ρ ] + O ( | x | − m − 1 ) .

The case 0 < ρ < 1 is somehow peculiar because the Mittag-Leffler function  E ρ ( − x ) is completely monotonic, that is for all n ∈ N, for all x ∈ R + , ( − 1 ) n E ρ ( n ) ( − x ) ≥ 0. It follows in particular that E ρ ( − x ) has no real zeroes for x ∈ R + when 0 < ρ < 1. Conversely, in the region 1 < ρ < 2, the Mittag-Leffler function E ρ ( − x ) shows a damped oscillation converging to zero as x → ∞, as can be seen in Figure 1.

The graph of the Mittag-Leffler function E ρ ( − x ) with parameter ρ = 1.9. Plotted with Mathematica [20]

The global behaviour of the Mittag-Leffler function E ρ ( − x ) is established in [12, Corollary 3.1]: ∀ x > 0 : | E ρ ( − x ) | ≤ M 1 + x , 0:\hspace{1em}|{E_{\rho }}(-x)|\le \frac{M}{1+x},\]]]> for some constant M depending only on ρ.

For any ρ > 00$]]>, the derivative of the Mittag-Leffler function E ρ ( − x ) satisfies d d x E ρ ( − x ) = − 1 ρ E ρ , ρ ( − x ) , where E ρ , ρ ( z ) is the two-parametric Mittag-Leffler function.

Recall that, by definition, the two-parametric Mittag-Leffler function is (see, for example, [12]), E ρ , τ ( z ) = ∑ k = 0 ∞ z k Γ ( ρ k + τ ) , z ∈ C , ℜ ( ρ ) > 0 , τ ∈ C . 0,\tau \in \mathbb{C}.\]]]> As for the original Mittag-Leffler function E ρ ( z ), also E ρ , ρ ( z ) is an entire function of order ω = 1 ρ . Moreover, it is proved in [17, (1.2.5)] that a similar asymptotic expansion to the one given in (10) for the Mittag-Leffler function holds for the two-parametric Mittag-Leffler function: E ρ , ρ ( − x ) = − ∑ k = 1 m ( − x ) − k Γ [ ρ − k ρ ] + O ( | x | − m − 1 ) . However, since 1 Γ [ − n ] = 0 for any integer n ≥ 1, we see that the previous estimate becomes, for any m ≥ 1, E ρ , ρ ( x ) = x − 1 + O ( | x | − m − 1 ) . We are thus led to the following global estimate for the behaviour of the two-parametric Mittag-Leffler function E ρ , ρ ( − x ): ∀ x > 0 : | E ρ , ρ ( − x ) | ≤ M 2 1 + x , 0:\hspace{1em}|{E_{\rho ,\rho }}(-x)|\le \frac{{M_{2}}}{1+x},\]]]> for some constant M 2 depending only on ρ.

 A generalized Wright–Fox type function

In this paper, we will encounter the following special function, G ρ μ ( z ) = ∑ k = 0 ∞ a k z k = ∑ k = 0 ∞ μ k Γ ( k ρ + 1 ) z k , where ( μ ) n = Γ ( μ + n ) Γ ( μ ) is the Pochhammer symbol. This function has some similarities to the Mittag-Leffler function, and in fact both are special cases of the generalized Wright–Fox type function, defined as p Ψ q ( a 1 , A 1 ) , … , ( a p , A p ) ( b 1 , B 1 ) , … , ( b q , B q ) ; z = ∑ k = 0 ∞ ∏ i = 1 p Γ ( a i + A i k ) ∏ j = 1 q Γ ( b j + B j k ) z k k ! , where the series converges, with a i , b j ∈ C and A i , B j ∈ R, see [14]. We recover both the two-parametric Mittag-Leffler function E ρ , τ ( z ) and G ρ μ ( z ) in p Ψ q with the following choices: E ρ , τ ( z ) = 1 Ψ 1 ( 1 , 1 ) ( τ , ρ ) ; z G ρ μ ( z ) = 2 Ψ 1 ( 1 , 1 ) , ( μ , 1 ) ( 1 , ρ ) ; z .

In the next lemma, we prove that G ρ μ ( z ) is well defined and that it is an entire function, as was E ρ .

Assume that ρ > 11$]]>. Then the function G ρ ( z ) is an entire function of the variable z ∈ C of order ω = 1 ρ − 1 .

Since G ρ μ is defined in terms of a power series, to prove that it is an entire function it is sufficient to study the radius of convergence. By the ratio test, it is sufficient to study the behaviour of R = lim n → ∞ a n a n + 1 = lim n → ∞ Γ ( μ ) Γ ( ρ n + ρ + 1 ) Γ ( μ + n + 1 ) Γ ( μ + n ) Γ ( μ ) Γ ( ρ n + 1 ) = 1 μ + n Γ ( ρ n + ρ + 1 ) Γ ( ρ n + 1 ) . Stirling’s series can be used to find the asymptotic behaviour of the last ratio [19]: Γ ( ρ n + 1 ) Γ ( ρ n + ρ + 1 ) = ( ρ n ) − ρ 1 + O ( n − 1 ) , therefore, by taking the limit R = lim n → ∞ a n a n + 1 = lim n → ∞ n ρ μ + n , we see that R is infinite provided ρ > 11$]]>.

By using again Stirling’s approximation, it is possible to compute the order ω of the entire function G ρ μ ( z ). Recall that ω is the infimum of all κ such that max | z | = r | G ρ μ ( z ) | < e r κ , eventually for r → ∞. For the entire function G ρ μ ( z ) = ∑ a k z k , the order is given by the formula ω = lim sup k → ∞ k log ( k ) log ( 1 / | a k | ) ; for large k, we have 1 a k ≈ e − k ρ ( k ρ ) k ρ + 1 / 2 e − k ( k ) k + 1 2 + μ , log ( 1 / a k ) ≈ ( ρ − 1 ) k log ( k ) + ( 1 − ρ ) k + ( k ρ + 1 / 2 ) log ( ρ ) − μ log ( k ) , which implies ω = 1 ρ − 1 .  □

Further properties of the function G ρ μ will be proved in the sequel. In particular, as a corollary of the computation in Lemma 7, we obtain the following asymptotic behaviour for large x: G ρ μ ( − x ) ≤ C 1 x + 1 x μ , x → ∞ , for μ ≠ 1. A plot of he restriction of G ρ μ to the negative real numbers can be seen in Figure 2.

The graph of the function G ρ μ ( − x ) with parameters μ = 4 and ρ = 1.9. Plotted with Mathematica [20]

 The empirical mean process

We consider in this section the solution of the fractional evolution equation (3) and the properties of the empirical mean process Y n ( t ) defined in (2).

The relation between Mittag-Leffler function and Riemann–Liouville fractional derivative is specified through the following relation [12, (3.7.44)] E ρ ( λ t ρ ) − λ ∫ 0 t E ρ ( λ s ρ ) g ρ ( t − s ) d s = 1 . By comparing it with the definition of the scalar resolvent function s α ( t ) in (4) we obtain the relation s α ( t ) = E ρ ( − α t ρ ) , t > 0 , for all α > 0 . 0,\hspace{2.5pt}\text{for all}\hspace{2.5pt}\alpha >0\text{.}\]]]>

For the sake of completeness, let us recall the behaviour of the Gamma distribution.

Assume that α is a random variable on the probability space ( Ω ¯ , F ¯ , P ¯ ) with a Gamma distribution of parameters μ and λ, i.e. α has the density f ( x ∣ μ , λ ) = 1 Γ ( μ ) λ e − λ x λ x μ − 1 x > 0 . 0.\]]]> Then the moments of α are given by E ¯ [ α n ] = 1 λ n μ n .

We apply this proposition to compute the expected value (with respect to the measure P ¯) of the random process s α ( t ).

Assume that α has a Gamma distribution of parameters μ and λ on the probability space ( Ω ¯ , F ¯ , P ¯ ) . Then the expected value of s α ( t ) satisfies E ¯ [ s α ( t ) ] = G ρ μ ( − t ρ / λ ) = ∑ k = 0 ∞ μ k Γ ( k ρ + 1 ) − t ρ λ k .

Let us consider the representation of the scalar resolvent kernel s α ( t ) in terms of the Mittag-Leffler function E ρ ( x ) given in (15) and the definition of the Mittag-Leffler function as an entire function given in (9): E ¯ [ s α ( t ) ] = E ¯ ∑ k = 0 ∞ ( − 1 ) k Γ ( k ρ + 1 ) α k t k ρ . By an application of Fubini’s theorem we obtain E ¯ [ s α ( t ) ] = ∑ k = 0 ∞ ( − 1 ) k Γ ( k ρ + 1 ) E ¯ α k t k ρ , and by (16) we obtain E ¯ [ s α ( t ) ] = ∑ k = 0 ∞ ( − 1 ) k Γ ( k ρ + 1 ) μ k t ρ λ k .  □

The limit of the empirical mean process

In this section, we will study the asymptotic behaviour of the empirical mean process Y n ( t ) = ∫ 0 t 1 n ∑ k = 1 n s α k ( t − τ ) d W ( τ ) . By the law of large numbers 1 n ∑ k = 1 n s α k ( t ) ⟶ E ¯ [ s α ( t ) ] , t > 0 , 0,\]]]> hence a natural candidate for the limit process of Y n ( t ) is given by Y ( t ) = ∫ 0 t G ρ μ − ( t − s ) ρ λ d W ( s ) .

The sequence Y n ( t ) converges in L 2 ( Ω , P ), uniformly in [ 0 , T ], to Y ( t ), P ¯-a.s.

In order to prove the convergence of Y n to Y, we first notice that f n ( t ) = 1 n ∑ k = 1 n s α k ( t ) satisfies, by estimate (11), | f n ( t ) | ≤ M ; therefore, the same applies to the function G ρ μ ( − t ρ / λ ), and Y ( t ) is a Gaussian process.

Moreover, the function s α ( t ) is continuous and bounded for α ∈ R + and t ∈ [ 0 , T ], for arbitrary T > 00$]]>, then the uniform law of large numbers [16, Lemma 2.4] implies that G ρ μ ( − t ρ / λ ) is continuous and bounded in t ∈ [ 0 , T ] and sup t ∈ [ 0 , T ] 1 n ∑ k = 1 n s α k ( t ) − E ¯ [ s α ( t ) ] → 0 P ¯ -a.s. It follows from the dominated convergence theorem that E sup t ∈ [ 0 , T ] | Y n ( t ) − Y ( t ) | 2 ≤ 4 ∫ 0 T 1 n ∑ k = 1 n s α k ( t ) − E ¯ [ s α ( t ) ] 2 d t → 0 P ¯ -a.s.  □

 Weak convergence of the empirical mean process in <inline-formula id="j_vmsta216_ineq_130"><alternatives><mml:math> <mml:mi mathvariant="italic">C</mml:mi> <mml:mo mathvariant="normal" fence="true" stretchy="false">(</mml:mo> <mml:mo fence="true" stretchy="false">[</mml:mo> <mml:mn>0</mml:mn> <mml:mo mathvariant="normal">,</mml:mo> <mml:mi mathvariant="italic">T</mml:mi> <mml:mo fence="true" stretchy="false">]</mml:mo> <mml:mo mathvariant="normal" fence="true" stretchy="false">)</mml:mo></mml:math><tex-math><![CDATA[$C([0,T])$]]></tex-math></alternatives></inline-formula>

We have proved before that Y n ( t ) and Y ( t ) are, P ¯-a.s., Gaussian, continuous stochastic processes on ( Ω , F , P ).

In the remaining part of this section, the whole construction is thought to hold P ¯-a.s., meaning that we consider α given and we are in the space ( Ω , F , P ). We aim to consider the P-weak convergence of Y n to Y in the space C ( [ 0 , T ] ): Y n ⇒ Y if and only if ϕ ( Y n ) → ϕ ( Y ) in distribution, for every bounded Borel measurable ϕ : C ( [ 0 , T ] ) → R .

The sequence Y n P-weakly converges to the process Y in the space C ( [ 0 , T ] ), P ¯-almost surely.

By [5, Theorems 12.3 and 8.1], weak convergence of Y n to Y follows if we prove that the following two conditions hold:

convergence of the finite dimensional distributions, and

the moment condition E [ | Y n ( t ) − Y n ( s ) | κ ] ≤ C ( κ , δ ) | t − s | 1 + δ , where C ( κ , δ ) depends on the constants κ and δ only, compare [5, formula (12.51)].

We shall prove each of the above-described conditions in the following two lemmas, which therefore together imply the thesis of the theorem.

For any k ∈ N, t 1 , … , t k ∈ [ 0 , T ], x 1 , … , x k ∈ R, it holds that P ( Y n ( t 1 ) ≤ x 1 , … , Y n ( t k ) ≤ x k ) converges to P ( Y ( t 1 ) ≤ x 1 , … , Y ( t k ) ≤ x k ) as n → ∞.

Both Y n and Y are Gaussian processes, so their finite dimensional distributions are fully determined by their mean and covariance. For any k ∈ N, t 1 , … , t k ∈ [ 0 , T ], the mean and covariance of Y n ( t 1 ) , … , Y n ( t k ) converge to mean and covariance of Y ( t 1 ) , … , Y ( t k ), thanks to Theorem 1 and the pointwise (in t) convergence of Y n ( t ) → Y ( t ).  □

The moment condition (17) holds, and the family of empirical mean processes { Y n } n is tight.

First, we notice that by exploiting the Gaussianity of Y n , it is sufficient to search an estimate for κ = 2. Further, we can reduce the problem to the estimate of a single X k , since E [ | Y n ( t ) − Y n ( s ) | 2 ] ≤ 1 n ∑ k = 1 n E [ | X k ( t ) − X k ( s ) | 2 ] . Next, we write (for simplicity, we drop the index k in the following computation) E [ | X ( t ) − X ( s ) | 2 ] = E ∫ 0 s [ s α ( t − u ) − s α ( s − u ) ] d W ( u ) + ∫ s t s α ( t − u ) d W ( u ) 2 ≤ 2 ∫ 0 s [ s α ( t − u ) − s α ( s − u ) ] 2 d u + 2 ∫ s t [ s α ( t − u ) ] 2 d u ≤ 2 ∫ 0 s [ s α ( t − s + u ) − s α ( u ) ] 2 d u + 2 ∫ 0 t − s [ s α ( u ) ] 2 d u . Recall that, by estimate (11), | s α ( t ) | = | E ρ ( − α t ρ ) | ≤ M, hence the second integral in the right hand side of previous estimate is bounded by 2 M 2 ( t − s ). As far as the first term is concerned, we get ∫ 0 s [ s α ( t − s + u ) − s α ( u ) ] 2 d u = ∫ 0 s [ E ρ ( − α ( t − s + u ) ρ ) − E ρ ( − α u ρ ) ] 2 d u = [ v = α 1 / ρ u ] α − 1 / ρ ∫ 0 α 1 / ρ s [ E ρ ( − ( α 1 / ρ ( t − s ) + v ) ρ ) − E ρ ( − v ρ ) ] 2 d v = α − 1 / ρ ∫ 0 α 1 / ρ s [ s 1 ( α 1 / ρ ( t − s ) + v ) − s 1 ( v ) ] 2 d v = α − 1 / ρ ∫ 0 α 1 / ρ s ∫ v v + α 1 / ρ ( t − s ) s 1 ′ ( x ) d x 2 d v . The derivative of the scalar resolvent kernel can be written in terms of the two-parametric Mittag-Leffler function, see formula (12) s 1 ′ ( x ) = d d x E ρ ( − x ρ ) = − x ρ − 1 E ρ , ρ ( − x ρ ) . Recalling the bound (13) we see that s 1 ′ ( x ) ≤ x ρ − 1 M 2 1 + x ρ ≤ M 3 1 + x , for some constant M 3 .

In conclusion, we get E [ | X ( t ) − X ( s ) | 2 ] ≤ [ 2 M 3 2 α 2 / ρ T 2 + 2 M 2 ] ( t − s ) . This implies E [ | Y n ( t ) − Y n ( s ) | 2 ] ≤ K n ( t − s ) , where the constant K n can be explicitly given as K n = 2 M 3 2 T 2 1 n ∑ k = 1 n ( α k ) 2 / ρ + 2 M 2 . Notice that, P ¯-almost surely, K n → n → ∞ K ¯ = 2 M 3 2 T 2 E ¯ [ α 2 / ρ ] + 2 M 2 . It follows that for some n 0 ∈ N, K n ≤ K ¯ + 1 for any n ≥ n 0 . Set C ( 2 , 0 ) = max n ∈ { 1 , … , n 0 } K n ∨ ( K ¯ + 1 ) . Since Y n ( t ) − Y n ( s ) is a centred Gaussian random variable, by standard arguments we conclude that E [ | Y n ( t ) − Y n ( s ) | 2 n ] ≤ C ( 2 n , n − 1 ) | t − s | n and C ( 2 n , n − 1 ) = ( 2 n − 1 ) ! ! C ( 2 , 0 ) n .  □

 Almost sure convergence of the empirical mean process in <inline-formula id="j_vmsta216_ineq_173"><alternatives><mml:math> <mml:mi mathvariant="italic">C</mml:mi> <mml:mo mathvariant="normal" fence="true" stretchy="false">(</mml:mo> <mml:mo fence="true" stretchy="false">[</mml:mo> <mml:mn>0</mml:mn> <mml:mo mathvariant="normal">,</mml:mo> <mml:mi mathvariant="italic">T</mml:mi> <mml:mo fence="true" stretchy="false">]</mml:mo> <mml:mo mathvariant="normal" fence="true" stretchy="false">)</mml:mo></mml:math><tex-math><![CDATA[$C([0,T])$]]></tex-math></alternatives></inline-formula>

In this section, we consider the convergence of the paths of Y n to Y in C ( [ 0 , T ] ) with respect to the uniform topology.

In our assumptions, we have that lim n → ∞ sup t ∈ [ 0 , T ] | Y n ( t ) − Y ( t ) | → 0 , almost surely with respect to P × P ¯.

Our proof will follow from a somehow broader result. Let f n be a sequence of smooth functions on [ 0 , T ] that converges (at least pointwise) to f. We shall see later which further assumptions are required. Then we compute ( f n ∗ W ˙ ) ( t ) − ( f ∗ W ˙ ) ( t ) = f n ( t − s ) W ( s ) | s = 0 s = t + ( W ∗ f ˙ n ) ( t ) − f ( t − s ) W ( s ) | s = 0 s = t − ( W ∗ f ˙ ) ( t ) = [ f n ( 0 ) − f ( 0 ) ] W ( t ) + ( W ∗ ( f n − f ) ˙ ) ( t ) . It follows that there exists a constant C finite P-a.s. such that sup t ∈ [ 0 , T ] ( f n ∗ W ˙ ) ( t ) − ( f ∗ W ˙ ) ( t ) ≤ C | f n ( 0 ) − f ( 0 ) | + ∫ 0 T f n ˙ ( s ) − f ˙ ( s ) d s . Therefore, conditions sufficient for the required convergence are as follows:

f n ( 0 ) → f ( 0 ),

f n ˙ ( s ) → f ˙ ( s ) pointwise in [ 0 , T ], and

sup n ∈ N sup s ∈ [ 0 , T ] | f n ˙ ( s ) | ≤ C.

 Then an application of Lebesgue’s dominated convergence theorem shows that the following lemma holds.

Under conditions 1. to 3., the stochastic convolution processes X n = ( f n ∗ W ˙ ) converge in C ( [ 0 , T ] ) (endowed with the uniform topology) to the process X = ( f ∗ W ˙ ).

We are now ready to prove the main result of this section.

With the notation f n ( t ) = 1 n ∑ k = 1 n s α k ( t ) , f ( t ) = G ρ μ ( − t ρ / λ ) , introduced above, we see that f n ˙ ( t ) = 1 n ∑ k = 1 n d d t s α k ( t ) , and by using the identity s α ( t ) = E ρ ( − α t ρ ) and formula (12) we compute d d t E ρ ( − α t ρ ) = − α t ρ − 1 E ρ , ρ ( − α t ρ ) , where E ρ , ρ ( x ) is the two-parametric Mittag-Leffler function; recalling the bound (13), we get d d t E ρ ( − α t ρ ) ≤ α t ρ − 1 M 2 1 + α t ρ , which implies the bound d d t E ρ ( − α t ρ ) ≤ M α 1 / ρ sup x > 0 x ρ − 1 1 + x ρ ︸ ≤ 1 . 0}{\sup }\frac{{x^{\rho -1}}}{1+{x^{\rho }}}}}.\]]]> Since α has a Gamma distribution of parameters μ and λ, we compute E ¯ [ α 1 / ρ ] = 1 Γ ( μ ) ∫ 0 ∞ x 1 / ρ λ e − λ x λ x μ − 1 d x = 1 λ 1 / ρ Γ ( μ ) ∫ 0 ∞ λ e − λ x λ x μ + 1 / ρ − 1 d x = Γ ( μ + 1 / ρ ) λ 1 / ρ Γ ( μ ) .

Therefore, we get sup t ∈ [ 0 , T ] f n ˙ ( t ) ≤ M 1 n ∑ k = 1 n α k 1 / ρ ; by the law of large numbers in the space ( Ω ¯ , P ¯ ) we obtain that 1 n ∑ k = 1 n α k 1 / ρ → E ¯ [ α 1 / ρ ] < ∞ , hence there exists a constant C = C ( ω ¯ ), finite P ¯-almost surely, such that sup n ∈ N sup t ∈ [ 0 , T ] f n ˙ ( t ) ≤ C .

By another application of the law of large numbers we get that the sequence f n ˙ ( t ) = 1 n ∑ k = 1 n d d t s α k ( t ) converges to G ρ , ρ μ ( t ) = E ¯ − α t ρ − 1 E ρ , ρ ( − α t ρ ) . Since G ρ μ ( t ρ / λ ) = E ¯ E ρ ( − α t ρ ) , taking the derivative inside the expectation we obtain that d d t G ρ μ ( t ρ / λ ) = G ρ , ρ μ ( t ) , that is the pointwise convergence required in condition 2.

Finally, since f n ( 0 ) = 1 = f ( 0 ) for every n ∈ N, we conclude by applying Lemma 6.  □

 Asymptotic behaviour and stationary solutions

In this section we discuss the asymptotic behaviour as t → ∞ of the stochastic process Y ( t ) defined in (7). We already know that Y ( t ) is a Gaussian process with zero mean and variance given by σ t 2 = ∫ 0 t G ρ μ ( − u ρ / λ ) 2 d u = ∫ 0 t E ¯ [ s α ( u ) ] 2 d u . We are interested in the convergence of the laws L ( Y ( t ) ) as t → ∞.

A standard assumption in this section will be the following: μ > 1 2 ρ . \frac{1}{2\rho }.\]]]>

By a direct computation, we see that the random variable α − 1 / ρ belongs to L 2 ( Ω ¯ , P ¯ ) if and only if μ − 2 ρ > 00$]]>, i.e. condition (18) holds.

Assume that α has a Gamma distribution with parameters μ > 1 2 ρ \frac{1}{2\rho }$]]> and λ > 00$]]>. There exists a Gaussian random variable η such that Y ( t ) converges in law to η as t → ∞.

We first extend the Wiener process W ( t ) to a process with time t ∈ R. Let W ˜ ( t ) be a standard Brownian motion on the space ( Ω , F , P ), independent of W ( t ), and set for t ≥ 0 W ( − t ) = W ˜ ( t ).

Now, for any − s < 0, define the stochastic convolution process Y − s ( t ) = ∫ − s t E ¯ [ s α ( t − u ) ] d W ( u ) . It is easily seen that L ( Y − s ( 0 ) ) = L ( Y ( s ) ), for s ≥ 0.

In next lemma we prove that { Y − s ( 0 ) } forms a Cauchy sequence in the space L 2 ( Ω , P ).

Under our assumptions, for any ε > 00$]]> there exists T > 00$]]> such that, for any t > s > Ts>T$]]>, E Y − t ( 0 ) − Y − s ( 0 ) 2 ≤ ε.

Recall that Y − t ( 0 ) = ∫ − t 0 E ¯ [ s α ( − u ) ] d W ( u ) is a centred Gaussian process, where s α ( t ) = E ρ ( − α t ρ ) satisfies the bound | s α ( t ) | ≤ M 1 + α t ρ . In order to prove the mean square convergence of Y − t ( 0 ) as t → ∞, we need to prove that the quantity E [ Y − t ( 0 ) − Y − s ( 0 ) 2 ] is bounded by a quantity that goes to 0 as T → ∞, for t > s > Ts>T$]]>.

We estimate E [ Y − t ( 0 ) − Y − s ( 0 ) 2 ] = ∫ − s − t E ¯ [ s α ( − u ) ] 2 d u ∫ s t E ¯ s α ( u ) 2 d u ∫ s t E ¯ M 1 + α u ρ 2 d u . Since α has a Gamma distribution with parameters μ and λ, we have E [ | Y − t ( 0 ) − Y − s ( 0 ) | 2 ] ≤ ∫ s t ∫ 0 ∞ M 1 + y u ρ λ ( λ y ) μ − 1 e − λ y d y 2 d u = [ z = λ y ] ∫ s t 1 Γ ( μ ) ∫ 0 ∞ M 1 + z u ρ λ z μ − 1 e − z d z 2 d u ≤ M 2 Γ ( μ ) 2 ∫ s t ∫ 0 1 1 1 + z u ρ λ z μ − 1 e − z d z + ∫ 1 ∞ 1 1 + z u ρ λ z μ − 1 e − z d z 2 d u ≤ M 2 Γ ( μ ) 2 ∫ s t ∫ 0 1 1 1 + z u ρ λ z μ − 1 d z + ∫ 1 ∞ 1 u ρ λ z μ − 1 e − z d z 2 d u ≤ M 2 Γ ( μ ) 2 ∫ s t F 2 , 1 ( 1 , μ , 1 + μ , − u ρ λ ) + Γ ( μ ) λ u ρ 2 d u , where F 2 , 1 is the hypergeometric function; notice that, for large u, we have F 2 , 1 ( 1 , μ , 1 + μ , − u ρ λ ) ∼ c λ , μ u − μ ρ , 0 < μ < 1 , c λ , 1 log ( u ) u , μ = 1 , c λ , μ u − ρ , μ > 1 , 1,\end{array}\right.\]]]> hence E [ Y − t ( 0 ) − Y − s ( 0 ) 2 ] ≤ C M , μ , λ , ρ ∫ s t u − 2 ρ 1 μ > 1 + log ( u ) 2 u 2 1 μ = 1 + u − 2 μ ρ 1 0 < μ < 1 d u 1}}+\frac{\log {(u)^{2}}}{{u^{2}}}{\mathbb{1}_{\mu =1}}+{u^{-2\mu \rho }}{\mathbb{1}_{0<\mu <1}}\right\}\hspace{0.1667em}\mathrm{d}u\right)\end{array}\]]]> and the right hand side is bounded, for every μ > 1 2 ρ \frac{1}{2\rho }$]]>, for every t > s > Ts>T$]]>, by C T 1 − 2 ρ ( μ ∧ 1 ) (for μ = 1, the bound is given by ( 2 + log ( T ) ) 2 T ); under the assumption μ > 1 2 ρ \frac{1}{2\rho }$]]>, the exponent is always negative, hence the quantity E [ Y − t ( 0 ) − Y − s ( 0 ) 2 ] can be taken arbitrarily small by choosing T large enough.  □

The result follows from the previous lemma. We have seen that { Y − t ( 0 ) } is a Cauchy sequence, hence it converges in L 2 ( Ω , P ) to some random variable η.

Now, a simple modification of the computations leading to Lemma 7 implies that sup t > 0 E | Y − t ( 0 ) | 2 ≤ C < ∞ , 0}{\sup }\mathbb{E}\left[|{Y_{-t}}(0){|^{2}}\right]\le C<\infty ,\]]]> and the stability of the Gaussian distribution with respect to the convergence in law implies that η has a Gaussian distribution as well, with variance σ 2 = sup t > 0 E | Y − t ( 0 ) | 2 = lim t → ∞ E | Y − t ( 0 ) | 2 . 0}{\sup }\mathbb{E}\left[|{Y_{-t}}(0){|^{2}}\right]=\underset{t\to \infty }{\lim }\mathbb{E}\left[|{Y_{-t}}(0){|^{2}}\right].\]]]> But, as already stated above, Y ( t ) ∼ Y − t ( 0 ), and we obtain the thesis.  □

Let us consider the following modification of our setting. Let W ( t ) be a two-side Wiener process, and define ξ k ( t ) = ∫ − ∞ t s α k ( t − s ) d W ( s ) , where, as opposite to (5), we are taking the integral on an infinite time interval ( − ∞ , t ).

The process ξ k ( t ) is a stationary Gaussian process with L ( ξ k ( t ) ) ∼ N ( 0 , ‖ s α ( · ) ‖ L 2 ( R + ) 2 ) , P ¯-almost surely.

The only point worth some care is the fact that the variance is finite. This follows from the computation in [6], see in particular Remark 3.2 therein, where it is proved that ‖ s α ( · ) ‖ L 2 ( R + ) 2 ≤ c 1 α 1 / ρ .  □

Assume that condition (18) holds. Then the process η ( t ) = ∫ − ∞ t G ρ μ ( − ( t − s ) ρ / λ ) d W ( s ) is a stationary Gaussian process with L ( η ( t ) ) = L ( η ) ∼ N ( 0 , σ 2 ), where η is the random variable defined in Theorem 4.

This follows from the computation leading to the proof of Lemma 7 (compare also to Remark 2 for the necessity of condition (18)).  □

Consider the empirical mean processes η n ( t ), defined by η n ( t ) = 1 n ∑ k = 1 n ξ k ( t ) , t ∈ R . Then, for every t ∈ R, the sequence η n ( t ) converges in L 2 ( Ω , P ) as n → ∞ to the random variable η ( t ) defined in (19), P ¯-almost surely.

Before we provide the proof, let us recall a few known results about convergence of a sequence of Gaussian random variables (see, e.g., [1]).

Assume that { X n } is a sequence of Gaussian random variables converging in law to a random variable X. Then X is also a Gaussian random variable.

Assume that X n ∼ N ( m n , σ n 2 ) for every n and that m n → m and σ n 2 → σ 2 as n → ∞. Then X n converges in law to X ∼ N ( m , σ 2 ).

Assume further that σ 2 = 0. Then the convergence holds also in probability and in L 2 , hence the convergence holds in L p for every p > 00$]]>.

By definition, the difference η ( t ) − η n ( t ) has a Gaussian distribution with zero mean and variance δ n 2 = ∫ − ∞ t G ρ μ ( − ( t − u ) ρ / λ ) − 1 n ∑ k = 1 n s α k ( t − u ) 2 d u = ∫ 0 ∞ E ¯ [ s α ( u ) ] − 1 n ∑ k = 1 n s α k ( u ) 2 d u . Taking into account the results in Lemma 8, the claim of the theorem follows by proving that lim n → ∞ ∫ 0 ∞ 1 n ∑ k = 1 n s α k ( t ) − E ¯ [ s α ( t ) ] 2 d t = 0 holds P ¯-a.s.

As opposed to Theorem 1, here we deal with an infinite horizon convergence problem. However, the uniform law of large numbers [16, Lemma 2.4] implies that for every m ∈ N there exists a negligible set N m ⊂ Ω ¯ such that the convergence sup t ∈ [ m , m + 1 ] 1 n ∑ k = 1 n s α k ( t ) − E ¯ [ s α ( t ) ] → 0 holds outside N m .

By setting N = ⋃ m N m we have P ¯ ( N ) = 0 and the pointwise convergence 1 n ∑ k = 1 n s α k ( t ) − E ¯ [ s α ( t ) ] → 0 for t ≥ 0, P ¯-a.s.

It remains to prove a dominated convergence theorem which holds P ¯-a.s. For this purpose, let us define f n ( t ) = 1 n ∑ k = 1 n s α k ( t ) , which satisfies, by estimate (11), | f n ( t ) | ≤ 1 n ∑ k = 1 n M 1 + α k t ρ . For simplicity, we split the estimate in two parts as follows: | f n ( t ) | ≤ M , t ∈ [ 0 , 1 ] , and, for some ε > 00$]]> to be chosen later, an explicit computation leads to | f n ( t ) | ≤ C ( M , ε , ρ ) 1 n ∑ k = 1 n α k − ( 1 + ε ) / 2 ρ t − ( 1 + ε ) / 2 , t ≥ 1 , where C ( M , ε , ρ ) is a constant depending only on the stated quantities. The sequence of random variables { α k − ( 1 + ε ) / 2 ρ } satisfies the law of large numbers provided that each of its term is integrable, i.e. under the condition 0 < ε < 2 μ ρ − 1 . This is a nonempty interval thanks to the inequality in (8).

Therefore, outside an event of zero measure we have, eventually in n, | f n ( t ) | ≤ f ( t ) : = C ( M , ε , ρ ) E ¯ [ α − ( 1 + ε ) / 2 ρ ] + 1 t − ( 1 + ε ) / 2 , t ≥ 1 , and the function f ( t ) is square integrable on R + .

In order to conclude the proof, it remains to apply the dominated convergence theorem to the integrals in (20). The sequence of functions 1 n ∑ k = 1 n s ( α k , t ) − E ¯ [ s ( α , t ) ] 2 converges pointwise to 0 and its dominated by the integrable function | f ( t ) | 2 + | G ρ μ ( − t ρ / λ ) | 2 (compare with (14)). Then Lebesgue’s dominated convergence theorem implies that the limit in (20) holds.  □

The next, and final, step is to prove that the weak convergence of the sequence of processes η n to η holds in the space C ( [ 0 , T ] ) of continuous functions on [ 0 , T ], for fixed T > 00$]]>.

The sequence η n P-weakly converges to the process η in the space C ( [ 0 , T ] ), P ¯-almost surely.

Notice that the proof is quite similar to that of Theorem 2: a consequence of the convergence of finite dimensional distributions and tightness of the laws. In the sequel we shall only emphasise the differences in computations.

The convergence of the finite dimensional distributions follows from the fact that all the processes involved are Gaussian and thus we have established L 2 and pointwise convergence of η n ( t ) to η ( t ) for all t, as in Lemma 4.

The next result is the analogue of Lemma 5 and deals with the tightness estimate corresponding to (17).

There exist constants δ , κ > 00$]]> such that for all s < t ∈ [ 0 , T ] E [ | η n ( t ) − η n ( s ) | κ ] ≤ C ( κ , δ ) | t − s | 1 + δ , where C ( κ , δ ) depends on the constants κ and δ only.

By an analysis of the proof of Lemma 5, we see that the key point is to estimate the quantity E [ | ξ ( t ) − ξ ( s ) | 2 ] = E ∫ − ∞ s [ s α ( t − u ) − s α ( s − u ) ] d W ( u ) + ∫ s t s α ( t − u ) d W ( u ) 2 ≤ 2 ∫ 0 ∞ [ s α ( t − s + u ) − s α ( u ) ] 2 d u + 2 ∫ 0 t − s [ s α ( u ) ] 2 d u and, in particular, the first term, where now the integration interval is  R + . As far as the first term is concerned, we get ∫ 0 ∞ [ s α ( t − s + u ) − s α ( u ) ] 2 d u = α − 1 / ρ ∫ 0 ∞ ∫ v v + α 1 / ρ ( t − s ) s 1 ′ ( x ) d x 2 d v . We recall the bound s 1 ′ ( x ) ≤ x ρ − 1 M 2 1 + x ρ ≤ M 3 1 + x , and divide the integral in two parts, for small v (say, v ∈ [ 0 , 1 ]) and for large v (say v ≥ 1). We thus estimate the right hand side of the previous estimate by R = α − 1 / ρ ∫ 0 1 M 3 α 1 / ρ ( t − s ) 2 d v + M 3 2 α − 1 / ρ ∫ 1 ∞ log ( v + α 1 / ρ ( t − s ) ) − log ( v ) 2 d v . Recalling that log ( 1 + x ) ≤ x for all x > 00$]]>, R ≤ M 3 2 α 1 / ρ ( t − s ) 2 + M 3 2 α − 1 / ρ ∫ 1 ∞ α 1 / ρ ( t − s ) v 2 d v ≤ 2 M 3 2 α 1 / ρ ( t − s ) 2 . In conclusion, we get E [ | ξ ( t ) − ξ ( s ) | 2 ] ≤ [ 2 M 3 2 α 2 / ρ T 2 + 2 M 2 ] ( t − s ), which further implies E [ | η n ( t ) − η n ( s ) | 2 ] ≤ K n ( t − s ) and, finally, we obtain E [ | η n ( t ) − η n ( s ) | 2 n ] ≤ C n | t − s | n .  □

 Acknowledgments

We would like to thank Prof. Luisa Beghin for a fruitful discussion and for pointing us towards Wright–Fox generalized functions, as well as the anonymous referees for their comments and remarks, which helped us in improving the paper.

 References Baldi, P.: Stochastic Calculus: An Introduction Through Theory and Exercises. Universitext. Springer, Cham (2017). MR3726894. https://doi.org/10.1007/978-3-319-62226-2 Barndorff-Nielsen, O.E., Schmiegel, J.: Brownian semistationary processes and volatility/intermittency. In: Advanced Financial Modelling, pp. 1–26. De Gruyter (2009). MR2648456. https://doi.org/10.1515/9783110213140.1 Barndorff-Nielsen, O.E., Schmiegel, J.: Ambit Processes; with Applications to Turbulence and Tumour Growth. In: Benth, F.E., Di Nunno, G., Lindstrøm, T., Øksendal, B., Zhang, T. (eds.) Stochastic Analysis and Applications. Abel Symposia, pp. 93–124. Springer, Berlin, Heidelberg (2007). MR2397785. https://doi.org/10.1007/978-3-540-70847-6_5 Barndorff-Nielsen, O.E., Benth, F.E., Veraart, A.E.D.: Ambit Processes and Stochastic Partial Differential Equations. In: Di Nunno, G., Øksendal, B. (eds.) Advanced Mathematical Methods for Finance, pp. 35–74. Springer, Berlin, Heidelberg (2011). MR2752540. https://doi.org/10.1007/978-3-642-18412-3_2 Billingsley, P.: Convergence of Probability Measures, 1st edn. John Wiley & Sons, Inc., New York-London-Sydney (1968). MR0233396 Bonaccorsi, S., Tubaro, L.: Mittag-Leffler’s Function and Stochastic Linear Volterra Equations of Convolution Type. Stoch. Anal. Appl. 21(1), 61–78 (2003). MR1954075. https://doi.org/10.1081/SAP-120017532 Clément, P., Da Prato, G.: Some results on stochastic convolutions arising in Volterra equations perturbed by noise. Atti Accad. Naz. Lincei, Rend. Lincei, Mat. Appl. 7(3), 147–153 (1996) MR1454409 Douissi, S., Es-Sebaiy, K., Tudor, C.A.: Hermite Ornstein–Uhlenbeck processes mixed with a Gamma distribution. Publ. Math. (Debr.) 96(1-2), 23–44 (2020). MR4062581. https://doi.org/10.5486/pmd.2020.8443 Erdélyi, A., Magnus, W., Oberhettinger, F., Tricomi, F.: Higher Trascendental Functions vol. III. McGraw-Hill, New York (1955) MR0717723 Es-Sebaiy, K., Tudor, C.A.: Fractional Ornstein-Uhlenbeck Processes Mixed with a Gamma Distribution. Fractals 23(03), 1550032 (2015). MR3375690. https://doi.org/10.1142/S0218348X15500322 Es-Sebaiy, K., Farah, F.-E., Hilbert, A.: Weyl multifractional Ornstein–Uhlenbeck processes mixed with a Gamma distribution. Probab. Math. Stat. 40(2), 269–295 (2020). MR4206415. https://doi.org/10.37190/0208-4147.40.2.5 Gorenflo, R., Kilbas, A.A., Mainardi, F., Rogosin, S.: Mittag-Leffler Functions, Related Topics and Applications. Springer Monographs in Mathematics. Springer, Berlin, Heidelberg (2020). MR4179587. https://doi.org/10.1007/978-3-662-61550-8 Igloi, E., Terdik, G.: Long-range Dependence trough Gamma-mixed Ornstein-Uhlenbeck Process. Electron. J. Probab. 4 (1999). MR1713649. https://doi.org/10.1214/EJP.v4-53 Luchko, Y.: The Wright function and its applications. In: Kochubei, A., Luchko, Y. (eds.) Basic Theory, pp. 241–268. De Gruyter (2019). MR3888404 Mittag-Leffler, G.: Sur la représentation analytique d’une branche uniforme d’une fonction monogène: Cinquième note. Acta Math. 29, 101–181 (1905). MR1555012. https://doi.org/10.1007/BF02403200 Newey, W.K., McFadden, D.: Large Sample Estimation and Hypothesis Testing. In: Engle, R.F., McFadden, D. (eds.) Handbook of Econometics vol. IV, pp. 2111–2245. Elsevier science B.V. (1994) MR1315971 Popov, A.Y., Sedletskii, A.M.: Distribution of roots of Mittag-Leffler functions. J. Math. Sci. 190(2), 209–409 (2013). MR2883249. https://doi.org/10.1007/s10958-013-1255-3 Riesz, M.: L’Intégrale de Riemann-Liouville et le probléme de Cauchy. Acta Math. 81, 222 (1949) MR0030102. https://doi.org/10.1007/BF02395016 Tricomi, F., Erdélyi, A.: The asymptotic expansion of a ratio of gamma functions. Pac. J. Math. 1(1), 133–142 (1951). MR0043948. https://doi.org/10.2140/pjm.1951.1.133 Mathematica, Version 12.2. Wolfram Research, Inc., Champaign, IL (2020)