Eric Grivel

Bastien Berthelot

Pierrick Legrand

A matrix approach to get the variance of the square of the fluctuation function of the DFA

Keywords: DFA, FA, statistical analysis, variance

The fluctuation analysis (FA) and the detrended fluctuation analysis (DFA) make it possible to estimate the Hurst exponent H, which characterizes the selfsimilarity of a signal. Both are based on the fact that the so-called fluctuation function, which can be seen as an approximation of the standard deviation of the process scaled in time by multiplying the time variable by a positive constant, depends on H. The main novelty of the paper is to provide the expression of the variance of the square of the fluctuation function, by using a matrix formulation.

We show that it depends on the correlation function of the signal under study when it is zero-mean and Gaussian. Illustrations are given when dealing with a zero-mean white Gaussian noise. Moving average processes and first-order autoregressive processes are also addressed.

Introduction

In signal processing, combining signal feature extraction and classification is commonly used even if deep-learning based approaches become more and more popular. There are various types of features that can be considered: the signal power in some frequency bands or the power ratio such as the so-called "LF/HF ratio". As an alternative, quantities such as the mean of the signal, its power, the zero-crossing rate but also the parameters of an a priori model representing the data can be used. Autoregressive moving average parameters (ARMA) and mel frequency cepstral coefficients (MFCC) are often chosen by the practitioners. Other types of approaches aim at quantifying the regularity of the signal, which can be seen as the propensity to chaos. These techniques can be based on entropy analysis such as Shannon entropy or the multi-scale entropy [START_REF] Gao | Multiscale entropy analysis of biological signals: A fundamental bi-scaling law[END_REF]. When dealing with monofractal signals, the fluctuation analysis (FA) and the detrended fluctuation analysis (DFA) initially proposed by Peng early in the 90ies [START_REF] Peng | Long-range correlations in nucleotide sequences[END_REF] [3] make it possible to estimate the Hurst coefficient, which characterizes the overall regularity of a signal.

In the FA, the temporal mean of the signal that has been integrated and low-pass filtered is computed whereas the DFA aims at addressing not only wide-sensestationary processes but also non-stationary ones. It is based on the following ideas: 1/ removing the trend of the centered integrated signal. A great deal of interest has been paid to these methods and various types of contributions were presented. It should be noted that as various authors working in different fields from mathematics to signal processing passing by biomedical applications are interested in these methods, some topics may have been addressed in a similar manner at the same period:

1/ due to the simplicity of execution, these methods have been applied in various fields and more particularly in biomedical applications to classify biosignals such as electroencephalograms [START_REF] Zappasodi | Fractal dimension of eeg activity senses neuronal impairment in acute stroke[END_REF], eye movements [START_REF] Kelty-Stephen | Fractal fluctuations in gaze speed visual search[END_REF] speech signals [START_REF] Pitsikalis | Analysis and classification of speech signals by generalized fractal dimension features[END_REF] in order to detect a pathology, a stress, etc.

2/ as the DFA is based on the estimation of the trend of the centered integrated signal, different variants have been proposed and differ by the way the trend is estimated. While the standard DFA is based on a discontinuous piecewise linear function, extensions to higher-order polynomials have been proposed, leading to the higher-order DFA. In addition, the adaptive fractal analysis (AFA) [START_REF] Riley | A tutorial introduction to adaptive fractal analysis[END_REF], the regularized DFA (RDFA) [START_REF] Berthelot | Regularized DFA to study the gaze position of an airline pilot[END_REF] and the continuous DFA (CDFA) [START_REF] Berthelot | Alternative ways to compare the detendred fluctuation analysis and its variants. application to visual tunneling detection[END_REF] were developed to avoid discontinuities in the global trend. Low-pass linear filter based methods such as the detrending moving average (DMA) [START_REF] Alessio | Second-order moving average and scaling of stochastic time series[END_REF] were also developed. Finally, methods bridging the gap between the DFA and the DMA have been proposed. The higher-order DMA [START_REF] Arianos | Self-similarity of higher-order moving averages[END_REF] can be interpreted as a solution based on the locally estimated scatterplot smoothing (LOESS) whereas the one based on the locally weighted scatterplot Smoothing (LOWESS) has been recently accepted for publication [START_REF] Berthelot | New Variants of DFA based on LOESS and LOWESS methods: generalization of the detrended moving average[END_REF]. See figure 1. for intance in [START_REF] Kantelhardt | Multifractal detrended fluctuation analysis of nonstationary time series[END_REF] as well as fast versions [START_REF] Tsujimoto | Fast algorithm of long-range cross-correlation analysis using savitzky-golay detrending filter and its application to biosignal analysis[END_REF].

