Hurst index estimation in stochastic differential equations driven by fractional Brownian motion

TL;DR

This paper develops new estimators for the Hurst index in stochastic differential equations driven by fractional Brownian motion, utilizing Malliavin calculus to analyze quadratic variations and their asymptotic behavior.

Contribution

It introduces a novel approach to estimate the Hurst index in SDEs driven by fractional Brownian motion using higher order increments and Malliavin calculus techniques.

Findings

Consistent estimators for the Hurst index are constructed.

Asymptotic properties of the estimators are analyzed.

Simulation results demonstrate estimator effectiveness.

Abstract

We consider the problem of Hurst index estimation for solutions of stochastic differential equations driven by an additive fractional Brownian motion. Using techniques of the Malliavin calculus, we analyze the asymptotic behavior of the quadratic variations of the solution, defined via higher order increments. Then we apply our results to construct and study estimators for the Hurst index.

Equations

$$X_t=x+{\int }_0^{t}f(s,X_s)ds+B_t^H,t\in T$$

$$\mathbb{E}[B_t^H{B}_s^H]=\frac{1}{2}(t^{2H}+s^{2H}-|t-s{|}^{2H})\text{for}s,t\in T.$$

$$\sum _{q=0}^p{a}_q=0.$$

$$\sum _{i=0}^p{a}_i{i}^k=0\text{ for }0\le k<M(a)\text{ and }\sum _{i=0}^p{a}_i{i}^{M(a)}\ne 0.$$

$$b_i=\sum _{k=0}^i{a}_k\text{mbox}fori=0,1,\dots ,p.$$

$${\Delta }_a{X}_j:=\sum _{i=0}^p{a}_i{X}_{(i+j)\Delta }$$

$${\Delta }_a{X}_j$$

$$=a_0(X_{j\Delta }-X_{(j+1)\Delta })+(a_0+a_1)X_{(j+1)\Delta }+a_2{X}_{(j+2)\Delta }+\cdots +a_p{X}_{(j+p)\Delta }=\cdots$$

$$=a_0(X_{j\Delta }-X_{(j+1)\Delta })+\cdots +\sum _{i=0}^{p-1}{a}_i(X_{(j+p-1)\Delta }-X_{(j+p)\Delta })+\underset{=0}{\underbrace{\sum _{i=0}^p{a}_i{X}_{(j+p)\Delta }}}$$

$$=\sum _{i=0}^{p-1}\sum _{k=0}^i{a}_k(X_{(i+j)\Delta }-X_{(i+j+1)\Delta })=\sum _{i=0}^{p-1}{b}_i(X_{(i+j)\Delta }-X_{(i+j+1)\Delta }).$$

$$\mathrm{Cov}({\Delta }_a{B}_i^H,{\Delta }_a{B}_j^H)=-\frac{{\Delta }^{2H}}{2}\sum _{k=0}^p\sum _{l=0}^p{a}_k{a}_l|i+k-j-l{|}^{2H}.$$

$$V(a,n,\Delta ,B^H)=\frac{1}{n}\sum _{j=1}^{n-p}\left(\frac{({\Delta }_a{B}_j^H{)}^2}{{\sigma }_{a,\Delta }}-1\right),$$

$${\sigma }_{a,\Delta }=\mathrm{Var}({\Delta }_a{B}_j^H)=-\frac{{\Delta }^{2H}}{2}\sum _{k=0}^p\sum _{l=0}^p{a}_k{a}_l|k-l{|}^{2H}.$$

$$\sqrt{n}V(a,n,\Delta ,B^H)\stackrel{(d)}{\to }N(0,{\sigma }_H),$$

$$\sqrt{n}(V(a^1,n,\Delta ,B^H),\dots ,V(a^N,n,\Delta ,B^H))\stackrel{(d)}{\to }N(0,{\Sigma }_H),$$

$$X_t=x+{\int }_0^{t}f(s,X_s)ds+B_t^H,$$

$$\Vert f{\Vert }_{\infty }+\Vert f^{\prime }{\Vert }_{\infty }\le M,$$

$$Y_t={\int }_0^{t}f(s,X_s)ds,\text{mbox}foreveryt\in T.$$

$$V(a,n,\Delta ,X)=\frac{1}{n}\sum _{j=1}^{n-p}\left(\frac{({\Delta }_a{X}_j{)}^2}{{\sigma }_{a,\Delta }}-1\right),$$

$$\sqrt{n}V(a,n,\Delta ,X)\stackrel{(d)}{\to }N(0,{\sigma }_H),$$

$$\sqrt{n}V(a,n,\Delta ,X)=\sqrt{n}V(a,n,\Delta ,B^H)+R_n$$

$$R_n=n^{-1/2}{\Delta }^{-2H}\sum _{r=0}^1{C}^r\sum _{i=0}^{n-p}({\Delta }_a{B}_i^H{)}^r({\Delta }_a{Y}_i{)}^{2-r}=:R_{n,0}+R_{n,1}$$

$$|\Delta_{a}Y_{i}|=\Bigl{|}\sum_{j=0}^{p-1}b_{j}\int_{(i+j)\Delta}^{(i+j+1)\Delta}f(s,X_{s})ds\Bigr{|}\lesssim\Delta,$$

$$\mathbb{E}[|R_{n,r}|]\lesssim n^{-1/2}{\Delta }^{-2H}{\Delta }^{2-r}\sum _{i=0}^{n-p}\mathbb{E}[|{\Delta }_a{B}_i^H{|}^r].$$

$$\mathbb{E}[|{\Delta }_a{B}_i^H|]=\frac{2\sqrt{{\sigma }_{a,\Delta }}}{\sqrt{2\pi }}\lesssim {\Delta }^H,$$

