VMSTA Modern Stochastics: Theory and Applications 2351-60542351-60462351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA234 10.15559/23-VMSTA234 Research Article Parameter estimation for fractional mixed fractional Brownian motion based on discrete observations https://orcid.org/0000-0001-7208-3130 RalchenkoKostiantynkostiantynralchenko@knu.ua https://orcid.org/0000-0003-0153-3047 YakovlievMykytamykyta.yakovliev@gmail.com∗ Taras Shevchenko National University of Kyiv, Kyiv, Ukraine Corresponding author. 2024 512202311112938202321120232112023 © 2024 The Author(s). Published by VTeX2024 Open access article under the CC BY license.

The object of investigation is the mixed fractional Brownian motion of the form X t = κ B t H 1 + σ B t H 2 , driven by two independent fractional Brownian motions B 1 H and B 2 H with Hurst parameters H 1 < H 2 . Strongly consistent estimators of unknown model parameters ( H 1 , H 2 , κ 2 , σ 2 ) ⊤ are constructed based on the equidistant observations of a trajectory. Joint asymptotic normality of these estimators is proved for 0 < H 1 < H 2 < 3 4 .

 Fractional Brownian motion mixed model strong consistency ergodic theorem asymptotic normality 60G22 62F10 62F12 Introduction

The paper is devoted to the statistical identification of the following stochastic process: X t = κ B t H 1 + σ B t H 2 , t ≥ 0 , where B H 1 and B H 2 are two independent fractional Brownian motions with Hurst indices H 1 and H 2 respectively. The process (1) is known in the literature as fractional mixed fractional Brownian motion [16]. We aim to estimate all four unknown parameters H 1 , H 2 , κ 2 , σ 2 based on equidistant observations { X k h , k = 0 , 1 , 2 , … }, h > 0. In order to get an identifiable model, we assume that 0 < H 1 < H 2 < 1, σ > 0, and κ > 0 throughout the article.

The application of fractional Brownian motion is highly widespread for modelling real-world processes that evolve over time and have the properties of self-similarity and long- or short-range dependence, see e.g., [4, 5, 11, 24–26, 32]. Mixed models like (1) containing fractional Brownian motions with different Hurst indices are important from the perspective of practical applications, in particular, to economics and financial mathematics; the motion of a lower Hurst index usually corresponds to a component responsible for changes over a short period of time, and the motion of a larger one corresponds to a long-term component. The particular case of the process (1) with H 1 = 1 / 2 is known as mixed fractional Brownian motion. It was firstly proposed by Cheridito [3] for better financial modeling, its properties can be found in [33]. There are already existing practical applications of this model [9, 27, 30, 32]. The model (1) with arbitrary H 1 and H 2 was considered in [8, 16]. Its application to the pricing of credit default swaps was studied in [10]. A more general model, containing several fractional Brownian motions with different Hurst indices can be found in [28].

Various methods for parameter estimation in the mixed Brownian–fractional Brownian model (i.e., in the case H 1 = 1 / 2) were proposed in [29] (the maximum likelihood method), [31] (Powell’s optimization algorithm), and [6] (approach based on power variations). The general case of H 1 , H 2 ∈ ( 0 , 1 ) is less studied. We can mention only the papers [18, 19], where the following mixed model with a trend was considered: X t = θ t + κ B t H 1 + σ B t H 2 , t ≥ 0 . The authors applied the maximum likelihood approach for the estimation of the unknown drift parameter θ in model (2) based on the continuous observations. They constructed an estimator, which involves a solution to the integral equation with weak polar kernel. An approach to parameter estimation based on the numerical solution to this equation was recently proposed in [21]. Estimation of the drift parameter by discrete observations was discussed in [20].

In Refs. [2, 7, 14, 17] the parameter estimation for model (2) with H 1 = 1 / 2 was studied. The simultaneous estimation of all four parameters of this model was studied in [7] by the maximum likelihood method, and in [14], where the approach of [6] was generalized and, moreover, ergodic-type estimators were proposed. Ref. [2] considered the estimation of the drift parameter θ assuming that H and σ are known and κ = 1. The drift parameter estimation for a more general model driven by a stationary Gaussian process was discussed in [17].

The goal of the present paper is to construct strongly consistent estimators and to study their joint asymptotic normality. We extend the approach of [14], based on the ergodic theorem and generalized method of moments. It is worth mentioning that asymptotic normality does not hold for H 2 > 3 4 , which is typical for models with fractional Brownian motion, see, e.g., [22]. The estimators remain strongly consistent but have non-Gaussian asymptotic distribution.

The paper is organized as follows. In Section 2 we construct strongly consistent estimators with the use of four statistics based on the various increments of the process X. Section 3 is devoted to asymptotic normality. First, we prove asymptotic normality of constructed additional statistics and evaluate their asymptotic covariance matrix. Thereafter, we obtain the main result on the asymptotic normality of our estimators by applying the delta method. The theoretical results are illustrated by numerical study, which is presented in Section 4. All technical proofs together with some auxiliary results are collected in Section 5.

We use the following notation. The symbol cov stands for the covariance of two random variables and for the covariance matrix of a random vector. In the paper, all the vectors are column ones. The superscript ⊤ denotes transposition. Weak convergence in distribution is denoted by → d .

 Construction of strongly consistent estimators

Let h > 0 be fixed. Assume that the process X t defined in (1) is observed at points t k = k h, k = 0 , 1 , 2 , … . Let us denote the increment Δ X k : = X ( k + 1 ) h − X k h = κ Δ B k H 1 + σ Δ B k H 2 . The estimation procedure will be based on the following result.

The Gaussian sequence Δ X k , k = 0 , 1 , … , is stationary and ergodic.

This lemma allows us to apply the ergodic theorem for the construction of estimators. Namely, if g : R l + 1 → R is a Borel function such that E | g ( Δ X 0 , Δ X 1 , … , Δ X l ) | < ∞, then 1 N ∑ k = 0 N − 1 g ( Δ X k , Δ X k + 1 , … , Δ X k + l ) → E g ( Δ X 0 , Δ X 1 , … , Δ X l ) a. s., as N → ∞. The main idea is to obtain four different convergences by choosing different functions g, and then to construct the estimators by solving the corresponding system of four equations. To this end, we introduce the following four statistics ξ N : = 1 N ∑ k = 0 N − 1 X ( k + 1 ) h − X k h 2 , η N : = 1 N ∑ k = 0 N − 1 X ( k + 2 ) h − X k h 2 , ζ N : = 1 N ∑ k = 0 N − 1 X ( k + 4 ) h − X k h 2 , ϕ N : = 1 N ∑ k = 0 N − 1 X ( k + 8 ) h − X k h 2 . Denote also ϑ = ( H 1 , H 2 , κ 2 , σ 2 ) ⊤ .

The following a.s. convergences hold as N → ∞:

ξ N → E ξ 0 = κ 2 h 2 H 1 + σ 2 h 2 H 2 = : f 1 ( ϑ ),

η N → E η 0 = κ 2 h 2 H 1 2 2 H 1 + σ 2 h 2 H 2 2 2 H 2 = : f 2 ( ϑ ),

ζ N → E ζ 0 = κ 2 h 2 H 1 2 4 H 1 + σ 2 h 2 H 2 2 4 H 2 = : f 3 ( ϑ ),

ϕ N → E ϕ 0 = κ 2 h 2 H 1 2 6 H 1 + σ 2 h 2 H 2 2 6 H 2 = : f 4 ( ϑ ).

Let us introduce the notation f ( ϑ ) : = f 1 ( ϑ ) , f 2 ( ϑ ) , f 3 ( ϑ ) , f 4 ( ϑ ) ⊤ , τ N = ( ξ N , η N , ζ N , ϕ N ) ⊤ , N ∈ N , τ 0 = ( E ξ 0 , E η 0 , E ζ 0 , E ϕ 0 ) ⊤ . Then f ( ϑ ) = τ 0 . Therefore, in order to obtain estimators for parameters ( H 1 , H 2 , κ 2 , σ 2 ) ⊤ we need to construct function q such that q ( τ 0 ) = f − 1 ( τ 0 ), and then required estimators would be defined as ( H ˆ N ( 1 ) , H ˆ N ( 2 ) , κ ˆ N 2 , σ ˆ N 2 ) ⊤ : = q ( τ N ). In other words, in order to construct strongly consistent estimators, we need to solve the estimation equation f ( ϑ ) = τ N with respect to unknown parameters. The solution to (7) can be found in explicit form, see Subsection 5.3 for details. This leads to the following result.

Let 0 < H 1 < H 2 < 1. The statistic ϑ ˆ : = ( H ˆ N ( 1 ) , H ˆ N ( 2 ) , κ ˆ N 2 , σ ˆ N 2 ) ⊤ , where H ˆ N ( 1 ) = 1 2 log 2 + η N ζ N − ξ N ϕ N + d N 2 ( η N 2 − ξ N ζ N ) , H ˆ N ( 2 ) = 1 2 log 2 + η N ζ N − ξ N ϕ N − d N 2 ( η N 2 − ξ N ζ N ) , κ ˆ N 2 = 2 ( η N 2 − ξ N ζ N ) η N ζ N − ξ N ϕ N + d N log 2 h ξ N d N + 2 η N 3 − 3 ξ N η N ζ N + ξ N 2 ϕ N 2 d N , σ ˆ N 2 = 2 ( η N 2 − ξ N ζ N ) η N ζ N − ξ N ϕ N − d N log 2 h ξ N d N − 2 η N 3 + 3 ξ N η N ζ N − ξ N 2 ϕ N 2 d N with d N : = ( ξ N 2 ϕ N 2 − 6 ξ N η N ζ N ϕ N − 3 η N 2 ζ N 2 + 4 η N 3 ϕ N + 4 ξ N ζ N 3 ) + 1 / 2 is a strongly consistent estimator of ϑ = ( H 1 , H 2 , κ 2 , σ 2 ) ⊤ .

Here log 2 + x = log 2 x , if x > 0 , 0 , if x ≤ 0 , ( x ) + 1 / 2 = x , if x > 0 , 0 , if x ≤ 0 .

The constructed estimators can be generalized for mixed model with trend (2) by introducing the estimator X N h N = 1 N ∑ k = 0 N − 1 ( X ( k + 1 ) h − X k h ) to define the estimator of θ. By subtracting X N h N from previously constructed statistics (4) we can obtain a system of equations similar to (7).

 Asymptotic normality

We will start with a study of the asymptotic properties of the vector τ N = ( ϕ N , ξ N , η N , ζ N ). The previously defined stationary Gaussian sequence Δ X k has the autocovariance function ρ ˜ ( i ) = cov ( Δ X 0 , Δ X i ) = κ 2 h 2 H 1 ρ ( i , H 1 ) + σ 2 h 2 H 2 ρ ( i , H 2 ) , i ∈ Z , where ρ ( i , H ) : = 1 2 | i + 1 | 2 H − 2 | i | 2 H + | i − 1 | 2 H denotes the autocovariance function of the stationary sequence { B k + 1 H − B k H , k ≥ 0 }, which is known as fractional Gaussian noise.