4/ mathematical analyzes were also led to better understand the behaviors of the DFA and its variants. In [START_REF] Taqqu | Estimators for long range dependence: an empirical study[END_REF][START_REF] Taqqu | On estimating the intensity of long range dependence in finite and infinite variance time series[END_REF][START_REF] Kantelhardt | Detecting long-range correlations with detrended fluctuation analysis[END_REF][START_REF] Crato | Statistical properties of detrended fluctuation analysis[END_REF][START_REF] Kiyono | Establishing a direct connection between detrended fluctuation analysis and fourier analysis[END_REF][START_REF] Höll | Detrended fluctuation analysis and the difference between external drifts and intrinsic diffusion-like nonstationarity[END_REF][START_REF] Kiyono | Theory and applications of detrending-operation-based fractalscaling analysis[END_REF][START_REF] Løvsletten | Consistency of detrended fluctuation analysis[END_REF][START_REF] Höll | Theoretical foundation of detrending methods for fluctuation analysis such as detrended fluctuation analysis and detrending moving average[END_REF][START_REF] Sikora | Probabilistic properties of detrended fluctuation analysis for gaussian processes[END_REF][START_REF] Grivel | Dfa-based abacuses providing the hurst exponent estimate for short-memory processes[END_REF], the authors often aim at analyzing some properties of the square of the fluctuation function. Thus in [START_REF] Bardet | Asymptotic properties of the detrended fluctuation analysis of long range dependence processes[END_REF], Bardet et al. study the time series defined by the square of the fluctuation function computed for one part of size N . Among the results presented, they give an asymptotic 1 expression of its statistical mean, its variance but also its covariance function when a fractional Gaussian Noise is studied. These expressions depend on N and H. They also show that the expectation of the square of the fluctuation function of the sum of two independent zero-mean processes is the sum of the square of 1 i.e. when N becomes large.

the fluctuation function of the processes. In [START_REF] Höll | The relationship between the detrendend fluctuation analysis and the autocorrelation function of a signal[END_REF], Höll et al. propose to express the square of the fluctuation function from an estimation of the correlation function while in [START_REF] Höll | Theoretical foundation of detrending methods for fluctuation analysis such as detrended fluctuation analysis and detrending moving average[END_REF] the problem is presented differently, but the ideas are quite similar. Thus, the signal under study is assumed to be the sum of two terms: one being deterministic and the second being stochastic. After defining the so-called path as the integrated version of the stochastic part of the signal, Höll et al. analyze what they called the mean squared displacement (MSD) of the path, which is the correlation function of the path for a lag equal2 to 0. They express it from the correlation function of the stochastic part and from the Hurst exponent H for different types of random parts. Instead of presenting the various steps of the approaches as a set of equations to be implemented, other authors have introduced a matrix formulation. Thus, Lovsletten in [START_REF] Løvsletten | Consistency of detrended fluctuation analysis[END_REF] and our team [START_REF] Grivel | Dfa-based abacuses providing the hurst exponent estimate for short-memory processes[END_REF][START_REF] Berthelot | 2D Fourier Transform Based Analysis Comparing the DFA with the DMA[END_REF] independently took advantage of the matrix form to conduct statistical analysis or to give physical interpretations of some steps of the DFA. Their approach slightly differ: the mean is not necessarily removed in the algorithms; the detailed expressions of the matrices for each variant can be given or not; some simplifications are made or not. In [START_REF] Berthelot | 2D Fourier Transform Based Analysis Comparing the DFA with the DMA[END_REF], when the signal is wide sense stationary (w.s.s.), the matrix formulation is a good way 1/ to express the statistical mean of the square of the fluctuation function from the theoretical correlation function of the signal3 2/ to represent the processing chain from the integration step to the detrending step by an "equivalent" linear filter whose impulse response depends on the way to estimate the trend. The differences between the DFA and its variants can hence be highlighted by comparing the frequency responses of the equivalent filters. The non-stationary case is also addressed. In [START_REF] Grivel | Dfa-based abacuses providing the hurst exponent estimate for short-memory processes[END_REF], abacuses are derived to link the parameters of the processes modeling the signal under study, the number of samples available and the value provided by the DFA when the statistical mean of the square of the fluctuation function is considered.

In [START_REF] Sikora | Probabilistic properties of detrended fluctuation analysis for gaussian processes[END_REF], the authors focus their attention on the distribution of the fluctuation function of the DFA when a centered Gaussian process is studied. Starting from the definition of the square of the fluctuation function as the sum of the square of elementary differences between the integrated signal and the trend, they propose to express each elementary error as a quadratic form and more particularly as the product of one "difference vector" storing the samples of difference between the signal and its trend in each frame multiplied by its transpose. Then, the authors propose to derive the statistical mean and the variance of each elementary error.

They express them in terms of the eigenvalues of the correlation matrix of the difference vector. They apply it to a Gaussian white noise and a fractionally integrated noise. Moreover, the specificity of the DFA has been also studied in [START_REF] Maraun | Tempting long-memory on the interpretation of dfa results[END_REF] and [START_REF] Morariu | Detrended fluctuation analysis of autoregressive processes[END_REF] where an empirical analysis for autoregressive (AR) processes is led. The analysis of the fluctuation function of the DFA when dealing with a first-order AR process was also addressed in [START_REF] Höll | The fluctuation function of the detrended fluctuation analysis -investigation on the ar(1) process[END_REF].

Remark: a more exhaustive state of the art is provided in [START_REF] Berthelot | Alternative ways to compare the detendred fluctuation analysis and its variants. application to visual tunneling detection[END_REF] and [START_REF] Grivel | Dfa-based abacuses providing the hurst exponent estimate for short-memory processes[END_REF]. As they have been recently published, the reader may refer to both papers.

In this paper, based on a matrix formulation, we derive the expression of the variance of the square of the fluctuation function for any process under study.

It appears that it depends on the correlation function of the signal when the latter is zero-mean. The reasoning is first illustrated with a zero-mean white Gaussian noise for which an analytical expression of the variance depending on the parameters of the DFA and the noise variance is provided. The impact of the order of the variant is presented. A comparison with the statistical properties of the square of the fluctuation function obtained with the FA is also done. Then, the case of the q th -order moving average (MA) process and a first-order AR process is addressed, by presenting the methodology to follow. The illustrations are provided for short-memory processes, but similar results could be obtained with long-memory processes.

