Statistical inference for generalized Ornstein-Uhlenbeck processes

In this paper, we consider the problem of statistical inference for generalized Ornstein-Uhlenbeck processes of the type \[ X_{t} = e^{-\xi_{t}} \left( X_{0} + \int_{0}^{t} e^{\xi_{u-}} d u \right), \] where \(\xi_s\) is a L{\'e}vy process. Our primal goal is to estimate the characteristics of the L\'evy process \(\xi\) from the low-frequency observations of the process \(X\). We present a novel approach towards estimating the L{\'e}vy triplet of \(\xi,\) which is based on the Mellin transform technique. It is shown that the resulting estimates attain optimal minimax convergence rates. The suggested algorithms are illustrated by numerical simulations.


Introduction
Let (ξ t ) t≥0 be a Lévy process with a Lévy triplet µ, σ 2 , ν . The main object of our study is the so-called generalized Ornstein-Uhlenbeck (GOU) process defined as The GOU processes have recently got much attention in the literature. A comprehensive study of the GOU processes and an extended list of references can be found in the theses of Behme [2], where, in particular, it is shown that X t satisfies the following SDE: The popularity of GOU processes is related to the fact they appear to be useful in several applications. For instance, the process (1) determines the volatility process in the COGARCH (COntinious Generalized AutoRegressive Conditionally Heteroscedastic) model introduced in Klüppelberg et al. [18]. One important result from the theory of GOU processes is that, under some conditions, the process (1) is stationary with invariant stationary distribution given by the distribution of the following exponential functional of ξ : In fact, the properties of the functional (2) have been widely studied in the literature and we refer to the survey by Bertoin and Yor [7] for a theoretical background of the exponential functionals. In particular, it is known that the Mellin transform of the density π of exponential functional, satisfies the following recursive formula where φ(z) is a Laplace exponent of the process ξ, i.e., φ(z) := − log E e −zξ1 , and complex z is taken from the strip Υ := z ∈ R : 0 < Re(z) < θ with θ := sup x ≥ 0 : E[e −xξ1 ] ≤ 1 . (4) The recursive formula (3) first appeared for real z in the paper by Maulik and Zwart [25]. The validity of (3) for complex z was recently shown by Kuznetsov, Pardo and Savov [20]. If ξ t is a subordinator, the parameter θ is equal to infinity. Let us note that the functional A ∞ appeared in such application areas as finance (see, e.g. the monograph by Yor [33]), carousel systems (see Litvak and Adan [23], Litvak and van Zwet [24]), self-similar fragmentations (see Bertoin and Yor [7]), and information transmission problems (especially TCP/IP protocol, see Guillemin, Robert and Zwart [13]). For the detailed discussion on the physical interpretations, we refer to Comtet, Monthus and Yor [10] and the dissertation by Monthus [26].
In this paper, we mainly focus on the case when ξ is a subordinator with finite Lévy measure. In terms of the Lévy triplet (µ, σ 2 , ν), this means that µ > 0, σ = 0, ν(IR − ) = 0 and moreover λ := ν(IR + ) < ∞. Suppose that the process (1) is observed at equidistant time points 0 = t 0 < t 1 < . . . < t n . Since under some mild assumptions the process is stationary and the invariant distribution is given by the distribution of the exponential functional A ∞ (see Fasen, [12]), we assume that X t0 , . . . , X tn are also distributed as A ∞ . Our main goal is statistical inference on the Lévy triplet (µ, σ 2 , ν) based on the observations X t0 , . . . , X tn . More precisely, we will pursue the following two aims: (1) estimation of the drift term µ and the intensity parameter λ; (2) estimation of the Lévy measure ν.
To the best of our knowledge, the statistical inference for GOU processes of the form (1) from their low-frequency observations has not been yet studied in the literature. In fact the resulting statistical problem is quite challenging and needs careful treatment. Indeed, the only connection between the stationary distribution of a GOU process, which can be estimated from the data, and the parameters of the underlying Lévy process is given by the recurrent relation (3) which is rather implicit. The main idea of our procedure for estimating the parameters of the process ξ can be described as follows. First, by making use of (3), we estimate the Laplace exponent φ(z) at the points z = u • + iv ∈ Υ, where u • > 0 is fixed and v varies on the equidistant grid between εV n and V n (with ε > 0 and V n → ∞ as n → ∞) Afterwards, we use the representation (v) → 0 as v → ∞ by the Riemann-Lebesgue lemma, upon taking real and imaginary parts of the left and right hand sides of (5), we are able to consequently estimate the parameters µ and λ. With no doubt, the second aim, a complete recovering of the Lévy measure ν, is the most difficult task. Since the estimates of the parameters µ and λ are already obtained, we can estimate by (5) the Fourier transform The last step of this procedure, the estimation of the Lévy measure ν, is based on the regularised inverse Fourier transform formula. The above estimation algorithm bears some similarity to the spectral estimation algorithm introduced by Belomestny and Reiss [4], [5]. Let us also mention that the problem of statistical inference for Lévy processes (or some their generalizations) observed at low frequency was the subject of many studies, see, e.g.
The paper is organized as follows. In the next section, we formulate our main assumptions and give some examples. In Section 3, the main estimation algorithm is presented and discussed in details. Next, we analyze the convergence rates of the proposed algorithms in Section 4 and provide some numerical examples in Section 5. The proofs of our theoretical results are collected in Section 6.

Main setup
In this article, we study the class of subordinators with finite Lévy measures as possible choice for the Lévy process (ξ t ). In terms of the Lévy triplet (µ, λ, ν), this means that A detailed discussion of the subordination theory as well as various examples of such processes (Gamma, Poisson, tempered stable, inverse Gaussian, Meixner processes, etc.), are given in [1], [6], [11], [29], [30]. Note that in the case of subordinators, the truncation function in the Lévy-Khinchine formula can be omitted, that is, the characteristic exponent of ξ is equal to Later on, we also need the Laplace exponent of ξ, which is defined as Under the assumption (6) the Laplace exponent φ(·) is given by Let us summarise the main properties of the functional A ∞ = ∞ 0 e −ξt dt in this case.

Estimation of the Lévy triplet
In the sequel, we suppose that we are given by the observations X t0 , X t1 , . . . , X tn of the process (1) at the equidistant time points 0 = t 0 < t 1 < . . . < t n with t j = j · ∆ for some ∆ > 0. Assuming that the process X t is stationary (see [2], [12]), we get that the random variables X k := X t k , k = 1, . . . , n, have all the same distribution, which coincides with the distribution of A ∞ . The first step of our estimation procedure consists in the estimation the Laplace exponent φ(z) for z = u • + iv, where u • > 0 is fixed and v varies. An estimator of φ(z) can be obtained from the recursive formula (3) for the Mellin transform of π. In [9], this formula is proved for real positive z such that φ(z) > 0 and M(z + 1) < ∞. The case of complex z is considered in [20], where one can also find some generalizations of the formula (3) to the integrals with respect to Brownian motion with drift. In particular, applying Theorem 2 from [20], we get that (3) holds for any z ∈ Υ. In the situation when (ξ t ) is a subordinator, the set Υ coincides with the positive half-plane (equivalently, the parameter θ is equal to infinity) due to and then define an estimate of the Laplace exponent φ(z) by If the sequence X 1 , . . . , X n has some mixing properties, then we can expect that Y n (z) → φ(z) in probability.

Estimation of λ and µ
Under our assumptions, the Laplace exponent of the Lévy process (ξ t ) can be represented in the form: where λ := I R+ ν(dx) andν(dx) := e −u • x ν(dx). The general idea of the procedure described below is to estimate the Laplace exponent φ(·) at the points z = u • + iv with v ∈ R and then use the relation (13) for the estimation of parameters. By the Riemann-Lebesque lemma, F[ν](−v) → 0 as v → +∞ (see, e.g., [17]) and we conclude from (13) that φ(u • + iv) is approximately (at least for large v) a linear function in v with the slope µ and the intercept term λ. This observation suggests that a properly weighted least-squares approach can be applied to estimate µ and λ. Let V n be a sequence of positive real numbers and w(·) be a nonnegative weight function supported on [0, 1]. Define a scaled weight function w n (v) = V −1 n w(v/V n ) and introduce the estimators of the parameters λ and µ as the solution of the following optimization problem: with Y n (z) defined in (12). The above optimisation problem admits an explicit solution given by Taking into account the definition of the weight function w n (·), we get also some equivalent representations of the estimators µ n and λ n µ n = arg min In practice, we need to replace the above integrals by sums. To this end, let the numbers α 1 , . . . , α M constitute an equidistant grid on the set [ε, 1] for some ε > 0. We estimate the Mellin transform M(z) for all z ∈ {u • + iα m V n , m = 1, . . . , M } and z ∈ {u • − 1 + iα m V n , m = 1, . . . , M } and so get the estimates of the Laplace exponent at the discrete points z = u • + iα m V n (see above). Now we define an estimate of the parameter µ viâ Afterwards, we estimate the parameter λ bŷ The whole algorithm is described below.
Denote v m,n := α m V n . Algorithm:

Estimation of the Lévy measure ν
As a result of Algorithm 1, we obtain the estimates µ n and λ n of the parameters µ and λ, respectively. Based on (13), we first define an estimate for the Fourier transform ofν viâ Next we estimate the measure ν by a regularised Fourier inversion formula where K is a regularizing symmetric kernel supported on [−1, 1]. Note that with a slight abuse of notation, we use ν also for the density of the Lévy measure, and ν n for an estimate of this density. In what follows, we also use the notation ν n = e −u • x ν n . The formal description of the algorithm is given below.
Denote v m,n := α m V n . Algorithm: 1-2 The first two steps coincide with ones of Algorithm 1.
Remark 3.1. It is a worth mentioning that the estimation Algorithms 1 and 2 can be applied to a more general situation when where the process τ t is a difference between two subordinators, i.e., τ t = τ + t +τ − t , and τ + and τ − are the processes of finite variation with Lévy measures ν + and ν − concentrated on IR + and IR − , respectively. In fact, in this case, the formula (13) still holds with Therefore, the consequent estimation of µ, λ and the Fourier transform of the measure e −u • x ν(dx), as well as the estimation of ν are still possible.

Convergence
In order to analyse the convergence properties of the estimates µ n , λ n and ν n we need to further specify the class of Lévy processes (ξ t ).
Definiton 4.1. For s ∈ N ∪ {0} and R > 0, let G(s, R) denote the set of all Lévy triplets (µ, 0, ν), such that ν is supported on R + and max ν(R + ), Note that if (21) holds, thenν is s-times (weakly) differentiable with It turns out that the convergence rates of the estimates µ n , λ n and ν n crucially depend on the asymptotic behaviour of the Mellin transform of A ∞ . In order to specify this behaviour, let us fix some u 0 > 0 and introduce two classes of probability densities: where α, β ∈ IR, L > 0 and for any density p, M[p] stands for the Mellin transform of p. Before we formulate the main convergence results, let us look at some examples.
with N, m j ∈ N, ρ j > 0, g jk > 0. First note that the assumption (6) obviously holds. Let us now check (21). We can apply the well-known Erdélyi lemma to derive for any exponentially decaying and smooth function f on IR + , and some complex c 1 depending on f . Therefore, we conclude that where c 2 > 0 depends on u • . Hence for any s < k * − 1, the condition (21) holds for some R > 0. Furthermore, taking into account the asymptotic behaviour of the Gamma function (see, e.g., formula 8.328 from [14]): we derive where ζ 1 , . . . , ζ K are the roots of the equation see [19]. Therefore, for any u • > 1/2, we conclude that π ∈ E(π/2, L) with any L > 0.
Example 4.3. Next, we provide an example of a Lévy process ξ t with A ∞ = ∞ 0 e −ξt dt having a density from P(β, L). Consider a subordinator T with drift µ > 0 and the Lévy density The exponential functional A ∞ of the process (ξ t ) has a density of the form with some C 1 > 0, see [9]. In other words, A ∞ has the same distribution as ξ/µ, where the r.v. ξ has the Beta distribution with parameters α = b + 1 and β = a/µ = λ/µ. The Mellin transform of the function π(x) in the half-plane Re(s) > −α is hence given by .
Using (25), we conclude that the Mellin transform of A ∞ has a polynomial decay in this case. More precisely, as |v| → ∞ and therefore π ∈ P(λ/µ, L).
Let us now formulate the main result concerning the convergence of the estimates µ n and λ n .
Theorem 4.4 (upper bounds for µ n and λ n ). Let (ξ t ) be a Lévy process with a triplet from G(s, R). Suppose that the sequence X 0 , X 1 , . . . , X n is αmixing and strictly stationary. Denote the α-mixing coefficients of the sequence X 0 , X 1 , . . . , X n by α(s).
Then the quadratic risks of the estimates µ n and λ n , under the choice V n = n 1/(2β+2s+3) , satisfy the following asymptotic relations then the choice leads to the rates E |µ n − µ| 2 log −2(s+2) (n), Proof. Proof is given in Section 6.1.
In a similar way, we can establish the upper bounds for the risk ofν n . In the theorem formulated below, the quality of the estimateν n is measured in terms of the mean integrated squared error (MISE) with some A > 0.
(ii) If the density of A ∞ belongs to the class E(α, L) and then under the choice Proof. Proof is given in Section 6.2.
The next theorem shows that the rates obtained in the previous theorem are optimal up to a logarithmic factor.  (29) where the infimum is taken over all possible estimatesν n of the functionν based on i.i.d. sample X 1 , . . . , X n from the distribution π T of A ∞ := ∞ 0 e −ξt dt such that the Lévy triplet T of (ξ t ) belongs to G(s, R).
Proof. Proof is given in Section 6.3.
An important condition of Theorems 4.4 and 4.5 is (26), which means that the sequence X 0 , X 1 , . . . , X n is exponentially α-mixing. Since β-mixing coefficient between two sigma-algebras is larger than or equal to the corresponding α-mixing coefficient, it is sufficient to show that X 0 , X 1 , . . . , X n is an exponentially β-mixing sequence (see Section 1.1 from [8]). For the case of the GOU processes (1), the latter question was addressed in [12]. The sufficient conditions for exponential β-mixing given in [12] are: 1. the distribution of A ∞ has a Pareto-like asymptotic behaviour, that is, with some α > 0 and C > 0; 2. there exist A > 0, B > A and h > 0 such that ψ(A) = 0, ψ(B) < ∞ with ψ given in (7), and As it is proved in [22], both conditions are guaranteed by the positiveness of µ and the existence of a positive zero of the function ψ(·). We refer also to [21] for some further results in this direction.

Simulation study
Example 1. Consider the subordinator τ t with the Lévy density Note that in this case, λ = I R+ ν(u)du = a. Define a Lévy process where W t is a Brownian motion. The Laplace exponent of ξ t is given by In [9], it is shown that the exponential functional A ∞ = ∞ 0 e −ξt dt is finite for any µ and σ, and moreover the density function π of A ∞ satisfies the following differential equation Some special cases are considered below: 1. In the case µ = 0, σ = 0 (pure jump process), this equation has a solution and therefore A ∞ 1, a), where G(α, β) is a Gamma distribution with shape parameter α and rate β. 2. If µ > 0, σ = 0 (pure jump process with drift), then In this situation A ∞ 3. In the case µ = 0, σ = 0, the equation (33) also allows for the closed form solution. Assuming for simplicity σ 2 /2 = 1, µ = −(b + 1), we get the solution of (33) in the following form: where we denote by I µ the modified Bessel function of the first kind, µ = a + 1/4, and the constant C is later chosen to guarantee the condition ∞ 0 π 3 (x)dx = 1.
For our numerical study, we assume that the data are generated from the distribution of (2), where the process (ξ t ) is defined by (31) with µ = 1.8, σ = 0, and the subordinator τ t in the form (30) with a = 0.7, b = 0.2. A sample from the distribution of the integral A ∞ can be simulated from the corresponding Beta-distribution, see (35). In the first step, we estimate the Mellin transform M(z) for z = u + iv with u = u • = 29 and u = u • + 1 = 30 and v lying on the equidistant grid between −30 and 30. Next, we estimate the Laplace exponent of ξ by the formula (12). Figure 1 graphically compares the proposed estimator of the Laplace exponent φ(u • + iv) with its theoretical values (µ + a/(b + u • + iv)) · (u • + iv) . Estimates for the parameters µ and λ = a are given in (15) and (17), respectively. The boxplots of this estimates based on 25 simulation runs are presented on Figure 2.

