VMSTA Modern Stochastics: Theory and Applications 2351-60542351-60462351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA267 10.15559/24-VMSTA267 Research Article Statistical inference for nth-order mixed fractional Brownian motion with polynomial drift https://orcid.org/0000-0003-0336-9092 El OmariMohamedelomari.m@ucd.ac.maa ESIM, Polydisciplinary Faculty of Sidi Bennour, Chouaib Doukkali University, Av. of Faculties, El Jadida, 24000, Morocco 2025 1911202412216918727520242110202429102024 © 2025 The Author(s). Published by VTeX2025 Open access article under the CC BY license.

The mixed model with polynomial drift of the form X ( t ) = θ P ( t ) + α W ( t ) + σ B H n ( t ) is studied, where B H n is the nth-order fractional Brownian motion with Hurst index H ∈ ( n − 1 , n ) and n ≥ 2, independent of the Wiener process W. The polynomial function P is known, with degree d ( P ) ∈ [ 1 , n ). Based on discrete observations and using the ergodic theorem estimates of H, α 2 and σ 2 are given. Finally, a continuous time maximum likelihood estimator of θ is provided. Both strong consistency and asymptotic normality of the proposed estimators are established.

 nth-order fractional Brownian motion maximum likelihood estimator ergodicity consistency asymptotic normality 60G18 60G22 62F10 62F12 Introduction

Since the introduction of nth-order fractional Brownian motions (n-fBm’s) in [26], less attention is devoted to them from the statistical point of view. Only a few papers tackled the problem of estimation for models driven by n-fBm’s, like estimating the Hurst parameter H ∈ ( n − 1 , n ), , [26] or estimating the drift with possibly other diffusion parameters [5, 10, 11]. The importance of n-fBm’s lies in themselves as an extension of the fractional Brownian motion. They can capture a wide class of 1 / f δ -nonstationary signals with the range δ ∈ ( 1 , ∞ ), and their smoothness is controlled by the order n. As self-similar processes they arise in image processing for modeling texture having multiscale patterns such as natural scenes [17, 25], bone texture radiographs [21], or rough surfaces [33] and in Nile River data [3, p. 84]. For other new application areas of n-fBm’s, we refer the reader to [14, 15]. For theoretical analysis of the n-fBm, we refer to [9, 10, 16, 30]. Recently, some extensions of the n-fBm are proposed and studied in [9, 10]. In the present work, we consider the continuous time model X ( t ) = θ P ( t ) + α W ( t ) + σ B H n ( t ) , t ≥ 0 , where W is a Wiener process independent of the n-fBm B H n with the known order n ≥ 2 and unknown Hurst parameter H ∈ ( n − 1 , n ). We assume that P is a known polynomial function of degree d ( P ) ∈ [ 1 , n ), thus the results we provide here can be viewed as complementary to those in [11]. Clearly, the approach of [11] used to estimate θ and σ relies on the crucial condition d ( P ) ≥ n which is violated here. Therefore it is not applicable. Instead, we apply the method of moments and ergodic theorem to estimate α 2 , σ 2 and H. Finally, we follow [4, 23] to estimate the drift parameter θ.

The model (1) is one of the class of mixed models that have gained much attention from many researchers. The first initiative was taken by Cheridito [6]. In his work the so-called mixed fBm (linear combination of two independent Wiener processes and fBm) was introduced, and shown to be a good tool to describe the fluctuations of financial assets in the presence of both arbitrage-free and long-range dependence. These processes do not exhibit the self-similarity property, but instead the so-called mixed self-similarity arises. It was shown in [6] that the mixed fBm is equivalent to a Wiener process if and only if H ∈ ( 3 / 4 , 1 ] ∪ { 1 / 2 }. We think this may be restrictive for the range of Hurst index H. Instead, we suggest the use of the nth-order mixed fBm shown to be equivalent to a Wiener process [9, Theorem 2.2] with a wide range of H. The statistical study for mixtures of fBm’s has attracted many researchers (e.g., [1, 7, 8, 18, 32, 34] among others). The procedures of estimation used in these references rely on one of the following methods: least squares method, mixed power-variations, and maximum likelihood approach which may be combined with some numerical optimization methods. In [18, 27, 28], the authors used the ergodic theorem to estimate all the parameters at once by the generalized method of moments. This approach seems very practical, and one can carefully adapt these procedures to solve the problem of estimation for the more general model: X ( t ) = θ t + ∑ k = 1 p α k B H k k , t ≥ 0, where B H k k , k = 1 , 2 , … , p, are independent fBm’s with different Hurst parameters.

The rest of the paper is organized as follows. In Section 2 we present basic properties of the n-fBm. Section 3 is devoted to our main results. Namely, we discuss the asymptotic behavior of the estimators of H, α 2 and σ 2 in Subsection 3.1, while Subsection 3.2 is devoted to examine the estimator of θ. In Section 4 we gather our concluding remarks on the performance of the proposed estimators. Finally, some auxiliary results with their proofs are relegated to Appendix A. On a stochastic basis ( Ω , F , ( F ) t , P ) we are concerned with convergence in law and almost surely, that are denoted as ⟶ L and ⟶ P − a.s , respectively. If ξ, ζ are two random variables, then E θ ( ξ ), V ar θ ( ξ ), and ξ ≜ θ ζ will stand for the mean, the variance, the covariance and equality in law under P θ , but if P is the probability of interest we omit θ. We systematically use the notations a ∧ b : = min ( a , b ), a ∨ b : = max ( a , b ), ( a ) + = max ( a , 0 ) for any , F t X denotes the natural filtration of the process X, while the transpose of a matrix M is denoted by M ′ . When taking increments of a function g ( t ) at t = t 0 , we simplify notations by setting Δ h g ( t 0 ) = Δ h g ( t ) | t = t 0 .

 <italic>n</italic>th-order fractional Brownian motion

The n-fBm (hereafter B H n ( t ), H ∈ ( n − 1 , n ), n ≥ 1 is integer) was first introduced in [26]. It is formally defined as a zero-mean Gaussian process starting at zero with the integral representation B H n ( t ) = 1 Γ ( H + 1 / 2 ) ∫ − ∞ 0 ( t − s ) H − 1 / 2 − ∑ j = 0 n − 1 H − 1 / 2 j ( − s ) H − 1 / 2 − j t j d B ( s ) + 1 Γ ( H + 1 / 2 ) ∫ 0 t ( t − s ) H − 1 / 2 d B ( s ) , for all t ≥ 0 , where B ( t ) is a two-sided standard Brownian motion, Γ ( x ) denotes the gamma function and α j = α ( α − 1 ) ⋯ ( α − ( j − 1 ) ) j ! , α 0 = 1 (by convention) . We retrieve the fBm in the case n = 1, as equation (2) reduces to the integral representation of the fBm given in [22]. The process B H n ( t ) satisfies the following properties for which the proofs can be found in [9, 10, 26].