$$\mathbb{E}[|R_{n,0}|]\lesssim {\Delta }^{-2H}{n}^{-1/2}{\Delta }^{2}n=n^{1/2-\alpha (2-2H)}$$

$$\mathbb{E}[|R_{n,1}|]\lesssim {\Delta }^{-2H}{n}^{-1/2}\Delta n{\Delta }^H=n^{1/2-\alpha (1-H)},$$

$$\sqrt{n}V(a,n,\Delta ,X)\stackrel{(d)}{\to }N(0,{\sigma }_H),$$

$$D_s{X}_t={\int }_s^t{f}^{\prime }(r,X_r)D_s{X}_{r}dr+1$$

$$D_s{X}_t=e^{{\int }_s^t{f}^{\prime }(r,X_r)dr},$$

$$e^{-tM}\le D_s{X}_t\le e^{tM}$$

$$\mathbb{E}[|{⟨D_{\cdot }(Y_u-Y_v),1_{[s,t]}(\cdot )⟩}_{\mathcal{H}}|]\lesssim (t-s)(u-v),$$

Full text

2010 AMS Classification Numbers: 60G15, 60H05, 60G18.

Key Words and Phrases: Hurst index estimation; Stochastic differential equation; fractional Brownian motion; Quadratic variation; Malliavin calculus; Central limit theorem.

1 Introduction

Fractional Brownian motion (fBm) is a widely used generalization of the standard Brownian motion that incorporates long-range dependence while preserving the self-similarity and Gaussianity. These features make it a popular stochastic processes in the mathematical modelling of various complex systems from financial applications to surface growth (e.g. [2, 1]). The need to detect and analyze such systems has led to the development of many, now classical, statistical instruments, including wavelet analysis, R/S estimators, generalized variations etc. and we refer to the monographs [3, 9, 18] for further reference.

Fractional Brownian motion is homogeneous in space. Many physical systems, however, are subject to an external force and are confined to certain locations with high probability. It may therefore be appropriate to consider stochastic differential equations (SDE) as a generalization of the fBm model. fBm then takes the role of a random forcing which heavily influences the underlying dynamics. Thus the identification of this random forcing is of particular interest.

In our work, we will consider an SDE with additive fractional Brownian noise, i.e.

[TABLE]

where $T$ is an interval in $[0,\infty)$, $x\in\mathbb{R}$ is the initial state and $f$ is a measurable deterministic function satisfying suitable assumptions. Stochastic models as (1) appear in mathematical finance (for example, in [12] the log-volatility is assumed to satisfy an SDE of the form (1)). Other potential applications come from climatology, stochastic processes as in (1) (with Gaussian or non-Gaussian noise) are utilized as models for problems related to the earth’s energy balance, see [10, 8].

In fBm-related models particular interest lies in the estimation of the Hurst index (see Section 2 below) since it constitutes the characteristic parameter of the driving forcing.

The purpose of our endeavour is to estimate the Hurst index $H$ in (1) based on the discrete observations of the process $(X_{t})_{t\in T}$. Here we will deploy a well-known method, based on the quadratic variations of $X$. This method is known to work well for self-similar processes (see e.g. [18] and references therein).

We will show that it can be also applied to (1), although $X$ is not a self-similar process. The method exploits the fact that the absolutely continuous integral component of (1) does not affect the roughness captured by the quadratic variation of $X$.

We will define the sequence $V(a,n,\Delta,B^{H})$ (see equation (7) below) of the so-called quadratic $a$-variations, defined in terms of higher order increments of $X$ over a filter $a$. We give the limit behavior in distribution of this sequence making use of the Malliavin calculus to handle correlations and combine it with already known results concerning the variations of the fBm. Interestingly, in this setting the case $H<\frac{1}{2}$ is easier to handle than the case of $H>\frac{1}{2}$. This is due to the fact that, when studying the limit of the sequence (7), one needs to take into account the correlations between the increments of the fBm, but also the correlations between the increments of the fBm and the increments of the Lebesque integral in (1). If $H<\frac{1}{2}$ these joint correlations are always dominated by those of the fBm, which is not the case for $H>\frac{1}{2}$, when a supplementary assumption is needed on the mesh $\Delta$ in (7).

By a standard procedure, we then construct a quadratic variation estimator for the Hurst index of the model (1). We prove its consistency and its asymptotic normality, using the limit behavior of the quadratic $a$-variations.

Our work is structured as follows. In Section 2 we introduce the objects of study, fBm, filters $a$ and the corresponding quadratic $a$-variations of a process. We also quote the underlying results on the behavior of the $a$-variations of fBm. Section 3 studies the properties and the $a$-variations of the solution to (1), via the techniques of the Malliavin calculus. Section 4 is devoted to the parameter estimation of the Hurst index from discrete observations of $X$. In Section 5 we present simulations and discuss statistical properties of the derived estimators. In the Appendix we give a short review of the definitions and results from Malliavin calculus relevant to this work.

2 Preliminaries: Fractional Brownian motion and its variations

A fractional Brownian motion $(B^{H}_{t})_{t\in T}$ with Hurst index $H\in(0,\,1)$ is a centered Gaussian process on the interval $T$ (in this work we consider $T=[0,\,1]$ or $T=\mathbb{R}^{+}$) with covariance function

[TABLE]

From this definition it is apparent that for $H=\tfrac{1}{2}$ we have that $\mathbb{E}[B^{H}_{t}B^{H}_{s}]=\min(s,t)$ and the fBm with $H=\tfrac{1}{2}$ is the standard Brownian motion.

For a natural number $p\in\mathbb{N}$ let us consider a $(p+1)$-tuple $a=(a_{0},\dots,a_{p})$ of real numbers with zero sum, i.e.

