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

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-changingHurst exponentH(t) and a variance parameterC in amultifractional Brownianmotionmodel in the presence of white noise.Themethod 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.


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 longrange 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.
The time-changing Hurst exponent () characterizes the path regularity of process  at time  since sample paths near  with small (), close to 0, are space filling and highly irregular, while paths with large (), close to 1, are very smooth.The variance constant  determines the energy level of the process.Alternatively, a spectral representation of mBm is given by where  is a constant scale (variance) parameter and  the standard Brownian [7].Several approaches were proposed to estimate timechanging Hurst exponent () and variance  from sample paths of mBm signals.Benassi et al. [8] investigated estimation of a continuously differentiable () without the direct estimation of .A local version of quadratic variations was used in several researches to estimate the constant Hurst exponent [11][12][13].Recently Fhima et al. [14] adopted the increment ratio statistic method for () estimation only.For an overview of estimating constant (), the reader is also referred to Beran [15] including various statistical methods or Bardet and Bertrand [16] concentrating on wavelet approaches.Estimation of both () and variance parameter  has received little attention from the statistics community while  is mostly treated as a nuisance parameter.When a signal is modeled with mBm, the estimation of () can be improved by the accurate estimation of  from covariance structures involving both () and .For that purpose, the application of a local version of quadratic variations for estimating () and  in mBm was discussed in Coeurjolly [17], in which  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 ((), ) given a path of () 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  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 estimate () and  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 () and  simultaneously together with a noise level in a regression setup specified by the empirical moments.
This paper is organized as follows.Section 2 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.