B H n ( t ) is self-similar with exponent H, i.e. for any c > 0, the processes B H n ( c t ) and c H B H n ( t ) have the same finite dimensional distributions.

B H n ( t ) admits derivatives up to order n − 1 starting at zero and the ( n − 1 )th derivative coincides with the fBm, that is, d n − 1 d t n − 1 ( B H n ( t ) ) = B H − ( n − 1 ) 1 ( t ).

B H n ( t ) exhibits a long-range dependence with stationarity of increments achieved at order n, which means that the increments Δ s ( k ) B H n ( t ), s > 0, k ≤ n − 1, are nonstationary and Δ s ( n ) B H n ( t ) is a stationary process. Here Δ l ( k ) g ( x ) denotes the increments of order k of the function g ( x ) with the explicit form Δ l ( k ) g ( x ) = ∑ j = 0 k ( − 1 ) k − j k j g ( x + j l ) . If g is a multivariate function, then we distinguish the increments with respect to variables as Δ l , x j ( k ) g ( x 1 , … , x m ), j = 1 , … , m. For a definition and detailed discussion on the long-range dependence, we refer the reader to [12, Subsection 2.2].

The covariance function of the process B H n ( t ) is given by G H , n ( t , s ) = ( − 1 ) n C H n 2 | t − s | 2 H − ∑ j = 0 n − 1 ( − 1 ) j 2 H j [ ( t s ) j | s | 2 H + ( s t ) j | t | 2 H ] , where C H n is a nonnegative constant defined recursively by C H 1 = 1 / ( Γ ( 2 H + 1 ) sin ( π H ) ) and for n ≥ 2 C H n = C H − ( n − 1 ) 1 ( 2 H ) ( 2 H − 1 ) ⋯ ( 2 H − ( 2 n − 3 ) ) = 1 Γ ( 2 H + 1 ) | sin ( π H ) | . In particular, V ar ( B H n ( t ) ) = C H n 2 H − 1 n − 1 | t | 2 H .

For any n ≥ 2 the process B H n ( t ) is a semimartingale with finite variation.

B H n ( t ) is only a Markov process in the case n = 1 and H = 1 / 2.

B H n ( t ) can be extended to an α-order fBm U H α ( t ) (see [10]) defined as U H α ( t ) = 1 Γ ( α + 1 ) ∫ 0 t ( t − s ) α d B H ( s ) , H ∈ ( 0 , 1 ) , α ∈ ( − 1 , ∞ ) . whenever this integral exists. Here B H denotes a one-sided fBm. In the case α = 0, we retrieve the fBm B H . If α = n − 1, then U H α ( t ) coincides with the n-fBm with the Hurst parameter H ′ = H + ( n − 1 ).

The main results Estimating the noise parameters <italic>H</italic>, <inline-formula id="j_vmsta267_ineq_084"><alternatives><mml:math> <mml:msup> <mml:mrow> <mml:mi mathvariant="italic">α</mml:mi> </mml:mrow> <mml:mrow> <mml:mn>2</mml:mn> </mml:mrow> </mml:msup></mml:math><tex-math><![CDATA[${\alpha ^{2}}$]]></tex-math></alternatives></inline-formula> and <inline-formula id="j_vmsta267_ineq_085"><alternatives><mml:math> <mml:msup> <mml:mrow> <mml:mi mathvariant="italic">σ</mml:mi> </mml:mrow> <mml:mrow> <mml:mn>2</mml:mn> </mml:mrow> </mml:msup></mml:math><tex-math><![CDATA[${\sigma ^{2}}$]]></tex-math></alternatives></inline-formula>

Fix h , l > 0 and consider the process X defined by (1). After taking the increments of order n we can get a stationary sequence Z k , k = 1 , 2 , …. In fact, we have Z h ( t ) : = Δ h ( n ) X ( t ) = α Δ h ( n − 1 ) Δ h W ( t ) + σ Δ h ( n ) B H n ( t ) = α ∑ j = 0 n − 1 d j ( n − 1 ) Y h , j ( t ) + σ Δ h ( n ) B H n ( t ) , where Y h , j ( t ) = W ( t + h ( j + 1 ) ) − W ( t + h j ) and d j ( n ) = ( − 1 ) n − j n j . Now, we define our sequence { Z k ( l ) } k ≥ 0 as Z k ( l ) = Z l h ( t ) | t = k h , k = 1 , 2 , …. This sequence is stationary, and one can observe that Y l h , j ( t ), j = 0 , 1 , … , ( n − 1 ), with fixed t are independent, and independent of Δ l h ( n ) B H n ( t ). Therefore, E [ ( Z 0 ( l ) ) 2 ] = α 2 ∑ j = 0 n − 1 ( d j ( n − 1 ) ) 2 E | W ( l h ( j + 1 ) ) − W ( l h j ) | 2 + σ 2 E | Δ l h ( n ) B H n ( 0 ) | 2 = l h α 2 ∑ j = 0 n − 1 ( d j ( n − 1 ) ) 2 + ( l h ) 2 H σ 2 ( − 1 ) n C H n 2 ∑ j = − n n ( − 1 ) j 2 n n + j | j | 2 H = : Ψ l , n ( H , α 2 , σ 2 ) . Let us introduce some statistics which will be essential for establishing the asymptotic behavior of our estimators. For fixed h > 0 and l = 1 , 2 , 4 we define J l , N = 1 N ∑ k = 0 N − 1 ( ∑ j = 0 n d j ( n ) X ( h ( k + l j ) ) ) 2 .

For a fixed h , l > 0 we have J l , N ⟶ P -a.s. Ψ l , n ( H , α 2 , σ 2 ) , as N → ∞ .

The statistic J l , N defined in (7) can be rewritten as J l , N = 1 N ∑ k = 0 N − 1 ( Z k ( l ) ) 2 . If { Z k ( l ) } k ≥ 0 is ergodic, then Lemma 1 follows immediately. Since { Z k ( l ) } k ≥ 0 is a centered stationary Gaussian sequence, it suffices to show that its autocovariance function vanishes as the time lag tends to infinity. Indeed, we have E ( Z 0 ( l ) Z k ( l ) ) = E ( α ∑ j = 0 n − 1 d j ( n − 1 ) Y l h , j ( 0 ) + σ Δ l h ( n ) B H n ( 0 ) ) × ( α ∑ j = 0 n − 1 d j ( n − 1 ) Y l h , j ( k h ) + σ Δ l h ( n ) B H n ( k h ) ) = α 2 ∑ r 1 , r 2 = 0 n − 1 d r 1 ( n − 1 ) d r 2 ( n − 1 ) E ( Y l h , r 1 ( 0 ) Y l h , r 2 ( k h ) ) + σ 2 E ( Δ l h ( n ) B H n ( 0 ) Δ l h ( n ) B H n ( k h ) ) = h α 2 ∑ r 1 , r 2 = 0 n − 1 d r 1 ( n − 1 ) d r 2 ( n − 1 ) { ( l ∧ ( k + l ( r 2 − r 1 ) + l ) ) + − ( l ∧ ( k + l ( r 2 − r 1 ) ) ) + } + σ 2 G H , n ( n ) ( k h ) ∼ σ 2 D H l 2 n h 2 H k 2 H − 2 n ⟶ 0 , as k ↑ ∞ ( because H ∈ ( n − 1 , n ) ) , where D H is some nonnegative constant independent of k. Here G H , n ( n ) ( τ ) denotes the autocovariance function of Δ l h ( n ) B H n with time lag τ. Formula (8) is justified by equations (26)–(29) in [26].  □

