Testing of Multifractional Brownian Motion

Fractional Brownian motion (FBM) is a generalization of the classical Brownian motion. Most of its statistical properties are characterized by the self-similarity (Hurst) index 0<H<1. In nature one often observes changes in the dynamics of a system over time. For example, this is true in single-particle tracking experiments where a transient behavior is revealed. The stationarity of increments of FBM restricts substantially its applicability to model such phenomena. Several generalizations of FBM have been proposed in the literature. One of these is called multifractional Brownian motion (MFBM) where the Hurst index becomes a function of time. In this paper, we introduce a rigorous statistical test on MFBM based on its covariance function. We consider three examples of the functions of the Hurst parameter: linear, logistic, and periodic. We study the power of the test for alternatives being MFBMs with different linear, logistic, and periodic Hurst exponent functions by utilizing Monte Carlo simulations. We also analyze mean-squared displacement (MSD) for the three cases of MFBM by comparing the ensemble average MSD and ensemble average time average MSD, which is related to the notion of ergodicity breaking. We believe that the presented results will be helpful in the analysis of various anomalous diffusion phenomena.


Introduction
Over the last decades, massive advances in single-particle tracking (SPT), partially based on superresolution microscopy of fluorescently tagged tracers, or fluorescence correlation spectroscopy allow experimentalists to obtain insight into the motion of submicron tracer particles or even single molecules in complex environments, such as living biological cells, down to nanometer precision and at submillisecond time resolution [1,2].
Such transient behavior has, for instance, been observed for green fluorescent proteins in cells or for the motion of lipid molecules in protein-crowded bilayer membranes [25,32]. Fractional Brownian motion (FBM) introduced by Kolmogorov in 1940 and rediscovered by Mandelbrot and van Ness in 1968 [33][34][35] is a generalization of the classical Brownian motion (BM). Most of its statistical properties are characterized by the self-similarity (Hurst) index 0 < H < 1. FBM is H-self-similar, namely for every c > 0 we have B H (ct) D = c H B H (t) in the sense of all finite dimensional distributions, and has stationary increments. It is the only Gaussian process satisfying these properties. FBM is the overdamped description for viscoelastic motion and thus intimately connected to the fractional Langevin processes, an attractive framework for many physical systems [36], for instance, of lipid molecules in bilayer membranes [22,25,37]. The second moment of FBM reads EB 2 H (t) = σ 2 t 2H , where EB 2 H (1) = σ 2 > 0. As a consequence, for H < 1/2 we obtain subdiffusive dynamics with persistent motion, whereas for H > 1/2, the process is superdiffusive and antipersistent. Since FBM is the classical model for power-law dependence a number of statistical tests have been already introduced for this process in the literature. Let us mention here the tests based on the autocovariance function (ACVF), MSD, and detrending moving average statistics [38][39][40].
FBM has stationary increments that do not allow us to model processes whose regularity of paths and "memory depth" change in time [41]. Several generalizations of FBM have been proposed recently. One of these, called multifractional Brownian motion (MFBM), was proposed by Peltier and Véhel [42] with time-varying Hurst exponent H(t) which is a Hölder function. The variance at time t of the MFBM B H(t) (t) is given by Var(B H(t) (t)) = σ 2 t 2H(t) [43]. The time-varying Hurst exponent H(t) characterizes the path regularity of the process at time t: 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 σ 2 determines the "energy level" of the process. This natural extension of FBM results in some loss of some of FBMs basic properties, in particular, the increments of MFBM are non-stationary and the process is no longer self-similar.
Other, similar generalizations are limited to a piecewise constant H [44] but, what is important from a data analysis point of view, is that they lead to continuous Gaussian processes with stationary increments. Let us also mention an idea involving an appropriate class of covariance functions. Ryvkina [45] uses such covariance functions to define Gaussian processes to extend FBM and MFBM to a class of fractional Brownian motions with a variable Hurst parameter parameterized by a set of all measurable functions with values in (1/2, 1), and different from MFBMs. However, from a biological data point of view, such a range for H values is not practical since it only corresponds to a superdiffusive (long-range dependent) case.
MFBMs have become popular as flexible models in describing real-life signals of high-frequency features in geoscience, microeconomics, and turbulence, to name a few [43]. They are closely related to the notion of transient diffusion dynamics observed in biological experiments. The article is structured as follows. In Section 2, the MFBM is defined and its basic properties are presented. We also recall formulas for the ensemble average MSD, time average MSD, and present three Hurst exponent functions that will be analyzed in the sequel. In Section 3, the main results are presented. We introduce a statistical test on MFBM based on its ACVF which is presented as a quadratic form. Next, the power of the test is studied for the three cases corresponding to different Hurst exponent functions. We show the areas where the test is very strong in distinguishing between the processes and the cases when it fails in this respect. Finally, Section 4 summarizes and concludes our work.