Let 0 < H 1 < H 2 < 3 4 . The vector τ N defined by (6) is asymptotically normal, namely, N ( τ N − τ 0 ) = N ξ N − E ξ 0 η N − E η 0 ζ N − E ζ 0 ϕ N − E ϕ 0 → d N ( 0 → , Σ ˜ ) with the asymptotic covariance matrix Σ ˜, which can be presented explicitly as Σ ˜ = Σ ˜ 11 Σ ˜ 12 Σ ˜ 13 Σ ˜ 14 Σ ˜ 12 Σ ˜ 22 Σ ˜ 23 Σ ˜ 24 Σ ˜ 13 Σ ˜ 23 Σ ˜ 33 Σ ˜ 34 Σ ˜ 14 Σ ˜ 24 Σ ˜ 34 Σ ˜ 44 , where Σ ˜ 11 = 2 ∑ i = − ∞ ∞ ρ ˜ ( i ) 2 , Σ ˜ 12 = ∑ i = − ∞ ∞ ρ ˜ ( i ) ( 4 ρ ˜ ( i ) + 4 ρ ˜ ( i + 1 ) ) , Σ ˜ 22 = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 12 ρ ˜ ( i ) + 16 ρ ˜ ( i + 1 ) + 4 ρ ˜ ( i + 2 ) ) , Σ ˜ 13 = ∑ i = − ∞ ∞ ρ ˜ ( i ) ( 8 ρ ˜ ( i ) + 12 ρ ˜ ( i + 1 ) + 8 ρ ˜ ( i + 2 ) + 4 ρ ˜ ( i + 3 ) ) , Σ ˜ 23 = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 28 ρ ˜ ( i ) + 48 ρ ˜ ( i + 1 ) + 32 ρ ˜ ( i + 2 ) + 16 ρ ˜ ( i + 3 ) + 4 ρ ˜ ( i + 4 ) ) , Σ ˜ 33 = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 88 ρ ˜ ( i ) + 160 ρ ˜ ( i + 1 ) + 124 ρ ˜ ( i + 2 ) + 80 ρ ˜ ( i + 3 ) + 40 ρ ˜ ( i + 4 ) + 16 ρ ˜ ( i + 5 ) + 4 ρ ˜ ( i + 6 ) ) , Σ ˜ 14 = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 16 ρ ˜ ( i ) + 28 ρ ˜ ( i + 1 ) + 24 ρ ˜ ( i + 2 ) + 20 ρ ˜ ( i + 3 ) + 16 ρ ˜ ( i + 4 ) + 12 ρ ˜ ( i + 5 ) + 8 ρ ˜ ( i + 6 ) + 4 ρ ˜ ( i + 7 ) ) , Σ ˜ 24 = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 60 ρ ˜ ( i ) + 112 ρ ˜ ( i + 1 ) + 96 ρ ˜ ( i + 2 ) + 80 ρ ˜ ( i + 3 ) + 64 ρ ˜ ( i + 4 ) + 48 ρ ˜ ( i + 5 ) + 32 ρ ˜ ( i + 6 ) + 16 ρ ˜ ( i + 7 ) + 4 ρ ˜ ( i + 8 ) ) , Σ ˜ 34 = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 216 ρ ˜ ( i ) + 416 ρ ˜ ( i + 1 ) + 376 ρ ˜ ( i + 2 ) + 320 ρ ˜ ( i + 3 ) + 256 ρ ˜ ( i + 4 ) + 192 ρ ˜ ( i + 5 ) + 132 ρ ˜ ( i + 6 ) + 80 ρ ˜ ( i + 7 ) + 40 ρ ˜ ( i + 8 ) + 16 ρ ˜ ( i + 9 ) + 4 ρ ˜ ( i + 10 ) ) , Σ ˜ 44 = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 688 ρ ˜ ( i ) + 1344 ρ ˜ ( i + 1 ) + 1260 ρ ˜ ( i + 2 ) + 1136 ρ ˜ ( i + 3 ) + 984 ρ ˜ ( i + 4 ) + 816 ρ ˜ ( i + 5 ) + 644 ρ ˜ ( i + 6 ) + 480 ρ ˜ ( i + 7 ) + 336 ρ ˜ ( i + 8 ) + 224 ρ ˜ ( i + 9 ) + 140 ρ ˜ ( i + 10 ) + 80 ρ ˜ ( i + 11 ) + 40 ρ ˜ ( i + 12 ) + 16 ρ ˜ ( i + 13 ) + 4 ρ ˜ ( i + 14 ) ) .

The condition H < 3 4 is essential since otherwise the series in Theorem 2 do not converge. Moreover, for H > 3 4 , the asymptotic normality of τ N does not hold; using the results of [15], it is possible to establish the convergence of the suitably normalized vector τ N − τ 0 to a special non-Gaussian distribution.

Now we are ready to provide the main result on asymptotic normality of our estimator (8).

Let 0 < H 1 < H 2 < 3 4 . The estimator ϑ ˆ N is asymptotically normal, namely, N ϑ ˆ N − ϑ = N H ˆ N ( 1 ) − H 1 H ˆ N ( 2 ) − H 2 κ ˆ N 2 − κ 2 σ ˆ N 2 − σ 2 → d N 0 → , Σ 0 with the asymptotic covariance matrix Σ 0 that can be found by the formula Σ 0 = ( f ′ ( ϑ ) ) − 1 Σ ˜ f ′ ( ϑ ) − ⊤ , where Σ ˜ is defined in Theorem 2 and f ′ ( ϑ ) = 2 κ 2 h 2 H 1 log h 2 σ 2 h 2 H 2 log h h 2 H 1 h 2 H 2 2 κ 2 ( 2 h ) 2 H 1 log ( 2 h ) 2 σ 2 ( 2 h ) 2 H 2 log ( 2 h ) ( 2 h ) 2 H 1 ( 2 h ) 2 H 2 2 κ 2 ( 4 h ) 2 H 1 log ( 4 h ) 2 σ 2 ( 4 h ) 2 H 2 log ( 4 h ) ( 4 h ) 2 H 1 ( 4 h ) 2 H 2 2 κ 2 ( 8 h ) 2 H 1 log ( 8 h ) 2 σ 2 ( 8 h ) 2 H 2 log ( 8 h ) ( 8 h ) 2 H 1 ( 8 h ) 2 H 2 .

Simulation study

In this section, we present results of estimators performance in numerical simulations. For each generated trajectory, we will estimate asymptotic covariance matrices defined in Theorem 2 and Theorem 3 by using values of estimators (9)–(12). For each set of parameters with 0 < H 1 < H 2 < 1, we generate 1000 trajectories of the process X and calculate the empirical means, empirical standard deviations of the estimates, and percentage of iterations with constant d N defined in (13) being equal to zero. Additionally, for 0 < H 1 < H 2 < 3 4 , we calculate the average approximate theoretical standard deviation (ASD) – the square root of the estimated asymptotic variance divided by N ( ( Σ ˆ N 0 ) i i / N ) – and the coverage probability (CP) for α = 5 %, based on the estimator of the asymptotic covariance matrix. We are computing these statistics only for H 1 < H 2 < 3 4 , excluding interval ( 3 4 , 1 ) where the series (17) will not converge as per Lemma 3. Therefore, the asymptotic covariance matrix and corresponding statistics (ASD and CP) will be undefined. Moreover, when any of the estimates H ˆ N will be out of boundaries ( 0 , 3 4 ), or the calculated constant d N defined in (13) is zero (e.g., in cases where discriminant D defined in the proof of Theorem 1 by (30) is nonpositive), then ASD and CP will also be considered undefined.

To numerically approximate limiting covariances, series alike (17) are divided into the sum of two convergent series ∑ i = − ∞ + ∞ ρ ˜ ( i + α ) ρ ˜ ( i + β ) = ∑ i = 0 + ∞ ρ ˜ ( i + α ) ρ ˜ ( i + β ) + ∑ i = 0 + ∞ ρ ˜ ( i + 1 − α ) ρ ˜ ( i + 1 − β ) , where α , β ∈ Z.

Then, series ∑ i = 0 + ∞ ρ ˜ ( i + α ) ρ ˜ ( i + β ) are approximated with precision δ = 10 − 3 by evaluating the sum S m = ∑ i = 0 m ρ ˜ ( i + α ) ρ ˜ ( i + β ), where m is the smallest number such that | S m − S m − 1 | < 10 − 3 . Obtained approximations will be used to compute the asymptotic covariance matrix defined in Theorem 2.

For H ˆ N ( 1 ) , H ˆ N ( 2 ) , κ ˆ N 2 , σ ˆ N 2 , Σ ˆ N 0 and fixed time step h = 1 we vary the horizon T = h 2 n for n ∈ { 8 , 10 , 12 , 14 , 16 , 18 , 20 }. For all simulations, the values σ = 1 and κ = 1 were used.

The estimator H ˆ N ( 1 ) with σ 2 = 1, κ 2 = 1 ( h = 1)

N
H 1 H 2 2 8 2 10 2 12 2 14 2 16 2 18 2 20
0.1 0.3 Mean 0.0722 0.0434 −0.0230 −0.1009 −0.0844 0.0482 0.0842
S.dev. 0.5041 0.4727 0.5063 0.5160 0.4459 0.1629 0.0791
ASD 0.7648 0.7039 1.1629 0.3871 0.2011 0.1239 0.0812
CP% 100.00 100.00 100.00 92.18 75.61 75.58 87.91
0.1 0.5 Mean 0.0184 −0.0455 −0.0009 0.0674 0.0953 0.0988 0.0991
S.dev. 0.5618 0.5945 0.3733 0.1710 0.0824 0.0414 0.0205
ASD 1.1951 0.5894 0.2696 0.1457 0.0829 0.0439 0.0222
CP% 100.00 100.00 99.30 89.82 93.76 96.85 96.50
0.1 0.7 Mean −0.0408 −0.0071 0.0773 0.0945 0.1002 0.0998 0.1005
S.dev. 0.5931 0.3828 0.1513 0.0715 0.0342 0.0179 0.0090
ASD 0.7428 0.3223 0.1516 0.0742 0.0375 0.0188 0.0094
CP% 100.00 100.00 100.00 99.30 97.09 95.90 95.40
0.1 0.9 Mean 0.0131 0.0725 0.0886 0.0956 0.0972 0.0986 0.0993
S.dev. 0.3842 0.1554 0.0707 0.0350 0.0187 0.0113 0.0075
0.3 0.5 Mean 0.2096 0.1797 0.1161 0.0665 0.1731 0.2658 0.2901
S.dev. 0.5181 0.5738 0.4782 0.5254 0.3445 0.1209 0.0564
ASD 1.1101 1.0150 0.7507 0.4497 0.2337 0.1124 0.0570
CP% 100.00 100.00 100.00 100.00 92.01 87.11 90.80
0.3 0.7 Mean 0.1748 0.1524 0.2330 0.2764 0.2931 0.2989 0.3001
S.dev. 0.5138 0.4661 0.2630 0.1165 0.0563 0.0263 0.0137
ASD 1.4395 0.5732 0.2682 0.1248 0.0566 0.0277 0.0138
CP% 100.00 100.00 100.00 100.00 98.32 96.30 95.60
0.3 0.9 Mean 0.1863 0.2636 0.2863 0.2934 0.2969 0.2983 0.2988
S.dev. 0.4268 0.1859 0.0824 0.0423 0.0211 0.0124 0.0078
0.5 0.7 Mean 0.4213 0.3300 0.2774 0.3367 0.4563 0.4859 0.4974
S.dev. 0.6302 0.6258 0.5085 0.3908 0.1461 0.0763 0.0367
ASD 2.0105 1.2367 0.8386 0.4502 0.1964 0.0830 0.0366
CP% 100.00 100.00 100.00 100.00 99.84 99.09 98.09
0.5 0.9 Mean 0.3188 0.3701 0.4577 0.4878 0.4943 0.4975 0.4986
S.dev. 0.5170 0.4018 0.1477 0.0627 0.0299 0.0162 0.0095
0.7 0.9 Mean 0.5487 0.4918 0.5172 0.5854 0.6706 0.6898 0.6957
S.dev. 0.6543 0.6735 0.4569 0.2908 0.0914 0.0396 0.0193