The statistic Θ ˆ = ( H ˆ , α 2 ˆ , σ 2 ˆ ) defined by H ˆ = 1 2 log 2 + ( J 4 , N − 2 J 2 , N J 2 , N − 2 J 1 , N ) , α 2 ˆ = J 2 , N − 2 2 H ˆ J 1 , N h ( 2 − 2 2 H ˆ ) ∑ j = 0 n − 1 ( d j ( n − 1 ) ) 2 , σ 2 ˆ = 2 ( − 1 ) n Γ ( 2 H ˆ + 1 ) | sin ( π H ˆ ) | ( J 2 , N − 2 J 1 , N ) h 2 H ˆ ( 2 2 H ˆ − 2 ) ∑ j = − n n ( − 1 ) j ( 2 n n + j ) | j | 2 H ˆ is a strongly consistent estimator of Θ = ( H , α 2 , σ 2 ), as N → ∞, where log 2 + ( x ) = log 2 ( x ) , if x > 0 , log 2 + ( x ) = 0 , if x ≤ 0 .

The strong consistency follows readily by Lemma 1 combined with the continuous mapping theorem. Let us now justify the construction of our estimators (9)–(11). According to Lemma 1 we have the following almost sure convergences, as N → ∞, J 1 , N ⟶ Ψ 1 , n ( H , α 2 , σ 2 ) = h α 2 V + h 2 H σ 2 U H , J 2 , N ⟶ Ψ 2 , n ( H , α 2 , σ 2 ) = 2 h α 2 V + 2 2 H h 2 H σ 2 U H , J 4 , N ⟶ Ψ 4 , n ( H , α 2 , σ 2 ) = 4 h α 2 V + 4 2 H h 2 H σ 2 U H , where V = ∑ j = 0 n − 1 ( d j ( n − 1 ) ) 2 and U H = ( − 1 ) n C H n 2 ∑ j = − n n ( − 1 ) j ( 2 n n + j ) | j | 2 H . It is not hard to see that J 4 , N − 2 J 2 , N ⟶ 2 2 H ( 2 2 H − 2 ) h 2 H σ 2 U H , J 2 , N − 2 J 1 , N ⟶ ( 2 2 H − 2 ) h 2 H σ 2 U H , which in turn implies ( J 4 , N − 2 J 2 , N ) / ( J 2 , N − 2 J 1 , N ) ⟶ 2 2 H , as N → ∞. It is then natural to define H ˆ by (9). Substituting H by H ˆ in (12) and using similar procedures one can justify (10)–(11).  □

The estimators given in (10)–(11) depend on H ˆ. As a result, our estimators should be computed in the following order: H ˆ, α 2 ˆ and σ 2 ˆ.

The following notations are needed to establish the asymptotic normality of the estimator Θ ˆ = ( H ˆ , α 2 ˆ , σ 2 ˆ ) ′  : Θ = ( H , α 2 , σ 2 ) ′ , Ψ ( · ) = ( Ψ 1 , n ( · ) , Ψ 2 , n ( · ) , Ψ 4 , n ( · ) ) ′ , Ψ l , n ( H , α 2 , σ 2 ) = l α 2 V + l 2 H σ 2 U H , V = h ∑ j = 0 n − 1 ( d j ( n − 1 ) ) 2 , U H = ( − 1 ) n h 2 H f H 2 Γ ( 2 H + 1 ) | sin ( π H ) | , f H = ∑ j = − n n ( − 1 ) j 2 n n + j | j | 2 H .

Consider the vector J N = ( J 1 , N , J 2 , N , J 4 , N ) ′ , where J l , N is defined by (7). Then for any H ∈ ( n − 1 , n − 1 / 4 ) we have N ( J N − Ψ ( Θ ) ) ⟶ L N ( 0 , Σ ) , as N → ∞ , where

First, we recall that J l , N = 1 N ∑ k = 0 N − 1 ( Z k ( l ) ) 2 , Z k ( l ) = Δ l h ( n ) X ( k h ). Let and set F N ( λ ) = ∑ j = 1 3 λ j { N ( J 2 j − 1 , N − Ψ 2 j − 1 , n ( Θ ) ) }. Straightforward computations lead to F N ( λ ) = ∑ j = 1 3 λ j { N ( 1 N ∑ k = 0 N − 1 ( Z k ( 2 j − 1 ) ) 2 − Ψ 2 j − 1 , n ( Θ ) ) } = 1 N ∑ k = 0 N − 1 ∑ j = 1 3 λ j { ( Z k ( 2 j − 1 ) ) 2 − E ( Z k ( 2 j − 1 ) ) 2 } = : 1 N ∑ k = 0 N − 1 { f ( Z k ) − E f ( Z k ) } , Z k = ( Z k ( 1 ) , Z k ( 2 ) , Z k ( 4 ) ) , where f ( x , y , z ) = λ 1 x 2 + λ 2 y 2 + λ 3 z 2 . According to the Cramér–Wold device, the vector N ( J N − Ψ ( Θ ) ) is asymptotically normal if and only if F N ( λ ) is asymptotically normal for each . Since f is of Hermite rank 2 (see Lemma 3 in the Appendix), and because of stationarity of the sequence { Z k } k = 0 ∞ , we need only to examine the asymptotics of the covariance functions R i , j ( k ) : = E ( Z m ( 2 i − 1 ) Z m + k ( 2 j − 1 ) ) = E ( Z 0 ( 2 i − 1 ) Z k ( 2 j − 1 ) ) , 1 ≤ i , j ≤ 3 . In fact, we shall prove that ∑ k = − ∞ ∞ R i , j ( k ) 2 < ∞ for each i, j. Set l = 2 i − 1 , p = 2 j − 1 . We have R i , j ( k ) = E ( Δ l h ( n ) X ( 0 ) Δ p h ( n ) X ( k h ) ) = α 2 ∑ r 1 r 2 = 0 n − 1 d r 1 ( n − 1 ) d r 2 ( n − 1 ) E ( Y l h , r 1 ( n ) ( 0 ) Y p h , r 2 ( n ) ( k h ) ) + σ 2 E ( Δ l h ( n ) B H n ( 0 ) Δ p h ( n ) B H n ( k h ) ) = h α 2 ∑ r 1 r 2 = 0 n − 1 d r 1 ( n − 1 ) d r 2 ( n − 1 ) { ( l ∧ ( k + p r 2 − l r 1 + p ) ) + − ( l ∧ ( k + p r 2 − l r 1 ) ) + } + σ 2 E ( Δ l h ( n ) B H n ( 0 ) Δ p h ( n ) B H n ( k h ) ) ∼ C H ( n ) ( l , h ) h 2 H − 2 n k 2 H − 2 n , as k → ∞ . In (18)–(19) we used the stationarity of increments of W and the fact that Y l h , r ( n ) are independent of B H n , while (20) is justified by Lemma 4 in the Appendix. It follows that ∑ k = 1 ∞ R i , j ( k ) 2 ∼ C H ( n ) ( l , h ) 2 h 4 H − 4 n ∑ k = 1 ∞ k 4 H − 4 n < ∞ , provided that H < n − 1 / 4. Similarly, we have ∑ k = − ∞ − 1 R i , j ( k ) 2 < ∞. Hence, by [2, Theorem 4] one has F N ( λ ) ⟶ L N ( 0 , σ 2 ( λ ) ), as N → ∞, with