[TABLE]

Such a $(p+1)$-tuple $a=(a_{0},\dots,a_{p})$ will be called a filter of length $p+1$. The order of a filter $a$, denoted by $M(a)$, is defined to be the order of the first non-zero moment of the $(p+1)$-tuple $a$, i.e.

[TABLE]

We will frequently need the partial sum of the components of $a$ and use the notation

[TABLE]

Since by definition a filter has zero sum, clearly for any filter $a$ we have $M(a)\geq 1$. For instance, $a=(a_{0},a_{1})=(-1,1)$ is a filter of order 1 and of length 2, while $a=(a_{0},a_{1},a_{2})=(1,-2,1)$ is a filter of order 2 and of length 3.

For a process $(X_{t})_{t\in T}$, a filter $a$ of length $p+1$ and a fixed mesh size $\Delta>0$ we consider the sum

[TABLE]

with $j\in\mathbb{N}$ such that $(p+j)\Delta\in T$. For example, if $a=(a_{0},a_{1})=(-1,1)$, then $\Delta_{a}X_{j}=X_{(j+1)\Delta}-X_{j\Delta}$ while if $a=(a_{0},a_{1},a_{2})=(1,-2,1)$ then $\Delta_{a}X_{j}=\Delta_{j\Delta}-2\Delta_{(j+1)\Delta}+X_{(j+2)\Delta}$.

We will refer to $\Delta_{a}X_{j}$ as the increment of the process $X$ at time $j$ over the filter $a$. Note that due to the zero sum condition on $a$ we can rewrite $\Delta_{a}X_{j}$ in the following way:

[TABLE]

We will refer to this form as differences representation. The correlation of increments of the fBm over a filter plays an important role in our calculations. From (2) we can see that for a zero sum vector $a\in\mathbb{R}^{(p+1)}$ and $i,\,j\in\mathbb{N}$ such that $\Delta_{a}B^{H}_{i}$, $\Delta_{a}B^{H}_{j}$ are well defined one has

[TABLE]

For a filter $a\in\mathbb{R}^{p+1}$ and a fractional Brownian motion $(B^{H}_{t})_{t\in T}$ we define its (normalized) quadratic $a$-variation in the following way:

[TABLE]

where $n\in\mathbb{N}$ with $n-p>0$ is the number of discrete observations on a grid of mesh size $\Delta>0$ (that may depend on $n$) and

[TABLE]

In what follows, we will assume that $\Delta=n^{-\alpha}$ for some $\alpha>0$ whenever the quadratic $a$-variation is considered. Observe that the choice $a=(-1,\,1)$ and $\alpha=1$ would yield the formula for the usual normalized quadratic variation of the fBm (see e.g. [18]).

Let us recall the main result concerning the asymptotic behavior of the quadratic $a$-variation of fBm (see [6] and [11]). It states that the sequence (7) converges to zero almost surely for any filter $a$ and if $M(a)>H+\frac{1}{4}$ then it satisfies a central limit theorem.

Theorem 1**.**

Let $(B^{H}_{t})_{t\in T}$ be a fractional Brownian motion with Hurst index $H\in(0,1)$ and let $a$ be a filter of order $M(a)$ and of length $p+1$ with $p\geq 1$. Then

$V(a,\,n,\,\Delta,\,B^{H})\stackrel{{\scriptstyle\text{a.s.}}}{{\to}}0$* as $n\to\infty.$* 2. 2.

If $M(a)>H+\frac{1}{4}$, then

[TABLE]

where $\sigma_{H}>0$ is an explicit constant.

In the construction of the estimators in Section 4 we need a multidimensional version of this statement which can also be found in [11] and [6]:

Theorem 2**.**

Let $(B^{H}_{t})_{t\in T}$ be a fractional Brownian motion with Hurst index $H\in(0,1)$. For $i=1,\dots,N$ let $a^{i}$ be a filter of order $M(a^{i})>H+\frac{1}{4}$. Then

[TABLE]

*where $\Sigma_{H}$ is a positive definite matrix. *

Notice that, in the case of a filter of order 1, we have the well-known restriction $H<\frac{3}{4}$ for the CLT of the sequence $V(a,\,n,\,\Delta,\,B^{H})$ to hold. This restriction can be lifted by taking a filter of order greater than one, i.e. by considering higher-order increments of the fBm.

3 Stochastic differential equations driven by the fractional Brownian motion

This section introduces the main object of our study: stochastic differential equations driven by the fractional Brownian motion. We will then study how the central limit theorems of the previous section carry over to the $a$-variation of the solution to the SDE.

3.1 Stochastic differential equations with fBm

For $t\in T\subseteq\mathbb{R}^{+}$ we consider the SDE

[TABLE]

where $(B^{H}_{t})_{t\in T}$ is a fractional Brownian motion with Hurst index $H\in(0,1)$ and $x$ is a real number. We assume that $f\in C^{0,\,1}_{b}(T\times\mathbb{R})$ (i.e., continuous in the time component and continuously differentiable in the space component as well as bounded, together with his partial derivative with respect to the second variable, in both components) with

[TABLE]

for a fixed constant $M>0$. We denote by $f^{\prime}$ the derivative of $f$ with respect to its second variable while $\|\cdot\|\infty$ stands for the infinity norm.

In the course of the paper we will be interested in two cases:

• (S1)

: the SDE (9) with $T=[0,\,1]$,

• (S2)

: the SDE (9) with $T=\mathbb{R}^{+}$ for $H\leq\frac{1}{2}$.

The existence of pathwise solutions to (9) in both cases can be seen by considering the process $X-B^{H}$ and thus reducing (S) to an ordinary differential equation. We refer to [15], [4] or [13], among others, for the existence and uniqueness of the solution under assumption (10).