Model and Methods
Let us start with a definition of the MFBM.
In general, MFBM has non-stationary increments. Its increment process Y(t) where Corr is the correlation function, i.e., Corr( [46]. Furthermore, the function H(t) can be pointwise interpreted as a local self-similarity parameter, i.e., where B H is a fractional Brownian motion with index H ≡ H(u), s(u) is a scaling function [47] and the convergence is on the space of continuous functions endowed with the topology of the uniform convergence on compact sets.

Mean-Squared Displacement
Let us now recall different estimators of the MSD for a sample of n trajectories X 1 , X 2 , . . . X n , each with N observations, that is, X i consists of X i (t 1 ), X i (t 2 ), . . . , X i (t N ) equally spaced in time, i = 1, 2, . . . , n. Ensemble average MSD (EAMSD) is defined as follows: where m = 1, . . . , N and ∆t = t 2 − t 1 . Time average MSD (TAMSD) is defined for each trajectory as Finally, we consider ensemble and time average MSD (EATAMSD) which is an average of TAMSDs: Physicists often observe systems where EA and EATAMSD are different. Such behavior is called weak ergodicity breaking. In mathematics, the notion of ergodicity is restricted to stationary processes. Since the increments of MFBM lack stationarity, when analyzing the results, one has to be exceedingly meticulous.

Three Cases of the Hurst Exponent Function
Following [48], we consider three basic families of the function H(t), namely for some time horizon T > 0. Furthermore, in the sequel we consider only case with parameter σ 2 = 1 in (2). Such functions are continuous and as a consequence satisfy the Hölder condition, so in order to MFBM be properly defined we only require H(t) ∈ (0, 1), for all t ∈ [0, T] [42]. Such choice of considered functions can be interpreted as follows. In the linear case, MFBM can switch steadily from short-to long-range dependence or vice versa, whereas in the logistic case such change is quite rapid and it happens between two levels. The latter case closely resembles instantaneous change in dependence or jump-type regime switching (such cases would lead to non-Hölder function). An alternative function to the logistic which is also considered is the literature is the arctan function [49]. Finally, the periodic case represents a situation where such changes are gradual and repetitive.
In the paper, we focus on the following special cases with specified parameters: We choose those specific parameters so that all of the cases have a similar "average" behavior of the function H(t), i.e., its mean is close to 0. 45 is directly related to the variance of the model at time τ, i.e., EAMSD(τ) = Var(X(τ)), thus, from (2), it should behave like τ 2H(τ) . For the linear case, see Figure 1, since H (1) (t) increases steadily from 0.3 to 0.6 we can see trajectories exhibit more variability for bigger times. This is also related to the EAMSD dynamics.
In addition, such a model exhibits weak ergodicity breaking behavior (i.e., lack of equality between EATAMSD and EAMSD), which can be inferred from the bottom panel.
Next, for the logistic case, see Figure 2, H (2) (t) increases quite rapidly from 0.3 to 0.6 near t = 500. As a consequence, we can see a switch in the behavior of simulated trajectories: for times t < 500 they exhibit far less variability than for t > 500. Intuitively, for times t < 500 trajectories locally exhibit short-range dependence, whereas for t > 500 they locally exhibit long-range dependence. Again, we can see weak ergodicity breaking on the bottom panel.   Finally, for the periodic case, see Figure 3, H (3) (t) varies between 0.3 and 0.6. We can clearly see two different regimes of behavior: for times when H (3) (t) is bigger, trajectories are smoother and generally have larger values, in contrast to times when H (3) (t) is smaller. Despite lack of stationarity, here, EAMSD almost always lies in the confidence region of EATAMSD, which could suggest there is no weak ergodicity breaking.

Results
In applications, it is crucial to be able to check whether a stochastic model describes empirical data well. Despite dedicated identification methods for the MFBM [50][51][52][53], to the best of the authors' knowledge, there is no rigorous statistical test designed for such process. Here, we propose an approach using a simple test statistic which also contains useful information about the process itself.

Test
For the testing purposes, we follow an approach based on the ACVF which was introduced by Balcerek and Burnecki [38]. ACVF is a very popular statistic and it is also one of the simplest quadratic forms. For a random sample X N = {X(1), X(2), . . . , X(N)} and τ ∈ {1, 2, · · · , N − 1}, it is defined as follows: Here, we only consider a version of ACVF without subtracting the sample mean as it does not influence performance of tests based on this statistic for a centered process [38] and it makes the formulas much simpler.
Let us now introduce a matrix and I is the indicator. To summarize, the matrix A(τ) is either diagonal (for τ = 0) with elements 1 N on diagonal or Toeplitz, with only two nonzero subdiagonals (starting at (1 + τ)th row and (1 + τ)th column) with elements 1 2 1 N−τ . The statistic ACVF N can be now expressed as a quadratic form (as shown in [38]) as a generalized χ 2 distribution, that is where Z i s are i.i.d standard normal variables (so Z 2 i has a χ 2 1 distribution) and λ k (τ) are eigenvalues of the matrix Σ N (τ) = Σ 1/2 A(τ)Σ 1/2 with Σ being the (theoretical) autocovariance matrix of our trajectory X N . It is important to note that this result is true regardless of whether the considered model is stationary or not.
Let us now formulate a test for checking whether a random sample X N comes from the MFBM with function H : [0, T] → (0, 1), where T is the time horizon: We will use ACVF N as a test statistic with its distribution given by Equation (9) to calculate critical regions of such test for a given significance level. Naturally, eigenvalues λ i (τ) depend on the matrix A(τ) as well as on the matrix Σ. Elements of Σ are given by ACVF (2) and are calculated using the function H(t) from the null hypothesis.

Three Power Case Studies
The power of the test is the probability to reject the null hypothesis when the alternative is true. The power is an important characteristic of any statistical test. We consider the following null hypotheses, which correspond to the examples presented in Figures 1-3. In our studies, for all considered cases, we calculate the power of the test by using Monte Carlo simulations. We assume that the significance level is equal to 5%. In our Monte Carlo simulations we consider the time horizon T = 1000 and equally spaced time points t = 1, 2, . . . , 1000. For each set of parameters from the alternative hypothesis, we simulate n = 1000 trajectories, calculate test statistic (7), and check if the null hypothesis is rejected at 5% significance level. In the test statistic, we consider only τ = 1 since other choices of τ lead to worse results. Finally, we estimate the power of this test for each considered case by calculating the fraction of rejected null hypotheses.
We present the results in the form of power functions with arguments being the parameters of the function H(t) from the alternative hypothesis. For all of the cases, we considered the alternative coming from the same family of functions as the function H(t) in the null hypothesis, i.e., linear alternative for linear null, etc.
First, let us consider testing MFBM with the linear Hurst exponent function. We can see the power function related to that case in Figure 4. The left panel presents the power function with respect to parameters a and b from the alternative hypothesis H 1 : H(t) = at + b, the right panel the corresponding heat map. We can see "layers" (regions on the heat map with the same color) for which our test has a very similar power. For example, the deep blue region on the right panel corresponds to the processes indistinguishable from the null hypothesis process. We believe that the shape of the region is related to the construction of our test, namely our test statistic ACVF N takes into account all addends X(t)X(t + τ) with the same weight, thus it is not that relevant whether t is big or not. In the case of MFBM, for which neither the process nor its increments are stationary, it might be an important factor. As a consequence, we can see that the test has a difficulty in distinguishing between MFBMs with increasing and decreasing Hurst exponent functions if their means are similar. However, this conjecture is not very precise, namely the mean case H ≡ 0.45, which matches the alternative hypothesis with a = 0, b = 0.45, yields a much higher power of the test than the significance level. We also note that for parameters from the null hypothesis: a = 0.3, b = 0.3 power of such test is approximately equal to 5%, which is the assumed significance level. On the heat map, white regions represent areas for which the process MFBM is not well-defined, i.e., H(t) / ∈ (0, 1) for some t ∈ [0, 1000]. Let us now consider the second case, that is MFBM with the logistic function H(t). We can observe the power function in Figure 5. Left panel presents the power function with respect to parameters b and c from the alternative hypothesis H 1 : the right panel the corresponding heat map. The parameter b is related to the local self-similarity parameter for t < 500, and c for t > 500.
Similarly to the case with the linear function null hypothesis, here we can observe "layers" in which parameters b and c are almost symmetric (e.g., the null hypothesis case where b = 0.3 and c = 0.6 is closely related to the case b = 0.6 and c = 0.3). Moreover, we can see that the power function seems to be quite high in the cases when a tested sample has the b parameter close to the value 0.3 from the null hypothesis, but c is far from 0.6, or vice versa.  Lastly, let us consider the case of MFBM with the periodic function H(t). We can observe the power function in Figure 6. The left panel presents the power function with respect to parameters a and b from the alternative hypothesis H 1 : H(t) = a sin 4π t T + b, the right panel the corresponding heat map. Parameter b is related to the "mean" behavior of the function H(t), whereas the parameter a corresponds to its amplitude. Again, similarly to the two previous null hypotheses, we can observe "layers" of similar power values. Those layers are symmetric with respect to a = 0. This means that for alternatives with opposite parameters a the test seems to return the same power. Let us note that this is not intuitive, namely such opposite as are related to completely different local behaviors of the self-similarity parameter. On the heat map, white regions represent areas for which process MFBM is not well defined, i.e., H(t) / ∈ (0, 1) for some t ∈ [0, 1000]. Finally, we would like to emphasize that the introduced test requires MFBM parameters to be fixed (we test if the data follow MFBM with fixed parameters). In practice, when analyzing empirical data the parameters are often estimated. In the literature, methods for estimation of the Hurst exponent function H(t) in the MFBM framework have been already introduced [50][51][52] and later combined to improve both the goodness of fit and the computational speed of the algorithm [53].

Discussion and Conclusions
For power-law anomalous diffusion of the form (1) with constant anomalous diffusion exponent α a number of models exist, including continuous-time random walks, fractional Langevin equation motion, FBM, or scaled Brownian motion [8]. These models all have different physical properties such as the PDF or their ergodic and aging properties [8].
In this paper, we concentrated on MFBM which is a generalization of FBM for Hölder continuous functions H(t) that allows the Hurst exponent to vary in time. The time-varying Hurst exponent has an impact on both the statistical properties of the process and trajectory characteristics. MFBM helps to model phenomena whose regularity of paths and anomalous diffusion exponent change in time. The process has no longer stationary increments and it is not self-similar but the variance scales in a natural way as t 2H(t) .
Following the idea of testing FBM based on the ACVF statistic [38], in this paper, we introduced a rigorous statistical test on MFBM with the ACVF statistic presented as a quadratic form. We derived the distribution of the statistic which is the generalized χ 2 . In order to study the efficiency of the test, we took into consideration three possible classes of the Hurst exponent function, namely linear, logistic, and periodic. For those cases, we conducted power studies with the help of Monte Carlo simulations. As alternatives, we considered MFBMs within the same class of the Hurst exponent function but with different parameters.
We found ranges of the parameters where the test is more sensitive to differences and ranges where it fails to distinguish between the models. It appears that for the linear Hurst exponent function the test is most sensitive to changes in the mean of the function. If the means are similar then the test often fails, even if the functions have completely different patterns, namely, one is increasing and the other decreasing. The latter observation may sound like a serious objection for using the test, but, in practice, an experimentalist knows whether the anomalous diffusion exponent increases or decreases in time. In the logistic case, the situation is different, namely the mean does not matter much as for the linear case. Now, the test is most sensitive to deviations from the true values of the two levels (c and b) with the exception that replacing c with b does not change the power of the test (so, again, it does not matter if the function increases or decreases). For the periodic case, we have again a different situation. The test is sensitive to the changes of the amplitude of the sine function and the value of the free term but it does not detect a sign of the parameter related to the magnitude.
Finally, we note that we checked the behavior of the test for other sets of parameters of the null hypotheses and different sample lengths, and the conclusions were similar. We only found that the range of possible H values from the null hypothesis has an influence on the width of the acceptance regions (the wider the range the wider the acceptance region, which is reasonable). We present some of the additional tests' power simulation studies in Appendix A. Figure A1 presents different cases of the null hypothesis for the linear case, Figure A2 for the logistic case, and Figure A3 for the periodic case. Tests for the linear case were performed for length N = 1000, whereas logistic and periodic cases were studied for length N = 200.
In sum, we introduced a rigorous statistical test for MFBM based on ACVF statistic presented as a quadratic form. We highlighted the weak and strong points of the test. Improving the efficiency of the test will be a subject of our future studies. We believe that the obtained results can help to understand the mechanisms underlying various anomalous diffusion phenomena.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analysis, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A
In Figures A1-A3, we present power functions of tests related to different null hypotheses. In Figure A1, we consider  In Figure A2, we consider  In Figure A3, we consider H 0 : H(t) = 0.1 sin 4π t 1000 + 0.8 (top panel), so a case in which function H varies periodically between 0.7 and 0.9, i.e., in the strong superdiffusion regime; H 0 : H(t) = 0.2 sin 4π t 1000 + 0.7 (middle panel), so a case in which function H varies periodically between 0.5 and 0.9, i.e., in the superdiffusion regime; and H 0 : H(t) = 0.5 sin 4π t 1000 + 0.4 (bottom panel), so a case in which function H varies periodically between 0.1 and 0.9, i.e., in the whole spectrum of anomalous diffusion.