where Σ is defined by (17). In the fourth equality above we used the fact that for any centered Gaussian vector ( ξ , ζ ), we have . It is not hard to see that E e i λ ′ { N ( J N − Ψ ( Θ ) ) } = E e i F N ( λ ) ⟶ e − 1 2 λ ′ Σ λ , i = − 1 , which completes the proof of (16).  □

Now, we are ready to apply the Delta method with the use of (16) to get the joint asymptotic normality of the estimators H ˆ, α 2 ˆ and σ 2 ˆ.

Note that in both Theorem 2 and Theorem 3 we restrict ourselves to the range H ∈ ( n − 1 , n − 1 / 4 ) so that the series (21) converges. This is a crucial condition needed in establishing the asymptotic normality.

Consider Θ ˆ, the estimator of Θ = ( H , α 2 , σ 2 ) ′ defined by (9)–(11), and let H ∈ ( n − 1 , n − 1 / 4 ). Then we have N ( Θ ˆ − Θ ) ⟶ L N ( 0 , ( D Ψ − 1 ) Σ ( D Ψ − 1 ) ′ ) where D Ψ − 1 denotes the inverse of the Jacobian matrix D Ψ = ∂ H Ψ 1 , n ( Θ ) ∂ α 2 Ψ 1 , n ( Θ ) ∂ σ 2 Ψ 1 , n ( Θ ) ∂ H Ψ 2 , n ( Θ ) ∂ α 2 Ψ 2 , n ( Θ ) ∂ σ 2 Ψ 2 , n ( Θ ) ∂ H Ψ 4 , n ( Θ ) ∂ α 2 Ψ 4 , n ( Θ ) ∂ σ 2 Ψ 4 , n ( Θ ) = σ 2 U H ′ V U H σ 2 2 2 H ( 2 log ( 2 ) U H + U H ′ ) 2 V 2 2 H U H σ 2 4 2 H ( 2 log ( 4 ) U H + U H ′ ) 4 V 4 2 H U H , where U H ′ = ∂ H U H with V, U H and Ψ l , n ( Θ ), l > 0, being defined by (13)–(15).

First, observe that Θ ˆ is constructed by solving the equation J N = Ψ ( Θ ˆ ), that is, Θ ˆ = Ψ − 1 ( J N ). To get the asymptotic normality of Θ ˆ, it suffices to apply the Delta method to J N with the function Ψ − 1 . If D Ψ − 1 , the Jacobian matrix of Ψ − 1 , exists then by [31, Theorem 3.1] we have N ( Θ ˆ − Θ ) ⟶ L D Ψ − 1 N ( 0 , Σ ) = N ( 0 , ( D Ψ − 1 ) Σ ( D Ψ − 1 ) ′ ). Note that D Ψ − 1 exists if D Ψ is invertible. Evaluating the determinant of D Ψ gives | D Ψ | = σ 2 U H ′ V U H σ 2 2 2 H ( 2 log ( 2 ) U H + U H ′ ) 2 V 2 2 H U H σ 2 4 2 H ( 2 log ( 4 ) U H + U H ′ ) 4 V 4 2 H U H = V σ 2 U H U H ′ 1 1 2 2 H ( 2 log ( 2 ) U H + U H ′ ) 2 2 2 H 4 2 H ( 4 log ( 2 ) U H + U H ′ ) 4 4 2 H = V σ 2 U H U H ′ 1 1 log ( 2 ) 2 2 H + 1 U H + ( 2 2 H − 2 ) U H ′ 0 ( 2 2 H − 2 ) log ( 2 ) 2 4 H + 2 U H + ( 2 4 H − 4 ) U H ′ 0 ( 2 4 H − 4 ) = − V σ 2 U H ( 2 2 H − 2 ) log ( 2 ) 2 2 H + 1 U H + ( 2 2 H − 2 ) U H ′ 1 log ( 2 ) 2 4 H + 2 U H + ( 2 4 H − 4 ) U H ′ ( 2 2 H + 2 ) = − V σ 2 U H ( 2 2 H − 2 ) { log ( 2 ) 2 2 H + 1 ( 2 2 H + 2 ) U H − log ( 2 ) 2 4 H + 2 U H } = V σ 2 U H 2 log ( 2 ) 2 2 H + 1 ( 2 2 H − 2 ) 2 > 0 . Hence, D Ψ is invertible and the proof is complete.  □

 Estimating the drift parameter <italic>θ</italic>

The approach used in Subsection 3.1 to estimate H, α 2 and σ 2 is suitable for the model (1) with no restrictions on d ( P ), the degree of P. Albeit, the drift parameter θ should be estimated differently, because taking higher increments cancels out the drift term. In the case where d ( P ) ≥ n, we may readily estimate θ by θ ˜ t : = X ( t ) P ( t ) , which satisfies the usual properties. But if d ( P ) < n this estimator becomes inconsistent under P. Assume that the parameters H, α 2 and σ 2 are known (or previously estimated) with α > 0, and let P 0 be the probability measure under which W is a Wiener process independent of B H n . Then the following results are obtained.

The maximum likelihood estimator (MLE) of θ is given by θ ˆ T = α M T ⟨ M ⟩ T , where M t = E 0 ( ∫ 0 t P ′ ( s ) d W ( s ) ∣ F t X ) is a Gaussian martingale with quadratic variation ⟨ M ⟩ t = E 0 ( M t 2 ).