The estimator H ˆ N ( 2 ) with σ 2 = 1, κ 2 = 1 ( h = 1)

N
H 1 H 2 2 8 2 10 2 12 2 14 2 16 2 18 2 20
0.1 0.3 Mean 0.2296 0.2845 0.3834 0.5157 0.4571 0.3634 0.3182
S.dev. 0.5205 0.6239 0.5643 0.5944 0.3578 0.1543 0.0594
ASD 2.2888 2.0419 1.4188 0.8483 0.4334 0.1937 0.0791
CP% 100.00 100.00 100.00 100.00 100.00 100.00 99.88
0.1 0.5 Mean 0.5437 0.6554 0.6072 0.5302 0.5093 0.5024 0.5003
S.dev. 0.6078 0.5234 0.2668 0.0985 0.0411 0.0193 0.0094
ASD 0.9381 0.5332 0.2632 0.1158 0.0481 0.0218 0.0107
CP% 100.00 100.00 100.00 100.00 99.66 98.37 96.80
0.1 0.7 Mean 0.7912 0.7303 0.7069 0.7018 0.7006 0.7000 0.7001
S.dev. 0.4123 0.1341 0.0554 0.0254 0.0126 0.0066 0.0032
ASD 0.3193 0.1257 0.0557 0.0265 0.0139 0.0064 0.0032
CP% 100.00 100.00 100.00 99.77 95.98 94.50 94.70
0.1 0.9 Mean 0.8818 0.8825 0.8882 0.8930 0.8959 0.8979 0.8991
S.dev. 0.1599 0.0601 0.0381 0.0272 0.0201 0.0151 0.0121
0.3 0.5 Mean 0.4179 0.4999 0.5838 0.6819 0.5997 0.5329 0.5097
S.dev. 0.6199 0.6409 0.5791 0.4784 0.2349 0.0921 0.0418
ASD 0.9471 0.7822 0.6198 0.3870 0.1737 0.0861 0.0412
CP% 100.00 100.00 100.00 100.00 97.81 90.36 90.90
0.3 0.7 Mean 0.7636 0.8383 0.7481 0.7094 0.7023 0.7007 0.7004
S.dev. 0.5831 0.4365 0.1412 0.0574 0.0271 0.0128 0.0066
ASD 0.4241 0.1793 0.0914 0.0466 0.0256 0.0130 0.0065
CP% 100.00 100.00 99.38 90.74 93.60 95.30 94.00
0.3 0.9 Mean 0.9267 0.8926 0.8898 0.8936 0.8954 0.8975 0.8981
S.dev. 0.2765 0.0827 0.0440 0.0301 0.0213 0.0166 0.0130
0.5 0.7 Mean 0.5823 0.6867 0.8179 0.8324 0.7494 0.7181 0.7056
S.dev. 0.6097 0.6288 0.5302 0.3556 0.1229 0.0581 0.0265
ASD 0.5033 0.2277 0.1997 0.1055 0.0672 0.0422 0.0245
CP% 100.00 100.00 99.54 82.13 78.62 85.38 92.58
0.5 0.9 Mean 0.8978 0.9350 0.8980 0.8936 0.8939 0.8965 0.8978
S.dev. 0.5267 0.2147 0.0793 0.0413 0.0274 0.0198 0.0150
0.7 0.9 Mean 0.7878 0.8862 0.9672 0.9309 0.9021 0.8963 0.8965
S.dev. 0.6431 0.5628 0.4037 0.1334 0.0578 0.0316 0.0213

The estimator κ ˆ N 2 with σ 2 = 1, κ 2 = 1 ( h = 1)

N
H 1 H 2 2 8 2 10 2 12 2 14 2 16 2 18 2 20
0.1 0.3 Mean NaN NaN NaN Inf 1.0591 1.0732 1.0460
S.dev. NaN NaN NaN NaN 0.7876 0.6474 0.4604
ASD 4.3622 4.9405 5.0730 2.8475 1.5443 0.9266 0.6165
CP% 100.00 100.00 93.08 80.66 70.73 74.27 89.53
0.1 0.5 Mean NaN Inf 1.1104 1.0594 1.0340 1.0098 1.0006
S.dev. NaN NaN 0.6161 0.4314 0.2528 0.1333 0.0672
ASD 4.4678 2.4422 1.0822 0.5453 0.2928 0.1484 0.0743
CP% 100.00 100.00 99.30 92.92 96.15 97.97 96.20
0.1 0.7 Mean Inf 1.0400 1.0183 1.0054 1.0035 1.0004 1.0010
S.dev. NaN 0.4440 0.2538 0.1260 0.0623 0.0330 0.0164
ASD 1.8115 0.6744 0.2911 0.1395 0.0693 0.0346 0.0173
CP% 100.00 100.00 100.00 100.00 97.79 95.70 95.70
0.1 0.9 Mean 1.0054 0.9872 0.9878 0.9922 0.9950 0.9972 0.9987
S.dev. 0.3366 0.1663 0.0823 0.0437 0.0264 0.0177 0.0132
0.3 0.5 Mean NaN NaN Inf 1.0451 1.0379 1.0466 1.0232
S.dev. NaN NaN NaN 0.8444 0.7341 0.5658 0.3768
ASD 4.8428 4.2202 3.2605 2.1396 1.0325 0.6043 0.3862
CP% 100.00 100.00 100.00 96.51 83.23 77.25 84.00
0.3 0.7 Mean NaN 1.1084 1.0696 1.0255 1.0106 1.0026 1.0017
S.dev. NaN NaN 0.5169 0.3305 0.1797 0.0877 0.0455
ASD 3.5137 1.3194 0.6535 0.3286 0.1770 0.0917 0.0461
CP% 100.00 100.00 0.9979 0.9074 0.9496 0.9570 0.9460
0.3 0.9 Mean Inf 0.9970 0.9867 0.9887 0.9923 0.9955 0.9968
S.dev. NaN 0.2908 0.1463 0.0802 0.0436 0.0285 0.0201
0.5 0.7 Mean NaN NaN Inf 1.0557 1.0674 1.0502 1.0270
S.dev. NaN NaN NaN 0.7737 0.6163 0.4567 0.2703
ASD 4.6051 2.2181 2.1464 1.2319 0.7646 0.4713 0.2722
CP% 100.00 100.00 92.17 77.05 78.14 86.29 93.96
0.5 0.9 Mean NaN 1.0008 0.9750 0.9800 0.9814 0.9905 0.9940
S.dev. NaN 0.5558 0.3666 0.1997 0.1057 0.0621 0.0402
0.7 0.9 Mean NaN NaN −Inf 0.9524 0.9568 0.9572 0.9690
S.dev. NaN NaN NaN 0.6566 0.4605 0.2777 0.1643

The estimator σ ˆ N 2 with σ 2 = 1, κ 2 = 1 ( h = 1)

N
H 1 H 2 2 8 2 10 2 12 2 14 2 16 2 18 2 20
0.1 0.3 Mean NaN NaN NaN Inf 0.9764 0.9499 0.9781
S.dev. NaN NaN NaN NaN 0.7882 0.6480 0.4606
ASD 4.3572 4.9387 5.0724 2.8472 1.5442 0.9265 0.6165
CP% 100.00 100.00 93.08 80.25 70.95 74.27 89.53
0.1 0.5 Mean NaN −Inf 0.8916 0.9408 0.9656 0.9902 0.9995
S.dev. NaN NaN 0.6168 0.4320 0.2527 0.1334 0.0672
ASD 4.4647 2.4412 1.0818 0.5452 0.2928 0.1484 0.0748
CP% 100.00 100.00 99.54 93.07 96.15 97.87 96.40
0.1 0.7 Mean −Inf 0.9609 0.9795 0.9941 0.9959 0.9994 0.9989
S.dev. NaN 0.4474 0.2513 0.1245 0.0617 0.0326 0.0162
ASD 1.8071 0.6716 0.2897 0.1388 0.0690 0.0344 0.0172
CP% 100.00 100.00 100.00 100.00 97.49 95.90 96.20
0.1 0.9 Mean 1.0170 1.0253 1.0273 1.0157 1.0099 1.0064 1.0067
S.dev. 0.5639 0.3708 0.2810 0.2096 0.1452 0.1039 0.0784
0.3 0.5 Mean NaN NaN −Inf 0.9554 0.9623 0.9536 0.9769
S.dev. NaN NaN NaN 0.8435 0.7340 0.5662 0.3768
ASD 4.8407 4.2199 3.2605 2.1396 1.0325 0.6043 0.3862
CP% 100.00 100.00 100.00 96.80 83.23 77.03 84.00
0.3 0.7 Mean NaN −Inf 0.9327 0.9886 0.9978 0.9984 0.9982
S.dev. NaN NaN 0.5156 0.3295 0.1787 0.0872 0.0453
ASD 3.5101 1.3180 0.6527 0.3283 0.1768 0.0915 0.0460
CP% 100.00 100.00 100.00 90.61 95.07 95.50 94.80
0.3 0.9 Mean −Inf 1.0177 1.0233 1.0235 1.0092 1.0092 1.0044
S.dev. NaN 0.4173 0.2704 0.1891 0.1256 0.1035 0.0803
0.5 0.7 Mean NaN NaN −Inf 0.9450 0.9323 0.9497 0.9730
S.dev. NaN NaN NaN 0.7735 0.6158 0.4568 0.2702
ASD 4.5998 2.2161 2.1461 1.2318 0.7645 0.4712 0.2721
CP% 100.00 100.00 91.71 77.05 78.30 86.42 94.06
0.5 0.9 Mean NaN 1.0025 1.0274 1.0277 1.0211 1.0122 1.0075
S.dev. NaN 0.6212 0.4051 0.2374 0.1476 0.1002 0.0726
0.7 0.9 Mean NaN NaN Inf 1.0399 1.0375 1.0414 1.0288
S.dev. NaN NaN NaN 0.6359 0.4416 0.2525 0.1337

Results of parameters estimations provided in Tables 1–4 indicate that constructed estimators perform better, the higher the difference is between H 1 and H 2 . That impact is especially significant for estimators of variances κ 2 and σ 2 since the denominator in (11) and (12) tends to 0 when the difference between two Hurst parameters tends to zero. Also, the smaller this difference is, the higher the percentage of iterations with constant d N being equal to zero, which is presented in Table 5. Additionally, we also observe that constructed estimators perform better for higher values of H 1 and H 2 . This also might be caused by the difference between H 1 and H 2 since most of the values in numerators and denominators of constructed estimators are highly dependent on the corresponding difference, as was shown in the proof of Theorem 1. And higher values of H 1 and H 2 might result in a more accurate estimate of the difference H ˆ N ( 2 ) − H ˆ N ( 1 ) in terms of the ratio between the estimated value and its deviation, and therefore – a lower percentage of iterations with the estimated difference being negative.