Example 2. Consider the compound Poisson process
where q ∈ (0, 1) is fixed, N t is a Poisson process with intensity λ and η k are i.i.d. random variables with a distribution L. The integral A ∞ admits the representation where T n is the jump time of N, i.e., T n = inf {t : N t = n}, and S n = n k=1 η k . Note that if η k take only positive values, then ξ t is a subordinator. For the overview of the properties of the integral A ∞ in the particular case η k ≡ 1 (that is, ξ t is a Poisson process up to a constant), we refer to [7].
Fix some positive α and consider the case when L is the standard normal distribution truncated on the interval (α, +∞). The density function of L is given by where p(·) and F (·) are the density and the distribution functions of the standard Normal distribution. In this case, the Laplace exponent of ξ t is equal to where the function F (·) can be calculated for complex arguments from the error function: In this example, we aim to estimate the Lévy measure of the process (ξ t ), which is given by For our numerical study, we take q = 0.5, α = 0.1, and λ = 1. First, we estimate the Laplace exponent by (12). The quality of the corresponding estimate at the complex points z = u • + iv with u • = 1 and v ∈ [−5, 5] can be visually seen in Figure 3. Next, we proceed with the estimation of the Fourier transform of the measureν(x) := e −u • x ν(x) by applying (18). For the last step of the Algorithm 2, i.e. the reconstruction of the Lévy measure by (19), we follow [3] and use the so-called flat-top kernel, which is defined as follows: The quality of the resulting estimateν n is shown in Figure 4. 6. Proofs