We proceed as in [4, 23] and construct the right likelihood function. In the first stage we rewrite (1) as X ( t ) = α W ‾ ( t ) + σ B H n ( t ) with W ( t ) = W ‾ ( t ) − θ α P ( t ), where W ‾ is a Wiener process under P θ . Here, P θ is the probability measure under which W ‾ and B H n are independent. Clearly, P θ ≪ P 0 , where P 0 corresponds to the case when θ = 0. The corresponding Radon–Nikodym derivative is given by d P θ d P 0 = exp ( θ α ∫ 0 T P ′ ( s ) d W ‾ ( s ) − θ 2 2 α 2 ∫ 0 T P ′ ( s ) 2 d s ) and is justified by the Girsanov theorem. But it is not the likelihood of (1) unless F t W ‾ = F t X , otherwise, we shall consider its conditional expectation L T ( X , θ ) = E 0 ( d P θ d P 0 ∣ F T X ) = E 0 { exp ( θ α ∫ 0 T P ′ ( s ) d W ‾ ( s ) − θ 2 2 α 2 ∫ 0 T P ′ ( s ) 2 d s ) ∣ F T X } = E 0 { exp ( θ α ∫ 0 T P ′ ( s ) d W ( s ) − θ 2 2 α 2 ∫ 0 T P ′ ( s ) 2 d s ) ∣ F T X } = exp ( − θ 2 2 α 2 ⟨ ζ ⟩ T ) E 0 { exp ( θ α ζ T ) ∣ F T X } = exp ( θ α M T − θ 2 2 α 2 ( ⟨ ζ ⟩ T − V T ) ) , where ζ t = ∫ 0 t P ′ ( s ) d W ( s ) with its quadratic variation ⟨ ζ ⟩ t , M T = E 0 ( ζ T ∣ F T X ) and V T = E 0 ( ( M T − ζ T ) 2 ∣ F T X ). In (23) we used the fact that W ‾ coincides with W under P 0 , while (24) is due to the Gaussian nature of ζ T given F t X with mean M T and variance V T . It is worth to mention that L T ( X , θ ) is the Radon–Nikodym derivative of μ θ with respect to μ 0 , where μ θ is the probability induced by X on the space of continuous functions with the usual supremum topology under probability P θ .

Since M is an F t X -measurable Gaussian martingale, then its quadratic variation ⟨ M ⟩ is deterministic and one has ⟨ M ⟩ t = E 0 ( M t 2 ) = E 0 ( M t E 0 ( ζ t ∣ F t X ) ) = E 0 ( M t ζ t ) = 1 2 E 0 { M t 2 + ζ t 2 − ( M t − ζ t ) 2 } = E 0 ( ζ t 2 ) − E 0 ( M t − ζ t ) 2 = ⟨ ζ ⟩ t − E 0 ( M t − ζ t ) 2 . As a result, we have the martingale property: E 0 { M t 2 − [ ⟨ ζ ⟩ t − E 0 ( M t − ζ t ) 2 ] ∣ F s X } = M s 2 − [ ⟨ ζ ⟩ s − E 0 ( M s − ζ s ) 2 ] , ∀ t > s . On the other hand, for all t > s E 0 { M t 2 − [ ⟨ ζ ⟩ t − V t ] ∣ F s X } = E 0 { E 0 ( ζ t 2 − ⟨ ζ ⟩ t ∣ F t X ) ∣ F s X } = E 0 { E 0 ( ζ t 2 − ⟨ ζ ⟩ t ∣ F s ) ∣ F s X } , ( F s X ⊂ F t X and F s X ⊂ F s ) = E 0 ( ζ s 2 − ⟨ ζ ⟩ s ∣ F s X ) = E 0 ( ζ s 2 ∣ F s X ) − ⟨ ζ ⟩ s = M s 2 − [ ⟨ ζ ⟩ s + M s 2 − E 0 ( ζ s 2 ∣ F s X ) ] = M s 2 − [ ⟨ ζ ⟩ s − V s ] . We used V t = E 0 ( ζ t 2 ∣ F t X ) − M t 2 in the first equality above, while the third equality follows by the martingale property of ( ζ t , F t ). Combining (26)–(27) with (25) we can assert that ⟨ M ⟩ t = ⟨ ζ ⟩ t − E 0 ( M t − ζ t ) 2 = ⟨ ζ ⟩ t − V t . Hence, (24) becomes L T ( X , θ ) = exp ( θ α M T − θ 2 2 α 2 ⟨ M ⟩ T ) .

This expression of the likelihood at hand allows us to define the MLE of θ as θ ˆ T = α M T ⟨ M ⟩ T . Hence, the proof is complete.  □

The estimator θ ˆ T defined by (22) can be rewritten as θ ˆ T = ∫ 0 T K ( T , s ) d X ( s ) ∫ 0 T K ( T , s ) P ′ ( s ) d s , where the kernel K ( · , · ) is defined by K ( T , s ) = 1 T g T ( s T ), for all s ∈ [ 0 , T ] with g T ( s ) = T α P ′ ( T s ) + ∑ j ≥ 0 a j ψ j ( s ) , for all s ∈ [ 0 , T ] , and a j = − σ 2 T 2 H γ j α ( α 2 + σ 2 T 2 H − 1 γ j ) ∫ 0 1 P ′ ( T u ) ψ j ( u ) d u . Here ( γ j , ψ j ) j = 0 ∞ denotes the system of eigenvalues and eigenvectors of the integral operator with the kernel G H − 1 , n − 1 defined by formula (3).

Under P θ , the model (1) can be rewritten as X ( t ) = σ ∫ 0 t B H − 1 n − 1 ( s ) d s + ∫ 0 t α d W ‾ ( s ), for all t ∈ [ 0 , T ]. Set ζ ‾ T = ∫ 0 T P ′ ( s ) d W ‾ ( s ). It is not hard to see that ( ζ ‾ T , X ( t ) ) forms a Gaussian system, thus we can assert that (e.g., [20, Lemma 10.1]) for each t ∈ [ 0 , T ] there is a function s ↦ K ( t , s ) defined on [ 0 , t ] such that E θ ( ζ ‾ T ∣ F t X ) = ∫ 0 t K ( t , s ) d X ( s ) , P θ -a.s. In particular, M T = E 0 ( ζ T ∣ F T X ) = ∫ 0 T K ( T , s ) d X ( s ) , P 0 -a.s. Note that we used the fact that W ‾ = W under P 0 . We know also that ⟨ M ⟩ T = E 0 ( M T 2 ) = E 0 ( M T ζ T ) = E 0 { ∫ 0 T P ′ ( s ) d W ( s ) ( α ∫ 0 T K ( T , s ) d W ( s ) + σ ∫ 0 T K ( T , s ) d B H n ( s ) ) } = α ∫ 0 T K ( T , s ) P ′ ( s ) d s . It follows that θ ˆ T = α M T ⟨ M ⟩ T = ∫ 0 T K ( T , s ) d X ( s ) ∫ 0 T K ( T , s ) P ′ ( s ) d s .