The percentage (%) of iterations with constant d N being equal to zero with σ 2 = 1, κ 2 = 1 ( h = 1)

N
H 1 H 2 2 8 2 10 2 12 2 14 2 16 2 18 2 20
0.1 0.3 2.9 1.3 0.3 0 0 0 0
0.1 0.5 1.1 0.1 0 0 0 0 0
0.1 0.7 0.2 0 0 0 0 0 0
0.1 0.9 0 0 0 0 0 0 0
0.3 0.5 3.8 1.5 0.1 0 0 0 0
0.3 0.7 0.9 0 0 0 0 0 0
0.3 0.9 0.1 0 0 0 0 0 0
0.5 0.7 3.5 2.2 0.1 0 0 0 0
0.5 0.9 0.6 0 0 0 0 0 0
0.7 0.9 4.3 1.1 0.4 0 0 0 0

The coverage probability (%) for estimators ( ξ N , η N , ζ N , ϕ N ) ⊤ with σ 2 = 1, κ 2 = 1 ( h = 1)

N
H 1 H 2 2 8 2 10 2 12 2 14 2 16 2 18 2 20
0.1 0.3 ξ N 94.29 92.19 95.38 94.24 94.68 92.25 93.14
η N 100.0 98.44 99.23 99.59 98.67 99.27 98.95
ζ N 100.0 100.0 100.0 100.0 100.0 100.0 99.88
ϕ N 100.0 100.0 100.0 100.0 100.0 100.0 100.0
0.1 0.5 ξ N 95.06 94.55 94.20 94.25 93.54 94.41 95.10
η N 98.77 97.52 97.45 99.12 97.73 98.27 98.50
ζ N 97.53 100.0 98.84 99.41 98.75 99.19 99.20
ϕ N 100.0 100.0 99.54 99.26 99.55 99.19 99.40
0.1 0.7 ξ N 95.83 95.11 96.33 95.24 94.48 94.00 94.40
η N 93.75 95.49 97.10 94.66 94.68 94.10 93.70
ζ N 90.62 95.11 93.82 91.88 91.77 90.00 91.10
ϕ N 88.54 92.48 91.31 89.91 90.26 87.50 88.60
0.3 0.5 ξ N 95.00 95.28 95.21 95.64 94.36 95.77 94.80
η N 95.00 96.23 97.01 96.51 97.34 97.29 96.50
ζ N 96.67 98.11 98.20 97.67 99.06 98.70 97.80
ϕ N 100.0 99.06 97.01 98.84 98.90 98.81 97.90
0.3 0.7 ξ N 96.15 94.21 94.24 92.70 93.07 94.00 93.70
η N 98.72 93.82 94.65 92.18 91.71 91.80 91.40
ζ N 96.15 93.05 91.98 90.61 90.35 89.40 89.00
ϕ N 93.59 93.44 89.51 87.22 88.88 86.90 85.70
0.5 0.7 ξ N 89.36 96.15 95.85 91.06 94.18 92.56 94.27
η N 87.23 97.69 92.17 87.68 91.04 89.16 90.88
ζ N 89.36 96.92 87.56 84.30 88.52 86.81 87.91
ϕ N 91.49 97.69 85.25 82.37 86.95 84.60 86.00

The percentage (%) of iterations with undefined asymptotic covariance matrix Σ ˆ N 0 with σ 2 = 1, κ 2 = 1 ( h = 1)

N
H 1 H 2 2 8 2 10 2 12 2 14 2 16 2 18 2 20
0.1 0.3 96.50 93.60 87.00 75.70 54.90 31.60 14.00
0.1 0.5 91.90 79.80 56.90 32.20 11.80 1.60 0
0.1 0.7 90.40 73.40 48.20 13.80 0.40 0 0
0.3 0.5 94.00 89.40 83.30 65.60 36.20 7.70 0
0.3 0.7 92.20 74.10 51.40 23.30 4.70 0 0
0.5 0.7 95.30 87.00 78.30 58.60 36.40 23.40 5.70

We observe that the higher the difference is between the true values of H 1 and H 2 , the closer the ASD is to the empirical standard deviation of the estimator and the closer the coverage probability to the theoretical 95 %. At the same time, results similar to those described above are observed in the performance of our estimates for asymptotic variance of (4) constructed in Theorem 2, which is shown in Table 6 by CP estimates. This table also shows that for higher values of H 1 and H 2 , the accuracy of (4) is the lower, the bigger the increment step (e.g., ϕ N is less accurate than ζ N , which is less accurate than η N , which is less accurate than ξ N ). Unfortunately, the percentage of iterations with undefined asymptotic covariance matrix Σ ˆ N 0 is high but decreases for bigger samples, which is shown in Table 7, due to estimated H ˆ N ( 1 ) , H ˆ N ( 2 ) violating the conditions of Theorem 3. Also, the percentage of iterations with the undefined Σ ˆ N 0 is the lower, the higher the difference between true H 1 and  H 2 .

Additionally, by conducting provided simulations, we observed that the closer parameters H 1 and H 2 to 0.75, the lower the convergence rate in the numerical approximation of the limiting covariance. To study this effect, we have conducted simulations for the series ∑ i = 0 + ∞ ρ ˜ ( i ) 2 . We estimate (16) by finding the smallest number m ∈ N such that for S m = ∑ i = 0 m ρ ˜ ( i ) 2 the inequality | S m − S m − 1 | ≤ δ holds for some predefined δ.

This estimation was conducted for δ ∈ { 10 − 3 , 10 − 4 , 10 − 5 , 10 − 6 , 10 − 7 , 10 − 8 , 10 − 9 } for σ = 1, κ = 1, H 1 = 0.1 and H 2 ∈ { 0.3 , 0.5 , 0.7 }. Corresponding estimation results are presented in Table 8. As one can see, in order to estimate S m with the precision δ = 10 − 9 for H 2 = 0.3, only a few hundred iterations are needed, while for H 2 = 0.7, it requires dozens of millions of iterations. Moreover, lowering the precision from 10 − 4 to 10 − 3 has almost no impact on the estimated value S m for low H 2 , but for high H 2 , it could result in around three times difference.

The convergence of the series (16) with σ 2 = 1, κ 2 = 1 and H 1 = 0.1 ( h = 1)

N
H 2 10 − 3 10 − 4 10 − 5 10 − 6 2 − 7 10 − 8 10 − 9
0.3 m 4 8 16 35 76 168 377
S m 4.45362 4.45427 4.45445 4.45451 4.45452 4.45452 4.45453
0.5 m 3 4 7 12 22 42 78
S m 4.18198 4.18203 4.18206 4.18207 4.18208 4.18208 4.18208
0.7 m 5.9 · 10 6 4.9 · 10 7 7.2 · 10 7 7.2 · 10 7 7.2 · 10 7 7.2 · 10 7 7.2 · 10 7
S m 9981.86 35818.4 45796.6 45796.6 45796.6 45796.6 45796.6
Proofs Auxiliary properties of covariance functions

The following statement provides some important properties of sequences ρ ˜ ( i ) and ρ ( i , H ), i ∈ Z, defined by (14) and (15), respectively.

For any H ∈ ( 0 , 1 ), ρ ( i , H ) ∼ H ( 2 H − 1 ) | i | 2 H − 2 as i → ∞.

For any H ∈ ( 0 , 3 4 ), ∑ i = − ∞ ∞ ρ 2 ( i , H ) < ∞.

For any 0 < H 1 < H 2 < 3 4 and for all α , β ∈ Z, ∑ i = − ∞ + ∞ ρ ˜ ( i + α ) ρ ˜ ( i + β ) < ∞ .

First two statements are well known, see, e.g., [22]. In order to prove the third one, we observe that ∑ i = − ∞ ∞ ρ ˜ 2 ( i ) ≤ 2 κ 4 h 4 H 1 ∑ i = − ∞ ∞ ρ 2 ( i , H 1 ) + 2 σ 4 h 4 H 2 ∑ i = − ∞ ∞ ρ 2 ( i , H 2 ) < ∞ , by the second statement. Then, using the Cauchy–Schwarz inequality, we get for all α , β ∈ Z, ∑ i = − ∞ + ∞ ρ ˜ ( i + α ) ρ ˜ ( i + β ) ≤ ∑ i = − ∞ + ∞ ρ ˜ 2 ( i + α ) · ∑ i = − ∞ + ∞ ρ ˜ 2 ( i + β ) = ∑ i = − ∞ + ∞ ρ ˜ 2 ( i ) < ∞ .  □

 Proofs of ergodic results

Since B H 1 and B H 2 are independent centered Gaussian processes with stationary increments, we see that { Δ X k , k ≥ 0 } is a stationary Gaussian sequence with E Δ X k = 0. Therefore, in order to establish the ergodicity of Δ X k , it suffices to prove that its autocovariance function ρ ˜ ( k ) vanishes at infinity. In turn, this fact follows immediately from the first statement of Lemma 3: ρ ˜ ( k ) = κ 2 h 2 H 1 ρ ( k , H 1 ) + σ 2 h 2 H 2 ρ ( k , H 2 ) ∼ κ 2 h 2 H 1 H 1 ( 2 H 1 − 1 ) k ( 2 H 1 − 2 ) + σ 2 h 2 H 2 H 2 ( 2 H 2 − 1 ) k ( 2 H 2 − 2 ) → 0 , as k → ∞.  □

( i ) Taking g ( x 0 ) = x 0 2 in (3), we obtain that ξ N → E X h 2 = E κ B h H 1 + σ B h H 2 2 = κ 2 E ( B h H 1 ) 2 + σ 2 E ( B h H 2 ) 2 = κ 2 h 2 H 1 + σ 2 h 2 H 2 .

( ii ) Applying g ( x 0 , x 1 ) = ( x 0 + x 1 ) 2 to (3), we arrive at η N → E X 2 h 2 = E κ B 2 h H 1 + σ B 2 h H 2 2 = κ 2 E ( B 2 h H 1 ) 2 + σ 2 E ( B 2 h H 2 ) 2 = κ 2 h 2 H 1 2 2 H 1 + σ 2 h 2 H 2 2 2 H 2 .

Statements ( iii ) and ( iv ) can be proved similarly; we need to choose g ( x 0 , x 1 , x 2 , x 3 ) = ( x 0 + x 1 + x 2 + x 3 ) 2 and g ( x 0 , … , x 7 ) = ( ∑ i = 0 7 x i ) 2 for (3), respectively.  □

 Construction of estimators. Proof of Theorem <xref rid="j_vmsta234_stat_003">1</xref>