Clearly, solutions to (S2) also solve (S1), if restricted to the unit interval. Therefore statements made for (S2) also apply to (S1). We make this distinction mainly for notational reasons: in (S1) we consider only $\alpha\geq 1$ (where $\Delta=n^{-\alpha}$), such that we do not leave the unit interval as the number of observations ($n$) grows. In (S2) $\alpha$ is allowed to be less than one, and the observation window $n\Delta$ grows with $n$.

We further denote by $Y$ the absolutely continuous component of the equation (9), i.e.

[TABLE]

If $X$ is the solution to the SDE (9) (in either the setting (S1) or (S2)) and $a\in\mathbb{R}^{p+1}$ is a filter of length $p+1$ we define the quadratic $a$-variation of $X$ to be

[TABLE]

where the number of observations $n\in\mathbb{N}$ satisfies $n-p>0$ and the mesh size is chosen to satisfy $\Delta=n^{-\alpha}$ for some $\alpha>0$. Here $\sigma_{a,\,\Delta}=Var(\Delta_{a}B^{H}_{j})$ is again the variance of the increment of $B^{H}$ along the filter $a$.

3.2 Central limit theorems for the quadratic $a$-variation of SDE

This section comprises the core results of this work, combining the limit results of the fBm from Section 2 with estimates for the solution of SDE obtained via Malliavin calculus. We obtain central limit theorems for the SDE (9) providing the basis for the estimation of $H$ to be discussed in Section 4. The results are presented in the Theorems 3 and 4 where the cases $H<\frac{1}{2}$ and $H\in[\frac{1}{2},\,1)$ are treated separately. The proof of Theorem 3 consists of direct estimates and does not require any assumptions on the order $M(a)$ of the filter $a$.

If $H\geq\frac{1}{2}$ these estimates are not sufficient and we rely on the Malliavin Calculus to provide a finer analysis of the correlation. The actual proof of Theorem 4 is postponed to Section 3.4 after some auxiliary results.

First we will consider the case $H<\frac{1}{2}$. Notice that in this case we have $M(a)>H+\frac{1}{4}$ for any filter $a$.

Theorem 3**.**

Assume $H<\frac{1}{2}$, let $X$ be a solution of the SDE (S2) and let a be a filter with $p+1$ components. For $\alpha>\frac{1}{2(1-H)}$ (with $\Delta=n^{-\alpha}$) we have

[TABLE]

with $\sigma_{H}>0$ from Theorem 1.

Proof.

By a binomial expansion we can write:

[TABLE]

with

[TABLE]

with some constants $C^{0}$, $C^{1}$. Due to Slutsky’s lemma and part (b) of Theorem 1 it is enough to show that $\mathbb{E}[|R_{n,\,r}|]\to 0$ for $r=0,\,1$. Using the differences representation (5) we obtain with $b_{j}$ given by (4)

[TABLE]

since $f$ is uniformly bounded by (10). The notation $f\lesssim g$ means that there exists a strictly positive constant $c$ such that $f\leq cg$. Thus, we deduce for $R_{n,\,r},r=0,1$:

[TABLE]

Moreover, since for $i\in\{0,\dots n\}$, $\Delta_{a}B^{H}_{i}$ is normally distributed, $\mathbb{E}[|\Delta_{a}B^{H}_{i}|]$ can be calculated directly:

[TABLE]

since by (8), $\sigma_{a,\,\Delta}=-\Delta^{2H}\underset{i>j}{\sum\limits_{i,\,j=0}^{p}}a_{i}a_{j}(i-j)^{2H}.$ Hence, we have

[TABLE]

as well as

[TABLE]

both of which converge to zero because of our assumption on $\alpha$. ∎

Remark 1**.**

*Note that for $H<1/2$ one can choose $\alpha=1$, thus considering the usual equidistant partition of the interval $[0,\,1]$. But it is also possible to chose a mesh $\Delta$ less that $\frac{1}{n}$.

In the second theorem the case H\in\big{[}\frac{1}{2},\,1\big{)} is being considered. For this case we will assume $\alpha\geq 1$, which means that the observations do not leave the interval $[0,\,1]$. For $H\geq\frac{1}{2}$ the bounds on the correlation terms in the above proof fail to converge for $\alpha=1$, therefore, we need more elaborate techniques to handle the remainder in Section 3.4. The assumption $\alpha\geq 1$ is purely technical: it permits the use of Lemma 1 in Section 3.3, the core estimate in our Malliavin calculus approach.

Theorem 4**.**

Let $(B^{H}_{t})_{t\in[0,\,1]}$ be a fractional Brownian motion with Hurst index H\in\big{[}\frac{1}{2},\,1\big{)} and a filter $a$ satisfying $M(a)>H+\frac{1}{4}$. If the process $(X_{t})_{t\in[0,\,1]}$ solves (S1) and if $\alpha\geq 1$, $\alpha>\max\left(\frac{2H-1}{2-2H},\,\frac{1}{4-4H}\right)$ (with $\Delta=n^{-\alpha}$), we have

[TABLE]

where $\sigma_{H}$ is the constant from Theorem 1.

Remark 2**.**

All in all, for $H\in\left(\frac{1}{2},\,\frac{3}{4}\right)$ we have the convergence conditions $\alpha>\frac{1}{4-4H}$ and $\alpha>\frac{2H-1}{2-2H}$, both of which are satisfied for $\alpha\geq 1$. When $H>\frac{3}{4}$, we need to chose $\alpha>1$.

Remark 3**.**

