Modern Stochastics: Theory and Applications

The paper is devoted to the statistical identification of the following stochastic process: 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}}$, ${\kappa ^{2}}$, ${\sigma ^{2}}$ based on equidistant observations $\{{X_{kh}},k=0,1,2,\dots \hspace{0.1667em}\}$, $h\gt 0$. In order to get an identifiable model, we assume that $0\lt {H_{1}}\lt {H_{2}}\lt 1$, $\sigma \gt 0$, and $\kappa \gt 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}}\in (0,1)$ is less studied. We can mention only the papers [18, 19], where the following mixed model with a trend was considered: 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 $\kappa =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}}\gt \frac{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 $\stackrel{\text{d}}{\to }$.

Let $h\gt 0$ be fixed. Assume that the process ${X_{t}}$ defined in (1) is observed at points ${t_{k}}=kh$, $k=0,1,2,\dots \hspace{0.1667em}$. Let us denote the increment The estimation procedure will be based on the following result.

This lemma allows us to apply the ergodic theorem for the construction of estimators. Namely, if $g:{\mathbb{R}^{l+1}}\to \mathbb{R}$ is a Borel function such that $\operatorname{\mathsf{E}}|g(\Delta {X_{0}},\Delta {X_{1}},\dots ,\Delta {X_{l}})|\lt \infty $, then a. s., as $N\to \infty $. 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 Denote also $\vartheta ={({H_{1}},{H_{2}},{\kappa ^{2}},{\sigma ^{2}})^{\top }}$.

Let us introduce the notation Then $f(\vartheta )={\tau _{0}}$. Therefore, in order to obtain estimators for parameters ${({H_{1}},{H_{2}},{\kappa ^{2}},{\sigma ^{2}})^{\top }}$ we need to construct function q such that $q({\tau _{0}})={f^{-1}}({\tau _{0}})$, and then required estimators would be defined as ${({\widehat{H}_{N}^{(1)}},{\widehat{H}_{N}^{(2)}},{\hat{\kappa }_{N}^{2}},{\hat{\sigma }_{N}^{2}})^{\top }}:=q({\tau _{N}})$. In other words, in order to construct strongly consistent estimators, we need to solve the estimation equation 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.

We will start with a study of the asymptotic properties of the vector ${\tau _{N}}=({\phi _{N}},{\xi _{N}},{\eta _{N}},{\zeta _{N}})$. The previously defined stationary Gaussian sequence $\Delta {X_{k}}$ has the autocovariance function where denotes the autocovariance function of the stationary sequence $\{{B_{k+1}^{H}}-{B_{k}^{H}},k\ge 0\}$, which is known as fractional Gaussian noise.

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

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\lt {H_{1}}\lt {H_{2}}\lt 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\lt {H_{1}}\lt {H_{2}}\lt \frac{3}{4}$, we calculate the average approximate theoretical standard deviation (ASD) – the square root of the estimated asymptotic variance divided by N ($\sqrt{{({\hat{\Sigma }_{N}^{0}})_{ii}}/N}$) – and the coverage probability (CP) for $\alpha =5\% $, based on the estimator of the asymptotic covariance matrix. We are computing these statistics only for ${H_{1}}\lt {H_{2}}\lt \frac{3}{4}$, excluding interval $(\frac{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 ${\widehat{H}_{N}}$ will be out of boundaries $(0,\frac{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 where $\alpha ,\beta \in \mathbb{Z}$.

Then, series ${\textstyle\sum _{i=0}^{+\infty }}\tilde{\rho }(i+\alpha )\tilde{\rho }(i+\beta )$ are approximated with precision $\delta ={10^{-3}}$ by evaluating the sum ${S_{m}}={\textstyle\sum _{i=0}^{m}}\tilde{\rho }(i+\alpha )\tilde{\rho }(i+\beta )$, where m is the smallest number such that $|{S_{m}}-{S_{m-1}}|\lt {10^{-3}}$. Obtained approximations will be used to compute the asymptotic covariance matrix defined in Theorem 2.

For ${\widehat{H}_{N}^{(1)}}$, ${\widehat{H}_{N}^{(2)}}$, ${\hat{\kappa }_{N}^{2}}$, ${\hat{\sigma }_{N}^{2}}$, ${\hat{\Sigma }_{N}^{0}}$ and fixed time step $h=1$ we vary the horizon $T=h{2^{n}}$ for $n\in \{8,10,12,14,16,18,20\}$. For all simulations, the values $\sigma =1$ and $\kappa =1$ were used.

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 ${\kappa ^{2}}$ and ${\sigma ^{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 ${\widehat{H}_{N}^{(2)}}-{\widehat{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.

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., ${\phi _{N}}$ is less accurate than ${\zeta _{N}}$, which is less accurate than ${\eta _{N}}$, which is less accurate than ${\xi _{N}}$). Unfortunately, the percentage of iterations with undefined asymptotic covariance matrix ${\hat{\Sigma }_{N}^{0}}$ is high but decreases for bigger samples, which is shown in Table 7, due to estimated ${\widehat{H}_{N}^{(1)}}$, ${\widehat{H}_{N}^{(2)}}$ violating the conditions of Theorem 3. Also, the percentage of iterations with the undefined ${\hat{\Sigma }_{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 We estimate (16) by finding the smallest number $m\in \mathbb{N}$ such that for ${S_{m}}\hspace{-0.1667em}=\hspace{-0.1667em}{\textstyle\sum _{i=0}^{m}}\tilde{\rho }{(i)^{2}}$ the inequality $|{S_{m}}-{S_{m-1}}|\le \delta $ holds for some predefined δ.

This estimation was conducted for $\delta \in \{{10^{-3}},{10^{-4}},{10^{-5}},{10^{-6}},{10^{-7}},{10^{-8}}$, ${10^{-9}}\}$ for $\sigma =1$, $\kappa =1$, ${H_{1}}=0.1$ and ${H_{2}}\in \{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 $\delta ={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.