Let us briefly explain, how to solve the estimating equation (7). Denote u = κ 2 h 2 H 1 , v = σ 2 h 2 H 2 , x = 2 2 H 1 , y = 2 2 H 2 , ξ = ξ N , η = η N , ζ = ζ N , ϕ = ϕ N . Then (7) takes the form ξ = u + v , η = u x + v y , ζ = u x 2 + v y 2 , ϕ = u x 3 + v y 3 . Using equations (18) and (19) we can express the variables u and v in terms of x and y as follows: v = η − ξ x y − x , u = η − ξ y x − y . Combining (22) with (20) we get ζ = x 2 · η − ξ y x − y + y 2 · η − ξ x y − x = 1 y − x − η x 2 + ξ y x 2 + η y 2 − ξ x y 2 = 1 y − x − ξ ( y − x ) x y + η ( y 2 − x 2 ) = − ξ x y + η ( y + x ) . Similarly, from (21) and (22) we derive ϕ = x 3 · η − ξ y x − y + y 3 · η − ξ x y − x = 1 y − x − ξ ( y 2 − x 2 ) x y + η ( y 3 − x 3 ) = − ξ ( y + x ) x y + η ( y 2 + x y + x 2 ) . Multiplying (23) by x + y results in ζ ( x + y ) = − ξ ( x + y ) x y + η ( y 2 + 2 x y + x 2 ) . Finally, by subtracting (24) from (25) we get ζ ( x + y ) − ϕ = η x y, whence x = ϕ − ζ y ζ − η y . At the same time, we can express x via y using (23): x = ζ − η y η − ξ y . From equations (26) and (27) we obtain ϕ − ζ y ζ − η y = ζ − η y η − ξ y , whence we get the following quadratic equation for y y 2 η 2 − ξ ζ − y η ζ − ξ ϕ + ζ 2 − η ϕ = 0 . It is easy to see that x satisfies the same quadratic equation, in view of symmetry. Therefore, x and y ( y > x) are two roots of (28), which are given by η ζ − ξ ϕ ± D 2 ( η 2 − ξ ζ ) , where D = ( η ζ − ξ ϕ ) 2 − 4 ( η 2 − ξ ζ ) ( ζ 2 − η ϕ ) = ξ 2 ϕ 2 − 6 ξ η ζ ϕ − 3 η 2 ζ 2 + 4 η 3 ϕ + 4 ξ ζ 3 . It is not hard to derive the following relations from equations (18)–(21): η ζ − ξ ϕ = ( u x + v y ) ( u x 2 + v y 2 ) − ( u + v ) ( u x 3 + v y 3 ) = − u v ( y − x ) 2 ( y + x ) ; η 2 − ξ ζ = ( u x + v y ) 2 − ( u + v ) ( u x 2 + v y 2 ) = − u v ( y − x ) 2 ; ζ 2 − η ϕ = ( u x 2 + v y 2 ) 2 − ( u x + v y ) ( u x 3 + v y 3 ) = − u v x y ( y − x ) 2 . Thus D = ( η ζ − ξ ϕ ) 2 − 4 ( η 2 − ξ ζ ) ( ζ 2 − η ϕ ) = u 2 v 2 ( y − x ) 4 ( y + x ) 2 − 4 u 2 v 2 x y ( y − x ) 4 = u 2 v 2 ( y − x ) 6 > 0 . Therefore, the discriminant is well defined and the corresponding equation has two different roots. From equation (29) η ζ − ξ ϕ ± D 2 ( η 2 − ξ ζ ) = − u v ( y − x ) 2 ( y + x ) ± u v ( y − x ) 3 − 2 u v ( y − x ) 2 = − 1 2 · ( − y − x ± ( y − x ) ) , which provides x = η ζ − ξ ϕ + D 2 ( η 2 − ξ ζ ) , y = η ζ − ξ ϕ − D 2 ( η 2 − ξ ζ ) . Recall that x = 2 2 H 1 and y = 2 2 H 2 , so we get H 1 = 1 2 log 2 η ζ − ξ ϕ + D 2 ( η 2 − ξ ζ ) , H 2 = 1 2 log 2 η ζ − ξ ϕ − D 2 ( η 2 − ξ ζ ) . Also, from equation (22) we obtain v = η − ξ x y − x = 2 η 3 − 2 ξ η ζ − ξ η ζ + ξ 2 ϕ − ξ D − 2 D = ξ D − 2 η 3 + 3 ξ η ζ − ξ 2 ϕ 2 D . Taking into account that σ 2 h 2 H 2 = v and a log b = b log a we get σ 2 = 2 ( η 2 − ξ ζ ) η ζ − ξ ϕ − D log 2 h ξ D − 2 η 3 + 3 ξ η ζ − ξ 2 ϕ 2 D . And similarly, we get κ 2 = 2 ( η 2 − ξ ζ ) η ζ − ξ ϕ + D log 2 h ξ D + 2 η 3 − 3 ξ η ζ + ξ 2 ϕ 2 D .

Equations (31)–(34) are providing us with inverse function q = f − 1 , where f is defined by (5), such that for f ( ϑ ) = τ 0 we get q ( τ 0 ) = ϑ. Therefore, by substituting strongly consistent estimator τ N instead of τ 0 into equations (31)–(34) we obtain strongly consistent estimator ϑ ˆ N = ( H ˆ N ( 1 ) , H ˆ N ( 2 ) , κ ˆ N 2 , σ ˆ N 2 ) ⊤ for ϑ = ( H 1 , H 2 , κ 2 , σ 2 ) ⊤ defined by (9)–(12).  □

 Proof of Theorem <xref rid="j_vmsta234_stat_005">2</xref>

The proof consists of two parts: in the first part we compute the asymptotic covariance matrix Σ ˜, while the second part contains the proof of asymptotic normality.

Part 1: Identification of the asymptotic covariance matrix.

Firstly, we will find the explicit form of covariance matrix Σ ˜ by evaluating convergence limits as N → ∞ for the following variances and covariances: N Var ( ξ N ), N Var ( η N ), N Var ( ζ N ), N Var ( ϕ N ), N cov ( ξ N , η N ), N cov ( ξ N , ζ N ), N cov ( ξ N , ϕ N ), N cov ( η N , ζ N ), N cov ( η N , ϕ N ), N cov ( ζ N , ϕ N ).

1.1 Evaluation of convergence limit for N Var ( ξ N ). Using Isserlis’ theorem,1

Isserlis’ theorem [12]: if ( X 1 , X 2 , X 3 , X 4 ) is a zero-mean multivariate normal random vector, then E ( X 1 X 2 X 3 X 4 ) = E X 1 X 2 E X 3 X 4 + E X 1 X 3 E X 2 X 4 + E X 1 X 4 E X 2 X 3 . In particular, cov ( X 1 2 , X 2 2 ) = 2 cov ( X 1 , X 2 ) 2 .

and stationarity of the sequence Δ X k , we can write cov ( Δ X k ) 2 , ( Δ X j ) 2 = 2 cov ( Δ X k ) , ( Δ X j ) 2 = 2 ρ ˜ ( k − j ) 2 . Therefore N Var ( ξ N ) = 1 N Var ∑ k = 0 N − 1 Δ X k 2 = 1 N ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov ( Δ X k ) 2 , ( Δ X j ) 2 = 2 N ∑ k = 0 N − 1 ∑ j = 0 N − 1 ρ ˜ ( k − j ) 2 . Further, by rearranging sums, we get N Var ( ξ N ) = 2 N ∑ k = 0 N − 1 ∑ i = k − N + 1 k ρ ˜ ( i ) 2 = 2 N ∑ i = − N + 1 0 ∑ k = 0 N − 1 + i ρ ˜ ( i ) 2 + 2 N ∑ i = 1 N − 1 ∑ k = i N − 1 ρ ˜ ( i ) 2 = 2 ∑ i = − N + 1 N − 1 1 − i N ρ ˜ ( i ) 2 → 2 ∑ i = − ∞ ∞ ρ ˜ ( i ) 2 , as N → ∞ . where the passage to the limit can be justified by the dominated convergence theorem due to Lemma 1.

1.2. Evaluation of convergence limit for N cov ( ξ N , η N ). Write N cov ( ξ N , η N ) = 1 N cov ∑ k = 0 N − 1 ( Δ X k ) 2 , ∑ j = 0 N − 1 ( Δ X j + Δ X j + 1 ) 2 = 1 N cov ∑ k = 0 N − 1 ( Δ X k ) 2 , ∑ j = 0 N − 1 ∑ a , b = 0 1 Δ X j + a Δ X j + b = 1 N ∑ a , b = 0 1 ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov ( Δ X k ) 2 , Δ X j + a Δ X j + b . By Isserlis’ theorem, cov ( Δ X k ) 2 , Δ X j + a Δ X j + b = E ( Δ X k ) 2 Δ X j + a Δ X j + b − E ( Δ X k ) 2 E Δ X j + a Δ X j + b = E Δ X k Δ X j + a E Δ X k Δ X j + b + E Δ X k Δ X j + b E Δ X k Δ X j + a = 2 cov ( Δ X k , Δ X j + a ) cov ( Δ X k , Δ X j + b ) = 2 ρ ˜ ( j + a − k ) ρ ˜ ( j + b − k ) . Therefore arguing as in (35), we obtain N cov ( ξ N , η N ) = ∑ a , b = 0 1 2 N ∑ j = 0 N − 1 ∑ k = 0 N − 1 ρ ˜ ( j + a − k ) ρ ˜ ( j + b − k ) = ∑ a , b = 0 1 2 N ∑ j = 0 N − 1 ∑ i = j − N + 1 j ρ ˜ ( i ) ρ ˜ ( i + b − a ) = 2 ∑ a , b = 0 1 ∑ i = − N + 1 N − 1 1 − i N ρ ˜ ( i ) ρ ˜ ( i + b − a ) → 2 ∑ a , b = 0 1 ∑ i = − ∞ ∞ ρ ˜ ( i ) ρ ˜ ( i + b − a ) = ∑ i = − ∞ ∞ ρ ˜ ( i ) ( 4 ρ ˜ ( i ) + 4 ρ ˜ ( i + 1 ) ) , as N → ∞, where the last series converges according to Lemma 3.

1.3. Evaluation of convergence limit for N Var ( η N ). We have N Var ( η N ) = 1 N ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov ∑ a , b = 0 1 Δ X k + a Δ X k + b , ∑ c , d = 0 1 Δ X j + c Δ X j + d = ∑ a , b , c , d = 0 1 1 N ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov ( Δ X k + a Δ X k + b , Δ X j + c Δ X j + d ) .

By Isserlis’ theorem, cov ( Δ X k + a Δ X k + b , Δ X j + c Δ X j + d ) = E Δ X k + a Δ X k + b Δ X j + c Δ X j + d − E Δ X k + a Δ X k + b E Δ X j + c Δ X j + d = E Δ X k + a Δ X k + b E Δ X j + c Δ X j + d + E Δ X k + a Δ X j + c E Δ X k + b Δ X j + d + E Δ X k + a Δ X j + d E Δ X j + c Δ X k + b − E Δ X k + a Δ X k + b E Δ X j + c Δ X j + d = E Δ X k + a Δ X j + c E Δ X k + b Δ X j + d + E Δ X k + a Δ X j + d E Δ X j + c Δ X k + b = cov ( Δ k + a , Δ X j + c ) cov ( Δ X k + b , Δ X j + d ) + cov ( Δ X k + a , Δ X j + d ) cov ( Δ X k + b , Δ X j + c ) = ρ ˜ ( k − j + c − a ) ρ ˜ ( k − j + d − b ) + ρ ˜ ( k − j + c − b ) ρ ˜ ( k − j + d − a ) = ∑ γ = 0 1 ρ ˜ ( k − j + c − a + γ ( a − b ) ) ρ ˜ ( k − j + d − b + γ ( b − a ) ) .