The proofs of the Theorems 3 and 4 were based upon the demonstration that the differences $\sqrt{n}(V(a,\,n,\,\Delta,\,X)-V(a,\,n,\,\Delta,\,B^{H}))$ converge to zero in $L^{1}$. This implies that both theorems can be generalised to their multivariate versions, again by means of Slutsky’s lemma and Theorem 2.

3.3 Auxiliary results from Malliavin calculus

We postpone to the Appendix the technical definitions and results from Malliavin calculus. We refer to [4] for the fact that the solution $X$ of (S1) belongs to the Sobolev-space $\mathbb{D}^{1,\,2}$, implying that the process $Y$ defined above is also Malliavin differentiable (compare also the proofs of Lemma 5.1 and Lemma 5.3 in [13] for a generalisation with respect to the time dependence of the drift). Moreover, the Malliavin derivative of $X$ (with respect to the fBm $B^{H}$) satisfies $D_{s}X_{t}=0$ for $s>t$ and

[TABLE]

for $s\leq t$. This linear equation has the explicit solution

[TABLE]

and similarly to the equation considered in [4] this provides the bound

[TABLE]

almost everywhere for $s\leq t$, where $M$ is the bound on $\|f^{\prime}\|_{\infty}$ defined above in (10).

Moreover, we need the following technical result, which analyzes the correlation between the increments of the processes $B^{H}$ and $Y$. Below, $\mathcal{H}$ is the canonical Hilbert space of the fBm (see Section 6).

Lemma 1**.**

Let $H\geq\frac{1}{2}$ and $Y$ be defined by (11) for the SDE (S1). Then

for every $s,\,t,\,u,\,v\in[0,\,1]$ with $s<t$, $u<v$, we have

[TABLE] 2. 2.

for every $s_{i},\,t_{i},\,u_{i},\,v_{i}\in[0,\,1]$ with $s_{i}<t_{i}$, $u_{i}<v_{i}$, $i=1,\,2$.

[TABLE]

Proof.

First note that for $s\in[0,\,1]$

[TABLE]

This can be seen by approximating the integrals using the fact that the integrands are bounded. Since the derivative operator is closed, it carries over to the limit. The chain rule yields

[TABLE]

which results in the bound

[TABLE]

Now we will consider the case $H=\frac{1}{2}$, in which $B^{H}$ is the standard Brownian motion. Since the space $\mathcal{H}$ coincides with the space $L^{2}(\Omega)$, we have

[TABLE]

applying in the last step the above bound (13) on $|D_{\beta}(Y_{u}-Y_{v})|$. Part (a) is proven for $H=\frac{1}{2}$. Part (b) is obtained similarly by the same arguments and bounds.

For the case $H>\frac{1}{2}$ we first calculate (recall that the norm in $|\mathcal{H}|$ is defined in (21))

[TABLE]

consequently, $D_{s}(Y_{u}-Y_{v})\in|\mathcal{H}|$, which enables us to use the scalar product representation given in (22). We get for the case (a), by (13)

[TABLE]

Moreover, it holds that

[TABLE]

since for $0<a<b<1$ and $H>\frac{1}{2}$

[TABLE]

This proves part (a) of the lemma for $H>\frac{1}{2}$. The proof of part (b) follows analogously, since

[TABLE]

∎

3.4 Proof of Theorem 4

The correlation estimates of Lemma 1 allow us to prove the theorem.

Proof.

For both cases we consider the sum $R_{n,\,0}+R_{n,\,1},$ defined as in the proof of Theorem 3. For the summand $R_{n,\,0}$ we get similarly to the previous theorem by the boundedness of $f$

[TABLE]

which converges to zero for $\alpha>\frac{1}{2(2-2H)}$. We rewrite $R_{n,\,1}$ with the help of the differences representation (5). Then

[TABLE]

ehere $\delta$ denotes the Skorohod integral. This decomposition follows directly from the integration-by-parts formula (23) applied to $u=\sum_{j=0}^{n-p}F_{j}h_{j}$ for

[TABLE]

Notice that (see (5)) $h_{j}=D_{\cdotp}\Delta_{a}B^{H}_{j}$. We estimate the first summand in (14). We can write

[TABLE]

for any $j\in{0,\dots,n-p}$, where the last inequality follows from Lemma 1, point 1. Hence we obtain

[TABLE]

which goes to zero if $\alpha>\frac{1}{4-4H}$. Let us turn to the first component $n^{-1/2}\Delta^{-2H}\delta(u)$ and show that it vanishes in $L^{2}$ (implying convergence in probability).

[TABLE]

which follows from a direct application of (24) . An expansion of the differences representation (5) for $\Delta_{a}Y_{j},\Delta_{a}Y_{k}$ and the definition of $h_{j},h_{k}$ gives

[TABLE]

With the use of Lemma 1 applied to each expectation this is easily dominated by

[TABLE]

Analogously we obtain

[TABLE]

where we used the fact that the expectations are again bounded by $M^{2}\Delta^{2}$.

Observe that for fixed $i$, $l$, $\mu$, $\nu$ we will make the sum in (16) larger if we estimate

[TABLE]

where $\eta=\max(\mu,\,\nu)$.

Now the inner product satisfies

[TABLE]

We note that for $H=\frac{1}{2}$ the indicator functions are orthogonal. On the diagonal we have for $i=j$

[TABLE]

For $H>\frac{1}{2}$ also the off-diagonal elements contribute. For $r:=i-j$ large we have $(r+1)^{2H}-2r^{2H}+(r-1)^{2H}\sim 2H(2H-1)r^{2H-2}$ since it uniformly approximates the second derivative of the function $x\mapsto x^{2H}$ at $x=r$. Indeed, a Taylor expansion shows that for $r>1$

[TABLE]

We can now rearrange the off-diagonal part of the sum in (18) and use this approximation:

[TABLE]

