# Identification of Nonstandard Multifractional Brownian Motions under White Noise by Multiscale Local Variations of Its Sample Paths

Authors: Kwang-Il Ahn; Kichun Lee
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The Hurst exponent and variance are two quantities that often characterize real-life, high-frequency observations. Such real-life signals are generally measured under noise environments. We develop a multiscale statistical method for simultaneous estimation of a time-changing Hurst exponentH(t) and a variance parameter C in a multifractional Brownian motion model in the presence of white noise. The method is based on the asymptotic behavior of the local variation of its sample paths which applies to coarse scales of the sample paths. This work provides stable and simultaneous estimators of both parameters when independent white noise is present. We also discuss the accuracy of the simultaneous estimators compared with a few selected methods and the stability of computations with regard to adapted wavelet filters.


---

## Body

## 1. Introduction

Fractional Brownian motion (fBm) has been commonly used to characterize a wide range of complex signals in natural phenomena that exhibit self-similarity and long-range dependence since the pioneering work of Mandelbrot and Van Ness [1]. Examples of such complex signals in time are abundant in medicine, economics, and geoscience, to list a few. The fBm model is characterized by two parameters of the regularity level and the variance level of a signal. The regularity attribute, also called the Hurst exponent, expresses the strength of statistical similarity at many different frequencies, and the variance attribute describes an order of energy magnitude.To model path regularity varying with time, multifractional Brownian motion (mBm) has been proposed as a generalization of fractional Brownian motion (fBm). The theory and applications of both fBm and mBm models have attracted the interests of researchers in numerous problems of, for example, sea level fluctuations [2], currency exchange rates [3], and network traffic [4–6]. To model mBm, Lévy-Véhel and Peltier [7] proposed a mean average approach, and Benassi et al. [8] introduced a spectral approach. Lim and Muniandy [9, 10] also proposed a mBm model based on the fBm defined by the Riemann-Liouville type of fractional integral. The proposed models represent mBm as a Gaussian process W(t) with a covariance function involving Hurst exponent H by a function of time, H(t), and variance parameter C. Specifically, a Gaussian process (W(t))t≥0 is called mBm with Hurst function H(t) and its variance level (scaling factor) C if its covariance function is represented as
(1)𝔼[W(t)W(s)]=C22g(H(t),H(s))×{|t|H(t)+H(s)+|s|H(t)+H(s)-|t-s|H(t)+H(s)}
for H∈ℂη([0,1]); t,s∈[0,1]; g(x,y)=K(2x)K(2y)/K(x+y); and K(α)=Γ(α+1)sin(απ/2)/π, 0<α<2. The process is well defined, or square-integrable, if function H(t) is the Hölderian of order 0<η≤1 on [0,1]. Clearly, the process W(t) is not weakly stationary since the covariance function 𝔼[W(t)W(s)] does not depend on t-s only. From (1), we have 𝔼[W(t)2]=C2tH(t), and consequently, Var[W(1)]=C2. In this sense, C is called the variance level of the process.The time-changing Hurst exponentH(t) characterizes the path regularity of process W at time t since sample paths near t with small H(t), close to 0, are space filling and highly irregular, while paths with large H(t), close to 1, are very smooth. The variance constant C determines the energy level of the process. Alternatively, a spectral representation of mBm is given by
(2)W(t)=CK(2H(t))2∫ℝ‍eitλ-1|λ|H(t)+1/2dB(λ),
where C is a constant scale (variance) parameter and B the standard Brownian [7].Several approaches were proposed to estimate time-changing Hurst exponentH(t) and variance C from sample paths of mBm signals. Benassi et al. [8] investigated estimation of a continuously differentiable H(t) without the direct estimation of C. A local version of quadratic variations was used in several researches to estimate the constant Hurst exponent [11–13]. Recently Fhima et al. [14] adopted the increment ratio statistic method for H(t) estimation only. For an overview of estimating constant H(t), the reader is also referred to Beran [15] including various statistical methods or Bardet and Bertrand [16] concentrating on wavelet approaches. Estimation of both H(t) and variance parameter C has received little attention from the statistics community while C is mostly treated as a nuisance parameter. When a signal is modeled with mBm, the estimation of H(t) can be improved by the accurate estimation of C from covariance structures involving both H(t) and C. For that purpose, the application of a local version of quadratic variations for estimating H(t) and C in mBm was discussed in Coeurjolly [17], in which C was, however, locally estimated in each sample path. Moreover, the existence of noise in mBm signals has not been dealt with in the literature though real-life signals are commonly measured under noise environments.The main objective of this paper is to develop a stable and accurate estimation procedure for unknown parameters(H(t),C) given a path of W(t) in the presence of independent white noise. Previous approaches by Coeurjolly [17] relied on local sample paths in the absence of white noise that resulted in estimators of C sensitive to the sampled paths. It is widely accepted that noise occurs from a variety of sources such as measurement devices.In this paper, we assume that mBm signals are contaminated by a moderate amount of noise. We extend the quadratic variations method to estimateH(t) and C simultaneously for mBm by applying dilated high-pass filters to all sampled paths (all subsample paths from a given sample path) and aggregating all local conditions from the previous filtering step. This method includes filtering all sampled paths with a dilated filter possessing a sufficient number of vanishing moments to capture regularity conditions at associated coarse scales and generating stationary filtered signals. The method further calculates empirical moments of the filtered signals and then estimates H(t) and C simultaneously together with a noise level in a regression setup specified by the empirical moments.This paper is organized as follows. Section2 introduces local variations in a mBm setting, discussing the procedures and justification for the simultaneous estimators of unknown parameters. Section 3 discusses numerical simulations and computational issues with adapted wavelet filters. The appendix presents proofs for the propositions in the preceding sections.