The remainder of the paper is organized as follows: after introducing some notations in section 2, section 3 recalls the main steps of the FA and the DFA.

The matrix form is used at each step. In section 4, some comments are done and serve as a basis to derive the variance of the square of the fluctuation function of the DFA or its higher-order variant. The white noise case is first considered.

Then, and more generally, a way to study random processes is proposed.

Notations

Let us first define some notations that will be useful to present the main steps of the FA and the DFA. They are the same as the ones used in [START_REF] Berthelot | Alternative ways to compare the detendred fluctuation analysis and its variants. application to visual tunneling detection[END_REF] or [START_REF] Grivel | Dfa-based abacuses providing the hurst exponent estimate for short-memory processes[END_REF].

• 0 j×k is a matrix of size j × k filled with 0s and 1 j×k is a matrix of size j × k filled with 1s.

• J j = I j -1 j 1 j×j , with I j the identity matrix.

• H M = M -1 r=0 diag(1 1×M -r , -r
) is a lower triangular matrix filled with 1s. H N is similarly defined, except that M is replaced by N .

• T r(X, r) is the sum of the elements of the r th diagonal of the matrix X.

• Y and Y int are two column vectors of size M storing {y(n)} n=1,...,M and {y int (n)} n=1,...,M , being respectively the samples of the signal and the integrated centered ones. Given the notations introduced above, one has:

Yint = [yint(1), yint(2), ..., yint(M )] T = HM JM Y (1)
Introducing C j,k = [0 j×k I j 0 j×(M -(j+k)) ] a matrix of size (j, M ), the first LN elements of the vector Y int denoted as Y int (1 : LN ) can be expressed as follows:

Yint(1 : LN ) = [yint(1), yint(2), ..., yint(LN )] T = CLN,0Yint = (1)
CLN,0HM JM Y (2)

• Y l = Y ((l -1)N + 1 : lN ) is the l th part of size N of the vector Y .

Presentation of the FA and the DFA

A time-continuous process y is said to be self-similar with parameter H if and only if y(N t) Therefore, the standard deviation of y(N t) is equal to the standard deviation of y(t), up to a multiplicative factor equal to N H .

In [START_REF] Peng | Long-range correlations in nucleotide sequences[END_REF], Peng suggested estimating H on discrete-time signals with the FA, operating with the following two steps:

Step 1. the signal y is integrated.

Step 2. By denoting < . > the temporal mean, the square of the fluctuation function F 2 0 0 0 (N ) is computed for different values of the lag N .

F 2 0 0 0 (N ) = i+N j=i+1 y(j) 2 = i+N j=i+1 y 2 (j) + 2 N -1 r=1 i+N -r j=i+1 y(j)y(j + r) (3)
Due to the integration step, F 2 0 0 0 (N ) is proportional to N 2(H+1) . Therefore, log(F 2 0 0 0 (N )) is represented as an affine function of log(N ) to estimate H in the least-square (LS) sense.

Since the FA was sensitive to the existence of a trend, the DFA was developed by Peng in [START_REF] Peng | Mosaic organization of DNA nucleotides[END_REF]. Using a matrix formalism on M consecutive samples {y(m)} m=1,...,M of the signal, the standard DFA is defined by the following steps:

Step 1. The profile is computed as follows:

y int (m) = m i=1 (y(i) -µ y ) where µ y = 1 M M m=1 y(m) is the mean of y.
Step 2. Splitting the profile into L non-overlapping parts of length N , i.e. {y int,l (n)} l=1,...,L with 4 n ∈ [[1; N ]], the d th -degree-polynomial based trend of the l th part is given by: all the local trends can be estimated in the LS sense for d > 0 as follows:

T loc,l,d = A l,d θ l,d ∀l ∈ [[1; L]] (4)
Θ1 1 1,d = (A T 1 1 1,d A 1 1 1,d ) -1 A T 1 1 1,d C LN,0 Y int (5)
Using ( 1) and ( 5), the global trend vector T 1 1 1,d can be expressed by:

T 1 1 1,d = A 1 1 1,d Θ1 1 1,d = A 1 1 1,d (A T 1 1 1,d A 1 1 1,d ) -1 A T 1 1 1,d C LN,0 H M J M Y (6)
Step 3. The residual vector R 1 1 1,d , also known as detendred profile, is:

R 1 1 1,d = C LN,0 H M J M Y -T 1 1 1,d = B 1 1 1,d Y (7)
Then, given

Γ 1 1 1,d = 1 LN B T 1 1 1,d B 1 1 1,d , the power of the residual F 2 1 1 1,d (N )
, which corresponds to the square of the fluctuation function, can be deduced as follows:

F 2 1 1 1,d (N ) = T r Γ 1 1 1,d Y Y T (8)
In a recent paper [START_REF] Grivel | Dfa-based abacuses providing the hurst exponent estimate for short-memory processes[END_REF], we analyzed the structure of the matrix Γ 1 1 1,d and pointed out the fact that it was a block-diagonal matrix, where each block of dimension N × N is equal to γ 1 1 1,d defined by:

γ 1 1 1,d = 1 N H T N A 1 1 1,d T A 1 1 1,d H N (9)
where

A 1 1 1,d = I N -A 1,d (A T 1,d A 1,d ) -1 A T 1,d .
Taking into account the structure of the matrix Γ 1 1 1,d and its link with γ 1 1 1,d , the square of the fluctuation function F 2 1 1 1,d (N ) given in ( 8) can be rewritten as follows:

F 2 1 1 1,d (N ) = 1 L L l=1 Y T l γ 1 1 1,d Y l (10)
Remark: as we will see in the rest of the paper, ( 8) and ( 10) can be used to derive the statistical mean of F 2 1 1 1,d (N ) whereas [START_REF] Alessio | Second-order moving average and scaling of stochastic time series[END_REF] will be rather considered to get the variance.

Statistical analysis

Given the above brief presentation, after recalling the expression of the statistical mean, let us derive the variance of the square fluctuation function when the DFA (or the higher-order DFA) is used.

About the statistical mean of the square of the fluctuation function

When the signal is w.s.s., the expectation of the square of the fluctuation function can be expressed the same way. Indeed, using (3), one has:

E[F 2 0 0 0 (N )] = N -1 r=1-N T r(Γ 0 0 0 , r)R y,y (r) (11)
where R y,y (τ ) is the autocorrelation function of the signal y and Γ 0 0 0 is a square matrix of size N whose every element is equal to 1.

Similarly, as initially shown in [START_REF] Berthelot | Alternative ways to compare the detendred fluctuation analysis and its variants. application to visual tunneling detection[END_REF], taking the statistical mean of (8) and due to the link between Γ 1 1 1,d and γ 1 1 1,d , one gets:

E[F 2 1 1 1,d (N )] = N -1 r=-N +1 T r(Γ 1 1 1,d , r)R y,y (r) = N -1 r=-N +1 T r(γ 1 1 1,d , r)R y,y (r) (12)
Remark: The above result is not new. See [START_REF] Berthelot | Alternative ways to compare the detendred fluctuation analysis and its variants. application to visual tunneling detection[END_REF] and [START_REF] Grivel | Dfa-based abacuses providing the hurst exponent estimate for short-memory processes[END_REF]. It is recalled in the paper because it will be useful to compute the variance of the square of the fluctuation function.

Remark 2: The same result could be obtained using the new expression [START_REF] Alessio | Second-order moving average and scaling of stochastic time series[END_REF] we introduce in this paper. Indeed, one has:

E[F 2 1 1 1,d (N )] = 1 L L l=1 E[Y T l γ 1 1 1,d Y l ] = 1 L L l=1 E[T r(γ 1 1 1,d Y l Y T l )] (13) = 1 L L l=1 T r(γ 1 1 1,d E[Y l Y T l ]) When y is w.s.s., E[Y l Y T l ] = R Y,N
for any l = 1, ..., L. This is the correlation matrix of a vector storing N consecutive samples of the process y. It is a Toeplitz matrix. This leads to:

E[F 2 1 1 1,d (N )] = 1 L L l=1 T r(γ 1 1 1,d R Y,N ) = T r(γ 1 1 1,d R Y,N ) (14)
=

N -1 r=-N +1 T r(γ 1 1 1,d , r)R y,y (r)
Let us now address the variance of the square of the fluctuation function of the DFA, which is the main novelty of the paper.

About the variance of the square of the fluctuation function

To get the variance of the square of the fluctuation function, F 4 1 1 1,d (N ) needs to be first expressed because E[F 4 1 1 1,d (N )] must be computed. It satisfies:

F 4 1 1 1,d (N ) = 1 L 2 L l=1 Y T l γ 1 1 1,d Y l 2 (15) = 1 L 2 L l=1 Y T l γ 1 1 1,d Y l Y T l γ 1 1 1,d Y l + 2 L 2 L-1 l=1 L m=l+1 Y T l γ 1 1 1,d Y l Y T m γ 1 1 1,d Y m
The next step of the reasoning is to express F 4 1 1 1,d (N ) using a matrix form. Therefore, let us introduce the following two column vectors of length N (N +1) 2 :

Y T l = y 2 ((l -1)N + 1) ... y 2 (lN ) (16)
y((l -1)N + 1)y((l -1)N + 2) ... y((l -1)N + 1)y(lN ) y((l -1)N + 2)y((l -1)N + 3) ... y((l -1)N + 2)y(lN ) ... y(lN -1)y(lN )

The above vector Y T l is defined from the squares of the N samples of the vector Y l and the N 2 -N 2 cross-products that can be defined from two different samples of Y l .

γ T 1 1 1,d = γ 1 1 1,d (1, 1) ... γ 1 1 1,d (N, N ) 2γ 1 1 1,d (1, 2) ... 2γ 1 1 1,d (1, N ) (17) 2γ 1 1 1,d (2, 3) ... 2γ 1 1 1,d (2, N ) ... 2γ 1 1 1,d (N -1, N )
where γ 1 1 1,d (i, j) is the element of the matrix γ 1 1 1,d located at the i th row and the j th column.

Remark: there are other ways to define the vector Y T l and consequently γ T 1 1 1,d . Thus, another definition could be:

Y T l = y 2 ((l -1)N + 1) ... ... ... y((l -1)N + 1)y(lN ) (18)
y 2 ((l -1)N + 2) ... ... y((l -1)N + 2)y(lN ) y 2 ((l -1)N + 3) ... y((l -1)N + 3)y(lN ) ...

y 2 (lN )
In this case, the vector Y l consists of N blocks whose dimension decreases from becomes: As we will see, depending on the random process under study, the second definition can be more convenient than the first one.

N
γ T 1 1 1,d = γ 1 1 1,d (1, 1) ... γ 1 1 1,d (1, N ) 2γ 1 1 1,d (1, 2) ... 2γ 1 1 1,d (1, N ) (19)
At this stage and whatever the definition chosen for Y l (i.e. ( 16) or ( 18)), let us rewrite the scalar Y T l γ 1 1 1,d Y l as follows:

Y T l γ 1 1 1,d Y l = Y T l γ 1 1 1,d = γ T 1 1 1,d Y l (20)
This hence leads to:

   Y T l γ 1 1 1,d Y l Y T l γ 1 1 1,d Y l = Y T l γ 1 1 1,d γ T 1 1 1,d Y l = T r(γ 1 1 1,d γ T 1 1 1,d Y l Y T l ) E[Y T l γ 1 1 1,d Y l Y T l γ 1 1 1,d Y l ] = T r(γ 1 1 1,d γ T 1 1 1,d E[Y l Y T l ]) (21)
where E[Y l Y T l ] is the correlation matrix of the vector Y l defined from the cross moments of different orders of the process y.

Similarly, one has:

   Y T l γ 1 1 1,d Y l Y T m γ 1 1 1,d Y m = Y T l γ 1 1 1,d γ T 1 1 1,d Y m = T r(γ 1 1 1,d γ T 1 1 1,d Y m Y T l ) E[Y T l γ 1 1 1,d Y l Y T l γ 1 1 1,d Y m ] = T r(γ 1 1 1,d γ T 1 1 1,d E[Y m Y T l ]) (22)
where E[Y m Y T l ] is the cross-correlation matrix of the vectors Y m and Y l . When the signal is w.s.s., R Y is the correlation matrix of Y l for any value of l = 1, ..., L. Combining ( 15), ( 21) and ( 22) leads to:

E[F 4 1 1 1,d (N )] = 1 L 2 L l=1 T r(γ 1 1 1,d γ T 1 1 1,d E[Y l Y T l ]) + 2 L 2 L-1 l=1 L m=l+1 T r(γ 1 1 1,d γ T 1 1 1,d E[Y m Y T l ]) (23) = 1 L T r(γ 1 1 1,d γ T 1 1 1,d RY ) + 2 L 2 L-1 l=1 L m=l+1 T r(γ 1 1 1,d γ T 1 1 1,d E[Y m Y T l ])
Therefore, using ( 14) and ( 23), the variance of the square of the fluctuation function is given by:

V ar([F 2 1 1 1,d (N )]) = E[F 4 1 1 1,d (N )] -E 2 [F 2 1 1 1,d (N )] (24) = 1 L T r(γ 1 1 1,d γ T 1 1 1,d R Y ) + 2 L 2 L-1 l=1 L m=l+1 T r(γ 1 1 1,d γ T 1 1 1,d E[Y m Y T l ]) -T r 2 (γ 1 1 1,d R Y,N )
The method we present provides an expression of the variance using a matrix form. In terms of methodology, no approximation is done whatever the value of N and no assumption is made on the number of samples available or on the ergodicity of the random process under study. This expression can be considered for any w.s.s. process provided that R Y,N , E[Y m Y T l ] and R Y are available. In the next section, let us apply this approach on different processes.

Applications

We first propose to apply our analysis on a toy example which is a zero-mean white Gaussian noise. For this case, the means and the variances of the square of the fluctuation functions of the FA and the DFA are compared. Note that we provide analytical expressions of these quantities, which depend on L, N , the elements of the matrix Γ 1 1 1,d and the variance of the white noise. Then, we focus our attention on the way to proceed to obtain the variance of the square of fluctuation when applying the DFA or its higher-order variant on MA processes or on a first-order AR process.

Case of a white noise