To complete the proof we shall justify equations (29)–(30). First, observe that under P 0 we have d X ( t ) = α d W ‾ ( t ) + σ d B H n ( t ) = α d W ( t ) + σ B H − 1 n − 1 ( t ) d t . Second, by definition M T must satisfy the equation E 0 ( X ( s ) ζ T ) = E 0 ( X ( s ) M T ), for all s ∈ [ 0 , T ], which becomes, after some computations, α ∫ 0 s P ′ ( u ) d u = α 2 ∫ 0 s K ( T , u ) d u + σ 2 ∫ 0 T ∫ 0 s K ( T , u ) G H − 1 , n − 1 ( u , v ) d v d u . Differentiating this equation with respect to the variable s yields α P ′ ( s ) = α 2 K ( T , s ) + σ 2 ∫ 0 T K ( T , u ) G H − 1 , n − 1 ( u , s ) d u for all s ∈ [ 0 , T ]. Set g T ( s ) = T K ( T , T s ). After the change of variables v = u / T and substituting s by T s we get φ T ( s ) = g T ( s ) + λ 2 ∫ 0 1 g T ( u ) G ‾ T ( s , u ) d u ( λ = σ / α ) , G ‾ T ( s , u ) = T 2 H − 1 G H − 1 , n − 1 ( u , s ) , φ T ( s ) = T α P ′ ( T s ) , for all s ∈ [ 0 , 1 ]. We will construct a solution g T ( · ) to this equation. As the kernel G H − 1 , n − 1 is continuous symmetric and positive on [ 0 , 1 ] × [ 0 , 1 ], it follows by Mercer’s theorem (e.g., [13, p. 136]) that for all s , u ∈ [ 0 , 1 ] G ‾ T ( s , u ) = T 2 H − 1 ∑ j ≥ 0 γ j ψ j ( s ) ψ j ( u ) , with ( γ j , ψ j ) j = 0 ∞ being the system of eigenvalues and eigenvectors of the integral operator defined on L 2 ( [ 0 , 1 ] ) by G ( f ) ( t ) : = ∫ 0 1 G H − 1 , n − 1 ( t , s ) f ( s ) d s , where ( ψ j ) j = 0 ∞ forms an orthonormal basis of L 2 ( [ 0 , 1 ] ) and { γ j } j = 0 ∞ ↓ 0. A solution to (33) of the form y ( s ) = φ T ( s ) + ∑ j ≥ 0 a j ψ j ( s ) exists if and only if ∑ j ≥ 0 ψ j ( s ) [ a j + λ 2 T 2 H − 1 γ j ( ∫ 0 1 φ T ( u ) ψ j ( u ) d u + ∑ l ≥ 0 a l ∫ 0 1 ψ l ( u ) ψ j ( u ) d u ) ] = 0 . Identifying the coefficients we get, for each j ≥ 0, a j = − λ 2 T 2 H − 1 γ j 1 + λ 2 T 2 H − 1 γ j ∫ 0 1 φ T ( u ) ψ j ( u ) d u , and the proof is then complete.  □

The estimator θ ˆ T defined by (28)–(30) is unbiased with variance vanishing at infinity, that is, E θ ( θ ˆ T ) = θ and V a r θ ( θ ˆ T ) ⟶ 0, as T → ∞. Furthermore, ∫ 0 T K ( T , s ) P ′ ( s ) d s ( θ ˆ T − θ ) ≜ N ( 0 , α ) .

Taking expectations in both sides of (32) and observing that ⟨ M ⟩ T is deterministic we obtain E θ ( θ ˆ T ) = E θ ( α M T ⟨ M ⟩ T ) = α ⟨ M ⟩ T E 0 ( M T L T ( X , θ ) ) = α ⟨ M ⟩ T E 0 ( M T e θ α M T − θ 2 2 α 2 ⟨ M ⟩ T ) = α ⟨ M ⟩ T e − θ 2 2 α 2 ⟨ M ⟩ T E 0 ( α d d θ ( e θ α M T ) ) = α 2 ⟨ M ⟩ T e − θ 2 2 α 2 ⟨ M ⟩ T d d θ { e θ 2 2 α 2 ⟨ M ⟩ T E 0 ( L T ( X , θ ) ) } = θ . Similarly, we have E θ ( θ ˆ T 2 ) = α 2 ⟨ M ⟩ T 2 E 0 ( M T 2 L T ( X , θ ) ) = α 2 ⟨ M ⟩ T 2 E 0 ( M T 2 e θ α M T − θ 2 2 α 2 ⟨ M ⟩ T ) = α 2 ⟨ M ⟩ T 2 e − θ 2 2 α 2 ⟨ M ⟩ T E 0 ( α 2 d 2 d θ 2 ( e θ α M T ) ) = α 4 ⟨ M ⟩ T 2 e − θ 2 2 α 2 ⟨ M ⟩ T d 2 d θ 2 { e θ 2 2 α 2 ⟨ M ⟩ T E 0 ( L T ( X , θ ) ) } = α 2 ⟨ M ⟩ T + θ 2 . Hence, V ar ( θ ˆ T ) = α 2 / ⟨ M ⟩ T and by Gaussianity of θ ˆ T the last statement follows. There remains to show that ⟨ M ⟩ T ⟶ ∞ as T → ∞. Let d be the degree of P. By virtue of the self-similarity of W and the change of variables u = s / T one has M T = E 0 ( ζ T ∣ F T X ) = E 0 ( ∫ 0 T P ′ ( s ) d W ( s ) ∣ F T X ) = T 1 / 2 E 0 ( ∫ 0 1 P ′ ( T u ) d W ( u ) ∣ F T X ) = T 1 / 2 ∑ j = 1 d ( j a j ) T j − 1 E 0 ( ∫ 0 1 u j − 1 d W ( u ) ∣ F T X ) = : T 1 / 2 ∑ j = 1 d ( j a j ) T j − 1 A j , T . Let { T k } k = 1 ∞ ↑ ∞ be a sequence of nonnegative numbers and set ξ k = T k 1 / 2 − d M T k . Clearly, F ∞ X = σ ( ⋃ t ≥ 0 F t X ) = σ ( ⋃ k = 1 ∞ F T k X ) and F T 1 X ⊆ F T 2 X ⊆ ⋯ F. On the other hand, E 0 | ∫ 0 1 u j − 1 d W ( u ) | ≤ ∫ 0 1 u 2 j − 2 d u ≤ 1, for all j = 1 , 2 , … , d. Therefore, for each j we have almost surely (e.g., [29, p. 510]) A j , T k ⟶ A j , ∞ as k ↑ ∞. So we deduce that ξ k 2 ∼ ( d a d A d , ∞ ) 2 , as k ↑ ∞. One can readily verify that sup k ≥ 1 E 0 ( ξ k 4 ) < ∞, which guarantees the uniform integrability of { ξ k 2 } k . Thus (e.g., [29, Theorem 4]) T k 1 − 2 d ⟨ M ⟩ T k = E 0 ( ξ k 2 ) ⟶ E 0 ( d a d A d , ∞ ) 2 . Finally, by the arbitrariness of { T k } k we conclude that ⟨ M ⟩ T ∼ E 0 ( d a d A d , ∞ ) 2 T 2 d − 1 , as T → ∞ with d ∈ [ 1 , n ).  □