Therefore N Var ( η N ) = ∑ a , b , c , d , γ = 0 1 1 N ∑ k = 0 N − 1 ∑ j = 0 N − 1 ρ ˜ ( k − j + c − a + γ ( a − b ) ) × ρ ˜ ( k − j + d − b + γ ( b − a ) ) = ∑ a , b , c , d , γ = 0 1 ∑ i = − ( N − 1 ) N − 1 1 − | i | N ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) → ∑ a , b , c , d , γ = 0 1 ∑ i = − ∞ ∞ ρ ˜ ( i + c − a + γ ( a − b ) ) ρ ˜ ( i + d − b + γ ( b − a ) ) = ∑ i = − ∞ ∞ ρ ˜ ( i ) ∑ a , b , c , d , γ = 0 1 ρ ˜ ( i + d − c + a − b + 2 γ ( b − a ) ) , as N → ∞. After simplifications, we arrive at lim N → ∞ N Var ( η N ) = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 12 ρ ˜ ( i ) + 16 ρ ˜ ( i + 1 ) + 4 ρ ˜ ( i + 2 ) ) .

1.4. Evaluation of convergence limit for N cov ( ξ N , ζ N ).

We rewrite N cov ( ξ N , ζ N ) in the form N cov ( ξ N , ζ N ) = 1 N cov ∑ k = 0 N − 1 ( Δ X k ) 2 , ∑ j = 0 N − 1 ∑ a , b = 0 3 Δ X j + a Δ X j + b = 1 N ∑ a , b = 0 3 ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov ( Δ X k ) 2 , Δ X j + a Δ X j + b . By Isserlis’ theorem, N cov ( ξ N , ζ N ) = ∑ a , b = 0 3 2 N ∑ j = 0 N − 1 ∑ k = 0 N − 1 ρ ˜ ( j + a − k ) ρ ˜ ( j + b − k ) = ∑ a , b = 0 3 2 N ∑ j = 0 N − 1 ∑ i = j − N + 1 j ρ ˜ ( i ) ρ ˜ ( i + b − a ) = 2 ∑ a , b = 0 3 ∑ i = − N + 1 N − 1 1 − i N ρ ˜ ( i ) ρ ˜ ( i + b − a ) → 2 ∑ a , b = 0 3 ∑ i = − ∞ ∞ ρ ˜ ( i ) ρ ˜ ( i + b − a ) = ∑ i = − ∞ ∞ ρ ˜ ( i ) 8 ρ ˜ ( i ) + 12 ρ ˜ ( i + 1 ) + 8 ρ ˜ ( i + 2 ) + 4 ρ ˜ ( i + 3 ) , as N → ∞, where the last series converges according to Lemma 3.

1.5. Evaluation of convergence limit for N cov ( η N , ζ N ). We have: N cov ( η N , ζ N ) = 1 N cov ∑ k = 0 N − 1 ∑ a , b = 0 1 Δ X j + a Δ X j + b , ∑ j = 0 N − 1 ∑ c , d = 0 3 Δ X j + c Δ X j + d = 1 N ∑ a , b = 0 1 ∑ c , d = 0 3 ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov Δ X j + a Δ X j + b , Δ X j + c Δ X j + d . By Isserlis’ theorem, N cov ( η N , ζ N ) = ∑ a , b , γ = 0 1 ∑ c , d , = 0 3 1 N ∑ k = 0 N − 1 ∑ j = 0 N − 1 ρ ˜ ( k − j + c − a + γ ( a − b ) ) × ρ ˜ ( k − j + d − b + γ ( b − a ) ) = ∑ a , b , γ = 0 1 ∑ c , d , = 0 3 ∑ i = − ( N − 1 ) N − 1 1 − | i | N ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) → ∑ a , b , γ = 0 1 ∑ c , d , = 0 3 ∑ i = − ∞ ∞ ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) = ∑ i = − ∞ ∞ ρ ˜ ( i ) ∑ a , b , γ = 0 1 ∑ c , d , = 0 3 ρ ˜ ( i + d − c + a − b + 2 γ ( b − a ) ) , as N → ∞. After simplifications, we arrive at lim N → ∞ N cov ( η N , ζ N ) = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 28 ρ ˜ ( i ) + 48 ρ ˜ ( i + 1 ) + 32 ρ ˜ ( i + 2 ) + 16 ρ ˜ ( i + 3 ) + 4 ρ ˜ ( i + 4 ) ) .

1.6. Evaluation of convergence limit for N Var ( ζ N ). We have: N Var ( ζ N ) ) = 1 N cov ∑ k = 0 N − 1 ∑ a , b = 0 3 Δ X j + a Δ X j + b , ∑ j = 0 N − 1 ∑ c , d = 0 3 Δ X j + c Δ X j + d = 1 N ∑ a , b , c , d = 0 3 ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov Δ X j + a Δ X j + b , Δ X j + c Δ X j + d . By Isserlis’ theorem, N Var ( ζ N ) = ∑ γ = 0 1 ∑ a , b , c , d = 0 3 1 N ∑ k = 0 N − 1 ∑ j = 0 N − 1 ρ ˜ ( k − j + c − a + γ ( a − b ) ) × ρ ˜ ( k − j + d − b + γ ( b − a ) ) = ∑ γ = 0 1 ∑ a , b , c , d = 0 3 ∑ i = − ( N − 1 ) N − 1 1 − | i | N ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) → ∑ γ = 0 1 ∑ a , b , c , d = 0 3 ∑ i = − ∞ ∞ ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) = ∑ i = − ∞ ∞ ρ ˜ ( i ) ∑ γ = 0 1 ∑ a , b , c , d = 0 3 ρ ˜ ( i + d − c + a − b + 2 γ ( b − a ) ) , as N → ∞. After simplifications, we arrive at lim N → ∞ N Var ( ζ N ) ) = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 88 ρ ˜ ( i ) + 160 ρ ˜ ( i + 1 ) + 124 ρ ˜ ( i + 2 ) + 80 ρ ˜ ( i + 3 ) + 40 ρ ˜ ( i + 4 ) + 16 ρ ˜ ( i + 5 ) + 4 ρ ˜ ( i + 6 ) ) .

1.7. Evaluation of convergence limit for N cov ( ξ N , ϕ N ). We have: N cov ( ξ N , ϕ N ) = 1 N cov ∑ k = 0 N − 1 ( Δ X k ) 2 , ∑ j = 0 N − 1 ∑ a , b = 0 7 Δ X j + a Δ X j + b = 1 N ∑ a , b = 0 7 ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov ( Δ X k ) 2 , Δ X j + a Δ X j + b . By Isserlis’ theorem, N cov ( ξ N , ϕ N ) = ∑ a , b = 0 7 2 N ∑ j = 0 N − 1 ∑ k = 0 N − 1 ρ ˜ ( j + a − k ) ρ ˜ ( j + b − k ) = ∑ a , b = 0 7 2 N ∑ j = 0 N − 1 ∑ i = j − N + 1 j ρ ˜ ( i ) ρ ˜ ( i + b − a ) = 2 ∑ a , b = 0 7 ∑ i = − N + 1 N − 1 1 − i N ρ ˜ ( i ) ρ ˜ ( i + b − a ) → 2 ∑ a , b = 0 7 ∑ i = − ∞ ∞ ρ ˜ ( i ) ρ ˜ ( i + b − a ) , as N → ∞. After simplifications, we arrive at lim N → ∞ N cov ( ξ N , ϕ N ) = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 16 ρ ˜ ( i ) + 28 ρ ˜ ( i + 1 ) + 24 ρ ˜ ( i + 2 ) + 20 ρ ˜ ( i + 3 ) + 16 ρ ˜ ( i + 4 ) + 12 ρ ˜ ( i + 5 ) + 8 ρ ˜ ( i + 6 ) + 4 ρ ˜ ( i + 7 ) ) .

1.8. Evaluation of convergence limit for N cov ( η N , ϕ N ). We have: N cov ( η N , ϕ N ) = 1 N cov ∑ k = 0 N − 1 ∑ a , b = 0 1 Δ X j + a Δ X j + b , ∑ j = 0 N − 1 ∑ c , d = 0 7 Δ X j + c Δ X j + d = 1 N ∑ a , b = 0 1 ∑ c , d = 0 7 ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov Δ X j + a Δ X j + b , Δ X j + c Δ X j + d . By Isserlis’ theorem, N cov ( η N , ϕ N ) = ∑ a , b , γ = 0 1 ∑ c , d = 0 7 1 N ∑ k = 0 N − 1 ∑ j = 0 N − 1 ρ ˜ ( k − j + c − a + γ ( a − b ) ) × ρ ˜ ( k − j + d − b + γ ( b − a ) ) = ∑ a , b , γ = 0 1 ∑ c , d = 0 7 ∑ i = − ( N − 1 ) N − 1 1 − | i | N ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) → ∑ a , b , γ = 0 1 ∑ c , d = 0 7 ∑ i = − ∞ ∞ ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) = ∑ i = − ∞ ∞ ρ ˜ ( i ) ∑ a , b , γ = 0 1 ∑ c , d = 0 7 ρ ˜ ( i + d − c + a − b + 2 γ ( b − a ) ) , as N → ∞. After simplifications, we arrive at lim N → ∞ N cov ( η N , ϕ N ) = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 60 ρ ˜ ( i ) + 112 ρ ˜ ( i + 1 ) + 96 ρ ˜ ( i + 2 ) + 80 ρ ˜ ( i + 3 ) + 64 ρ ˜ ( i + 4 ) + 48 ρ ˜ ( i + 5 ) + 32 ρ ˜ ( i + 6 ) + 16 ρ ˜ ( i + 7 ) + 4 ρ ˜ ( i + 8 ) ) .

1.9. Evaluation of convergence limit for N cov ( ζ N , ϕ N ). We have: N cov ( ζ N , ϕ N ) = 1 N cov ∑ k = 0 N − 1 ∑ a , b = 0 3 Δ X j + a Δ X j + b , ∑ j = 0 N − 1 ∑ c , d = 0 7 Δ X j + c Δ X j + d = 1 N ∑ a , b = 0 3 ∑ c , d = 0 7 ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov Δ X j + a Δ X j + b , Δ X j + c Δ X j + d . By Isserlis’ theorem, N cov ( ζ N , ϕ N ) = ∑ γ = 0 1 ∑ a , b = 0 3 ∑ c , d = 0 7 1 N ∑ k = 0 N − 1 ∑ j = 0 N − 1 ρ ˜ ( k − j + c − a + γ ( a − b ) ) × ρ ˜ ( k − j + d − b + γ ( b − a ) ) = ∑ γ = 0 1 ∑ a , b = 0 3 ∑ c , d = 0 7 ∑ i = − ( N − 1 ) N − 1 1 − | i | N ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) → ∑ γ = 0 1 ∑ a , b = 0 3 ∑ c , d = 0 7 ∑ i = − ∞ ∞ ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) = ∑ i = − ∞ ∞ ρ ˜ ( i ) ∑ γ = 0 1 ∑ a , b = 0 3 ∑ c , d = 0 7 ρ ˜ ( i + d − c + a − b + 2 γ ( b − a ) ) , as N → ∞. After simplifications, we arrive at lim N → ∞ N cov ( ζ N , ϕ N ) = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 216 ρ ˜ ( i ) + 416 ρ ˜ ( i + 1 ) + 376 ρ ˜ ( i + 2 ) + 320 ρ ˜ ( i + 3 ) + 256 ρ ˜ ( i + 4 ) + 192 ρ ˜ ( i + 5 ) + 132 ρ ˜ ( i + 6 ) + 80 ρ ˜ ( i + 7 ) + 40 ρ ˜ ( i + 8 ) + 16 ρ ˜ ( i + 9 ) + 4 ρ ˜ ( i + 10 ) ) .