Let us consider a w.s.s. zero-mean white Gaussian noise with variance σ 2 m whose properties are recalled below: On the one hand, given the information provided above, one has for the FA:

   E[y 2 (k)] = σ 2 m E[y 4 (k)] = 3σ 4 m E[y 2 (k)]E[y 2 (l)] = σ 4 m with k = l (
E[F 2 0 0 0 (N )] = N σ 2 m ( 26
)
This confirms the statement presented above. Then, starting from (3), after development and simplification, one can show that the second-order moment F 2 0 0 0 (N ) for a white noise is given by E[F 4 0 0 0 (N )] = 3N 2 σ 4 m . This makes it possible to express the variance of the square of the fluctuation function when the FA is used:

V ar[F 2 0 0 0 (N )] = 3N 2 σ 4 m -N 2 σ 4 m = 2N 2 σ 4 m ( 27
)
On the other hand, one has for the DFA:

E[F 2 1 1 1,d (N )] = σ 2 m T r(Γ 1 1 1,d , 0) (28)
The study made in [START_REF] Höll | The relationship between the detrendend fluctuation analysis and the autocorrelation function of a signal[END_REF] and confirmed by numerical simulations in [START_REF] Grivel | Dfa-based abacuses providing the hurst exponent estimate for short-memory processes[END_REF] showed that:

E[F 2 1 1 1,1 (N )] = σ 2 m N 2 -4 15N ( 29
)
5 Let us recall that looking at the slope is useful when representing the logarithm of the fluctuation function with respect to log(N ).

This means that E[F 2 1 1 1,1 (N )] tends to become proportional to N provided that N takes large values.

Let us now address V ar[F 2 1 1 1,d (N )]. More particularly, let us express the three terms of the sum defining V ar([F 2 1 1 1,d (N )]) in [START_REF] Sikora | Probabilistic properties of detrended fluctuation analysis for gaussian processes[END_REF]. In this case, Y l and Y m are independent.

E[Y l Y T l ] = E[Y l ]E[Y T l ].
In addition, R Y,N = σ 2 m I N . Thus, using [START_REF] Løvsletten | Consistency of detrended fluctuation analysis[END_REF], one has:

T r(γ 1 1 1,d γ T 1 1 1,d E[Y m Y T l ]) = T r(γ T 1 1 1,d E[Y m Y T l ]γ 1 1 1,d ) = γ T 1 1 1,d E[Y m Y T l ]γ 1 1 1,d (30)
= indep. γ T 1 1 1,d E[Y m ]E[Y T l ]γ 1 1 1,d = (20)
E[Y T m γ 1 1 1,d Y m ]E[Y T l γ 1 1 1,d Y l ] = T r 2 (γ 1 1 1,d R Y,N ) = σ 4 m T r 2 (γ 1 1 1,d )
Note that the trace is removed in the second equality presented in [START_REF] Morariu | Detrended fluctuation analysis of autoregressive processes[END_REF] because

the quantity γ T 1 1 1,d E[Y m Y T l ]γ 1 1 1
,d is a scalar. Using the above result, the expression of the variance becomes:

V ar([F 2 1 1 1,d (N )]) = 1 L T r(γ 1 1 1,d γ T 1 1 1,d R Y ) + 2 L 2 L 2 -L 2 σ 4 m T r 2 (γ 1 1 1,d ) -σ 4 m T r 2 (γ 1 1 1,d ) (31)
= 1 L T r(γ 1 1 1,d γ T 1 1 1,d R Y ) - σ 4 m L T r 2 (γ 1 1 1,d )
The next step is to express T r(γ

1 1 1,d γ T 1 1 1,d R Y ).
First of all, let us give the expression of R Y based on the first definition (16) of Y l . It is a block diagonal matrix.

The first block of size N × N is defined by:

R Y (1 : N, 1 : N ) =         3σ 4 m σ 4 m . . . σ 4 m σ 4 m . . . . . . . . . . . . . . . . . . σ 4 m σ 4 m . . . σ 4 m 3σ 4 m         (32)
The N -1 other blocks are square and of dimensions l that decreases from N -1 to 1. They are equal to σ 4 m I l . Therefore, one has:

1 L T r(γ 1 1 1,d γ T 1 1 1,d R Y ) = σ 4 m L 3 N k=1 γ 2 1 1 1,d (k, k) + 2 N -1 k=1 N l=k+1 γ 1 1 1,d (k, k)γ 1 1 1,d (l, l) (33) + 4 N -1 k=1 N l=k+1 γ 2 1 1 1,d (k, l)
To end up, let us express T r 2 (γ 1 1 1,d ):

T r 2 (γ 1 1 1,d ) = N k=1 (γ 1 1 1,d (k, k)) 2 + 2 N -1 k=1 N l=k+1 γ 1 1 1,d (k, k)γ 1 1 1,d (l, l) (34)
Therefore, the variance of the square of the fluctuation function for the DFA -in its standard version where d = 1 or in its higher-order variant-can be rewritten as follows:

V ar[F 2 1 1 1,d (N )] = 2σ 4 m L N k=1 (γ 1 1 1,d (k, k)) 2 + 2 N -1 k=1 N l=k+1 (γ 1 1 1,d (k, l)) 2 (35) 5.1.2. Illustration
In Fig. 2, 3, 4 and 5,

E[F 2 0 0 0 (N )]± V ar(F 2 0 0 0 (N )) and E[F 2 1 1 1,d (N )]± V ar(F 2 1 1 1,d (N )
) are represented as functions of N when M = N, 2N, 3N, 4N (or equivalently

L = 1, 2, 3, 4).
Concerning the variance of the square of the fluctuation function, it becomes smaller and smaller when L increases. This is consistent with the expression obtained in (35). One can also look at the influence of the order d. For instance, for N = 110, the larger the order, the smaller the variance.

A comparison between the values of log(E[F 2 0 0 0 (N )] + V ar(F 2 0 0 0 (N ))) and log(E[F 2 1 1 1,1 (N )] + V ar(F 2 1 1 1,1 (N ))) as functions of log(N ) is finally shown Fig. 6 to illustrate all the steps of the FA and the standard DFA.

Case of a Gaussian moving average process or a first-order autoregressive process

Let us now address the case of a w.s.s. real Gaussian q th -order moving average process, whose k th sample, denoted as y(k), is defined as follows: where u(k) is the k th sample of the driving process, assumed to be white, Gaussian, zero-mean with variance σ 2 m . b 0 = 1 and {b j } j=0,...,q denote the MA

y(k) = q j=0 b j u(k -j) (36)

parameters.

As done for the white noise, let us recall some statistical properties of the MA process. It is a zero-mean process since:

E(y(k)) = q j=0 b j E(u(k -j)) = 0 (37)
Its correlation function R yy (τ ) = R yy (-τ ) satisfies for τ ≥ 0:

E(y(k)y(k -τ )) = E(y(k)y(k + τ )) =      σ 2 u q j=τ b j b j-τ for τ ≤ q 0 otherwise (38)
As the MA process is zero-mean, its third-order moment6 M 3,y (τ 1 , τ 2 ) is equal to the third-order cumulant C 3,y (τ 1 , τ 2 ) of the process y. As the process is Gaussian, it is equal to 0. Similarly, its fourth-order moment M 4,y (τ 1 , τ 2 , τ 3 ) can be expressed from the fourth-order cumulant C 4,y (τ 1 , τ 2 , τ 4 ) of the process and the correlation function as follows:

M 4,y (τ 1 , τ 2 , τ 3 ) = E(y(k)y(k + τ 1 )y(k + τ 2 )y(k + τ 3 )) = C 4,y (τ 1 , τ 2 , τ 4 ) (39) + R yy (τ 1 )R yy (τ 2 -τ 3 ) + R yy (τ 2 )R yy (τ 3 -τ 1 ) + R yy (τ 3 )R yy (τ 1 -τ 2 )
As the process is Gaussian, the fourth-order cumulant is null and the fourthorder moment of the MA process reduces to:

M4,y(τ1, τ2, τ3) = Ryy(τ1)Ryy(τ2 -τ3) + Ryy(τ2)Ryy(τ3 -τ1) + Ryy(τ3)Ryy(τ1 -τ2) (40)
where the correlation function is defined from the MA parameters and the variance of the driving process in (38). Usually, the lags τ 1 , τ 2 and τ 3 are ordered, but this does change the result if this is not the case.

Let us now recall the case of a w.s.s. real Gaussian 1 st -order autoregressive process defined by: where a 1 is the AR parameter whose modulus is strictly smaller than 1.

y(k) = -a 1 y(k -1) + u(k) (41)
Due to the statistical property of the driving process u, E(y(k)) = 0. In addition, its correlation function satisfies:

R yy (τ ) = (-a 1 ) |τ | 1 -a 2 1 σ 2 u (42)
As the process is zero-mean and Gaussian, its third-order moment M 3,y (τ 1 , τ 2 )

is equal to 0. The fourth-order moment also satisfies (40) and hence can be expressed from the correlation function (42).

As the mean of the square of the fluctuation function has been already addressed in a recent paper [START_REF] Grivel | Dfa-based abacuses providing the hurst exponent estimate for short-memory processes[END_REF] and for the sake of space, let us focus on how to get its variance, without any approximation.

How to get the variance of the square of the fluctuation function?

Let us express the three terms of the sum defining V ar([F 2 1 1 1,d (N )]) in [START_REF] Sikora | Probabilistic properties of detrended fluctuation analysis for gaussian processes[END_REF] given the properties presented above. Firstly, T r 2 (γ 1 1 1,d R Y,N ) can be expressed from the correlation function. Moreover, 1 L T r(γ

1 1 1,d γ T 1 1 1,d R Y ) = 1 L T r(γ 1 1 1,d γ T 1 1 1,d E[Y l Y T l ]) and E[Y m Y T l ]
) depend on the correlation between vector Y m and Y l , m and l being able to equal or not. It is hence defined from the fourth-order moment given in (40) and consequently on products between different values of the correlation function. Thus, as an example, when l = m and N = 3, one has:

RY =              
M4,y(0, 0, 0) M4,y(1, 0, 0) M4,y(2, 0, 0) M4,y(0, 1, 1) M4,y(0, 1, 2) M4,y(0, 2, 2) M4,y(0, 0, 1) M4,y(1, 0, 1) M4,y(2, 0, 1) M4,y(1, 1, 1) M4,y(1, 1, 2) M4,y(1, 2, 2) M4,y(0, 0, 2) M4,y(1, 0, 2) M4,y(2, 0, 2) M4,y(2, 1, 1) M4,y(2, 1, 2) M4,y(2, 2, 2) M4,y(0, 1, 1) M4,y(1, 1, 1) M4,y(2, 1, 1) M4,y(0, 0, 0) M4,y(1, 0, 0) M4,y(0, 1, 1) M4,y(0, 1, 2) M4,y(1, 1, 2) M4,y(2, 1, 2) M4,y(0, 0, 1) M4,y(1, 0, 1) M4,y(1, 1, 1) M4,y(0, 2, 2) M4,y(1, 2, 2) M4,y(2, 2, 2) M4,y(0, 1, 1) M4,y(1, 1, 1) M4,y(0, 0, 0)

              (44)
The above matrix is symmetric. The block located on the main diagonal has even a Toeplitz structure as the order of the lags does not have an influence.

The reader will understand that, as done for the white noise, one could combine all the equations that have been presented in this section to get an analytical expression of the variance V ar[F 2 1 1 1,d (N )], but the latter would be quite long. The matrix form remains more convenient.

Conclusions

Evaluating the variance of the square of the fluctuation function appearing in the DFA is of interest as the estimation of the Hurst exponent is more reliable when the variance of the square of the fluctuation function is small. The matrix formulation we introduced in [START_REF] Berthelot | Alternative ways to compare the detendred fluctuation analysis and its variants. application to visual tunneling detection[END_REF] and used in [START_REF] Grivel | Dfa-based abacuses providing the hurst exponent estimate for short-memory processes[END_REF] brings in much to derive it as there is no approximation. Then, we chose to illustrate our reasoning with a white noise. When comparing the DFA and the FA, one can see that the variance is much smaller with the DFA. We also show how to deal with moving average process and autoregressive process. We can see that the variance depends on the correlation function of the process when the process under study is Gaussian and zero mean. Finally, the standard DFA was considered, but the same reasoning could be done for the higher-order DFA.

  2/ decimating the resulting signal by a factor N , estimating the variance of the N sequences that are obtained and averaging the estimated variances to get what the authors usually call the square of the fluctuation function 3/ combining the values of the fluctuation function computed for different values of N to deduce the Hurst exponent.

Figure 1 :

 1 Figure 1: State of the art: three families of the DFA and its variants

d=

  N H y(t), where d = means an equivalency in terms of finite joint distribution and H is the self-similarity parameter or the Hurst exponent.

  with θ l,d = a l,0 a l,1 ... a l,d T storing the (d+1) polynomial parameters and A l,d a N × (d + 1) matrix whose c th column is defined by {[(l -1)N + n] c-1 } n=1,...,N , with c = 1, ..., d + 1. By introducing Θ 1 1 1,d = θ T 1,d . . . θ T L,d T and the (LN × (d + 1)L) block diagonal matrix A 1 1 1,d defined from the set of matrices {A l,d } l=1,...,L , the parameters of 4 The first LN ≤ M samples of the profile are processed, with L the integer part of the ratio M N .