Upper bounds for the quadratic risks of µ n and λ n
The next proposition is the main technical result for this section.
Proposition 6.1. Let ξ t be a Lévy triplet from G(s, R). Suppose that the sequence X 0 , X 1 , . . . , X n of observations of the exponential functional A ∞ := lim T →∞ A T = ∞ 0 e −ξt dt is α-mixing and strictly stationary. Denote the mixing coefficients of the sequence X 0 , X 1 , . . . , X n by α(s).
Then for any p ∈ {0, 1, . . . , n} we have and we have Note that

Upper bounds for M ISE(ν n )
Proposition 6.2. Let the assumptions of the Proposition 6.1 be fulfilled and let the kernel K(·) satisfy the assumption (28). Then the mean integrated squared error of the estimatorν n (x) satisfies the following asymptotic relation Proof. Recall that By the Parsenval's identity, The treatment of J 1 is based on the observation that We get that As it was shown before, E |M n (u • + iv) − M(u • + iv)| 2 n −1 , see (38)-(40). Therefore, To complete the proof, it is sufficient to note that Proof of Theorem 4.5 (i) Recall that if π ∈ P(β, L), then see the proof of Theorem 4.4. Taking into account that J 1 n −1 V 2β+3 n , we arrive at Choosing V n = n 1/(2β+2s+3) , we get