1.10. Evaluation of convergence limit for N Var ( ϕ N ). We have: N Var ( ϕ N ) = 1 N cov ∑ k = 0 N − 1 ∑ a , b = 0 7 Δ X j + a Δ X j + b , ∑ j = 0 N − 1 ∑ c , d = 0 7 Δ X j + c Δ X j + d = 1 N ∑ a , b , c , d = 0 7 ∑ k = 0 N − 1 ∑ j = 0 N − 1 cov Δ X j + a Δ X j + b , Δ X j + c Δ X j + d . By Isserlis’ theorem, N Var ( ϕ N ) = ∑ γ = 0 1 ∑ a , b , c , d = 0 7 1 N ∑ k = 0 N − 1 ∑ j = 0 N − 1 ρ ˜ ( k − j + c − a + γ ( a − b ) ) × ρ ˜ ( k − j + d − b + γ ( b − a ) ) = ∑ γ = 0 1 ∑ a , b , c , d = 0 7 ∑ i = − ( N − 1 ) N − 1 1 − | i | N ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) → ∑ γ = 0 1 ∑ a , b , c , d = 0 7 ∑ i = − ∞ ∞ ρ ˜ ( i + c − a + γ ( a − b ) ) × ρ ˜ ( i + d − b + γ ( b − a ) ) = ∑ i = − ∞ ∞ ρ ˜ ( i ) ∑ γ = 0 1 ∑ a , b , c , d = 0 7 ρ ˜ ( i + d − c + a − b + 2 γ ( b − a ) ) , as N → ∞. After simplifications, we arrive at lim N → ∞ N Var ( ϕ N ) = ∑ i = − ∞ + ∞ ρ ˜ ( i ) ( 688 ρ ˜ ( i ) + 1344 ρ ˜ ( i + 1 ) + 1260 ρ ˜ ( i + 2 ) + 1136 ρ ˜ ( i + 3 ) + 984 ρ ˜ ( i + 4 ) + 816 ρ ˜ ( i + 5 ) + 644 ρ ˜ ( i + 6 ) + 480 ρ ˜ ( i + 7 ) + 336 ρ ˜ ( i + 8 ) + 224 ρ ˜ ( i + 9 ) + 140 ρ ˜ ( i + 10 ) + 80 ρ ˜ ( i + 11 ) + 40 ρ ˜ ( i + 12 ) + 16 ρ ˜ ( i + 13 ) + 4 ρ ˜ ( i + 14 ) ) .

Part 2: Proof of asymptotic normality.

Let us define Y k = ( Y k ( 1 ) , Y k ( 2 ) , Y k ( 3 ) , Y k ( 4 ) ) ⊤ by Y k ( 1 ) : = Δ X k , Y k ( 2 ) : = Δ X k + 1 , Y k ( 3 ) : = Δ X k + 2 + Δ X k + 3 Y k ( 4 ) : = Δ X k + 4 + Δ X k + 5 + Δ X k + 6 + Δ X k + 7 . Then ξ N = 1 N ∑ k = 0 N − 1 Y k ( 1 ) 2 , η N = 1 N ∑ k = 0 N − 1 Y k ( 1 ) + Y k ( 2 ) 2 , ζ N = 1 N ∑ k = 0 N − 1 Y k ( 1 ) + Y k ( 2 ) + Y k ( 3 ) 2 , ϕ N = 1 N ∑ k = 0 N − 1 Y k ( 1 ) + Y k ( 2 ) + Y k ( 3 ) + Y k ( 4 ) 2 .

We shall prove the convergence of vector ( ξ N , η N , ζ N , ϕ N ) ⊤ with the help of the Cramér–Wold device. Let the parameters α, β, γ, λ ∈ R be such that α 2 + β 2 + γ 2 + λ 2 ≠ 0. We introduce the function f ( y ) = α y 1 2 + β ( y 1 + y 2 ) 2 + γ ( y 1 + y 2 + y 3 ) 2 + λ ( y 1 + y 2 + y 3 + y 4 ) 2 , y = ( y 1 , y 2 , y 3 ) ⊤ ∈ R 3 , so that α ξ N + β η N + γ ζ N + λ ϕ N = 1 N ∑ k = 0 N − 1 f ( Y k ) . Thus, we need to prove that the sequence N ( α ( ξ N − E ξ 0 ) + β ( η N − E η 0 ) + γ ( ζ N − E ζ 0 ) + λ ( ϕ N − E ϕ 0 ) ) = 1 N ∑ k = 0 N − 1 ( f ( Y k ) − E f ( Y k ) ) converges to a normal distribution. This fact can be established by application of the Breuer–Major theorem for stationary vectors [1, Theorem 4]. In order to apply this theorem, it suffices to verify the following condition: ∑ j ∈ Z r ( p , l ) ( j ) q < ∞ , for all p , l ∈ { 1 , 2 , 3 } , where r ( p , l ) ( k ) : = cov Y 1 ( p ) , Y 1 + k ( l ) , k ∈ Z , and q is the Hermite rank2

The function f : R d → R is said to have Hermite rank equal to q with respect to a random vector X if (a)  E [ ( f ( X ) − E [ f ( X ) ] ) p m ( X ) ] = 0 for every polynomial p m (on R d ) of degree m ≤ q − 1; and (b) there exists a polynomial p q of degree q such that E [ ( f ( X ) − E [ f ( X ) ] ) p m ( X ) ] ≠ 0, see [23].

of f with respect to Y 1 .

It is not hard to see that this Hermite rank q ≥ 2. Indeed, note that f ( Y 1 ) is a second-order polynomial of zero-mean normally distributed random variables Y 1 ( 1 ) , Y 1 ( 2 ) , Y 1 ( 3 ) , Y 1 ( 4 ) in which only terms of the second order are present. Therefore, according to Isserlis’ theorem, we will have that the expected value of the product of f ( Y 1 ) and any of Y 1 ( t ) , t ∈ { 1 , 2 , 3 }, will be equal to zero. All the obtained terms will have the form E G 1 G 2 G 3 , where ( G 1 , G 2 , G 3 ) ⊤ is a zero-mean multivariate normal random vector and according to Isserlis’ theorem E G 1 G 2 G 3 = 0. Therefore, the expected value of the product of f ( Y 1 ) and any first-order polynomial of { Y 1 } (a linear combination of Y 1 ( 1 ) , Y 1 ( 2 ) , Y 1 ( 3 ) , Y 1 ( 4 ) ) will be equal to zero, and therefore the Hermite rank q of the function f with respect to Y 1 cannot be equal to 1.

Since q ≥ 2, we easily see that in order to verify the condition (38), it suffices to prove that ∑ j ∈ Z r ( p , l ) ( j ) 2 < ∞ , for all p , l ∈ { 1 , 2 , 3 } .

In turn, (39) follows from Lemma 1, since using (36) and the definition of ρ ˜, we can represent r ( p , l ) via ρ ˜ as follows: r ( 1 , 1 ) ( k ) = r ( 2 , 2 ) ( k ) = ρ ˜ ( k ) , r ( 1 , 2 ) ( k ) = ρ ˜ ( k + 1 ) , r ( 2 , 1 ) ( k ) = ρ ˜ ( k − 1 ) , r ( 1 , 3 ) ( k ) = ρ ˜ ( k ) + ρ ˜ ( k + 2 ) + ρ ˜ ( k + 3 ) , r ( 3 , 1 ) ( k ) = ρ ˜ ( k − 2 ) + ρ ˜ ( k − 3 ) , r ( 2 , 3 ) ( k ) = + ρ ˜ ( k + 1 ) + ρ ˜ ( k + 2 ) , r ( 3 , 2 ) ( k ) = ρ ˜ ( k − 1 ) + ρ ˜ ( k − 2 ) , r ( 3 , 3 ) ( k ) = ρ ˜ ( k + 1 ) + 2 ρ ˜ ( k ) + ρ ˜ ( k − 1 ) , r ( 1 , 4 ) ( k ) = ρ ˜ ( k + 4 ) + ρ ˜ ( k + 5 ) + ρ ˜ ( k + 6 ) + ρ ˜ ( k + 7 ) , r ( 4 , 1 ) ( k ) = ρ ˜ ( k − 4 ) + ρ ˜ ( k − 5 ) + ρ ˜ ( k − 6 ) + ρ ˜ ( k − 7 ) , r ( 2 , 4 ) ( k ) = ρ ˜ ( k + 3 ) + ρ ˜ ( k + 4 ) + ρ ˜ ( k + 5 ) + ρ ˜ ( k + 6 ) , r ( 4 , 2 ) ( k ) = ρ ˜ ( k − 3 ) + ρ ˜ ( k − 4 ) + ρ ˜ ( k − 5 ) + ρ ˜ ( k − 6 ) , r ( 3 , 4 ) ( k ) = ρ ˜ ( k + 2 ) + 2 ρ ˜ ( k + 2 ) + 2 ρ ˜ ( k + 3 ) + 2 ρ ˜ ( k + 4 ) + ρ ˜ ( k + 5 ) , r ( 4 , 3 ) ( k ) = ρ ˜ ( k − 2 ) + 2 ρ ˜ ( k − 2 ) + 2 ρ ˜ ( k − 3 ) + 2 ρ ˜ ( k − 4 ) + ρ ˜ ( k − 5 ) , r ( 4 , 4 ) ( k ) = ρ ˜ ( k − 3 ) + 2 ρ ˜ ( k − 2 ) + 3 ρ ˜ ( k − 1 ) + 4 ρ ˜ ( k ) + 3 ρ ˜ ( k + 1 ) + 2 ρ ˜ ( k + 2 ) + ρ ˜ ( k + 3 ) . Indeed, for example, the series ∑ k = − ∞ + ∞ r ( 1 , 3 ) ( k ) 2 = ∑ k = − ∞ + ∞ ( ρ ˜ ( k + 2 ) + ρ ˜ ( k + 3 ) ) 2 = ∑ k = − ∞ + ∞ ρ ˜ ( k + 2 ) 2 + 2 ∑ k = − ∞ + ∞ ρ ˜ ( k + 2 ) ρ ˜ ( k + 3 ) + ∑ k = − ∞ + ∞ ρ ˜ ( k + 3 ) 2 , converges as a finite sum of series convergent according to Lemma 1. Similarly, other series in (39) are also finite.

Therefore, the assumptions of the Breuer–Major theorem for stationary vectors are satisfied, whence the desired weak converge of (37) to a zero-mean normal distribution follows.  □

 Proof of Theorem <xref rid="j_vmsta234_stat_007">3</xref>