Recall now that for the sum $\sum_{k=1}^{n}\frac{1}{k^{s}}$ with $n\in\mathbb{N}$ and $s\in\mathbb{R}$ the following asymptotic equality stems from the Euler–Maclaurin formula.

[TABLE]

where $\zeta(s)$ is the Riemann zeta function. Hence,

[TABLE]

As a final step we verify that correctly scaled we have for the second summand in (14)

[TABLE]

Substituting $\Delta=n^{-\alpha}$ gives us the condition $\alpha>\frac{1}{4-4H}$ as above for the vanishing of the first term and the aditional condition $\alpha>\frac{2H-1}{2-2H}$ for the second. This bound combined with (15) and (14) will imply the conclusion. ∎

4 Estimation of the Hurst parameter

Here we will construct estimators for the Hurst parameter $H$ of the Brownian motion $B^{H}$ driving the SDE (S2) and derive their properties, using the previous results. In particular, we will transfer the almost sure convergence result from Theorem 1 to the quadratic $a$-variation of the solution of the SDE (S2) and apply the central limit theorems from the previous section.

We shall fix in the sequel a filter $a=(a_{0},\dots,a_{p})$ is a filter of order $M(a)$.

For a stochastic process $(Z_{t})_{t\in[0,\,1]}$ we define

[TABLE]

with $n$ such that $n-p>0$ and $\Delta=n^{-\alpha}$ for some $\alpha>0$. Then

[TABLE]

holds, and we deduce from Theorem 1 that

[TABLE]

The following theorem shows the same result for the solution $X$ of (S1).

Theorem 5**.**

Let $(B^{H}_{t})_{t\in[0,\,1]}$ be a fractional Brownian motion with Hurst parameter $H\in(0,1)$ and let $(X_{t})_{t\in[0,\,1]}$ be the solution of the SDE (S2) driven by the fBm $B^{H}$. Then

[TABLE]

Proof.

We have

[TABLE]

so it is enough to show that the last two summands divided by $\sigma_{a,\,\Delta}$ ($\sim\Delta^{2H}$) converge to zero almost surely. Recall that we have from (12) $|\Delta_{a}Y_{j}|\lesssim\Delta$. It follows that

[TABLE]

which goes to [math] for every $\alpha$ as $n$ tends to infinity.

To estimate $\Delta_{a}B^{H}_{j}$ we refer to [14] for the fact that almost all sample paths of a fractional Brownian motion with the Hurst index $H$ are Hölder continuous of order $H-\varepsilon$ for any $\varepsilon>0$ (it follows from the Kolmogorov’s continuity theorem, since $B^{H}$ is self-similar). We get for almost every trajectory, using the differences representation (5):

[TABLE]

Note that the constant $C$ can vary depending on the trajectory. Then we get

[TABLE]

We can ensure that it converges to zero by setting $\varepsilon:=\min\left(\frac{1-H}{2},\,\frac{H}{2}\right)$, and the claim follows. ∎

Note that for the special case $a=(-1,\,1)$ and $\Delta=n^{-1}$ we can deduce directly from this result that the ”standard” estimator $\hat{H}$ given by

[TABLE]

converges almost surely to $H$.

Remark 4**.**

Since the relationship between $V(a,\,n,\,\Delta,X)$ and $U(a,\,n,\,\Delta,X)$ is the same as for the process $B^{H}$, it follows from Theorem 5 that $V(a,\,n,\,\Delta,X)\stackrel{{\scriptstyle\text{a.s.}}}{{\to}}0$.

Now we are in a position to construct an estimator $\hat{H}_{1}$ in the way it was done in [11] (see also [6]) for processes with stationary increments. We know from Theorem 5 that $U(a,\,n,\,\Delta,X)$ asymptotically equals $\sigma_{a,\,\Delta}$, and a simple rearrangement argument gives us

[TABLE]

Let us consider for $i\in\{1,\dots,m\}$ a family of filters $a^{(i)}=(a_{0}^{(i)},\dots,a_{p_{i}}^{(i)})$. For each member the associated variance $\sigma_{a,\Delta}$ is a linear combination of the $(\Delta d)^{2H}$, $d\in\{1,\dots,P\}$, where $P=\max_{i}p_{i}$, the numbers $p_{i}$ being not necessarily different. So if we consider a matrix $A=(A_{ij})\in\mathbb{R}^{m\times P}$ defined by

[TABLE]

then $U_{n}:=(U(a^{(i)},\,n,\,\Delta,X))_{i=1,\dots,m}$ tends almost surely to $AD$, where $D$ denotes the vector $((\Delta d)^{2H})_{d=1,\dots,P}$.

We choose the $a^{(i)}$ such that $A$ is of full rank and estimate $D$ by linear regression (preserving the almost sure convergence property): $\hat{D}=(A^{T}A)^{-1}A^{T}U_{n}$. Now we can consider $\log|\hat{D}|$, which tends to $2H(\log(d\Delta))_{d=1,\dots,P}$, and estimate $H$ by another simple linear regression:

[TABLE]

If we know the observations to be equidistant but do not have access to the actual mesh size $\Delta$, we can estimate $D$ in the same way (because $U_{n}$ does not directly depend on $\Delta$) and then perform a regression with the intercept term $2H\log(\Delta)$. This allows the estimation

[TABLE]

To establish some properties of the estimators $\hat{H}_{1}$ and $\hat{H}_{2}$, we first recall that, as a consequence of Theorem 5, we have the convergence result $U_{n}\stackrel{{\scriptstyle\text{a.s.}}}{{\to}}AD$ and therefore

[TABLE]

