Herwig Wendt
email: herwig.wendt@ens-lyon.fr

Patrice Abry
email: patrice.abry@ens-lyon.fr

BOOTSTRAP TESTS FOR THE TIME CONSTANCY OF MULTIFRACTAL ATTRIBUTES

Keywords: Multifractal analysis, Non parametric bootstrap, Stationarity test, Time constancy test, Wavelet Leaders

On open and controversial issue in empirical data analysis is to decide whether scaling and multifractal properties observed in empirical data actually exist, or whether they are induced by intricate non stationarities. To contribute to answering this question, we propose a procedure aiming at testing the constancy along time of multifractal attributes estimated over adjacent non overlapping time windows. The procedure is based on non parametric bootstrap resampling and on wavelet Leader estimations for the multifractal parameters. It is shown, by means of numerical simulations on synthetic multifractal processes, that the proposed procedure is reliable and powerful for discriminating true scaling behavior against non stationarities. We end up with a practical procedure that can be applied to a single finite length observation of data with unknown statistical properties.

INTRODUCTION

Scale invariance. Scale invariance, or scaling, is a property that has been extensively observed in empirical data of very different nature, such as turbulence, network traffic or biomedical signals. In essence, scale invariance can be defined as the fact that moments of order q of some multiresolution quantities T X (a, t) (e.g., wavelet coefficients) are characterized by power law behaviors with respect to the analysis scale a. The exponents of such power law behaviors, labeled ζ(q), are referred to as the scaling exponents and are often involved in detection, identification or classification tasks. Empirical multifractal analysis, whose goal is to analyze scale invariance in real data and to provide measurements of the scaling attributes, has been successfully used in various applications and is currently becoming a standard tool in empirical time series analysis. Scale invariance and non stationarity. In the practical multifractal analysis of empirical data, there has been, and there still is, an important controversy: Do scaling actually exist in data, or are they rather the consequence of non stationarities that conspire to mimic scaling behavior? To contribute to answering this question, let us first clarify the issue. There exist two major classes of stochastic processes used to model scale invariance: Self-similar and multifractal processes. Both classes consist of non stationary processes, and there is hence no contradiction between scale invariance and non stationarity in that respect. The controversy between scale invariance and non stationarity can in fact be cast in the following three categories: First, scaling actually exist but a smooth trend (in the mean or variance, for example), hence a non stationarity, is superimposed to the data and is likely to impair the analysis; Second, scaling exist in data but their parameters exhibit some form of variability with respect to time, for instance due to a change in experimental conditions; Third, scaling are not present in data but a strong non stationary variability is confused with a scaling property. The first category has been addressed in a number of research papers (cf. [START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF] and the references therein) and will not be further considered here. The second an third categories are much more involved as a non stationary variability can correspond to many different realities. Nevertheless, their detection is of crucial practical importance, since the blind analysis of such time series is likely to produce misleading interpretations of scaling. Constancy along time of scaling exponents. The discrimination of true scaling against various forms of non stationary variability can be addressed with the following heuristic: When data posses true scaling properties, scaling exponents estimated over the entire time series or over non overlapping adjacent windowed time series are statistically consistent. Conversely, when scaling exponents obtained over non overlapping adjacent subsets of the data are not statistically consistent, this can only be the signature of some form of non stationarity, whatever its precise and a priori unknown nature. Therefore, the issue of testing scale invariance against non stationarity can be meaningfully recast into a test of time constancy of scaling exponents estimated over adjacent non overlapping subsets of the analyzed time series. This is precisely the intuition developed in [START_REF] Veitch | A statistical test for the time constancy of scaling exponents[END_REF], where a time constancy test is developed for the (wavelet coefficient based) estimation of the Hurst parameter of Gaussian self similar stationary increment (H-sssi) processes. Goals of the present contribution: bootstrap and multifractal. Our goal is to extend such an approach to testing for the time constancy of scaling exponents of multifractal process. This implies two major changes in methodology: First, the description of multifractal processes requires a whole collection of attributes, related to positive and negative statistical orders q, while the Gaussian nature of the self similar processes analyzed in [START_REF] Veitch | A statistical test for the time constancy of scaling exponents[END_REF] allows to concentrate on the second order (q = 2) only. Second, estimations are no longer based on wavelet coefficients but on wavelet Leaders [START_REF] Wendt | Bootstrap for empirical multifractal analysis with application to hydrodynamic turbulences[END_REF]. While the former are obtained by a linear transform of the data, the latter stem from a non linear transform. Moreover, multifractal processes are necessarily non Gaussian, heavy tailed and strongly correlated [START_REF] Riedi | Multifractal processes[END_REF]. For these reasons, the design of statistical tests is significantly more complicated for multifractal processes. In particular, the properties of the statistics underlying the test can no longer be obtained analytically, as opposed to the Gaussian H-sssi case. To address such changes in goals and methodology and to overcome the corresponding difficulties, a non parametric bootstrap based test procedure is proposed. Its performance are assessed by application to a large number of realizations of synthetic processes whose (multifractal and non stationarity) properties are known and controlled a priori. We end up with a practical and operational non parametric test procedure, that exhibits satisfactory statistical performance and that can be applied to a single observation of empirical data to assess the true existence of scaling.

