Asymptotic Normality of a Hurst Parameter Estimator Based on the Modified Allan Variance

Abstract

In order to estimate the memory parameter of Internet traffic data, it has been recently proposed a log-regression estimator based on the so-called modified Allan variance (MAVAR). Simulations have shown that this estimator achieves higher accuracy and better confidence when compared with other methods. In this paper we present a rigorous study of the MAVAR log-regression estimator. In particular, under the assumption that the signal process is a fractional Brownian motion, we prove that it is consistent and asymptotically normally distributed. Finally, we discuss its connection with the wavelets estimators.

1. Introduction

It is well known that different kinds of real data (hydrology, telecommunication networks, economics, and biology) display self-similarity and long-range dependence (LRD) on various time scales. By self-similarity we refer to the property that a dilated portion of a realization has the same statistical characterization as the original realization. This can be well represented by a self-similar random process with a given scaling exponent 𝐻 (Hurst parameter). The long-range dependence, also called long memory, emphasizes the long-range time correlation between past and future observations and it is thus commonly equated to an asymptotic power law behaviour of the spectral density at low frequencies or, equivalently, to an asymptotic power-law decrease of the autocovariance function, of a given stationary random process. In this situation, the memory parameter of the process is given by the exponent 𝑑 characterizing the power law of the spectral density. (For a review of historical and statistical aspects of the self-similarity and the long memory see [1].)

Though a self-similar process cannot be stationary (and thus nor LRD), these two properties are often related in the following sense. Under the hypothesis that a self-similar process has stationary (or weakly stationary) increments, the scaling parameter 𝐻 enters in the description of the spectral density and covariance function of the increments, providing an asymptotic power law with exponent 𝑑=2𝐻−1. Under this assumption, we can say that the self-similarity of the process reflects on the long-range dependence of its increments. The most paradigmatic example of this connection is provided by the fractional Brownian motion and by its increment process, the fractional Gaussian noise [2].

In this paper we will consider the problem of estimating the Hurst parameter 𝐻 of a self-similar process with weakly stationary increments. Among the different techniques introduced in the literature, we will focus on a method based on the log-regression of the so-called modified Allan variance (MAVAR). The MAVAR is a well-known time-domain quantity generalizing the classic Allan variance [3, 4], which has been proposed for the first time as a traffic analysis tool in [5]. In a series of paper [5–7], its performance has been evaluated by simulation and comparison with the real IP traffic. These works have pointed out the high accuracy of the method in estimating the parameter 𝐻 and have shown that it achieves better confidence if compared with the well-established wavelet log-diagram.

The aim of the present work is to substantiate and enrich these results from the theoretical point of view, studying the rate of convergence of the estimator toward the memory parameter. In particular, our goal is to provide the limit properties and the precise asymptotic normality of the MAVAR log-regression estimator in order to compute the related asymptotic confidence intervals. This will be reached under the assumption that the signal process is a fractional Brownian motion. Although this hypothesis may look restrictive (indeed this estimator is designed for a larger class of stochastic processes), the obtained results are a first step toward the mathematical study of the MAVAR log-regression estimator. To our knowledge, there are no similar results in the literature. The present paper also provides the theoretical foundations (mathematical details and proofs) for [8]. Indeed the formulas analytically obtained here have been implemented and numerically tested in [8] for different choices of the regression weights and it has been shown that the numerical evidences are in good agreement with the theoretical results proven here.

For a survey on Hurst parameter estimators of a fractional Brownian motion, we refer to [9, 10]. However we stress once again that the MAVAR-estimator is not specifically a target for the fractional Brownian motion, but it has been thought and successfully used for more general processes.

Our theorems can be viewed as a counterpart of the already established results concerning the asymptotic normality of the wavelet log-regression estimators [11–13]. Indeed, although the MAVAR can be related in some sense to a suitable Haar-type wavelets family (see [14] for the classical Allan variance), the MAVAR and wavelets log-regression estimators do not match as the regression runs on different parameters (see Section 5). Hence, we adopt a different argument which in turns allows us to avoid some technical troubles due to the poor regularity of the Haar-type functions.

The paper is organized as follows. In Section 2 we recall the properties of self-similarity and long-range dependence for stochastic processes and the definition of fractional Brownian motion; in Section 3 we introduce the MAVAR and its estimator, with their main properties; in Section 4 we state and prove the main results concerning the asymptotic normality of the estimator; in Section 5 we make some comments on the link between the MAVAR and the wavelet estimators and on the modified Hadamard variance, which is a generalization of the MAVAR; in the Appendix we recall some results used along the proof.

2. Self-Similarity and Long-Range Dependence

In the sequel we consider a centered real-valued stochastic process 𝑋={𝑋(𝑡),𝑡∈ℝ}, that can be interpreted as the signal process. Sometimes it is also useful to consider the 𝜏-increment of the process 𝑋, which is defined, for every 𝜏>0 and 𝑡∈ℝ, as 𝑌_𝜏(𝑡)=𝑋(𝑡+𝜏)−𝑋(𝑡)𝜏. (2.1)

In order to reproduce the behavior of the real data, it is commonly assumed that 𝑋 satisfies one of the two following properties: (i) self-similarity; (ii) long range dependence.(i)The self-similarity of a centered real-valued stochastic process 𝑋, starting from zero, refers to the existence of a parameter 𝐻∈(0,1), called Hurst index or Hurst parameter of the process, such that, for all 𝑎>0, it holds {𝑋(𝑡),𝑡∈ℝ}=^𝑑𝑎^−𝐻𝑋(𝑎𝑡),𝑡∈ℝ. (2.2) In this case we say that 𝑋 is an 𝐻-self-similar process. (ii)We first recall that a centered real-valued stochastic process 𝑋 is weakly stationary if it is square integrable and its autocovariance function, 𝐶_𝑋(𝑡,𝑠)∶=Cov(𝑋(𝑡),𝑋(𝑠)), is a translation invariant, namely, if 𝐶_𝑋(𝑡,𝑠)=𝐶_𝑋(𝑡+𝑟,𝑠+𝑟),∀𝑡,𝑠,𝑟∈ℝ. (2.3) In this case we also set 𝑅_𝑋(𝑡)∶=𝐶_𝑋(𝑡,0). If 𝑋 is a weakly stationary process, we say that it displays a long-range dependence, or long memory, if there exists 𝑑∈(0,1) such that the spectral density of the process, 𝑓_𝑋(𝜆), satisfies the condition 𝑓_𝑋(𝜆)∼𝑐_𝑓|𝜆|^−𝑑as𝜆→0, (2.4) for some finite constant 𝑐_𝑓≠0, where we write 𝑓(𝑥)∼𝑔(𝑥) as 𝑥→𝑥₀, if lim_𝑥→𝑥0(𝑓(𝑥)/𝑔(𝑥))=1. Due to the correspondence between the spectral density and the autocovariance function, given by 𝑅_𝑋(𝑡)=12𝜋_ℝ𝑒^𝑖𝑡𝜆𝑓_𝑋(𝜆)𝑑𝜆, (2.5) the long-range condition (2.4) can be often stated in terms of the autocovariance of the process as 𝑅_𝑋(𝑡)∼𝑐_𝑅|𝑡|^−𝛽as|𝑡|→+∞, (2.6) for 𝛽=(1−𝑑)∈(0,1) and some finite constant 𝑐_𝑅≠0.