Then, evoking the continuous mapping theorem, we obtain $\hat{H}_{1}\stackrel{{\scriptstyle\text{a.s.}}}{{\to}}H$, $\hat{H}_{2}\stackrel{{\scriptstyle\text{a.s.}}}{{\to}}H$. In other words, both estimators are strongly consistent.

The question of asymptotic normality depends upon whether Theorems 3 and 4 are applicable. Choosing filters of order 2 one can ensure that the filter condition is always satisfied, but other assumptions require more care. In particular, if there is no prior knowledge about the true value of $H$, one has to assume that the mesh size condition $\alpha>\frac{1}{4-4H}$ is satisfied. If a certain bound for $H$ is known, say, if $H\in\left(0,\,\frac{3}{4}\right)$, then under the assumption that $\alpha\geq 1$ (which includes the usual partition of $[0,\,1]$) we obtain asymptotic normality. The following theorem summarizes our observations.

Theorem 6**.**

Let $X$ be a solution of the SDE (S2) and let the mesh size $\Delta$ be chosen such that the assumptions of Theorems 3 and 4 are satisfied. Moreover, let $a^{(i)}$ be filters of respective lengths $p^{(i)}+1$ (for $i=1,\dots,m$) such that $M(a^{(i)})>H+\frac{1}{4}$. Then the sequences $\sqrt{n}(\hat{H}_{1}-H)$ and $\sqrt{n}(\hat{H}_{2}-H)$ converge weakly to normally distributed, centered random variables.

Proof.

This is an application of the delta method (mentioned, for example, in [7]) for the vector $V_{n}:=(V(a^{(i)},\,n,\,\Delta,X))_{i=1,\dots,m}$. Its almost sure convergence was shown in Remark 4 and the multivariate convergence in law of $\sqrt{n}V_{n}$ follows due to Remark 3. Since for the construction of $\hat{H}_{1}$, $\hat{H}_{2}$ it underwent only linear and logarithmic transformations, the obtained estimators are indeed asymptotically normal. ∎

5 Simulation study

Similarly to [11] and [6] two filters $(a^{(1)}_{i})_{i\in\{0,\dots,p\}}$ and $(a^{(2)}_{i})_{i\in\{0,\dots,2p\}}$ are considered, where $a^{(2)}$ is obtained by ”thinning” the filter $a^{(1)}$ (i.e., $a^{(2)}_{2k}:=a^{(1)}_{k}$ for $k\in\{0,\dots,p\}$ and zero otherwise). In this case the estimator $\hat{H}_{1}$ simplifies to

[TABLE]

Note that this estimator is independent of the scaling: its strong consistency has been shown for all $\alpha>0$, hence, it is enough to know that the observations are equidistant, and the mesh size $\Delta$ can then be chosen appropriately. It is, moreover, by construction independent of deterministic multiplicative scaling factors of the fBm involved, which allows us to include the (slightly more general) case

[TABLE]

with some unknown $\sigma>0$ and observed over an unknown interval in our simulation study. From these observations different settings for the simulations arise.

We simulate $100$ trajectories of a process $X$ defined by

[TABLE]

Some trajectories are displayed in Figure 1 for different values of $H$. We calculate the MSE for $1000$, $2000$, $4000$ and $8000$ observations in the following settings:

1. (S1)

$H=0.7$ on an interval $[0,\,1]$ with the ”standard” estimator $\hat{H}$, $\sigma=1$, 2. (S2)

$H=0.7$ on an interval $[0,\,1]$ with $\hat{H}_{1}$ for $a^{(1)}:=(-1,\,1)$, $\sigma=1$, 3. (S3)

$H=0.7$ on an interval $[0,\,1]$ with $\hat{H}_{1}$ for $a^{(1)}:=(-1,\,1)$, $\sigma=5$, 4. (S4)

$H=0.7$ on an interval $[0,\,10]$ with $\hat{H}_{1}$ for $a^{(1)}:=(-1,\,1)$, $\sigma=1$, 5. (S5)

$H=0.98$ on an interval $[0,\,1]$ with the ”standard” estimator $\hat{H}$, $\sigma=1$, 6. (S6)

$H=0.98$ on an interval $[0,\,1]$ with $\hat{H}_{1}$ for $a^{(1)}:=(1,\,-2,\,1)$, $\sigma=1$, 7. (S7)

$H=0.98$ on an interval $[0,\,10]$ with $\hat{H}_{1}$ for $a^{(1)}:=(1,\,-2,\,1)$, $\sigma=1$.

We obtain the following results:

[TABLE]

.

Indeed, one can see that the estimators perform well for a broad range of conditions. For $H<\frac{3}{4}$, in case of known $\sigma$ and on the unit interval the simplest estimator $\hat{H}$ seems to perform best. However, for large $H$ its convergence becomes slower compared to that of $\hat{H}_{1}$ for the order $2$ filter. It is a known result (see e.g. [5]) that for the filter $(-1,\,1)$ the asymptotic rates of convergence for the quadratic variation of an fBm are $\sqrt{n\log(n)}$ in the case $H=\frac{3}{4}$ and $n^{2H-1}$ for $H>\frac{3}{4}$. This constraint seems to transfer to the process $X$. Figures 1 and 2 provide a visualisation for the difference in the convergence rates.

6 Appendix: Definitions and standard results from Malliavin calculus

In the following we will talk about the main definitions of Malliavin calculus for the fractional Brownian motion case, following mainly [14]. Throughout this section we will consider $T=[0,\,1]$.

For a real separable Hilbert space $(\mathcal{H},\,\langle\cdotp,\,\cdotp\rangle_{\mathcal{H}})$ we call a stochastic process $W=\{W(h),\,h\in\mathcal{H}\}$ in a complete probability space $(\Omega,\,\mathcal{F},\,\mathbb{P})$ a Gaussian process on $\mathcal{H}$ if $W$ is a centred Gaussian family of random variables such that for all $f,\,g\in\mathcal{H}$