EMPIRICAL MULTIFRACTAL ANALYSIS

The theory of multifractal analysis is not recalled here and the reader is referred to e.g., [START_REF] Wendt | Bootstrap for empirical multifractal analysis with application to hydrodynamic turbulences[END_REF][START_REF] Riedi | Multifractal processes[END_REF][START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF]. We concentrate here only on multifractal parameter estimation and on empirical multifractal analysis, referred to as multifractal formalism. Multiresolution quantities. Let X(t) denote the time series to be analyzed. Let dX (j, k) = R R X(t) 2 -j ψ0(2 -j t -k) dt denote its discrete wavelet transform (DWT) coefficients. The mother-wavelet ψ0(t) consists of a reference pattern with a compact time support, characterized by its number of vanishing moments

N ψ ≥ 1: ∀k = 0, 1, . . . , N ψ -1, R R t k ψ0(t)dt ≡ 0 and R R t N ψ ψ0(t)dt = 0. Also, it is such that the collection {2 -j/2 ψ0(2 -j t -k), j ∈ Z, k ∈ Z} forms an orthonormal basis of L 2 (R).
Let us moreover define dyadic intervals as λ = λ j,k = ˆk2 j , (k + 1)2 j ´, and let 3λ denote the union of the interval λ with its 2 adjacent dyadic intervals:

3λ j,k = λ j,k-1 ∪ λ j,k ∪ λ j,k+1
. Following [START_REF] Wendt | Bootstrap for empirical multifractal analysis with application to hydrodynamic turbulences[END_REF][START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF], wavelet Leaders are defined as:

LX (j, k) ≡ L λ = sup λ ′ ⊂3λ |d X,λ ′ |.
In other words, the wavelet Leader LX (j, k) consists of the largest wavelet coefficient |dX (j ′ , k ′ )| at all finer scales 2 j ′ ≤ 2 j in a narrow time neighborhood. It has been shown theoretically that the multifractal formalism is more efficient when wavelet Leaders, rather than wavelet coefficients, are chosen as multiresolution quantities (cf. [START_REF] Wendt | Bootstrap for empirical multifractal analysis with application to hydrodynamic turbulences[END_REF][START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF]). Multifractal attributes. The existence of scaling in data and the measurements of the corresponding parameters are commonly based on structure functions. The structure functions mostly consist of the sample moment estimates of order q of the LX (j, k), or of the sample cumulant estimates of order p of the ln LX (j, k), labeled S(2 j , q) and C(2 j , p), respectively. Scaling properties are practically and operationally defined via the following equations:

S(2 j , q) = Fq2 jζ(q) , (1)
C(2 j , p) = c 0 p + cp ln 2 j . (2)
The scaling exponents ζ(q) and the so-called log-cumulants cp represent the multifractal attributes of X. They are related together as ζ(q) = P p≥1 cpq p /p! [START_REF] Castaing | Log-similarity for turbulent flows[END_REF]. Estimation procedures. Based on Eqs. ( 1) and ( 2), estimations of the ζ(q) and cp are obtained by weighted linear regressions, from scales 2 j 1 to 2 j 2 , of log 2 S(2 j , q) and of log 2 e • C(2 j , p) versus log 2 2 j = j. This has been extensively assessed elsewhere [START_REF] Wendt | Bootstrap for empirical multifractal analysis with application to hydrodynamic turbulences[END_REF] and is not further detailed here.

TIME CONSTANCY TEST

In [START_REF] Veitch | A statistical test for the time constancy of scaling exponents[END_REF], a uniformly most powerful invariant test for the time constancy of the Hurst parameter H of Gaussian H-sssi processes is devised and analyzed. The test is constructed from wavelet coefficient based estimates Ĥ(m) , obtained from adjacent non overlapping subsets X (m) of X, and relies on Gaussianity, independence and (analytically) known statistics of the estimates. Notably, the variance of the estimates is known a priori and independent of the true H. The test statistic reads:

TH = M X m=1 1 σ 2 (m) 0 B @ Ĥ(m) - P M n=1 Ĥ(m) σ 2 (m) P M n=1 1 σ 2 (m) 1 C A 2 .
(

) 3
Under the null hypothesis (H constant), its distribution is known exactly, which enables the formulation of the test.

To adapt the test to multifractal processes, we have to extend it to any multifractal attributes θ ∈ {ζ(q), cp}, whose estimations are based on wavelet Leaders [START_REF] Wendt | Bootstrap for empirical multifractal analysis with application to hydrodynamic turbulences[END_REF]. This induces two major difficulties: i) Variances σ 2 (m) for the θ(m) are no longer known a priori and are likely to depend on the parameter values; ii) The null distribution of the test statistics T θ is no longer known a priori.

BOOTSTRAP TEST

To overcome such severe difficulties in the test formulation and following investigations reported in [START_REF] Wendt | Bootstrap for empirical multifractal analysis with application to hydrodynamic turbulences[END_REF], we propose a non parametric bootstrap based test procedure [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF][START_REF] Zoubir | Bootstrap Techniques for Signal Processing[END_REF]. Let θ denote the multifractal attribute under test. From the time series X to be analyzed, M wavelet Leader based subset estimates θ(m) of θ are obtained from adjacent non overlapping subsets X (m) . Assessing the time constancy of θ then amounts to testing the hypothesis that the random variables { θ(m) , m = 1, . . . , M } have identical mean:

H0 : θ (1) = θ (2) = • • • = θ (M ) . (4)
The test makes use of a bootstrap procedure on the wavelet Leaders, as sketched in Table 1. It is fully specified in e.g. [START_REF] Wendt | Bootstrap for empirical multifractal analysis with application to hydrodynamic turbulences[END_REF] and not further detailed here. We only recall that it consists of a moving time-blocks bootstrap, accounting for dependence amongst wavelet Leaders. Bootstrap test statistic. The test statistic consists of a modified version of Eq. ( 3). It is based on bootstrap variance estimates for the unknown variances, and on the Graybill Deal estimator instead of the maximum likelihood estimator of the consensus mean:

T θ = M X m=1 1 σ2 * (m) 0 B @ θ(m) - P M n=1 θ(m) σ2 * (m) P M n=1 1 σ2 * (m) 1 C A 2 . (5)
First, the {LX (j, k)} are cut into M subsets {LX (m) (j, k)}, corresponding to the subsets X (m) . The subset estimates θ(m) are computed by applying Eqs. ( 1) and ( 2) to the {LX (m) (j, k)}.

Second, the variance estimates σ2 * (m) for each θ(m) are obtained by applying the bootstrap as in Table 1 •) . Table 2 (left) sketches this procedure. Bootstrap null distribution estimation. A bootstrap estimate of the distribution of the test statistic T θ under H0 is obtained from the empirical distribution of:

T * θ = M X m=1 1 σ2 * * (m) 0 B B @ θ * (m) - P M n=1 θ * (m) σ2 * * (m) P M n=1 1 σ2 * * (m) 1 C C A 2 . ( 6
)
First, the bootstrap resampling of Table 1 is applied to the complete set {LX (j, k)} of Leaders, yielding the B resamples {L * (b) 2). Let us emphasize that resampling from the complete set of Leaders, rather than from subsets, is a crucial issue, as it ensures that the θ * (m) all have the same conditional distributions and thus that T * θ reproduces the statistics of T θ under H0, shall X satisfy H0 or H1 (cf.

X (j, k)}, b = 1, • • • , B. Each of these resamples is then cut into M sub- sets {(L * (b) X (j, k)) (m) },
for b = 1, • • • , B for j = 1, • • • , j2
From {LX (j, 1),

X (j, k)}, b = 1, • • • , B from a set of Leaders {LX (j, k)}. {LX } cut Ö Ö * G G a{L * (b) X }aa cut * G G {L * * (b,•) X } cut {LX (m) } * G G estimate {(LX (m) ) * } estimate {(L * (b) X ) (m) } estimate {(L * * (b,•) X ) (m) } estimate θ(m) θ(m) * → σ * (m) u } θ * (b) (m) θ * * (b,•) (m) → σ * * (b) (m) s { T θ T * (b) θ b = 1, • • • , B
Table 2: Procedure for obtaining T θ (left) and T * θ (right) from the wavelet Leaders {LX (j, k)} of X. "cut", "estimate" and " * " stand for cutting a set into M subsets, computing estimates θ from Eqs. ( 1) and ( 2), and bootstrap resampling as in Table 1, respectively.

. This procedure is summarized in Table 2 (right). Bootstrap test. The test is now readily formulated as:

d θ = 1 if T θ > T *
θ,C and 0 otherwise, [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF] where the test critical value T * θ,C is the (1-α) quantile of the empirical distribution of T * θ , for a certain preset significance level α. The critical value of α for which the observed test statistic T θ would be regarded as just decisive against H0 is called the p-value p θ of T θ .

PERFORMANCE ASSESSMENT AND RESULTS

Monte-Carlo simulations. To evaluate the performance of the proposed test procedures, we apply them to a large number NMC of realizations of length N of a synthetic multifractal process with a priori known and controlled multifractal properties. We choose a well known and easy to simulate multifractal process called multifractal random walk (MRW) for which c1, c2 = 0; p ≥ 3 : cp ≡ 0. For a definition of this process, see e.g. [START_REF] Bacry | Multifractal random walk[END_REF]. The simulation parameters are set to NMC = 1000, N = 2 15 , B = B2 = 99 and α = 0.1. For multifractal attribute estimation, we use Daubechies wavelets with N ψ = 3. The regression range is set to j1 = 3 and j2(M ) = log 2 N -log 2 M -(2N ψ -1) (cf. [START_REF] Wendt | Bootstrap for empirical multifractal analysis with application to hydrodynamic turbulences[END_REF]). The multifractal parameters (specified below) are chosen to correspond to realistic situations observed in actual data (for instance, c2 ≈ -0.025 is a commonly accepted value in turbulence, cf. [START_REF] Wendt | Bootstrap for empirical multifractal analysis with application to hydrodynamic turbulences[END_REF][START_REF] Castaing | Log-similarity for turbulent flows[END_REF]). Performance assessment. The performance of the test procedures, given that a certain hypothesis H (•) is true, are assessed by their mean rejection rates and p-values,

dH (•) θ = b EN M C {d θ |H (•) } (8) pH (•) θ = b EN M C {p θ |H (•) }, (9)
where b EN M C stands for the mean over Monte Carlo simulations. For space reasons, we only report results for θ ∈ {c1, c2} and M = {2, 4} with equal subset lengths. Similar results can be obtained for cp, p ≥ 3 and ζ(q), other choices of M , and splitting into subsets of non equal length.

Performance under H0

The performance under H0 are studied on processes with constant multifractal attributes {c1, c2}. Table 3 summarizes results for two different sets of parameters: {c1, c2} = {0.75, -0.01} (left) and {c1, c2} = {0.8, -0.02} (right). We observe that the mean rejection rates dH 0 θ are close to the preset significance level α for both c1 and c2, and for both parameter settings. Furthermore, the mean p values pH 0 θ are close to 0.5, indicating a satisfactory null distribution estimation. Indeed, under H0, the p-value would be uniformly distributed between on 0 and 1 if the test was based on the exact null distribution of the test statistic. The results lead us to the conclusion that the empirical distribution of T * θ is a satisfactory approximation of the null distribution of T θ under H0, and that it is robust with respect to the precise values of the multifractal parameters. For sake of completeness, we note that results not reported here show that the slight discrepancies in the observed test sizes are mainly due to small differences between the variance of θ(m) (as measured form Monte-Carlo simulations) and its bootstrap estimation σ * 2 (m) .

Performance under H1

To study the power of the proposed tests, we need to define an alternative hypothesis. One could imagine many forms of non stationary processes, a number of them being likely to mimic scaling behaviors when analyzed blindly over the entire time series. Here, we study one of the simplest such alternatives: processes possess piecewise constant multifractal attributes. H1 is thus analyzed with an alternative consisting of the concatenation of two truly multifractal processes of equal length with different multifractal attributes {c

(i) 1 , c (i)
2 }i=1,2. Two cases are investigated. Non constant c1, constant c2. In the first case, which we denote H1(c1), we set c

≡ c

(2) 2 = c2, i.e., c2 is constant, while c1 is not. Thus, dc 1 assesses the power of the test for time constancy of c1. The parameters are set to c2 = -0.02 and c = 0.8), corresponding to values that are in practice considered as being very close, the test rejects the time constancy hypothesis for c1 with a probability above 0.6 (M = 2) and close to 0.5 (M = 4). Conversely, the mean test decisions dc 2 reproduce closely the preset significance level α and remain constant when c

1 varies, indicating that the time constancy test for c2 is not subject to cross-influence from changes in c1. We conclude, first, that the test for time constancy of c1 is powerful, and second, that the test for constancy of c2 closely reproduces the level α, independently of c