## 2. Multiscale Local Variations of Multifractional Brownian Motion

Let us consider a case in which a discretized sample path(W′) is given by
(3)W′(iN)=W(iN)+σε(iN),i=1,…,N,
in which ε(i/N) is independent white noise and σ is the noise level. Hurst function H(t), generated by W(i/N), is assumed to be Hölderian function of order 0<η≤1 on [0,1]. In addition, noise magnitude σ is assumed to be sufficiently small compared to the variance of mBm. The covariance function of (W′) is
(4)𝔼[W′(t)W′(s)]=C22g(H(t),H(s))×{|t|H(t)+H(s)+|s|H(t)+H(s)-|t-s|H(t)+H(s)}+σ21(t=s),
where 1(A) is an indicator of relation A and g(x,y)=K(2x)K(2y)/K(x+y) as defined above. From the above covariance function, we have 𝔼[W(t)2]=C2t2H(t)+σ2, and consequently, Var[W(1)]=C2+σ2. Noticeably the estimation of C is nontrivial because of the dependence structure from the covariance function; that is to say, the sample variance of a sample path does not lead to the direct expression of C alone but to an expression mixed with all unknown parameters. The entries in (4) generate covariance matrix Σ, which depends on unknown parameters θ:=(H(t),C,σ)∈RN+2. The covariance matrix consists of N(N+1)/2 parameters (due to symmetry) that can be organized into an N(N+1)/2×1 vector Γ(θ). Model (4) is locally identifiable almost everywhere if Jacobian matrix ∂Γ(θ)/∂θ′, which is N(N+1)/2×(N+2), has full column rank [18].In order to weaken the dependence inW′(t) in (3), a differencing filter a of length l+1 and order p>1 (the number of vanishing moments) is applied. Filter a is defined by its taps (a0,…,al) such that
(5)∑q=0l‍aqqi=0,i=0,…,p-1,∑q=0l‍aqqi≠0,i=p.
Let us also introduce a(m) based on filter a, defined by
(6)ai(m)={ai/m,imisaninteger,0,otherwise.
We observe that a(m), the filter a dilated m times, focuses on a resolution at a low frequency, corresponding to a coarse space, as m increases. For m=1, it captures the finest level of detail. For example, a(1)=a by definition, and for a second-order filter a:=(1,-2,1), a(2) becomes (1,0,-2,0,1). Furthermore, we can choose a as high-pass wavelet filters corresponding to orthogonal wavelets such as Daubechies and Symlet wavelets. A detailed discussion of wavelet filters can be found in Daubechies [19] and Vidakovic [20].Let (Va(m)′) be a process consisting of (W′) filtered by a(m), that is,
(7)Va(m)′(jN)=∑q=0ml‍aq(m)W′(j-qN),forj=ml+1,…,N.
For example, when a is (1,-2,1), the filter is of order 2, and (Va′) represents the second-order differences of (W′). The process (Va(m)) is defined similarly with (W) instead of (W′): (Va(m)) is a process consisting of (W) filtered by a(m). The filtering by a(m) breaks the dependence structure between observations. Specifically the process (Va(m)′) is stationary due to the vanishing moment property of filter a(m). To verify it, we need to introduce a sufficiently small neighborhood covering j. Let ν(t) be an index set of a neighborhood of t, defined as
(8)ν(t)={j∈ℤ∣l<j≤N,|jN-t|≤ϵ}
for a parameter ϵ>0. We set ϵ to be a function of N in such a way that ϵ→0, ϵN→∞, and ϵηlog(N)→0 as N→∞. In other words, for a sufficiently large N, the size of one neighbor becomes sufficiently small while maintaining the summation of the sizes of all the neighbors that are sufficiently large. In addition, it is possible to make ϵη converge to zero faster than log(N) grows. Then we derive the covariance of (Va(m)′) as follows.Proposition 1.
Letj1∈ν(t1), j2∈ν(t2). Then, the covariance of Va(m)′(j/N) in (7), 𝔼[Va(m)′(j1/N)Va(m)′(j2/N)], depends on j1-j2 as follows:
(9)𝔼[Va(m)′(j1N)Va(m)′(j2N)]=C2g(H(t1),H(t2))NH(t1)+H(t2)×πa(m),H(t1)/2+H(t2)/2(j1-j2)+σ2∑p-q=j1-j2‍ap(m)aq(m)+𝒪(ϵηlogN),
where
(10)πa(m),h(k)=-12∑p,q=0ml‍ap(m)aq(m)|k-p+q|2h.The above proposition states thatVa(m)′(j/N) is weakly stationary as Gaussian. Particularly for j1=j2=j, as N→∞, it simplifies to
(11)E[Va(m)′(jN)2]=C2(mN)2H(t)πa,H(t)(0)+σ2∑q‍aq2.
Observably the above proposition deals with two pointwise positions, j1 and j2, for two neighborhoods near t1 and t2, respectively. Thus an aggregate behavior of each neighborhood is analyzed via the following setup.Let us define the second empirical moment of the filtered signalVa(m)′ as follows:
(12)S′(t,a(m))=1|ν(t)|∑j∈ν(t)‍Va(m)′(jN)2,fort∈[0,1],
which represents the average squared energy of the a(m)-filtered signal in the neighborhood of t. We notice that S′(t,a(m)) is random because Va(m)′(j/N) is random and its expectation E[S′(t,a(m))] equals that of Va(m)′(j/N)2 because Va(m)′(j/N) is weakly stationary. That is to say,
(13)E[S′(t,a(m))]=E[Va(m)′(jN)2]=C2(mN)2H(t)πa,H(t)(0)+σ2∑q‍aq2.
Now, to relate S′(t,a(m)) to E[S′(t,a(m))] more specifically, we define a statistic V(t,a(m)), called the m-scale local variation of (W), as
(14)V(t,a(m))=1|ν(t)|∑j∈ν(t)‍{Va(m)(j/N)2𝔼[Va(m)(j/N)2]-1},
where |ν(t)| is the cardinal number of ν(t). The statistic V(t,a(m)) captures the amount of deviations of the a(m)-filtered signal from the standard normal distribution near t because Va(m) is normalized by its standard deviation, the square root of 𝔼[Va(m)(j/N)2]. The definition based on the second order can be extended to the kth order Hermite polynomial in the summation of (14): the second Hermite polynomial is defined by t2-1. In this paper, we use local variations based on the second Hermite polynomial as the minimum asymptotic variance estimators, as shown in Coeurjolly [13] for fBm settings. Next, we connect the m-scale local variation V(t,a(m)) with the empirical moment S′(t,a(m)) through the following relationship.Proposition 2.
Letm-scale local variation V(t,a(m)) and the empirical moment S′(t,a(m)) be defined by (14) and (12), respectively, given W′(t) as above and a of order >1. Then
(15)logS′(t,a(m))=log𝔼[S′(t,a(m))]+V(t,a(m))(1+o(1)),asN⟶∞.The proposition connects the empirical momentS′(t,a(m)) and the log of its expectation through m-scale local variation V(t,a(m)). Since the m-scale local variation converges almost surely to 0 and its distribution follows normal distribution asymptotically [17], the above relationship establishes a regression setup. We also note that a filter of an order of at least 2 ensures asymptotic normality for all the values of the function H(t). For a filter of order 1, this convergence is available if and only if 0<suptH(t)<3/4.Next the relationship between the log of the expectation of the empirical momentS′(t,a(m)) and the parameters of interest (H(t), C, and σ) is derived naturally in the light of Proposition 2 and (13). Thus for t1,…,td∈[0,1] we obtain a regression model for logS′(ti,a(m)) as N→∞:
(16)logS′(ti,a(m))~log(C2(mN)2H(ti)πa,H(ti)(0)+σ2∑q‍aq2),∀i,m,
which is nonlinear with respect to H(ti), C, and σ. In particular, when the noise level σ is considered to be zero, the regression model simplifies to, for all i and m,
(17)logS′(ti,a(m))~2logC+2H(ti)log(mN)+log(πa,H(ti)(0))=[12log(mN)]×[2logC+log(πa,H(ti)(0))H(ti)],
which turns out to be linear with respect to H(ti) with intercept 2logC, if log(πa,H(ti)(0)) is negligible. The above regression model possesses a computational advantage though ignoring the presence of noise.Whenσ is nonzero, the following least square estimator of (H,C,σ) is introduced:
(18)(H^,C^,σ^)=argmin(H,C,σ)×∑i=1d‍∑m=1M‍(logS(ti,a(m))-log𝔼[S(ti,a(m))])2.
The computation of the least-square estimator is feasible because, based on (16), for fixed C as Co and σ as σo, the computation of H is separable into each H(ti). In other words, a solution of H(ti) is given by
(19)H^(ti)=argminh×∑m=1M‍((Co2(mN)2hπa,h(0)+σo2∑q‍aq2)logS′(ti,a(m))-log(Co2(mN)2hπa,h(0)+σo2∑q‍aq2))2.
Numerical approaches such as the bisection method can be used for the above procedure, which is nonlinear in h. The bisection method achieves a desired precision level, τ, for H^(ti) with the number of iterations greater than log2τ-1. In other words, 10 iterations, for example, results in precision τ<0.001.


## 3. Simulations and Comparisons

We present here a simulation study of the performance of the approach suggested in this paper, denoted byS-K-var. Simulation is done with the “known truth” of Hurst function H(t) the controlled signal variance and the signal-to-noise (SNR) ratio. Test functions are shown in Figure 1 with the step function for H(t) in Figure 1(a) and the straight-line function in Figure 1(c). Their illustrations of W′(t) are shown in Figures 1(b) and 1(d), correspondingly. For the sake of comparison, we chose several popular methods such as the local spectra slope, which is summarized in Gao [21] and denoted by LSS, and k-variation of variance-uncorrected, denoted by K-var, and the k-variation of variance-corrected version in Coeurjolly [13], denoted by K-var-VC. The average mean squared error (AMSE) was used as a performance measure to capture the difference between true H(t) and estimated H^(t). To simulate a sample path from a fBm on [0,1], we used the method of Wood and Chan [22]. One can simulate a standard mBm W with covariance matrix CH(·) by generating Z~N(0,IN) and estimating W:=CH(·)1/2Z. This method is exact in theory and sufficiently fast for a reasonable sample size N.Tested Hurst functions are shown in (a) step-functionH(t) and (c) straight-line H(t); their illustrations of signals are shown in (b) and (d), correspondingly.
(a)
(b)
(c)
(d)In this section we will use the following notations regarding filters:Diff.i denoting the filter of differences of order i, Db.i denoting a Daubechies wavelet filter of order i, and Sym.i denoting a Symlet wavelet filter of order i. We generate 1,000 series of length N=4096 for step-function H(t) and N=1024 for straight-line function H(t). A simple difference filter [1-21] (Diff.2) was used for straight-line H(t), and Db.6 was used for step-function H(t). For the local spectra slope of LSS, the length of the subsignal was set to be 512, which is sufficient for its numerical stability, and the two levels, by which spectral slopes are calculated, were 3 and 7. The size of a neighborhood of t, ν(t) in (8), is set to be 50 for S-K-var, K-var, and K-var-VC.Illustrations of the estimators under no noise are shown in Figure2. Estimation by K-var-VC most accurately follows true H(t) among the tested methods. Estimation results by K-var, considering no scale parameter C, notably deviate from true H(t). We note that the distance between K-var estimation and true H(t) relates to C. Estimators by LSS are bumpy because it assumes that subsignals during its computation follow fBm without considering the variability of H(t). We also observe that K-var-VC is more unstable than S-K-var.Illustrations of the estimators: in panel (a), varianceC was 2; in panel (b), variance C was 4; the best match to the true H(t) was the method of S-K-var among the four.
(a)
(b)Regarding the estimation ofC, the comparison between K-var-VC and S-K-var is shown in Figure 3, in which empirical confidence intervals for true C=2 are shown with the upper panels for no noise and the lower panels for SNR 10. We sampled 1000 series of W′(t) with C=2 and straight-line H(t) under white noise of SNR 10. Consistently, the estimation results by S-K-var at the right panels are more accurate, and their confidence intervals are sharper than those by K-var-VC. We also note that the confidence intervals by S-K-var in Figures 3(b) and 3(d) are constant in time since S-K-var employs a global constant in regression model (16). A noise level of SNR 10 heavily worsened the estimation results by K-var-VC while those by S-K-var yielded a slight increase in the confidence intervals. Accurate estimation of variance level C by K-var-VC leads to accurate estimation of H(t), which will be demonstrated in the following tests.Estimation ofC by K-var-VC with empirical 95% confidence intervals in blue is shown in (a) for no noise and (c) for SNR 10 when true C=2. Similarly, estimation of C by S-K-var in red is shown in (b) for no noise and (d) for SNR 10. S-K-var yields more stable and shorter confidence intervals.
(a)
K-var-VC, no noise
(b)
S-K-var, no noise
(c)
K-var-VC, SNR 10
(d)
S-K-var, SNR 10We comparedS-K-var with K-var-VC and K-var in terms of AMSE in various settings. The method LSS was dropped due to obvious poor performance as is shown in Figure 2. We varied the levels of variance C from 1 to 4 and 10 under SNR levels of 10, 30, and the infinity. The number of sample path W′(t) was 1000 in each setting, and AMSE was computed for each of the methods. The results are shown in Table 1 for each step-function and straight-line H(t). We observe that our proposed method S-K-var consistently outperforms the other methods except for only a few settings. Overall, there was little difference between K-var-VC and K-var in performance. This experimental result is not surprising since S-K-var reflects the existence of white noise and globally includes variance constant C.Table 1
Average mean squared error for theS-K-var, K-var-VC, and K-var methods in various settings of variance constant C and SNR levels based on 1000 sample paths of W′(t) for each of step-function and straight-line H(t). The asterisk marks show that S-K-var outperformed the other methods.
Parameters
Average mean squared error
Step-functionH(t)
Straight-lineH(t)
S-K-var
K-var-VC
K-var
S-K-var
K-var-VC
K-var
C
=
1, SNR 10
31.54*
178.6
107.1
84.62
215.5
81.31
30
22.85*
89.64
64.59
25.48
147.8
23.53
∞
9.940*
26.31
48.16
.8242*
9.508
.9376
C
=
4, SNR 10
15.56*
212.5
82.29
46.16*
176.9
175.2
30
4.406*
152.2
39.91
36.11*
92.07
120.2
∞
.9706*
10.00
8.173
10.04*
24.73
93.21
C
=
10, SNR 10
61.80*
178.2
281.5
89.94*
216.5
223.9
30
55.97*
87.35
237.0
42.21*
150.9
144.6
∞
33.86*
26.08
193.1
1.425*
9.288
77.32The effects of adapted filters are summarized in Figure4. The experiments were done with straight-line H(t), variance C=2, and SNR = 7. We observe that the performance of S-K-var on the estimation of C does not change much depending on the filter it uses. However, we mention that the variance of AMSEs tends to increase according to the filter size.Performance according to different filters with straight-lineH(t); in panel (a), box plots for estimations of C=2; in (b), box plots for AMSE of H(t). Filter numbers represent the following: 1, Diff.1; 2, Diff.2; 3, Diff.3; 4, Diff.4; 5, Sym.4; 6, Sym.6; 7, Sym.8; 8, Db.2; 9, Db.3; and as well as 10, Db.4.
(a)
(b)


## 4. Conclusion

To conclude, we proposed the joint estimators of the time-changing Hurst exponentH(t) and its variance coefficient C for mBm under white noise. The proposed method is based on filtering sampled paths with dilated high-pass filters to derive regularity conditions at associated scales. The second empirical moment, average squared energy, of the filtered signals near a time position is connected to the theoretical expectation and used to establish a regression setup through the asymptotic distribution of multiscale local variation statistics. The effectiveness of the approach was verified through numerical experiments that compared it with that of several other approaches. Simulation results show that the proposed approach yields more precise and stable estimation of Hurst exponents and variance constants under noiseless or noised conditions.



---

*Source: 794130-2013-10-31.xml*