For any γ > 1 / ( 2 d ( P ) − 1 ) the estimator θ ˆ N γ is strongly consistent, where d ( P ) is the degree of the polynomial function P.

While proving Theorem 4 we have shown that E θ ( θ ˆ T − θ ) 2 = α 2 / ⟨ M ⟩ T and ⟨ M ⟩ T ∼ E 0 ( d a d A d , ∞ ) 2 T 2 d − 1 , as T → ∞, where A d , ∞ = E 0 ( ∫ 0 1 u d − 1 d W ( u ) ∣ F ∞ X ) and a d is the leading coefficient of the polynomial function P of degree d = d ( P ) ∈ [ 1 , n ). Let ϵ > 0 and set T = N γ . We have ∑ N ≥ 1 P θ ( | θ ˆ N γ − θ | > ϵ ) ≤ ϵ − 2 ∑ N ≥ 1 E θ ( θ ˆ N γ − θ ) 2 ≤ ∑ N ≥ 1 α 2 ϵ 2 ⟨ M ⟩ N γ ∼ C ∑ N ≥ 1 N γ ( 1 − 2 d ) < ∞ , provided that γ > 1 / ( 2 d − 1 ) and C is some nonnegative constant. By using the Borel–Cantelli lemma we can assert that Corollary 1 is well established.  □

 Concluding remarks

In the present work an nth-order mixed fractional Brownian motion with polynomial drift was studied. Based on discrete observations, we constructed estimators ( H ˆ , α 2 ˆ , σ 2 ˆ ) of the noise parameters ( H , α 2 , σ 2 ) and examined their asymptotic behavior (consistency and asymptotic normality) by taking higher order increments and using the ergodic theorem. We stress that taking higher order increments cancels out the drift. Thus, we estimated the drift parameter θ independently by using the maximum likelihood approach based on continuous obervations { X ( t ) : t ∈ [ 0 , T ] }. Based on limiting distributions results with the use of Prohorov’s theorem, the estimators of α 2 and θ (with T = N being the time horizon or number of observations) show the same performance as those given in [11]. They have the rate of convergence of order O P ( N − 1 / 2 ) and O P ( N 1 / 2 − d ( P ) ), respectively. Unlike [11] in which the estimator of σ 2 has the rate of convergence of order O P ( N 1 / 2 − H ), the rate of convergence established here equals O P ( N − 1 / 2 ) and is of great importance when H ∈ ( 0 , 1 ). The limiting distributions provided in [11] are nonnormal (except for the estimator of α 2 ), while those given here are all normal. This would simplify the task of finding the associated confidence intervals. An extension of this work would be to study a model driven by an infinite mixture of higher order fBm’s (see [9]).

 Appendix

In this section we provide some auxiliary results needed to establish Theorem 2 and Theorem 3. Lemma 2 is elementary, while Lemma 3 can be found in [24] or [19]. Let us recall the definition of the operator of increments Δ · , · . Given a multivariate function we define its increments of order k with step h (with respect to x j ) as Δ h , x j ( k ) g ( x ) = ∑ r = 0 k ( − 1 ) k − r k r g ( x 1 , … , x j − 1 , x j + r h , x j + 1 , … , x m ) , for all x = ( x 1 , … , x m ) and j = 1 , 2 , … , m. By convention Δ h , x j ( 0 ) = id (identity application) and Δ h , x j ( 1 ) = Δ h , x j .

The following statements hold true.

The operator Δ · , · is linear and commutes with the expectation.

For a given function we have Δ l , t Δ h , s g ( t , s ) = Δ h , s Δ l , t g ( t , s ) = Δ − l , s Δ h , s g ( t , s ) , if g ( t , s ) = G ( t − s ) .

For any and with p + q > j we have Δ l , t ( p ) Δ h , t ( q ) ( t j ) = 0

Let be a function, and let ξ = ( ξ ( 1 ) , … , ξ ( d ) ) be a Gaussian vector such that f ( ξ ) has finite second moment. We define the Hermite rank of f with respect to ξ as rank ( f ) : = inf { τ : ∃ a polynomial P of degree τ with E [ ( f ( ξ ) − E f ( ξ ) ) P ( ξ ) ] ≠ 0 } .

Let q > − 1 / 2, and let ξ = ( ξ ( 1 ) , … , ξ ( d ) ) be a centered Gaussian vector. Then the function defined by: f ( ξ ) = ∑ j = 1 d λ j | ξ ( j ) | q is of Hermite rank 2.

Let B H n be the n-fBm defined by (2). Then for s ≥ 0 fixed and l , h > 0 we have

where C H ( n ) ( l , h ) is a constant independent of t, s given as C H ( n ) ( l , h ) = ( − 1 ) n ( 2 H 2 n ) ∑ r , k = 0 n ( − 1 ) r + k ( n r ) ( n k ) ( k h − r l ) 2 n 2 Γ ( 2 H + 1 ) | sin ( π H ) | .

Let t > s ≥ 0 and set κ n = ( − 1 ) n C H n 2 , c H , j = ( − 1 ) j ( 2 H j ). By using Lemma 2 we obtain

Finally, we substitute κ n and Δ − l , s ( n ) Δ h , s ( n ) ( s 2 n ) by their values with the use of (5) to complete the proof.  □

 Acknowledgement

The author would like to thank the Editor and two anonymous reviewers for their careful reading and constructive comments which have led to improvements of the article.

 Declarations