≡ c

(2) 1 = c1, i.e., c1 is constant, while c2 is not. Therefore, dc 2 assesses the power of the test for time constancy of c2. The parameters are set to c1 = 0.75 and c

(1) 2 = {-0.11, -0.09, • • • , -0.01}, c (2) 2
= -0.01. Fig. 1 (bottom) shows test decisions dc 1 (solid red line) and dc 2 (dashed blue line) as a function of the step size c

p -c

(2) p and thus on the precise hypothesis H1. Moreover, the T * cp,C equal the critical values under H0 (given by the rightmost points). This shows that the empirical distribution of T * θ under H1 provides us with a robust and accurate null distribution estimation, as test design demands.

CONCLUSIONS AND PERSPECTIVES

We have devised a practical procedure for discriminating the existence of true scaling properties against non stationarities. It consists of a bootstrap based test for the constancy along time of wavelet Leader based multifractal parameter estimates. We have shown, by means of numerical simulations, that this bootstrap based test procedure is reliable and powerful. Notably, the empirical distribution of T * θ under H1 yields an accurate estimation of the null distribution, a central feature for relevant test design. Our procedure successfully addresses this nontrivial issue by combining a "split then bootstrap" for T θ and a "bootstrap then split" for T * θ . It has heavy computational cost (due to double bootstrap) but remains, to the best of our knowledge, the only procedure practically available. It can be applied to a single observation of real data with unknown statistical characteristics. The impact of choosing M remains to be discussed in terms of standard trade-off between type-I and type-II errors: A test with too small M may miss non stationarities, choosing M too large results in a lack of power due to poor estimations, hence the existence of an optimal M for a given but unknown alternative hypothesis. The procedure can be further extended to testing the constancy along time of the whole structure functions (S(2 j , q), C(2 j , p)) or to testing jointly the constancy of a vector of multifractal attributes. This test is being applied to a set of biomedical data.

  to each subset {LX (m) (j, k)}, yielding B resamples {(LX (m) (j, k)) * (b) }, b = 1, • • • , B. Then Eqs. (1) and (2) are used on each of these resamples to obtain the bootstrap subset estimates θ(m) * (b) , from which the sample variance is finally estimated: σ2 * (m) = d Var * θ(m) * (

  and the subset estimations θ * (b) (m) are computed from Eqs. (1) and (

Fig. 2 )

 2 . The variance estimates σ2 * * (b) (m) of θ * (b) (m) are obtained by first applying bootstrap to each {L * (b) X (j, k)}, giving the B2 double bootstrap resamples

  {L * * (b,b 2 ) X (j, k)}, b2 = 1, • • • , B2. Each of these double bootstrap resamples is in turn cut into M subsets {(L * * (b,b 2 ) X (j, k)) (m) }, enabling the computation of the double bootstrap subset estimates θ * * (b,b 2 )

Table 3 :

 3 Mean rejection rates dH 0 θ and p values pH 0 θ (H0, M = 2).

1 = 1 Fig. 1 :

 111 Fig. 1: Test decisions dc 1 (solid red) and dc 2 (dashed blue) under H1(c1) (top) and H1(c2) (bottom) for M = 2 (circles) and M = 4 (squares) as a function of c (1) p -c (2) p . The solid black line indicates the preset significance level α.

1 .

 1 Constant c1, non constant c2. In the second case, which we denote H1(c2), we set c

2 .

 2 the roles of dc 1 and dc 2 , conclusions are similar to those obtained under H1(c1): satisfactory power of the test for time constancy of c2, and insensitivity of the test on c1 with respect to level change c Null distribution estimation under H1. Fig.2shows bootstrap test critical values (as defined in Eq. (7)) T * c 1 ,C under H1(c1) (left) and T * c 2 ,C under H1(c2) (right). The circles and the bars correspond, respectively, to b EN M C T * θ,C and to 1.64 • d Std N M C T * θ,C . We observe that the T * cp,C do not depend on the step size c

Fig. 2 :

 2 Fig. 2: Bootstrap test critical values T * θ,C under H1 (mean value and 1.64σ T * θ,C bars, obtained through Monte Carlo simulations): Left, T * c 1 ,C under H1(c1); Right T * c 2 ,C under H1(c2) (M = 2).

Table 1 :

 1 Moving blocks bootstrap for obtaining B bootstrap resamples {L * (b)

			• • • , LX (j, nj)}
	random draw, with replacement, circular and
	overlapping blocks to form an unsorted collection
	n	L	* (b) X (j, 1), • • • , L	* (b) X (j, nj) o
	end
	end