[TABLE]

For a fractional Brownian motion $(B^{H}_{t})_{t\in T}$ with Hurst parameter $H\in(0,\,1)$ let $\mathcal{H}$ be the closure of the set of indicator functions with respect to the inner product

[TABLE]

Then $B^{H}$ is via $B^{H}(1_{[0,\,s]})=B^{H}_{s}$ by definition a Gaussian process on $\mathcal{H}$, and this is how the notation $\mathcal{H}$ is used in this work.

For $H>\frac{1}{2}$ the space $\mathcal{H}$ does not only contain functions (see [17]), so we consider the following definition:

[TABLE]

It has been shown in [16] that $|\mathcal{H}|$ is a subspace of $\mathcal{H}$, which yields a sufficient condition for functions to belong to the space $\mathcal{H}$.

We have a useful representation of the inner product of two functions from the space $|\mathcal{H}|$. That is, if $H>\frac{1}{2}$, then for $f,\,g\in|\mathcal{H}|$ (see e.g. [14])

[TABLE]

For $H=\frac{1}{2}$ (that is, if $B^{H}$ is the usual Brownian motion) the space $\mathcal{H}$ is identical with the space $L^{2}([0,\,1])$ endowed with the usual inner product.

For the Gaussian process $(B^{H}_{t})_{t\in T}$ on $\mathcal{H}$ let $\mathcal{S}$ denote the set of smooth random variables of the form

[TABLE]

where $C^{\infty}_{p}(\mathbb{R}^{n})$ denotes the space of infinitely continuously differentiable functions $g:\mathbb{R}^{n}\mapsto\mathbb{R}$ such that $g$ and all its partial derivatives grow at most polynomially.

Then we can define the Malliavin derivative $D$ as an operator which acts on functions $F=g(B^{H}(h_{1}),\dots,B^{H}(h_{n}))\in\mathcal{S}$ by

[TABLE]

which means that $DF$ is a random variable with values in $\mathcal{H}$. It does not depend on the choice of the function $g$ and of $h_{1},\dots,h_{n}\in\mathcal{H}$.

For the subset of $\mathcal{S}$, whose functions $g$ have compact support, we will write $\mathcal{S}_{0}$.

$D$ is clearly linear and it holds: $D(B^{H}_{t}-B^{H}_{s})=1_{[s,\,t]}(\cdotp)$.

Moreover, $D:\mathcal{S}\to L^{2}(\Omega;\,\mathcal{H})$ is closable (this has been shown in [14]) and we will denote by $\mathbb{D}^{1,\,2}$ the closure of $\mathcal{S}$ with respect to the norm

[TABLE]

For a fixed $h\in\mathcal{H}$ we define the operator $D^{h}$ on the set $\mathcal{S}$ by

[TABLE]

Since Malliavin derivatives are $\mathcal{H}$-valued random variables, one can identify them with stochastic processes if those values in $\mathcal{H}$ are functions. In this case we will write $D_{s}X$ for $DX(s)$.

We have a chain rule for Malliavin derivatives, which is stated in [14]: Let $f\in C(\mathbb{R}^{m};\mathbb{R})$ be a function with $\|\partial_{i}f\|_{\infty}<M_{i}<\infty$ for some $M_{i}>0$ ($i=1,\dots,m$) and let $X_{1},\dots X_{m}\in\mathbb{D}^{1,\,2}$. Then $f(X_{1},\dots,X_{m})\in\mathbb{D}^{1,\,2}$ and

[TABLE]

Let $\mathcal{S}_{\mathcal{H}}$ denote the set

[TABLE]

Then for $u=\sum_{i=1}^{n}F_{i}h_{i}$ we can define $Du:=\sum_{i=1}^{n}DF_{i}\otimes h_{i}$ and consider the norm

[TABLE]

Now we can, just as for the space $\mathcal{S}$, consider the closure of $\mathcal{S}_{\mathcal{H}}$ with respect to this norm and call it $\mathbb{D}^{1,\,2}(\mathcal{H})$.

The divergence operator, denoted by $\delta$, is the adjoint of the derivative operator $D$. Its domain $\operatorname{Dom}(\delta)$ is the set of random variables $u\in L^{2}(\Omega;\,\mathcal{H})$ such that

[TABLE]

for all $F\in\mathbb{D}^{1,\,2}$, where $c_{u}$ is a constant.

For $u\in\operatorname{Dom}(\delta)$ $\delta(u)$ satisfies the following characterising equation

[TABLE]

for all $F\in\mathbb{D}^{1,\,2}$.

A result from [14] provides a more direct characterisation for a certain elementary subset of $\operatorname{Dom}(\delta)$: let $u=\sum_{i=1}^{n}F_{i}h_{i}$ be an element of $\mathcal{S}_{\mathcal{H}}$. Then $u$ belongs to $\operatorname{Dom}(\delta)$ and we have

[TABLE]

We need the following identity for the expectation of $\delta(u)\delta(v)$, where $u$ and $v$ belong to the domain of $\delta$ and have a certain specific form. In [14] a statement for more general $u$ and $v$ is given, however, in the present work we only need this special case: for $F_{1},\,F_{2}\in\mathbb{D}^{1,\,2}$, $h_{1},\,h_{2}\in\mathcal{H}$ the variables $u:=F_{1}h_{1}$, $v:=F_{2}h_{2}$ belong to the space $\operatorname{Dom}(\delta)$ and

[TABLE]

Acknowledgement: The financial support of the Collaborative Research Center ’Statistical modeling of non-linear dynamic processes’ (SFB 823) is gratefully acknowledged.