In order to prove Theorem 3 we apply the delta method (see, e.g., [13, Theorem B.6]) for function g constructed in Theorem 1 and sequence τ N that is asymptotically normal by Theorem 1, with the inverse function theorem. In order to apply the delta method we need to calculate the derivation matrix g ′ and prove that it is nonsingular. But considering the complicated structure of this function, we will calculate the derivation matrix f ′ for the function f defined in (5) and apply the inverse function theorem at the point ϑ since g is an inverse function of f, and thus g ′ = ( f ′ ) − 1 . Therefore, we will need to find matrix f ′ ( ϑ ) and to show that it is nonsingular.

We will start with evaluating partial derivatives of function f 1 . Since f 1 ( H 1 , H 2 , κ 2 , σ 2 ) = κ 2 h 2 H 1 + σ 2 h 2 H 2 , we see that ∂ f 1 ∂ H 1 = 2 log h κ 2 h 2 H 1 , ∂ f 1 ∂ H 2 = 2 log h σ 2 h 2 H 2 , ∂ f 1 ∂ κ 2 = h 2 H 1 , ∂ f 1 ∂ σ 2 = h 2 H 2 .

For the function f 2 ( H 1 , H 2 , κ 2 , σ 2 ) = κ 2 h 2 H 1 2 2 H 1 + σ 2 h 2 H 2 2 2 H 2 = κ 2 e 2 H 1 log ( 2 h ) + σ 2 e 2 H 2 log ( 2 h ) , we get ∂ f 2 ∂ H 1 = 2 log ( 2 h ) κ 2 h 2 H 1 2 2 H 1 , ∂ f 2 ∂ H 2 = 2 log ( 2 h ) σ 2 h 2 H 2 2 2 H 2 , ∂ f 2 ∂ κ 2 = h 2 H 1 2 2 H 1 , ∂ f 2 ∂ σ 2 = h 2 H 2 2 2 H 2 .

Similarly, one can evaluate partial derivatives for the functions f 3 and f 4 . We arrive at the following derivative matrix: f ′ ( ϑ ) = 2 log ( h ) κ 2 h 2 H 1 2 log ( h ) σ 2 h 2 H 2 h 2 H 1 h 2 H 2 2 log ( 2 h ) κ 2 ( 2 h ) 2 H 1 2 log ( 2 h ) σ 2 ( 2 h ) 2 H 2 ( 2 h ) 2 H 1 ( 2 h ) 2 H 2 2 log ( 4 h ) κ 2 ( 4 h ) 2 H 1 2 log ( 4 h ) σ 2 ( 4 h ) 2 H 2 ( 4 h ) 2 H 1 ( 4 h ) 2 H 2 2 log ( 8 h ) κ 2 ( 8 h ) 2 H 1 2 log ( 8 h ) σ 2 ( 8 h ) 2 H 2 ( 8 h ) 2 H 1 ( 8 h ) 2 H 2 . Now, we need to compute its determinant, so that we could apply the delta method. We will have det ( f ′ ( ϑ ) ) = 4 κ 2 σ 2 h 4 H 1 + 4 H 2 log ( h ) log ( h ) 1 1 log ( 2 h ) 2 2 H 1 log ( 2 h ) 2 2 H 2 2 2 H 1 2 2 H 2 log ( 4 h ) 4 2 H 1 log ( 4 h ) 4 2 H 2 4 2 H 1 4 2 H 2 log ( 8 h ) 8 2 H 1 log ( 8 h ) 8 2 H 2 8 2 H 1 8 2 H 2 = 4 κ 2 σ 2 h 4 H 1 + 4 H 2 log 2 0 0 1 1 2 2 H 1 2 2 H 2 2 2 H 1 2 2 H 2 2 · 2 4 H 1 2 · 2 4 H 2 2 4 H 1 2 4 H 2 3 · 2 6 H 1 3 · 2 6 H 2 2 6 H 1 2 6 H 2 = 4 κ 2 σ 2 h 4 H 1 + 4 H 2 log 2 2 2 H 1 2 2 H 2 2 2 H 2 2 · 2 4 H 1 2 · 2 4 H 2 2 4 H 2 3 · 2 6 H 1 3 · 2 6 H 2 2 6 H 2 − 2 2 H 1 2 2 H 2 2 2 H 1 2 · 2 4 H 1 2 · 2 4 H 2 2 4 H 1 3 · 2 6 H 1 3 · 2 6 H 2 2 6 H 1 . Considering that 2 2 H 1 2 2 H 2 2 2 H 2 2 · 2 4 H 1 2 · 2 4 H 2 2 4 H 2 3 · 2 6 H 1 3 · 2 6 H 2 2 6 H 2 = − 2 2 H 1 2 10 H 2 + 4 · 2 4 H 1 2 8 H 2 − 3 · 2 6 H 1 2 6 H 2 , and 2 2 H 1 2 2 H 2 2 2 H 1 2 · 2 4 H 1 2 · 2 4 H 2 2 4 H 1 3 · 2 6 H 1 3 · 2 6 H 2 2 6 H 1 = 2 10 H 1 2 2 H 2 + 3 · 2 6 H 1 2 6 H 2 − 4 · 2 8 H 1 2 4 H 2 , we get det ( f ′ ( ϑ ) ) = − 4 κ 2 σ 2 h 6 H 1 + 6 H 2 log ( 2 ) 2 2 H 2 − 2 2 H 1 4 ≠ 0 . Therefore, the derivative matrix is nonsingular, and the corresponding derivative matrix of the inverse function g ′ ( τ 0 ) = ( f ′ ( ϑ ) ) − 1 is nonsingular by the inverse function theorem. Thus, the delta method can be applied, which implies the statement of the theorem and provides the formula for the asymptotic covariance matrix Σ 0 defined using the inverse derivative matrix f ′ ( ϑ ).  □

 References Arcones, M.A.: Limit theorems for nonlinear functionals of a stationary Gaussian sequence of vectors. Ann. Probab. 22(4), 2242–2274 (1994). MR1331224 Cai, C., Chigansky, P., Kleptsyna, M.: Mixed Gaussian processes: a filtering approach. Ann. Probab. 44(4), 3032–3075 (2016). MR3531685. https://doi.org/10.1214/15-AOP1041 Cheridito, P.: Mixed fractional Brownian motion. Bernoulli 7(6), 913–934 (2001). MR1873835. https://doi.org/10.2307/3318626 Dai, Q., Singleton, K.J.: Specification analysis of affine term structure models. J. Finance 55, 1943–1978 (2000) Ding, Z., Granger, C.W.J., Engle, R.F.: A long memory property of stock market returns and a new model. J. Empir. Finance 1(1), 83–106 (1993) 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-Nouty, C.: The fractional mixed fractional Brownian motion. Stat. Probab. Lett. 65(2), 111–120 (2003). MR2017255 Filatova, D.: Mixed fractional Brownian motion: Some related questions for computer network traffic modeling. In: 2008 International Conference on Signals and Electronic Systems, pp. 393–396 (2008) He, X., Chen, W.: The pricing of credit default swaps under a generalized mixed fractional Brownian motion. Physica A 404, 26–33 (2014). MR3188756 Henry, O.T.: Long memory in stock returns: Some international evidence. Appl. Financ. Econ. 12(10), 725–729 (2002) Isserlis, L.: On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables. Biometrika 12(1/2), 134–139 (1918) Kubilius, K., Mishura, Y., Ralchenko, K.: Parameter Estimation in Fractional Diffusion Models. Bocconi & Springer Series, vol. 8, p. 390. Bocconi University Press & Springer (2017). MR3752152. https://doi.org/10.1007/978-3-319-71030-3 Kukush, A., Lohvinenko, S., Mishura, Y., Ralchenko, K.: Two approaches to consistent estimation of parameters of mixed fractional Brownian motion with trend. Stat. Inference Stoch. Process. 25(1), 159–187 (2022). MR4419677. https://doi.org/10.1007/s11203-021-09252-6 Major, P.: Non-central limit theorem for non-linear functionals of vector valued Gaussian stationary random fields. arXiv:1901.04086 [math.PR] (2019) Miao, Y., Ren, W., Ren, Z.: On the fractional mixed fractional Brownian motion. Appl. Math. Sci. (Ruse) 2(33-36), 1729–1738 (2008). MR2443901 Mishura, Y., Ralchenko, K., Shklyar, S.: Maximum likelihood drift estimation for Gaussian process with stationary increments. Austrian J. Stat. 46(3–4), 67–78 (2017) Mishura, Y.: Maximum likelihood drift estimation for the mixing of two fractional brownian motions. In: Trends in Mathematics, pp. 263–280. Springer (2016). MR3708386. https://doi.org/10.1007/978-3-319-07245-6_14 Mishura, Y., Voronov, I.: Construction of maximum likelihood estimator in the mixed fractional–fractional Brownian motion model with double long-range dependence. Mod. Stoch. Theory Appl. 2(2), 147–164 (2015). MR3389587. https://doi.org/10.15559/15-vmsta28 Mishura, Y., Ralchenko, K., Shklyar, S.: Parameter estimation for Gaussian processes with application to the model with two independent fractional Brownian motions. In: Stochastic Processes and Applications. Springer Proc. Math. Stat., vol. 271, pp. 123–146. Springer (2018). MR3903341 Mishura, Y., Ralchenko, K., Zhelezniak, H.: Numerical approach to the drift parameter estimation in the model with two fractional brownian motions. Commun. Stat., Simul. Comput. 1–15 (2022) Nourdin, I.: Selected Aspects of Fractional Brownian Motion. Springer (2012). MR3076266. https://doi.org/10.1007/978-88-470-2823-4 Nourdin, I., Peccati, G., Podolskij, M.: Quantitative Breuer-Major theorems. Stoch. Process. Appl. 121(4), 793–812 (2011). MR2770907 Paxson, V., Floyd, S.: Wide area traffic: the failure of Poisson modeling. IEEE/ACM Trans. Netw. 3(3), 226–244 (1995) Rostek, S., Schöbel, R.: A note on the use of fractional Brownian motion for financial modeling. Econ. Model. 30, 30–35 (2013) Sadique, S., Silvapulle, P.: Long-term memory in stock market returns: International evidence. Int. J. Financ. Econ. 6(1), 59–67 (2001) Sun, L.: Pricing currency options in the mixed fractional Brownian motion. Physica A 392(16), 3441–3458 (2013). MR3069168 Thäle, C.: Further remarks on mixed fractional brownian motion. Appl. Math. Sci. 38, 1885–1901 (2009) 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 Xiao, W.-L., Zhang, W.-G., Zhang, X., Zhang, X.: Pricing model for equity warrants in a mixed fractional Brownian environment and its algorithm. Physica A 391(24), 6418–6431 (2012). MR2972795 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) Zhang, W.-G., Xiao, W.-L., He, C.-X.: Equity warrants pricing model under fractional Brownian motion and an empirical study. Expert Syst. Appl. 36(2), 3056–3065 (2009) Zili, M.: On the mixed fractional Brownian motion. J. Appl. Math. Stoch. Anal. 2006, Art. ID 32435 (2006), 9 pp. MR2253522. https://doi.org/10.1155/JAMSA/2006/32435