Lower bounds for M ISE
Proof of Theorem 4.6. The general idea of the proof is to apply Theorem 2.7 from [32]. This theorem yields that (29) holds, if there exists a parameterized set of Lévy triplets for some s ∈ N ∪ 0, R > 0, L > 0 and a set of parameters {θ (j) , j = 0, . . . , M } such that the following two properties hold.
Below we present a detailed proof for the polynomial case.

Presentation of the models.
Consider an exponential Lévy model A 0,∞ = ∞ 0 e −ξ0,s ds, where ξ 0,s is a Lévy subordinator with a triplet (1, 0, ν 0 ) and ν 0 (x) = abe −bx for some 0 < a ≤ 1, 0 < b < 1. It is clear that (1, 0, ν 0 ) ∈ G(0, R) for some R > 0 and the Laplace exponent of ξ 0,s is given by see Example 1 from Section 5. For the case of general classes G(s, R) with s > 0, we could take a Lévy density of the form ν 0 (x) = b 1+s x s e −bx /Γ(s + 1). Fix some L > 0 and let us construct now a parameterized set of Lévy triplets T θ = (1, 0, ν θ ), θ ∈ {0, 1} L , with Lévy measure ν θ defined by where δ > 0 small enough, ∆ θ (x) := g θ (x) + a(g θ exp(−b·)(x)) , 2. Distributional properties of the models. In this step, we perform some technical calculations, which will be used later. It holds where L[g θ ](z) is the Laplace transform of the function g 0 (·), which is equal to We see that where φ θ (·) is the Laplace exponent of a Lévy process ξ θ,s with the Lévy triplet T θ . Furthermore, the Laplace transform of g 0 is given by The Mellin transform of the density π θ corresponding to the Lévy model T θ satisfies the following functional equation we derive the following infinite product representation for the ratio M θ (z)/M 0 (z) (1 + δL[g θ ](z + k)).
Furthermore, it can be proved that for some absolute constant c > 0. Note that the random variables A θ,∞ = ∞ 0 e −ξ θ,s ds with ξ θ,s being a Lévy process with the triplet T θ , satisfies 0 < A θ,∞ < 1 a.s. Moreover the density p 0 of the r.v. A 0,∞ has the form and the Mellin transform M 0 (z) of A 0,∞ is given by see Example 4.3.