Step 4 .

 4 Steps 2 and 3 are repeated for different values of N . As F 2 1 1 1,d (N ) is proportional to N 2α with α = H + 1 due to the signal integration, log(F 2 1 1 1,d (N )) is plotted as a linear function of log(N ). The final step consists in searching the straight line fitting the log-log representation. Its slope is estimated in the LS sense. Depending on the number of samples that are available, the practitioner can select the values of N he can use to deduce the slope. Note that H is a non linear function of the signal due to the logarithm function.

  to 1, starting by the element equal to y 2 ((l -1)N + b l ) and ending by the element equal to y((l -1)N + b l )y(lN ) with b l = 1, ..., N . In this case, γ T 1 1 1,d

2γ 1 1 1

 1 ,d[START_REF] Peng | Long-range correlations in nucleotide sequences[END_REF][START_REF] Peng | Mosaic organization of DNA nucleotides[END_REF] ... 2γ 1 1 1,d (2, N ) ... 2γ 1 1 1,d (N, N )

25 ) 5 . 1 . 1 .

 25511 Statistical mean and variance of the square of the fluctuation function Since the Hurst exponent is equal to -1 2 , F 2 0 0 0 (N ) should be proportional to N 2(H+1) = N . Similarly, F 2 1 1 1,d (N ) should be proportional to N . Thus, for this particular random process, when the square of the fluctuation functions are represented as functions of N , a straight line passing by the origin should be obtained, no matter what the value of the slope can be 5 .

Figure 2 :

 2 Figure 2: Evolution of E[F 2 0 0 0 (N )] ± V ar(F 2 0 0 0 (N )), i.e. the statistical mean of the square of the fluctuation function in the FA case, plus or minus its standard deviation, when N increases

Figure 3 :

 3 Figure 3: Evolution of E[F 2 1 1 1,1 (N )] ± V ar(F 2 1 1 1,1 (N )), i.e. the statistical mean of the square of the fluctuation function in the 1st-order DFA case, plus or minus its standard deviation, when N increases. The different colors and symbols represent different values of the length M : M = N, 2N, 3N, 4N etc.

Figure 4 :

 4 Figure 4: Evolution of E[F 2 1 1 1,2 (N )] ± V ar(F 2 1 1 1,2 (N )), i.e. the statistical mean of the square of the fluctuation function in the 2nd-order DFA case, plus or minus its standard deviation, when N increases. The different colors and symbols represent different values of the length M : M = N, 2N, 3N, 4N etc.

Figure 5 :

 5 Figure 5: Evolution of E[F 2 1 1 1,3 (N )] ± V ar(F 2 1 1 1,3 (N )), i.e. the statistical mean of the square of the fluctuation function in the 3rd-order DFA case, plus or minus its standard deviation, when N increases. The different colors and symbols represent different values of the length M : M = N, 2N, 3N, 4N etc.

Figure 6 :

 6 Figure 6: Comparison between log(E[F 2 i i i (N )] + V ar(F 2 i i i (N ))) plus its standard deviation, when log(N ) increases. The different colors and symbols represent the different approaches (FA and DFA d = 1) and the values of the length M considered for the DFA: M = N, 2N, 3N .

  To simplify the way to present the result we propose to use the second definition[START_REF] Crato | Statistical properties of detrended fluctuation analysis[END_REF] of the vector Y l . In this case, at the c th row of the b th l block of the row vector Y T l , the element is defined by y((l -1)N + b l )y((l -1)N + b l + c -1), where b l = 1, ..., N and c = 1, ..., N -b l + 1. Similarly, Y m is defined byy((m -1)N + b m )y((m -1)N + b m + r -1), where b m = 1, ..., N and r = 1, ..., N -b m + 1. Therefore, E[y((m -1)N + b m )y((m -1)N + b m + r -1)y((l -1)N + b l )y((l -1)N + b l + c -1)] is given by: (r -1, (l -m)N + b l -bm, (l -m)N + b l -bm + c -1) if l > m M4,y(c -1, (m -l)N + bm -b l , (m -l)N + bm -b l + r -1) if l < m M4,y(c -1, bm -b l , bm -b l + r -1) if l = mand bm ≥ b l M4,y(r -1, b l -bm, b l -bm + c -1) if l = m and b l > bm (43) It should be noted there are different cases because the minimum index among (m -1)N + b m , (m -1)N + b m + r -1, (l -1)N + b l and (l -1)N + b l + c -1 has to be searched to define the three lags characterizing the fourth-order moment. Thus, when l > m, the smallest value is necessary (m -1)N + b m . When m > l, the smallest value is necessary (l -1)N + b l . When m = l, this depends on b l and b m . Combining (38) and (40) for the MA process and (42) and (40) for the AR process, one can express the element of the block matrix E[Y m Y T l ]) located at the r th row and c th column of the block of dimension of (N -b r ) × (N -b c ).

If the deterministic part was the trend, this would correspond to the statistical mean of the square of the fluctuation function.

It is consistent with the results obtained by Höll[START_REF] Höll | The relationship between the detrendend fluctuation analysis and the autocorrelation function of a signal[END_REF] when the number of the samples tends to infinity provided the signal is w.s.s. and ergodic.

The lags {τ i } i are assumed to be positive and usually ordered.