Conflicts of interests.  No potential conflict of interest was reported by the author.

 References Almani, H.M., Hosseini, S., Tahmasebi, M.: Fractional Brownian motion with two-variable Hurst exponent. J. Comput. Appl. Math. 388, 113262 (2021). MR4199791. https://doi.org/10.1016/j.cam.2020.113262 Arcones, M.A.: Limit theorems for nonlinear functionals of a stationary Gaussian sequence of vectors. Ann. Probab. 2242–2274 (1994). MR1331224 Beran, J.: Statistics for Long-Memory Processes. Chapman & Hall (1994). MR1304490 Cai, C., Chigansky, P., Kleptsyna, M.: The maximum likelihood drift estimator for mixed fractional Brownian motion (2012). Preprint, available online at arXiv:1208.6253. Chaouch, H., Maroufy, H., Omari, M.: Statistical inference for models driven by n-th order fractional Brownian motion. Theory Probab. Math. Stat. 108, 29–43 (2023). MR4588238. https://doi.org/10.1090/tpms/1185 Cheridito, P.: Mixed fractional Brownian motion. Bernoulli 7(6), 913–934 (2001). MR1873835. https://doi.org/10.2307/3318626 Dozzi, M., Mishura, Y., Shevchenko, G.: Asymptotic behavior of mixed power variations and statistical estimation in mixed models. Stat. Inference Stoch. Process. 18(2), 151–175 (2015). MR3348583. https://doi.org/10.1007/s11203-014-9106-5 Dufitinema, J., Pynnönen, S., Sottinen, T.: Maximum likelihood estimators from discrete data modeled by mixed fractional Brownian motion with application to the Nordic stock markets. Commun. Stat., Simul. Comput. 51(9), 5264–5287 (2022). MR4491681. https://doi.org/10.1080/03610918.2020.1764581 El Omari, M.: Mixtures of higher-order fractional Brownian motions. Commun. Stat., Theory Methods 52(12), 4200–4215 (2023). MR4574931. https://doi.org/10.1080/03610926.2021.1986541 El Omari, M.: An α-order fractional Brownian motion with Hurst index H ∈ ( 0 , 1 ) and α ∈ R + . Sankhya A 85(1), 572–599 (2023). MR4540809. https://doi.org/10.1007/s13171-021-00266-z El Omari, M.: Parameter estimation for nth-order mixed fractional Brownian motion with polynomial drift. J. Korean Stat. Soc. 52(2), 450–461 (2023). MR4594997. https://doi.org/10.1007/s42952-023-00209-4 El Omari, M.: On the Gaussian Volterra processes with power-type kernels. Stoch. Models 40(1), 152–165 (2024). MR4694530. https://doi.org/10.1080/15326349.2023.2212763 Gohberg, I., Goldberg, S.: Basic Operator Theory. Birkhäuser (1981). MR0632943 Guo, J., Chen, Y., Hu, B., Yan, L.-Q., Guo, Y., Liu, Y.: Fractional Gaussian fields for modeling and rendering of spatially-correlated media. ACM Trans. Graph. 38(4), 1–13 (2019) Gupta, A., Joshi, S.D., Prasad, S.: A new approach for estimation of statistically matched wavelet. IEEE Trans. Signal Process. 53(5), 1778–1793 (2005). MR2148604. https://doi.org/10.1109/TSP.2005.845470 Han, X.: A Gladyshev theorem for trifractional Brownian motion and n-th order fractional Brownian motion. Electron. Commun. Probab. 26, 1–12 (2021). MR4319988. https://doi.org/10.1214/21-ecp422 Jennane, R., Harba, R.: Fractional Brownian motion: A model for image texture. In: EUSIPCO, Signal Processing, vol. 3, pp. 1389–1392 (1994) Kukush, A., Lohvinenko, S., Mishura, Y., Ralchenko, K.: Two approaches to consistent estimation of parameters of mixed fractional Brownian motion with trend. In: Statistical Inference for Stochastic Processes, vol. 25, pp. 159–187 (2022). MR4419677. https://doi.org/10.1007/s11203-021-09252-6 Liang, W., Yiming, D.: Wavelet-based estimator for the Hurst parameters of fractional Brownian sheet. Acta Math. Sci. 37(1), 205–222 (2017). MR3581537. https://doi.org/10.1016/S0252-9602(16)30126-6 Liptsert, R.S., Shiryaev, A.N.: Statistics of Random Processes I. General Theory, vol. 394. Springer (2001). MR1800857 Lundahl, T., Ohley, W.J., Kay, S.M., Siffert, R.: Fractional Brownian motion: A maximum likelihood estimator and its application to image texture. IEEE Trans. Med. Imaging 5(3), 152–161 (1986). https://doi.org/10.1109/TMI.1986.4307764 Mandelbrot, B.B., Van Ness, J.W.: Fractional Brownian motions, fractional noises and applications. SIAM Rev. 10(4), 422–437 (1968). MR0242239. https://doi.org/10.1137/1010093 Mishura, Y.: Maximum likelihood drift estimation for the mixing of two fractional Brownian motions. In: Stochastic and Infinite Dimensional Analysis, pp. 263–280. Springer (2016). MR3708386. https://doi.org/10.1007/978-3-319-07245-6_14 Morales, C.J.: Wavelet-based multifractal spectra estimation: Statistical aspects and applications. Boston University (2002). MR2703146 Pentland, A.P.: Fractal-based description of natural scenes. IEEE Trans. Pattern Anal. Mach. Intell. 6, 661–674 (1984). https://doi.org/10.1109/TPAMI.1984.4767591 Perrin, E., Harba, R., Berzin-Joseph, C., Iribarren, I., Bonami, A.: nth-order fractional Brownian motion and fractional Gaussian noises. IEEE Trans. Signal Process. 49(5), 1049–1059 (2001). https://doi.org/10.1109/78.917808 Ralchenko, K., Yakovliev, M.: Asymptotic normality of parameter estimators for mixed fractional Brownian motion with trend. Austrian J. Stat. 52(SI), 127–148 (2023). https://doi.org/10.17713/ajs.v52iSI.1770 Ralchenko, K., Yakovliev, M.: Parameter estimation for fractional mixed fractional Brownian motion based on discrete observations. Mod. Stoch. Theory Appl. 11(1), 1–29 (2024). MR4692015 Shiryaev, A.N.: Probability-1, vol. 95. Springer (1996). MR1368405. https://doi.org/10.1007/978-1-4757-2539-1 Sottinen, T., Viitasaari, L.: Transfer principle for n-th order fractional Brownian motion with applications to prediction and equivalence in law. Theory Probab. Math. Stat. 98, 199–216 (2019). MR3824687. https://doi.org/10.1090/tpms/1071 Van der Vaart, A.W.: Asymptotic Statistics, vol. 3. Cambridge University Press (1998). MR1652247. https://doi.org/10.1017/CBO9780511802256 Xiao, W.-L., Zhang, W.-G., Zhang, X.-L.: Maximum-likelihood estimators in the mixed fractional Brownian motion. Statistics 45(1), 73–85 (2011). MR2772157. https://doi.org/10.1080/02331888.2010.541254 Zabat, M., Vayer-Besançon, M., Harba, R., Bonnamy, S., Van Damme, H.: Surface topography and mechanical properties of smectite films. In: Trends in Colloid and Interface Science XI, pp. 96–102 (1997). Zhang, P., Sun, Q., Xiao, W.-L.: Parameter identification in mixed Brownian–fractional Brownian motions using Powell’s optimization algorithm. Econ. Model. 40, 314–319 (2014). https://doi.org/10.1016/j.econmod.2014.04.026