Multiscale Local Variations of Multifractional Brownian Motion
Let us consider a case in which a discretized sample path (W  ) is given by in which (/) is independent white noise and  is the noise level.Hurst function (), generated by (/), 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 where 1() is an indicator of relation  and (, ) = √(2)(2)/( + ) as defined above.From the above covariance function, we have E[() 2 ] =  2  2() +  2 , and consequently, Var[(1)] =  2 +  2 .Noticeably the estimation of  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  alone but to an expression mixed with all unknown parameters.The entries in (4) generate covariance matrix Σ, which depends on unknown parameters  := ((), , ) ∈  +2 .The covariance matrix consists of (+ 1)/2 parameters (due to symmetry) that can be organized into an ( + 1)/2 × 1 vector Γ().Model ( 4) is locally identifiable almost everywhere if Jacobian matrix Γ()/  , which is ( + 1)/2 × ( + 2), has full column rank [18].In order to weaken the dependence in   () in (3), a differencing filter a of length  + 1 and order  > 1 (the number of vanishing moments) is applied.Filter a is defined by its taps ( 0 , . . .,   ) such that Let us also introduce a () based on filter a, defined by We observe that a () , the filter a dilated  times, focuses on a resolution at a low frequency, corresponding to a coarse space, as  increases.For  = 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 (V  a () ) be a process consisting of (W  ) filtered by a () , that is, For example, when a is (1, −2, 1), the filter is of order 2, and (  a ) represents the second-order differences of (  ).The process (V a () ) is defined similarly with (W) instead of (W  ): (V a () ) is a process consisting of (W) filtered by a () .The filtering by a () breaks the dependence structure between observations.Specifically the process (V  a () ) is stationary due to the vanishing moment property of filter a () .To verify it, we need to introduce a sufficiently small neighborhood covering .Let ]() be an index set of a neighborhood of , defined as for a parameter  > 0. We set  to be a function of  in such a way that  → 0,  → ∞, and   log() → 0 as  → ∞.In other words, for a sufficiently large , 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() grows.Then we derive the covariance of (V  a () ) as follows. where The above proposition states that   a () (/) is weakly stationary as Gaussian.Particularly for Observably the above proposition deals with two pointwise positions,  1 and  2 , for two neighborhoods near  1 and  2 , respectively.Thus an aggregate behavior of each neighborhood is analyzed via the following setup.
Let us define the second empirical moment of the filtered signal   a () as follows: which represents the average squared energy of the a ()filtered signal in the neighborhood of .We notice that   (, a () ) is random because   a () (/) is random and its expectation [  (, a () )] equals that of   a () (/) 2 because   a () (/) is weakly stationary.That is to say, Now, to relate   (, a () ) to [  (, a () )] more specifically, we define a statistic (, a () ), called the -scale local variation of (W), as where |]()| is the cardinal number of ]().The statistic (, a () ) captures the amount of deviations of the a ()filtered signal from the standard normal distribution near  because  a () is normalized by its standard deviation, the square root of E[ a () (/) 2 ].The definition based on the second order can be extended to the th order Hermite polynomial in the summation of ( 14): the second Hermite polynomial is defined by  2 − 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 scale local variation (, a () ) with the empirical moment   (, a () ) through the following relationship.
Mathematical Problems in Engineering Proposition 2. Let -scale local variation (, a () ) and the empirical moment   (, a () ) be defined by ( 14) and (12), respectively, given   () as above and a of order > 1.Then The proposition connects the empirical moment   (, a () ) and the log of its expectation through -scale local variation (, a () ).Since the -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 ().For a filter of order 1, this convergence is available if and only if 0 < sup  () < 3/4.
Next the relationship between the log of the expectation of the empirical moment   (, a () ) and the parameters of interest ((), , and ) is derived naturally in the light of Proposition 2 and (13).Thus for  1 , . . .,   ∈ [0, 1] we obtain a regression model for log   (  , a () ) as  → ∞: which is nonlinear with respect to (  ), , and .In particular, when the noise level  is considered to be zero, the regression model simplifies to, for all  and , log   (  , a () ) ∼ 2 log  + 2 (  ) log (   ) + log ( a,(  ) (0)) Figure 2: Illustrations of the estimators: in panel (a), variance  was 2; in panel (b), variance  was 4; the best match to the true () was the method of S-K-var among the four.
which turns out to be linear with respect to (  ) with intercept 2 log , if log( a,(  ) (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, , ) is introduced: The computation of the least-square estimator is feasible because, based on ( 16), for fixed  as   and  as   , the computation of H is separable into each (  ).In other words, a solution of (  ) is given by Numerical approaches such as the bisection method can be used for the above procedure, which is nonlinear in ℎ.The bisection method achieves a desired precision level, , for Ĥ(  ) with the number of iterations greater than log 2  −1 .In other words, 10 iterations, for example, results in precision  < 0.001.

Simulations and Comparisons
We present here a simulation study of the performance of the approach suggested in this paper, denoted by S-K-var.Simulation is done with the "known truth" of Hurst function () the controlled signal variance and the signal-to-noise (SNR) ratio.Test functions are shown in Figure 1 with the step function for () in Figure 1(a) and the straight-line function in Figure 1(c).Their illustrations of   () 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 -variation of variance-uncorrected, denoted by K-var, and the -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 () and estimated Ĥ().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  with covariance matrix  (⋅) by generating  ∼ (0,   ) and estimating  :=  1/2 (⋅) .This method is exact in theory and sufficiently fast for a reasonable sample size .
In this section we will use the following notations regarding filters: Diff.i denoting the filter of differences of order ,  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 , ]() 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 Figure 2. Estimation by K-var-VC most accurately follows true () among the tested methods.Estimation results by K-var, considering no scale parameter , notably deviate from true ().We note that the distance between K-var estimation and true () relates to .Estimators by LSS are bumpy because it assumes that subsignals during its computation follow fBm without considering the variability of ().We also observe that K-var-VC is more unstable than S-K-var.
Regarding the estimation of , the comparison between K-var-VC and S-K-var is shown in Figure 3, in which empirical confidence intervals for true  = 2 are shown with the upper panels for no noise and the lower panels for SNR 10.We sampled 1000 series of   () with  = 2 and straight-line () under white noise of SNR 10.Consistently, the estimation results by S-K-var at the right  Table 1: Average mean squared error for the S-K-var, K-var-VC, and K-var methods in various settings of variance constant  and SNR levels based on 1000 sample paths of   () for each of step-function and straight-line ().The asterisk marks show that S-K-var outperformed the other methods.

Parameters
Average mean squared error Step-function () 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  by K-var-VC leads to accurate estimation of (), which will be demonstrated in the following tests.We compared S-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  from 1 to 4 and 10 under SNR levels of 10, 30, and the infinity.The number of sample path   () 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 straightline ().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 .
The effects of adapted filters are summarized in Figure 4.The experiments were done with straight-line (), variance  = 2, and SNR = 7.We observe that the performance of S-K-var on the estimation of  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.

Conclusion
To conclude, we proposed the joint estimators of the timechanging Hurst exponent () and its variance coefficient  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.

Figure 1 :
Figure 1: Tested Hurst functions are shown in (a) step-function () and (c) straight-line (); their illustrations of signals are shown in (b) and (d), correspondingly.

Figure 3 :
Figure 3: Estimation of  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  = 2. Similarly, estimation of  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.