Notice that if 𝑋 is a self-similar process, then it obviously cannot be weakly stationary. On the other hand, assuming that 𝑋 starts from zero and it is as 𝐻-self-similar process with weakly stationary increments, that is, the quantity 𝔼𝑋𝜏₂+𝑠+𝑡−𝑋(𝑠+𝑡)𝑋𝜏₁+𝑠−𝑋(𝑠) (2.7) does not depend on 𝑠, it turns out that the autocovariance function is given by 𝐶_𝑋()=𝜎𝑠,𝑡^2𝐻2|𝑡|^2𝐻−|𝑡−𝑠|^2𝐻+|𝑠|^2𝐻, (2.8) with 𝜎²_𝐻∶=𝔼[𝑋²(1)], which is clearly not a translation invariant. Consequently, denoting by 𝑌_𝜏 its 𝜏-increment process (see (2.1)), the autocovariance function of 𝑌_𝜏 is such that ([15]) 𝑅_𝑌𝜏(𝑡)∼𝜎²_𝐻𝐻(2𝐻−1)|𝑡|^2𝐻−2as|𝑡|→+∞. (2.9) In particular, if 𝐻∈(1/2,1), the process 𝑌_𝜏 displays a long-range dependence in the sense of (2.6) with 𝛽=2−2𝐻. Under this assumption, we thus embrace the two main empirical properties of a wide collection of real data.

A basic example of the connection between self-similarity and long-range dependence is provided by the fractional Brownian motion 𝐵_𝐻={𝐵_𝐻(𝑡),𝑡∈ℝ} [2]. This is a centered Gaussian process, starting from zero, with autocovariance function given by (2.8), where [16] 𝜎²_𝐻=1Γ()().2𝐻+1sin𝜋𝐻 (2.10) It can be proven that 𝐵_𝐻 is a self-similar process with Hurst index 𝐻∈(0,1), which corresponds, for 𝐻=1/2, to the standard Brownian motion. Moreover, its increment process 𝐺_𝜏,𝐻(𝑡)=𝐵_𝐻(𝑡+𝜏)−𝐵_𝐻(𝑡)𝜏, (2.11) called fractional Gaussian noise, turns out to be a weakly stationary Gaussian process [2, 17].

In the next sections we will perform the analysis of the modified Allan variance and of the related estimator of the memory parameter.

3. The Modified Allan Variance

In this section we introduce and recall the main properties of the modified Allan variance (MAVAR) and of the log-regression estimator of the memory parameter based on it.

Suppose that 𝑋 is a centered real-valued stochastic process, starting from zero, with weakly stationary increments. Let 𝜏₀>0 be the “sampling period” and define the sequence of times {𝑡_𝑘}_𝑘≥1 taking 𝑡₁∈ℝ and setting 𝑡_𝑖−𝑡_𝑖−1=𝜏₀, that is, 𝑡_𝑖=𝑡₁+𝜏₀(𝑖−1).

Definition 3.1. For any fixed integer 𝑝≥1, the modified Allan variance (MAVAR) is defined [3] as 𝜎²_𝑝=𝜎²𝜏₀1,𝑝∶=2𝜏²₀𝑝²𝔼1𝑝^𝑝_𝑗=1𝑋𝑡_𝑗+2𝑝−2𝑋𝑡_𝑗+𝑝+𝑋𝑡_𝑗²=12𝜏²𝔼1𝑝^𝑝_𝑗=1𝑋𝑡_𝑗+2𝜏−2𝑋𝑡_𝑗+𝜏+𝑋𝑡_𝑗², (3.1) where we set 𝜏∶=𝜏₀𝑝. For 𝑝=1 we recover the well-known Allan variance.

Let us assume that a finite sample 𝑋₁,…,𝑋_𝑛 of the process 𝑋 is given, and that the observations are taken at times 𝑡₁,…,𝑡_𝑛. In other words we set 𝑋_𝑖=𝑋(𝑡_𝑖) for 𝑖=1,…,𝑛.

A standard estimator for the modified Allan variance (MAVAR estimator) is given by 𝜎²_𝑝(𝑛)=𝜎²𝜏₀1,𝑝,𝑛∶=2𝜏²₀𝑝⁴𝑛_𝑝^𝑛𝑝_𝑟=1^{𝑝+𝑟−1}_𝑗′=𝑟𝑋_{𝑗′+2𝑝}−2𝑋_𝑗′+𝑝+𝑋_𝑗′²=1𝑛_𝑝^𝑛𝑝−1_𝑘=01√2𝜏₀𝑝^2𝑝_𝑗=1𝑋_{𝑘+𝑗+2𝑝}−2𝑋_{𝑘+𝑗+𝑝}+𝑋_𝑘+𝑗², (3.2) for 𝑝∈{1,…,⌊𝑛/3⌋} and 𝑛_𝑝∶=𝑛−3𝑝+1.

For 𝑘∈ℤ, let us set 𝑑_𝑝,𝑘=𝑑𝜏₀1,𝑝,𝑘∶=√2𝜏₀𝑝^2𝑝_𝑗=1𝑋_{𝑘+𝑗+2𝑝}−2𝑋_{𝑘+𝑗+𝑝}+𝑋_𝑘+𝑗, (3.3) so that the process {𝑑_𝑝,𝑘}_𝑘 turns out to be weakly stationary for each fixed 𝑝, with 𝔼[𝑑_𝑝,𝑘]=0, and we can write 𝜎²_𝑝=𝔼𝑑²_𝑝,0=𝔼𝑑²_𝑝,𝑘,𝜎²_𝑝(𝑛)=1𝑛_𝑝^𝑛𝑝−1_𝑘=0𝑑²_𝑝,𝑘. (3.4)

3.1. Some Properties

Let us further assume that 𝑋 is an 𝐻-self-similar process (see (2.2)), with 𝐻∈(1/2,1). Applying the covariance formula (2.8), it holds 𝜎²_𝑝=𝔼𝑑²_𝑝,0=𝜎²_𝐻𝜏^2𝐻−2𝐾𝐻,𝑝=𝜎²_𝐻𝐾𝐻,𝑝𝜏^2−2𝐻, (3.5) with 𝜎²_𝐻=𝔼[𝑋(1)²] and 𝐾2𝐻,𝑝∶=𝑝1−2^2𝐻−2+12𝑝^2𝑝−1ℓ_ℓ=1ℎ=1𝑃_2𝐻ℎ𝑝, (3.6) where 𝑃_2𝐻 is the polynomial of degree 2𝐻 given by 𝑃_2𝐻(𝑥)∶=−6|𝑥|^2𝐻+4|1+𝑥|^2𝐻+4|1−𝑥|^2𝐻−|2+𝑥|^2𝐻−|2−𝑥|^2𝐻. (3.7)

Since we are interested in the limit for 𝑝→∞, we consider the approximation of the two finite sums in (3.6) by the corresponding double integral, namely, 1𝑝^2𝑝−1ℓ_ℓ=1ℎ=1𝑃_2𝐻ℎ𝑝=¹₀^𝑦₀𝑃_2𝐻(𝑥)𝑑𝑥𝑑𝑦+𝑂_𝐻𝑝⁻¹. (3.8) Computing the integral and inserting the result in (3.6), we get 𝐾𝐻,𝑝=𝐾(𝐻)+𝑂_𝐻𝑝⁻¹, (3.9) where 𝐾(𝐻)2∶=^2𝐻+4+2^2𝐻+3−3^2𝐻+2−152().2𝐻+1)(2𝐻+2 (3.10) From (3.5) and (3.9), we get 𝜎²_𝑝−𝜎²_𝐻𝜏^2𝐻−2𝐾(𝐻)≤𝜎²_𝐻𝜏^2𝐻−2𝑂_𝐻𝑝⁻¹. (3.11)

Under the above hypotheses on 𝑋, one can also prove that the process {𝑑_𝑝,𝑘}_𝑝,𝑘 satisfies the stationary condition Cov𝑑_𝑝,𝑘,𝑑_{𝑝−𝑢,𝑘′}=Cov𝑑_{𝑝,𝑘−𝑘′},𝑑_{𝑝−𝑢,0}for0≤𝑢<𝑝. (3.12) To verify this condition, we write explicitly the covariance as 𝔼𝑑_𝑝,𝑘𝑑_{𝑝−𝑢,𝑘′}=12𝜏²𝑝−𝑢^2𝑝_𝑗=1^𝑝−𝑢_𝑗′=1𝔼𝑋_{𝑘+𝑗+2𝑝}−𝑋_{𝑘+𝑗+𝑝}𝑋_{𝑘′+𝑗′+2(𝑝−𝑢)}−𝑋_{𝑘′+𝑗′+(𝑝−𝑢)}−𝔼𝑋_{𝑘+𝑗+2𝑝}−𝑋_{𝑘+𝑗+𝑝}𝑋_{𝑘′+𝑗′+(𝑝−𝑢)}−𝑋_{𝑘′+𝑗′}−𝔼𝑋_{𝑘+𝑗+𝑝}−𝑋_𝑘+𝑗𝑋_{𝑘′+𝑗′+2(𝑝−𝑢)}−𝑋_{𝑘′+𝑗′+(𝑝−𝑢)}+𝔼𝑋_{𝑘+𝑗+𝑝}−𝑋_𝑘+𝑗𝑋_{𝑘′+𝑗′+(𝑝−𝑢)}−𝑋_{𝑘′+𝑗′}. (3.13) Setting 𝑅𝑡;𝜏₁,𝜏₂∶=𝔼𝑋𝜏₂+𝑠+𝑡−𝑋(𝑠+𝑡)𝑋𝜏₁(𝑠),𝑟𝜏+𝑠−𝑋𝑞,ℎ,𝑝,𝑢∶=𝑅₀𝑞+ℎ+𝑢;𝜏₀𝑝−𝑢,𝜏₀𝑝−𝑅𝜏₀𝑞+ℎ+𝑝;𝜏₀𝑝−𝑢,𝜏₀𝑝−𝑅𝜏₀𝑞+ℎ−𝑝−𝑢;𝜏₀𝑝−𝑢,𝜏₀𝑝+𝑅𝜏₀𝑞+ℎ;𝜏₀𝑝−𝑢,𝜏₀𝑝, (3.14) we get 𝔼𝑑_𝑝,𝑘𝑑_{𝑝−𝑢,𝑘′}=12𝜏²1𝑝−𝑢^2𝑝_𝑗=1^𝑝−𝑢_𝑗′=1𝑟𝑘−𝑘^′,𝑗−𝑗^′,𝑝,𝑢. (3.15) This immediately provides the stationary condition (3.12).

To better understand the behavior of the covariance 𝔼[𝑑_𝑝,𝑘𝑑_{𝑝−𝑢,𝑘′}] as 𝑘−𝑘^′ varies, we apply the covariance formula (2.8) and get 𝑅𝑡;𝜏₁,𝜏₂=𝜎²_𝐻2𝑡+𝜏₂^2𝐻+𝑡−𝜏₁^2𝐻−𝑡+𝜏₂−𝜏₁^2𝐻−|𝑡|^2𝐻. (3.16) Thus, from (3.14), 𝑟𝑞,ℎ,𝑝,𝑢=𝜎²_𝐻𝜏^32𝐻_𝑖=1𝑔_𝐻^(𝑖)𝑞𝑝,ℎ𝑝,𝑢𝑝, (3.17) where 𝑔_𝐻⁽¹⁾𝑞𝑝,ℎ𝑝,𝑢𝑝∶=𝑞+ℎ𝑝+𝑢𝑝+1^2𝐻+𝑞+ℎ𝑝+𝑢𝑝−1^2𝐻−2𝑞+ℎ𝑝+𝑢𝑝^2𝐻,𝑔_𝐻⁽²⁾𝑞𝑝,ℎ𝑝,𝑢𝑝1∶=−2𝑞+ℎ𝑝+2^2𝐻+𝑞+ℎ𝑝^2𝐻−2𝑞+ℎ𝑝+1^2𝐻,𝑔_𝐻⁽³⁾𝑞𝑝,ℎ𝑝,𝑢𝑝1∶=−2𝑞+ℎ𝑝+2𝑢𝑝^2𝐻+𝑞+ℎ𝑝+2𝑢𝑝−2^2𝐻−2𝑞+ℎ𝑝+2𝑢𝑝−1^2𝐻. (3.18) Inserting (3.17) in (3.15), we obtain 𝔼𝑑_𝑝,𝑘𝑑_{𝑝−𝑢,𝑘′}=𝜎²_𝐻𝜏^2𝐻−2𝐺_𝐻𝑘−𝑘^′𝑝,𝑢𝑝,𝑝, (3.19) where 𝜏=𝜏₀𝑝 as before, and 𝐺_𝐻𝑞𝑝,𝑢𝑝1,𝑝∶=21𝑝−𝑢^2𝑝_𝑗=1^𝑝−𝑢_𝑗′³_=1𝑖=1𝑔_𝐻^(𝑖)𝑞𝑝,𝑗−𝑗^′𝑝,𝑢𝑝. (3.20) Now we set 𝑞=𝑚𝑝+𝑟 with 𝑚∈ℤ and 𝑟 an integer number in {0,…,𝑝−1} and we study the asymptotic behavior of 𝐺_𝐻(𝑚+(𝑟/𝑝),𝑢/𝑝,𝑝) for |𝑚|→+∞.

The limit relation |𝑚+1|^2𝐻+|𝑚−1|^2𝐻−2|𝑚|^2𝐻∼2𝐻(2𝐻−1)|𝑚|^2𝐻−2as|𝑚|→+∞, (3.21) that is, lim_{|𝑚|→+∞}|𝑚+1|^2𝐻+|𝑚−1|^2𝐻−2|𝑚|^2𝐻|𝑚|^2𝐻−2=2𝐻(2𝐻−1), (3.22) implies that 𝑔_𝐻⁽¹⁾𝑟𝑚+𝑝,ℎ𝑝,𝑢𝑝(𝑟∼2𝐻2𝐻−1)𝑚+𝑝+ℎ𝑝+𝑢𝑝^2𝐻−2|𝑚|𝑔as→+∞,_𝐻(2)𝑟𝑚+𝑝,ℎ𝑝,𝑢𝑝(𝑟∼−𝐻2𝐻−1)𝑚+𝑝+ℎ𝑝+1^2𝐻−2|𝑚|𝑔as→+∞,_𝐻(3)𝑟𝑚+𝑝𝑢,𝑥,𝑝(𝑟∼−𝐻2𝐻−1)𝑚+𝑝+ℎ𝑝+2𝑢𝑝−1^2𝐻−2as|𝑚|→+∞, (3.23) and so, for |𝑚|→+∞, ³_𝑖=1𝑔_𝐻^(𝑖)𝑟𝑚+𝑝,ℎ𝑝,𝑢𝑝∼−𝐻(2𝐻−1)𝑚+𝑟+ℎ𝑝+1^2𝐻−2+𝑚+𝑟+ℎ+2𝑢𝑝−1^2𝐻−2−2𝑚+𝑟+ℎ+𝑢𝑝^2𝐻−2()×=−𝐻2𝐻−1𝑚+𝑟+𝑢+ℎ𝑝+𝑝−𝑢𝑝^2𝐻−2+𝑚+𝑟+𝑢+ℎ𝑝−𝑝−𝑢𝑝^2𝐻−2−2𝑚+𝑟+𝑢+ℎ𝑝^2𝐻−2=−𝐻(2𝐻−1)𝑝−𝑢𝑝^2𝐻−2×𝑚𝑝+𝑝−𝑢𝑟+𝑢+ℎ𝑝−𝑢+1^2𝐻−2+𝑚𝑝+𝑝−𝑢𝑟+𝑢+ℎ𝑝−𝑢−1^2𝐻−2−2𝑚𝑝𝑝−𝑢𝑞+𝑢+ℎ𝑝−𝑢^2𝐻−2∼𝑐_𝐻𝑝−𝑢𝑝^2𝐻−2𝑚𝑝+𝑝−𝑢𝑟+𝑢+ℎ𝑝−𝑢^2𝐻−4=𝑐_𝐻𝑝−𝑢𝑝²𝑚+𝑟+𝑢+ℎ𝑝^2𝐻−4∼𝑐_𝐻𝑝−𝑢𝑝²|𝑚|^2𝐻−4, (3.24) where 𝑐_𝐻∶=(−4!/2)₄^2𝐻. We can conclude that 𝐺_𝐻𝑞𝑝,𝑢𝑝,𝑝=𝐺_𝐻𝑟𝑚+𝑝,𝑢𝑝∼𝑐,𝑝^𝐻2𝑢1−𝑝|𝑚|^2𝐻−4as|𝑚|→+∞. (3.25)

3.2. The MAVAR Log-Regression Estimator

Let 𝑛 be the sample size, that is, the number of the observations.

Definition 3.2. Let 𝑝,ℓ∈ℕ such that 1≤(1+ℓ)𝑝≤𝑝_max(𝑛)=⌊𝑛/3⌋, and let 𝑤=(𝑤₀,…,𝑤_ℓ) be a vector of weights satisfying the conditions ^ℓ_ℓ=0𝑤_ℓ=0,^ℓ_ℓ=0𝑤_ℓlog(1+ℓ)=1. (3.26) The MAVAR log-regression estimator associated to the weights 𝑤 is defined as 𝛼_𝑛𝑝,𝑤=𝛼_𝑛𝜏₀,𝑝,𝑤∶=^ℓ_ℓ=0𝑤_ℓlog𝜎²_()1+ℓ𝑝(𝑛). (3.27) Roughly speaking, the idea behind this definition is to use the approximation 𝜎²_𝑝(𝑛),…,𝜎²_(1+ℓ)𝑝(𝑛)=𝜎²_𝑝(𝑛),…,𝜎²_(1+ℓ)𝑝(𝑛), (3.28) in order to get, by (3.11) and (3.26), 𝛼_𝑛𝑝,𝑤=^ℓ_ℓ=0𝑤_ℓlog𝜎²_()1+ℓ𝑝=^ℓ_ℓ=0𝑤_ℓ𝛼log(1+ℓ)+𝛼log𝜏₀𝑝+log𝜎²_𝐻𝐾(𝐻)=𝛼, (3.29) where 𝛼∶=2𝐻−2. Thus, given the data 𝑋₁,…,𝑋_𝑛, the following procedure is used to estimate 𝐻: (i)compute the modified Allan variance by (3.2), for integer values (1+ℓ)𝑝, with 1≤(1+ℓ)𝑝≤𝑝_max(𝑛)=⌊𝑛/3⌋; (ii)compute the weighted MAVAR log-regression estimator by (3.27) in order to get an estimate 𝛼 of 𝛼; (iii)estimate 𝐻 by 𝐻=(𝛼+2)/2. In the sequel we will give, under suitable assumptions, two convergence results in order to justify these approximations and to get the rate of convergence of 𝛼_𝑛(𝑝,𝑤) toward 𝛼=2𝐻−2. Obviously, we need to take 𝑝=𝑝(𝑛)→+∞ as 𝑛→+∞ in order to reach jointly the above two approximations.

4. The Asymptotic Normality of the Estimator

Since now on we will always assume that 𝑋 is a fractional Brownian motion (with 𝐻∈(1/2,1)) so that the process {𝑑_𝑝,𝑘}_𝑝,𝑘 is also Gaussian. Under this assumption, and with the notation introduced before, we can state the following results.

Theorem 4.1. Let 𝑝=𝑝(𝑛) be a sequence of integers such that 𝑝(𝑛)→+∞,𝑛𝑝(𝑛)⁻¹→+∞. (4.1) Let ℓ be a given integer, 𝜎²_𝑛(𝑝,ℓ) the vector (𝜎²_𝑝(𝑛),𝜎²_2𝑝(𝑛),…,𝜎²_(1+ℓ)𝑝(𝑛)) and, analogously, set 𝜎²(𝑝,ℓ)=(𝜎²_𝑝,𝜎²_2𝑝,…,𝜎²_(1+ℓ)𝑝). Then, it holds 𝑛𝑝𝜏₀𝑝^2−2𝐻𝜎²_𝑛𝑝,ℓ−𝜎²𝑝,ℓ−−−^𝑑→_𝑛→∞𝒩(0,𝑊(𝐻)), (4.2) where 𝑊(𝐻) is a suitable symmetric matrix.

From this theorem, as an application of the 𝛿-method, we can state the following result.

Theorem 4.2. Let 𝛼_𝑛(𝑝,𝑤) be defined as in (3.27), for some finite integer ℓ and a weight-vector 𝑤 satisfying (3.26). If 𝑝=𝑝(𝑛) is a sequence of integers such that 𝑝(𝑛)→+∞,𝑛𝑝(𝑛)⁻¹→+∞,𝑛𝑝(𝑛)⁻³→0, (4.3) then 𝑛𝑝𝛼_𝑛𝑝,𝑤−𝛼−−−^𝑑→_𝑛→∞𝒩0,𝑤^𝑇_∗𝑉(𝐻)𝑤_∗, (4.4) where 𝛼=2𝐻−2, the column-vector 𝑤_∗ is such that [𝑤_∗]_ℓ∶=𝑤_ℓ(1+ℓ)^2−2𝐻, and 𝑉(𝐻)=(𝜎²_𝐻𝐾(𝐻))⁻²𝑊(𝐻).

Let us stress that from the above result, and due to the condition 𝑛𝑝(𝑛)⁻¹→+∞, it follows that the estimator 𝛼_𝑛(𝑝,𝑤) is consistent.

Before starting the proof of the above theorems, we need the following lemma.

Lemma 4.3. Let 𝑝=𝑝(𝑛) be a sequence of integers such that 𝑝(𝑛)→+∞,𝑛𝑝(𝑛)⁻¹→+∞. (4.5) For two given integers ℓ^′,ℓ, with 0≤ℓ^′≤ℓ, set 𝑝_ℓ′=(ℓ^′+1)𝑝 and 𝑝_ℓ=(ℓ+1)𝑝. Then 𝑛𝑝𝜏₀𝑝^4−4𝐻Cov𝜎²_𝑝ℓ(𝑛),𝜎²_𝑝ℓ′(𝑛)−−−−→_𝑛→+∞𝑊_ℓ′,ℓ(𝐻), (4.6) where 𝑊_ℓ′,ℓ(𝐻) is a finite quantity.

Proof. Since 𝑛/𝑝→+∞, without loss of generality we can assume that (1+ℓ)𝑝≤𝑝_max(𝑛)=⌊𝑛/3⌋ for each 𝑛. Recall the notation 𝑛_𝑝=𝑛−3𝑝+1, and set 𝑛_ℓ=𝑛_𝑝ℓ∼𝑛 and 𝑛_ℓ′=𝑛_𝑝ℓ′∼𝑛. From the definition of the MAVAR estimator and applying the Wick’s rule for jointly Gaussian random variables (see the Appendix), we get Cov𝜎²_𝑝ℓ(𝑛),𝜎²_𝑝ℓ′(𝑛)=1𝑛_ℓ𝑛_ℓ′^{𝑛ℓ−1𝑛}_𝑘=0^ℓ′−1_𝑘′=0Cov𝑑²_𝑝ℓ,𝑘,𝑑²_{𝑝ℓ′,𝑘′}=2𝑛_ℓ𝑛_ℓ′^{𝑛ℓ−1𝑛}_𝑘=0^ℓ′−1_𝑘′=0Cov𝑑_𝑝ℓ,𝑘,𝑑_{𝑝ℓ′,𝑘′}²=2𝜎⁴_𝐻(1+ℓ)^4𝐻−4𝜏₀𝑝^4𝐻−4𝑛_ℓ𝑛_ℓ′^{𝑛ℓ−1𝑛}_𝑘=0^ℓ′−1_𝑘′=0𝐺²_𝐻𝑘𝑘−^′𝑝_ℓ,𝑢_ℓ′ℓ𝑝_ℓ,𝑝_ℓ, (4.7) where 𝑢_ℓ′ℓ=𝑝_ℓ−𝑝_ℓ′. Since, by (3.20), the function 𝐺_𝐻(𝑞/𝑝_ℓ,𝑢_ℓ′ℓ/𝑝_ℓ,𝑝_ℓ) only depends on 𝑞/𝑝_ℓ and 𝑢_ℓ′ℓ/𝑝_ℓ=(ℓ−ℓ^′)/(1+ℓ) as 𝑛→+∞, we rewrite the last line as 2𝜎⁴_𝐻(1+ℓ)^4𝐻−4𝜏₀𝑝^4𝐻−4𝑛_ℓ𝑛_ℓ′^{𝑛ℓ−1𝑛}_𝑘=0^ℓ′−1_𝑘′=0𝐺²_𝐻𝑘𝑘−^′𝑝_ℓ,ℓ−ℓ^′.1+ℓ (4.8) This term, multiplied by (𝑛/𝑝)(𝜏₀𝑝)^4−4𝐻, is equal to 2𝜎⁴_𝐻(1+ℓ)^4𝐻−41𝑝𝑛𝑛_ℓ𝑛_ℓ′^{𝑛ℓ−1𝑛}_𝑘=0^ℓ′−1_𝑘′=0𝐺²_𝐻𝑘𝑘−^′𝑝_ℓ,ℓ−ℓ^′().1+ℓ (4.9) It is easy to see that this quantity converges (to a finite strictly positive limit) if and only if 1𝑝_ℓ^{[𝑛/𝑝ℓ]𝑝}_𝑚=0^ℓ−1_𝑟=0𝐺²_𝐻𝑟𝑚+𝑝_ℓ,ℓ−ℓ^′(1+ℓ) (4.10) converges. From (3.25) it holds 𝐺²_𝐻𝑟𝑚+𝑝_ℓ,ℓ−ℓ^′()∼𝑐1+ℓ4^2𝐻1−ℓ−ℓ^′1+ℓ²𝑚^4𝐻−8as𝑚→+∞. (4.11) Thus, the quantity (4.10) is controlled in the limit 𝑛→∞ by the sum ∑^∞_𝑚=1𝑚^4𝐻−8 that is convergent.

Proof of Theorem 4.1. As before, without loss of generality, we can assume that (1+ℓ)𝑝≤𝑝_max(𝑛) for each 𝑛. Moreover, set again 𝑛_ℓ=𝑛_𝑝ℓ and 𝑛_ℓ′=𝑛_𝑝ℓ′. For a given real vector 𝑣^𝑇=(𝑣₀,…,𝑣_ℓ), let us consider the random variable 𝑇_𝑛=𝑇(𝑝(𝑛),ℓ,𝑣) defined as a linear combination of the empirical variances 𝜎²_𝑝(𝑛),…,𝜎²_(1+ℓ)𝑝(𝑛) as follows: 𝑇_𝑛∶=^ℓ_ℓ=0𝑣_ℓ𝜎²_()1+ℓ𝑝(𝑛)=^ℓ_ℓ=0𝑣_ℓ𝑛_ℓ^𝑛ℓ−1_𝑘=0𝑑²_{(1+ℓ)𝑝,𝑘}. (4.12) In order to prove the convergence stated in Theorem 4.1, we have to show that the random variable 𝑛/𝑝(𝜏₀𝑝)^2−2𝐻(𝑇_𝑛−∑^ℓ_ℓ=0𝑣_ℓ𝜎²_(1+ℓ)𝑝) converges to the normal distribution with zero mean and variance 𝑣^𝑇𝑊(𝐻)𝑣. To this purpose, we note that 𝑛𝑝𝜏₀𝑝^2−2𝐻𝑇_𝑛−^ℓ_ℓ=0𝑣_ℓ𝜎²_()1+ℓ𝑝=𝑉^𝑇_𝑛𝐴_𝑛𝑉_𝑛−𝔼𝑉^𝑇_𝑛𝐴_𝑛𝑉_𝑛, (4.13) where 𝑉_𝑛 is the random vector with entries 𝑑_{(1+ℓ)𝑝,𝑘}, for 0≤ℓ≤ℓ and 0≤𝑘≤𝑛_ℓ−1, and 𝐴_𝑛 is the diagonal matrix with entries 𝐴_{𝑛((1+ℓ)𝑝,𝑘),((1+ℓ)𝑝,𝑘)}=𝑣_ℓ𝑛_ℓ𝑛𝑝𝜏₀𝑝^2−2𝐻=𝑂𝑝^3/2−2𝐻𝑛^−1/2. (4.14) By Lemma 4.3, if 𝑊(𝐻) is the symmetric matrix with [𝑊(𝐻)]_ℓ′,ℓ=𝑊_ℓ′,ℓ when 0≤ℓ^′≤ℓ≤ℓ, it holds Var𝑉^𝑇_𝑛𝐴_𝑛𝑉_𝑛=𝑛𝑝𝜏₀𝑝^4−4𝐻ℓ_ℓ=0^ℓ_ℓ′=0𝑣_ℓ𝑣_ℓ′Cov𝜎²_()1+ℓ𝑝(𝑛),𝜎²_𝑝1+ℓ′(𝑛)−−−−→_𝑛→+∞^ℓ_ℓ=0^ℓ_ℓ′=0𝑣_ℓ𝑣_ℓ′𝑊_ℓ′,ℓ(𝐻)=𝑣^𝑇𝑊(𝐻)𝑣, (4.15) therefore, condition (1) of Lemma A.2 is satisfied.
In order to verify condition (2) of Lemma A.2, let 𝐶_𝑛 denote the covariance matrix of the random vector 𝑉_𝑛, and let 𝜌[𝐶_𝑛] denote its spectral radius. By Lemma A.3, we have 𝜌𝐶_𝑛≤^ℓ_ℓ=0𝜌𝐶_𝑛(1+ℓ)𝑝, (4.16) where 𝐶_𝑛((1+ℓ)𝑝) is the covariance matrix of the subvector (𝑑_{(1+ℓ)𝑝,𝑘})_{𝑘=0,…,𝑛ℓ−1}. Applying the spectral radius estimate (A.3), and from equality (3.19), we then have 𝜌𝐶_𝑛(1+ℓ)𝑝≤𝜎²_𝐻𝜏₀(1+ℓ)𝑝^2𝐻−2𝐺_𝐻0,0,𝑝_ℓ+2^𝑛ℓ−1_𝑞=1𝐺_𝐻𝑞𝑝_ℓ,0,𝑝_ℓ=𝑂𝑝^2𝐻−1. (4.17) In order to conclude, it is enough to note that 𝜌𝐴_𝑛𝜌𝐶_𝑛=𝑂𝑝^3/2−2𝐻𝑛^−1/2𝑂𝑝^2𝐻−1=𝑂𝑝^1/2𝑛^−1/2−−−−→_𝑛→∞0. (4.18)

Proof of Theorem 4.2. By assumptions (4.3) on the sequence 𝑝=𝑝(𝑛), and in particular from the condition 𝑛/(𝑝)³→0, and inequality (3.11), it holds 𝑛𝑝𝜏₀𝑝^2−2𝐻𝜎²_(1+ℓ)𝑝−𝜎²_𝐻𝜏₀(1+ℓ)𝑝^2𝐻−2𝐾(𝐻)≤𝜎²_𝐻(1+ℓ)^2𝐻−2𝑂𝑛^1/2𝑝^−3/2−−−−→_𝑛→+∞0. (4.19) Thus, from Theorem 4.1, we get 𝑛𝑝𝜏₀𝑝^2−2𝐻𝜎²_𝑛𝑝,ℓ−𝜎−−−²_∗^𝑑→_𝑛→∞𝒩(0,𝑊(𝐻)), (4.20) where 𝜎²_∗ is the vector with elements 𝜎²_∗ℓ∶=𝜎²_𝐻(1+ℓ)^2𝐻−2𝐾(𝐻)for0≤ℓ≤ℓ. (4.21) Now we observe that if 𝑓(𝑥)∶=∑^ℓ_ℓ=0𝑤_ℓlog(𝑥_ℓ), then, by (3.26) and (3.27), we have 𝛼_𝑛𝑝,𝑤=𝑓𝜎²𝑝,ℓ=𝑓𝜏₀𝑝^2−2𝐻𝜎²𝑝,ℓ. (4.22) Moreover, 𝛼=𝑓(𝜎²_∗) and ∇𝑓(𝜎²_∗)=(𝜎²_𝐻𝐾(𝐻))⁻¹𝑤_∗. Therefore, by the application of the 𝛿-method, the convergence (4.20) entails 𝑛𝑝𝛼_𝑛𝑝,𝑤−𝛼−−−^𝑑→_𝑛→∞𝒩0,∇𝑓𝜎²_∗^𝑇𝑊(𝐻)∇𝑓𝜎²_∗, (4.23) where ∇𝑓𝜎²_∗^𝑇𝑊(𝐻)∇𝑓𝜎²_∗=𝑤^𝑇_∗𝑉(𝐻)𝑤_∗, (4.24) and thus concludes the proof.

Remark 4.4. By Lemma 4.3, and since by definition 𝑉_ℓ′,ℓ(𝐻)=(𝜎²_𝐻𝐾(𝐻))⁻²𝑊_ℓ′,ℓ(𝐻), an estimate of 𝑉_ℓ′,ℓ(𝐻), with 0≤ℓ^′≤ℓ≤ℓ, is given by 𝑛𝑝𝜏₀𝑝^4−4𝐻𝜎⁴_𝐻𝐾(𝐻)²Cov𝜎²_𝑝ℓ(𝑛),𝜎²_𝑝ℓ′(𝑛). (4.25) Therefore, setting 𝜌²_𝑛(𝑤,𝐻) equal to the corresponding estimate of 𝑤^𝑇_∗𝑉(𝐻)𝑤_∗ divided by 𝑛/𝑝, that is 𝜌²_𝑛𝑤=2,𝐻𝐾(𝐻)^2ℓ_ℓ=0^ℓ_ℓ′=01+ℓ∧ℓ^′1+ℓ∨ℓ^′2−2𝐻𝑤_ℓ𝑤_ℓ′𝑛_ℓ𝑛_ℓ′×^{𝑛′ℓ∨ℓ−1𝑛}_𝑘=0^{′ℓ∧ℓ−1}_𝑘′=0𝐺²_𝐻𝑘−𝑘^′𝑝1+ℓ∨ℓ^′,ℓ−ℓ^′1+ℓ∨ℓ^′,𝑝1+ℓ∨ℓ^′, (4.26) and from (4.6), we obtain the convergence 𝜌_𝑛⁻¹𝑤,𝐻𝛼_𝑛𝑝,𝑤−𝛼−−−^𝑑→_𝑛→∞𝒩(0,1), (4.27) which can be used to obtain an asymptotic confidence interval for the parameter 𝐻.

4.1. An Alternative Representation of the Covariance

It is well known that an FBM (with 𝑋(0)=0 and 𝔼[𝑋²(1)]=𝜎²_𝐻) has the following stochastic integral representation (see [2, 16]): 𝑋(𝑡)=1Γ(𝐻+1/2)_ℝ(𝑡−𝑠)^+𝐻−1/2−(𝑠⁻)^𝐻−1/2𝑑𝑊_𝑠. (4.28) Using this representation and Ito’s formula, we can obtain another representation of the covariance Cov(𝜎²_𝑝(𝑛),𝜎²_𝑝−𝑢(𝑛)) and, consequently, another formula for the coefficient 𝜌_𝑛(𝑤,𝐻).

Indeed, assuming without loss of generality 𝑡₁=𝜏₀ so that 𝑋_𝑖=𝑋(𝑡_𝑖)=𝑋(𝑖𝜏₀), by (3.3) and (4.28), we can write for 0≤𝑢<𝑝 and 𝑘∈ℤ𝑑_{𝑝−𝑢,𝑘}=𝜎_𝐻√2𝜏₀𝑝^𝐻−3/2𝑝²𝑝−𝑢^2ℝ𝛾_𝐻𝑘𝑝𝑠,𝑝,𝑢,𝜏₀𝑝𝑑𝑊_𝑠, (4.29) where 𝛾_𝐻𝑘𝑝𝑠,𝑝,𝑢,𝜏₀𝑝Γ1∶=𝐻+2𝜎_𝐻⁻¹×1𝑝^𝑝−𝑢_𝑗=1𝑘+𝑗𝑝+2𝑝−𝑢𝑝−𝑠𝜏₀𝑝^+𝐻−1/2−2𝑘+𝑗𝑝+𝑝−𝑢𝑝−𝑠𝜏₀𝑝^+𝐻−1/2+𝑘+𝑗𝑝−𝑠𝜏₀𝑝^+𝐻−1/2. (4.30) By Ito’s formula we get 𝑑²_{𝑝−𝑢,𝑘}=𝔼𝑑²_{𝑝−𝑢,𝑘}+𝜎²_𝐻𝜏₀𝑝^2𝐻−3𝑝⁴𝑝−𝑢^4ℝ𝑍𝑘𝑝,𝑝,𝑢,𝜈𝛾_𝐻𝑘𝑝𝜈,𝑝,𝑢,𝜏₀𝑝𝑑𝑊_𝜈, (4.31) where 𝑍𝑘𝑝,𝑝,𝑢,𝜈∶=^𝜈_−∞𝛾_𝐻𝑘𝑝𝑠,𝑝,𝑢,𝜏₀𝑝𝑑𝑊_𝑠. (4.32) Inserting this expression in (3.4), we obtain 𝜎²_𝑝−𝑢(𝑛)−𝜎²_𝑝−𝑢=𝜎²_𝐻𝜏₀𝑝^2𝐻−3𝑝⁴𝑝−𝑢^4ℝ1𝑛^𝑛_𝑝−𝑢^{𝑝−𝑢−1}_𝑘=0𝑍𝑘𝑝,𝑝,𝑢,𝜈𝛾_𝐻𝑘𝑝𝜈,𝑝,𝑢,𝜏₀𝑝𝑑𝑊_𝜈. (4.33) It follows, using Ito’s formula again and after the change of variables 𝑠=𝜏₀(𝑝𝑣+𝑘^′) and 𝜈=𝜏₀(𝑝𝑟+𝑘^′), Cov𝜎²_𝑝(𝑛),𝜎²_𝑝−𝑢(𝑛)=𝜎⁴_𝐻𝜏₀𝑝^4𝐻−4𝑝⁴𝑝−𝑢⁴1𝑛_𝑝𝑛_𝑝−𝑢×^{𝑛𝑝−1𝑛}_𝑘=0^{𝑝−𝑢−1}_𝑘′=0ℝ𝛾_𝐻𝑘−𝑘^′𝑝,𝑝,0,𝑟𝛾_𝐻×0,𝑝,𝑢,𝑟_𝑟−∞𝛾_𝐻𝑘−𝑘^′𝑝,𝑝,0,𝑣𝛾_𝐻0,𝑝,𝑢,𝑣𝑑𝑣𝑑𝑟. (4.34)

5. Some Comments

5.1. The Modified Allan Variance and the Wavelet Estimators

Suppose that 𝑋 is a self-similar process with weakly stationary increments and consider the generalized process 𝑌={𝑌(𝑡),𝑡∈ℝ} defined through the set of identities 𝑌;𝐼_{[𝑡1,𝑡2]}=^𝑡2_𝑡1𝑌(𝑡)𝑑𝑡=𝑋𝑡₂−𝑋𝑡₁,∀𝑡₁,𝑡₂∈ℝ. (5.1) In short, we write 𝑌=𝑋̇. From this definition and with the notation introduced in Section 3, we can rewrite the MAVAR as 𝜎²𝜏₀=1,𝑝2𝜏²₀𝑝²𝔼1𝑝^𝑝_𝑖=1^{𝑡𝑖+2𝑝}_{𝑡𝑖+𝑝}𝑌(𝑡)𝑑𝑡−^{𝑡𝑖+𝑝}_𝑡𝑖𝑌(𝑡)𝑑𝑡² (5.2) and its related estimator as 𝜎²𝜏₀=1,𝑝,𝑛2𝜏²₀𝑝²𝑛_𝑝^𝑛𝑝−1_𝑘=01𝑝^𝑝_𝑖=1^{𝑡𝑖+𝑘+2𝑝}_{𝑡𝑖+𝑘+𝑝}𝑌(𝑡)𝑑𝑡−^{𝑡𝑖+𝑘+𝑝}_{𝑡𝑖+𝑘}𝑌(𝑡)𝑑𝑡². (5.3) Now we claim that, for 𝑝 fixed, the quantity 𝑑𝜏₀1,𝑝,𝑘∶=√2𝑝²𝜏₀^𝑝_𝑖=1^{𝑡𝑖+𝑘+2𝑝𝑡}_{𝑖+𝑘+𝑝}𝑌(𝑡)𝑑𝑡−^{𝑡𝑖+𝑘+𝑝}_{𝑡𝑖+𝑘}𝑌(𝑡)𝑑𝑡 (5.4) recalls a family of discrete wavelet transforms of the process 𝑌, indexed by 𝜏₀ and 𝑘. To see that, let us fix 𝑗∈ℕ and 𝑘∈ℤ and set 𝜏₀=2^𝑗 and 𝑡₁=2^𝑗, so that 𝑡_𝑖+𝑘=2^𝑗(𝑖+𝑘), for all 𝑖∈ℕ. With this choice on the sequence of times, it is not difficult to construct a function 𝜓(𝑠) such that 𝑑_𝑘,𝑗∶=𝑑2^𝑗,𝑝,𝑘=𝑌;𝜓_𝑘,𝑗with𝜓_𝑘,𝑗(𝑠)∶=2^−𝑗𝜓2^−𝑗𝑠−𝑘. (5.5) An easy check shows that the function 𝜓(𝑠)1∶=𝑝√𝑝^𝑝_𝑖=1𝜓^𝑖(𝑠),𝜓^𝑖(𝑠)1∶=𝐼2𝑝^{[𝑖+𝑝,𝑖+2𝑝]}(𝑠)−𝐼_{[𝑖,𝑖+𝑝]}(𝑠), (5.6) satisfies (5.5). Notice also that the components 𝜓^𝑖, 𝑖=1,…𝑝, are suitably translated and renormalized Haar functions.

In the case 𝑝=1, corresponding to the classical Allan variance, the function 𝜓 is exactly given by the Haar mother wavelet, as was already pointed out in [14].

Though the wavelet representation could be convenient in many respects, the Haar mother wavelet does not satisfy one of the conditions which are usually required in order to study the convergence of the estimator (see condition (W2) in [18]). Moreover, there is a fundamental difference between the two methods: in the wavelet setting the log-regression is done over the scale parameter 𝜏₀ for 𝑝 fixed, while the MAVAR log-regression given in (3.27) is performed over 𝑝 and for 𝜏₀ fixed.

5.2. The Modified Hadamard Variance

Further generalizing the notion of the MAVAR, one can define the modified Hadamard variance (MHVAR). For fixed integers 𝑝,𝑄, and 𝜏₀∈ℝ, set 𝜎²𝜏₀1,𝑝,𝑄∶=𝑄!𝜏²₀𝑝²𝔼1𝑝^𝑝𝑄_𝑖=1𝑞=0𝑐𝑄,𝑞𝑋𝑡_{𝑖+𝑞𝑝}², (5.7) where 𝑐(𝑄,𝑞)∶=(−1)^𝑞(𝑄!/𝑞!(𝑄−𝑞)!). Notice that for 𝑄=2 we recover the modified Allan variance. The MHVAR is again a time-domain quantity which has been introduced in [7] for the analysis of the network traffic. A standard estimator for this variance is given by 𝜎²_𝑛𝜏₀1,𝑝,𝑄∶=𝑄!𝜏²₀𝑝⁴𝑛^𝑛_𝑄,𝑝^𝑄,𝑝_ℎ=1^𝑝+ℎ−1_𝑖′^𝑄_=ℎ𝑞=0𝑐𝑄,𝑞𝑋_{𝑖′+𝑞𝑝}²=1𝑛^𝑛_𝑄,𝑝^{𝑄,𝑝−1}_𝑘=01𝑄!𝜏₀𝑝^2𝑝𝑄_𝑖=1𝑞=0𝑐𝑄,𝑞𝑋_{𝑘+𝑖+𝑞𝑝}², (5.8) for 𝑝=1,…,[𝑛/(𝑄+1)] and 𝑛_𝑄,𝑝∶=𝑛−(𝑄+1)𝑝+1.

Similarly to the analysis performed for the MAVAR, let us set 𝑑_𝑝,𝑘=𝑑𝜏₀1,𝑝,𝑄,𝑘∶=𝑄!𝜏₀𝑝^2𝑝𝑄_𝑖=1𝑞=0𝑐𝑄,𝑞𝑋_{𝑘+𝑖+𝑞𝑝}, (5.9) so that we can write 𝜎²𝜏₀,𝑝,𝑄=𝔼𝑑²_𝑝,0,𝜎²_𝑛𝜏₀=1,𝑝,𝑄𝑛^𝑛_𝑄,𝑝^{𝑄,𝑝−1}_𝑘=0𝑑²_𝑝,𝑘. (5.10) This suggests that convergence results, similar to Theorems 4.1 and 4.2, can be achieved also for the MHVAR and its related log-regression estimator.

5.3. The Case of Stationary Processes

In applications, MAVAR and MHVAR are also used in order to estimate the memory parameter of long-range dependent processes. This general case is not included in our analysis (which is restricted to the fractional Brownian motion) and it requires a more involved investigation. To our knowledge, there are no theoretical results along this direction.

Appendix

In this appendix we recall the Wick’s rule for jointly Gaussian random variables and some facts used along the proofs.

Wick’s Rule
Let us consider a family {𝑍_𝑖} of jointly Gaussian random variables with zero mean. The Wick’s rule is a formula that provides an easy way to compute the quantity 𝔼[𝑍_Λ]∶=𝔼[∏_𝑖∈Λ𝑍_𝑖], for any index-set Λ (see, e.g., [19]).

Since the 𝑍_𝑖’s are zero-mean random variables, if Λ has odd cardinality, we trivially get 𝔼[𝑍_Λ]=0. We then assume that |Λ|=2𝑘, for some 𝑘≥1. To recall the Wick’s rule, it is convenient to introduce the following graph representation. To the given index-set Λ we associated a vertex-set 𝑉 indexed by the distinct elements of Λ, and to every vertex 𝑗∈𝑉 we attached as many half-edges as many times the index 𝑗 appears in Λ. In particular there is a biunivocal correspondence between the set of half-edges and Λ, while |𝑉|≤|Λ|. Gluing together two half-edges attached to vertices 𝑖 and 𝑗, we obtain the edge (𝑖,𝑗). Performing this operation recursively over all remaining half-edges, we end up with a graph 𝐺, with vertex set 𝑉(𝐺) and edge-set 𝐸(𝐺). Let 𝒢_Λ denote the set of graphs (possibly not distinguishable) obtained by performing this “gluing procedure” in all possible ways.

With this notation, and for all index-sets Λ with even cardinality, the Wick’s rule for a family {𝑍_𝑖} of jointly centered Gaussian random variables provides the identity 𝔼𝑍_Λ=_{𝐺∈𝒢Λ(𝑖,𝑗)∈𝐸(𝐺)}𝔼𝑍_𝑖𝑍_𝑗. (A.1)

Example A.1. Consider the quantity 𝔼[𝑍²_𝑘𝑍²_𝑗], for 𝑗≠𝑘. By the graphical representation, we take a vertex set 𝑉={𝑗,𝑘} with two half-edges attached to each vertex and perform the gluing operation on the half-edges. We then obtain the following three graphs (identified by their edges): 𝐺₁=(𝑗,𝑘)(𝑗,𝑘), 𝐺₂=(𝑗,𝑘)(𝑗,𝑘),and𝐺₃=(𝑗,𝑗)(𝑘,𝑘). Thus, from the Wick’s rule, we get the identity 𝔼𝑍²_𝑘𝑍²_𝑗=2𝔼𝑍_𝑗𝑍_𝑘²+𝔼𝑍²_𝑗𝔼𝑍²_𝑘, (A.2) that has been used throughout the paper.

Now we recall some facts used in the proof of Theorem 4.1.

Denote by 𝜌[𝐴] the spectral radius of a matrix 𝐴={𝑎_𝑖,𝑗}_{1≤𝑖,𝑗≤𝑛}, then 𝜌[𝐴]≤max^𝑛_{1≤𝑗≤𝑛𝑖=1}𝑎_𝑖,𝑗. (A.3) Moreover the following lemmas hold.

Lemma A.2 (see [13]). Let (𝑉_𝑛) be a sequence of centered Gaussian random vectors and denote by 𝐶_𝑛 the covariance matrix of 𝑉_𝑛. Let (𝐴_𝑛) be a sequence of deterministic symmetric matrices such that (1)lim_𝑛→+∞Var[𝑉^𝑇_𝑛𝐴_𝑛𝑉_𝑛]=𝜆²∈[0,+∞), (2)lim_𝑛→+∞𝜌[𝐴_𝑛]𝜌[𝐶_𝑛]=0. Then 𝑉^𝑇_𝑛𝐴_𝑛𝑉_𝑛−𝔼[𝑉^𝑇_𝑛𝐴_𝑛𝑉_𝑛] converges in distribution to the normal law 𝒩(0,𝜆²).

Lemma A.3 (see [11]). Let 𝑚≥2 be an integer and 𝐶 a 𝑚×𝑚 covariance matrix. Let 𝑟 be an integer such that 1≤𝑟≤𝑚−1. Denote by 𝐶₁ the top left submatrix with size 𝑟×𝑟 and by 𝐶₂ the bottom right submatrix with size (𝑚−𝑟)×(𝑚−𝑟), that is, 𝐶₁=𝐶_{𝑖,𝑗1≤𝑖,𝑗≤𝑟},𝐶₂=𝐶_{𝑖,𝑗𝑟+1≤𝑖,𝑗≤𝑚}. (A.4) Then 𝜌[𝐶]≤𝜌[𝐶₁]+𝜌[𝐶₂].

Acknowledgments

The authors are grateful to Stefano Bregni for having introduced them to this subject, proposing open questions and providing useful insights. They are also grateful to Marco Ferrari for a helpful remark on the first version of this work. This work was partially supported by GNAMPA (2011).