Biases in the Simulation and Analysis of Fractal Processes

Fractal processes have recently received a growing interest, especially in the domain of rehabilitation. More precisely, the evolution of fractality with aging and disease, suggesting a loss of complexity, has inspired a number of studies that tried, for example, to entrain patients with fractal rhythms. This kind of study requires relevant methods for generating fractal signals and for assessing the fractality of the series produced by participants. In the present work, we engaged a cross validation of three methods of generation and three methods of analysis. We generated exact fractal series with the Davies–Harte (DH) algorithm, the spectral synthesis method (SSM), and the ARFIMA simulation method. The series were analyzed by detrended fluctuation analysis (DFA), power spectral density (PSD) method, and ARFIMA modeling. Results show that some methods of generation present systematic biases: DH presented a strong bias toward white noise in fBm series close to the 1/f boundary and SSM produced series with a larger variability around the expected exponent, as compared with other methods. In contrast, ARFIMA simulations provided quite accurate series, without major bias. Concerning the methods of analysis, DFA tended to systematically underestimate fBm series. In contrast, PSD yielded overestimates for fBm series. With DFA, the variability of estimates tended to increase for fGn series as they approached the 1/f boundary and reached unacceptable levels for fBm series. The highest levels of variability were produced by PSD. Finally, ARFIMA methods generated the best series and provided the most accurate and less variable estimates.


Introduction
e repeated measurement of physiological or behavioral events (stride durations, heartbeat intervals, and intertap intervals) is typically characterized by a marked variability. For a long time, this variability has just been considered a random perturbation, without any functional signification. However, a number of authors, during the last two decades, showed that these biological series presented a typical longrange correlational structure over time and especially a statistical self-similar (fractal) pattern [1][2][3][4]. Fractal processes have recently received a growing interest, especially in the domain of rehabilitation. More precisely, the evolution of fractality with aging and disease, suggesting a loss of complexity [5], has inspired a number of studies that tried, for example, to entrain patients with fractal rhythms [6][7][8].
In this domain, authors are confronted with two main methodological problems: e first one concerns the evaluation of the level of long-range correlations in physiological series. A number of different methods have been proposed, and their respective qualities were systematically assessed in comparative studies [9][10][11][12]. e second one refers to the generation of exact fractal signals necessary of providing experimental devices (metronomes and virtual environments) with controlled long-range correlation properties. A number of methods have been proposed for simulating such series [13][14][15]. e assessment of estimation methods, on the one hand, and simulation methods, on the other hand, raises a typical problem of circularity, performances, and biases in the former being analyzed on the basis of the latter, and vice versa [11]. When a bias is identified, it remains difficult to attribute the problem to the method of simulation or to the method of assessment. In order to overcome this problem, we propose in the present paper a cross-validation study, combining three methods of generation and three methods of analysis.
We first propose a formal introduction of the three main domains of definition of long-range correlated processes: the fractional Brownian motion framework [16], the spectral domain [17], and the autoregressive fractionally integrated moving average (ARFIMA) processes [18].

e fBm/fGn Model.
e fractional Brownian motion (fBm) denoted as B H (t) is a mathematical model of continuous stochastic process introduced by Mandelbrot and Van Ness [16] as a generalization of the Brownian motion, where increments do not need to be independent. B H (t) is characterized by the Hurst parameter (H), which can take any real value within the interval ]0, 1[. is value gives information about the nature and the strength of the correlation between successive increments in the process. If H is below 0.5, B H (t) is underdiffusive, and its increments are anticorrelated. In contrast, if H is above 0.5, B H (t) is overdiffusive, and its increments are positively correlated. In the case that H equals 0.5, B 0.5 (t) is an ordinary Brownian motion (normal diffusion), and its increments are an uncorrelated Gaussian white noise.
fBm has some fundamental properties. e first one is that its variance grows as a power function of the length of the time interval observed, with an exponent 2H: (1) e second is that fBm is a fractal process, characterized by statistical self-similarity: By definition, B H (t) is fully described by its autocovariance function c B H (t, s): As previously indicated, the successive increments of B H (t) can be correlated and the Hurst exponent informs about the nature of this correlation, and thus, the derivative of B H (t) should be a stationary correlated noise. However, this derivative cannot be calculated because in theory Brownian motion describes an infinitely broken continuous trajectory. In other words, nondifferentiability is a fundamental property of B H (t). One can still estimate this derivative using a discrete version of B H (t) where increments are defined on a time interval m. is estimate is a discrete process called fractional Gaussian noise (fGn), denoted as G m (i): As for the fBm, the fGn is fully described by its autocovariance function, which is easily derived from equation (3) with m � 1: 2.2. e Spectral Model. Stochastic fractal processes can also be defined in the frequency domain, on the basis of a scaling law that relates power (i.e., squared amplitude) to frequency according to an inverse power function, with an exponent β [9,10]. For an fGn process with exponent H, the power spectrum can be expressed as follows [12,19]: and for the corresponding fBm process [12,20]: en the power spectrum of fGn/fBm processes has the general form: with β ∈ ]−1, 1[ for fGn and β ∈ ]1, 3[ for fBm. is suggests that fGn and fBm could be considered as a continuum, with both families being characterized in the time domain by a scaling exponent α, α � H for fGn and α � H + 1 for fBm. is assumption has been exploited by Peng et al. [4] in the conception of the detrended fluctuation analysis (see Section 3). e scaling exponent α is linearly related to the spectral exponent β over the whole fGn/fBm continuum by α � (β + 1/2) [21]. e case β � 1 defines the so-called "1/f noise," which represents the boundary between fGn and fBm processes.

e ARFIMA Model.
A third approach to long-range correlated processes is provided by the autoregressive fractionally integrated moving average (ARFIMA) models [18]. is approach is an extension of the ARIMA (for autoregressive, integrated, moving average) framework, introduced by Box et al. [22], which intended to represent a variety of short-term relationships in time series. ARIMA models are potentially composed of three components. e autoregressive component suggests that the current observation y(t) is determined by a weighted sum of the p previous observations, plus a random perturbation ε: e moving average component supposes that the current observation depends on the value of the random perturbations that affect the q preceding observations, plus its own specific perturbation: Computational and Mathematical Methods in Medicine Finally, the differencing parameter d indicates the number of differencing that should be applied to the series before modeling. An ARIMA model is the combination of these three components and can be designated by the respective orders of the three processes as (p, d, q). ARIMA modeling has been used either for generating time series with specified p, d, q parameters or for determining the best p, d, q combination for accounting for a given series [22,23].
ARIMA models could be more conveniently expressed using the so-called backshift operator, defined as e generic ARIMA (p, d, q) model can then be rewritten as where ϕ(B) and θ(B) are, respectively, the autoregressive and the moving average operators, represented as polynomials in the backshift operator: [24]. In the initial formulation of the model, the d parameter was an integer [22]. Granger and Joyeux [18] showed that it was possible to provide this model with long-range dependence properties by allowing the differencing parameter d to take on fractional values, thereby obtaining an ARFIMA model.
Here we focus on the most simple model ARFIMA(0, d, 0), which is supposed to only contain long-range correlations. Using the backshift operator notation, this model is expressed as follows [25]: with Granger and Joyeux [18] derive a filter A(B) from equations (13) and (14) and demonstrate that the process can be rewritten as . (16) where d is a measure of the intensity of long-range correlations in the series. Note, however, that ARFIMA models account only for stationary series, d being bounded within the interval ]−0.5; 0.5[. In other words, ARFIMA models remain limited to fGn series. d is related to the spectral exponent β and the scaling exponent α by the following linear equations: Each of these theoretical frameworks provided specific methods for generating fractal signals or for assessing the fractality of empirical series. In the present work, we engaged a cross validation of three methods of generation and three methods of analysis. We selected one simulation method and one estimation method in each previously presented domain of definition. Our rational is that biases that are revealed by the three analysis methods should be attributed to the generation method and conversely biases that appear whatever the generation method should originate in the analysis method.

Methods
In order to explore the whole fGn/fBm continuum, we first generated series from α � 0.1 to 0.9, by steps of 0.1, and from 1.01 to 1.9, by steps of 0.1. ese values were used in most previous similar studies [9,10,26]. Additionally, in order to analyze more closely the behavior of simulation and analysis methods around to the 1/f boundary, we generated a series from α � 0.91 to 1.09, by steps of 0.01. is range of exponents was rarely considered in the literature (for a noticeable exception, see [27]). However, this focus on the 1/f boundary seems of particular interest, because the problems of fGn/fBm classification are concentrated within this interval [10] and also because one could have some doubts about the hypothesis of continuity around the 1/f boundary [28].
We used three methods of generation: the Davies-Harte algorithm, the spectral synthesis method, and the ARFIMA simulation method. ese methods are detailed below. For each selected α value and with each method, we generated 120 series of 1024 data points. In this section, all methods are written in the discrete time and frequency domain; for reading convenience, we keep the variable t for discrete time domain and f for discrete frequency domain.

Davies-Harte Algorithm (DH).
We used the algorithm proposed by Davies and Harte [13], for generating fGn series of length N (N being a power of 2). As previously indicated, an fGn process is fully described by its autocovariance function (see equation (5)). en, one can deduce the exact spectral power S expected for this autocovariance function, from the discrete Fourier transform of the following sequences of covariance values c G1 defined by equation (5): where f � 0, 1, . . ., N − 1. It is important to check that S(f) ≥ 0 for all f. Negativity would indicate that the sequence is not valid.

Computational and Mathematical Methods in Medicine
Let W gn (f), where f � 0, . . ., N − 1, be a white Gaussian noise. e randomized spectral amplitudes, V(f ), are calculated according to the following equations: , , Finally, the first N elements of the discrete Fourier transform of V are used to compute the simulated series x(t): where t � 1, 2, . . ., N.
We first generated fGn series for H values ranging from 0.1 to 0.9, by steps of 0.1. In order to explore more precisely the performance of DFA close to the 1/f boundary, we also generated fGn series for H values ranging from 0.91 to 0.99, by steps of 0.01. A second set of fGn series was generated, for H values ranging from 0.1 to 0.9, by steps of 0.1, and was integrated for obtaining fBm series for each corresponding H value (i.e., for α values ranging from 1.1 to 1.9). Finally, we generated fGn series for H ranging from 0.01 to 0.09 and integrated them for obtaining fBm series close to the 1/f boundary (i.e., for α values ranging from 1.01 to 1.09).

Spectral Synthesis Method (SSM).
e spectral synthesis method is designed to produce fBm and fGn series x(t) based on the characteristics of the power spectral density S(f) of these signals. In other words, the idea of SSM is to generate a right kind of S(f) that gives rise to fBm or fGn with an exponent 0 < H < 1.
Since S(f ) is obtained as the square of the modulus of the Fourier transform of the signal and its proportional to f −β (see equation (8)), We need to generate a complex series X(f) with a modulus r(f) proportional to 1/f β/2 . r(f) is obtained as follows: where W gn is a white Gaussian noise and f is the frequency ranging from 1 to N/2.
To generate the X(f ), we also need to generate a random phase in radian: where W un is a white noise with uniform distribution. We can create the complex coefficients: ese coefficients are extended to the whole range of expected values (in respect of the Shannon theorem), a(f ) from (N/2) + 1 to N being equal to a(f) from (N/2) to 1 and b(f ) from (N/2) + 1 to N being equal to b(f ) from (N/2)to 1.
en, the complex series X(f ) are generated as follows: Finally, the inverse Fourier transform of this complex number is computed to obtain the time series: where t � 1, 2, . . ., N.

ARFIMA Simulation Method.
We used the Matlab code proposed by Fatichi [14], based on equation (12), and in the present case, as we limited ourselves to ARFIMA (0, d, 0) models, on equation (13). is code just bounds the summation in equation (15) to k � 100. We used this procedure for generating fGn series (d ∈ ]−0.5, 0.5[). In order to obtain the whole set of series we needed, we computed fBm series by cumulated summation. We now present the three estimation methods we used. Note that for a better readability, the exponents of the simulated series and the estimates are expressed or converted in α metrics. Estimates are denoted with a circumflex (α, β, d).

Detrended Fluctuation Analysis.
e DFA algorithm works as follows, for a series x(t) of length N where t � 1, 2, . . ., N. e series is first integrated, by computing for each t the accumulated departure from the mean of the whole series: is integrated series is then divided into k nonoverlapping intervals of length n. e last N − kn data points are excluded from analysis. Within each interval, a least squares line is fitted to the data. e series X(t) is then locally detrended by subtracting the theoretical values X th (t) given by the regression. For a given interval length n, the characteristic size of fluctuation for this integrated and detrended series is calculated by the following equation: 4 Computational and Mathematical Methods in Medicine (28) In the original algorithm, this computation is repeated over all possible interval lengths, for example, from n � 10 to n � N/2 [9]. In the present paper, we applied the evenly spaced averaged version of DFA [29], which significantly reduces the variability of estimates. is procedure consists in dividing the (log) abscissa into P bins of length (log 10 (n max /n min ))/P, starting from log 10 (n min ). e P bins are defined as follows: bin(p) � ⎡ ⎣ log 10 n min + p − 1 P log 10 n max n min ; where p � 1, 2, . . ., P. In the present analyses, we set n min � 10, n max � 512, and P � 18.
Within each bin p, the average interval length n(p) and the average fluctuation size F(n(p)) are computed. A power law is expected as e exponent estimate (α) is obtained as the slope of the double logarithmic plot of F(n(p)) as a function of n(p).

Power Spectral Density Analysis (PSD).
is method works on the basis of the periodogram obtained by the fast Fourier transform algorithm and exploits the power law given by equation (8). β is estimated by calculating the negative slope (−β) of the line relating log (S(f )) to log f.
In the present paper, we used the improved version proposed by Fougere [30] and modified by Eke et al. [10], designated as low PSD we . is method uses a combination of preprocessing operations: first, the mean of the series is subtracted from each value, and then, a parabolic window is applied-each value in the series is multiplied by the following function: where t � 1, 2, . . ., N. Secondly, a bridge detrending is performed by subtracting from the data the line connecting the first and last point of the series: where t � 1, 2, . . ., N.
e Fourier transform is applied on the modified series x w,l , and the fitting of β excludes the high-frequency power estimates (f > 1/8 of maximal frequency). is method was proven by Eke et al. [10] to provide more reliable estimates of the spectral index. For allowing a direct comparison between methods, β was then converted into α metrics by a simple linear transform (α � (β + 1)/2).

ARFIMA Modeling.
e d parameter was estimated by fitting an ARFIMA(0, d, 0) model to the series using Whittle approximation of the maximum likelihood estimator. We used the ARFIMA(p, d, q) estimator package for Matlab proposed by Inzelt [31] on the Matlab central file exchange platform.
As previously explained, ARFIMA modeling holds only for fGn series and d is bounded within the interval ]−0.5, 0.5[. For fBm series, Diebolt and Guiraud [32] proposed to apply ARFIMA modeling to the corresponding fGn (obtained by differentiation) and then to estimate the theoretical fractional parameter of the fBm series by adding 1 to the d value obtained from the fGn. is strategy, however, requires an a priori assessment of the fGn/fBm classification of series. Some solutions, based on the preliminary application of methods working indifferently on fGn and fBm, such as DFA or PSD, have been proposed but yielded important percentages of misclassification around the 1/f boundary [9,10].
In order to apply ARFIMA modeling indifferently on fGn and fBm series, as DFA and PSD, we used the following strategy: the algorithm consists in finding the best Whittle approximate of the maximum likelihood estimator by constrained optimization. According to ARFIMA properties, the output parameter is bounded, yielding the parameter value d � 0.4999 if the algorithm did not converge on the upper bound. When d o , obtained from the original series, was equal to 0.4999, the series was considered as fBm and the algorithm was applied on the differentiated time series in order to obtain d diff estimate. In this nonstationary case, d � d diff + 1 else d � d o . is parameter estimate d was then converted into α metrics by a simple linear transform (α � (2d + 1)/2).

Results
We present in Figure 1 the relationships between the true α and the mean α values, for each simulation and estimation method. If we consider the global shape of results around the identity line, it seems obvious that estimation methods performed better when series were generated by the method belonging to the same domain (see the three graphs ranging along the descending diagonal). is was already observed by Eke et al. [11]. e first row depicts the results obtained by the three estimation methods with the series generated by the DH algorithm. A common bias appears in the three graphs, with a strong underestimation of α in fBm series, close to the 1/f boundary.
is reveals a clear disruption in the original fGn/fBm model. is phenomenon was not observed when series were generated with the two other methods.
Computational and Mathematical Methods in Medicine DFA tends to underestimate α for fBm series, especially when series were generated by SSM or ARFIMA modeling. Conversely, PSD tends to overestimate α in fBm series, especially when series were generated by the DH algorithm or by ARFIMA modeling. Figure 2 represents a focus of the previous results on the 1/f boundary (i.e., from α � 0.9 to α � 1.1). e first row confirms the clear disruption, around the 1/f boundary, for series generated with the DH algorithm. Additionally, it reveals a slight overestimation in fGn series, in the vicinity of 1/f noise, for PSD and ARFIMA modeling. In contrast, series generated by SSM and ARFIMA simulation present a clear continuity around the 1/f boundary, whatever the method of estimation. DFA tends to overestimate α for series generated by spectral synthesis, but not for series simulated with ARFIMA. PSD tends to overestimate α in this particular range. Finally, ARFIMA modeling slightly overestimates α for fBm series generated by SSM.
Finally, we present in Figure 3 the variability (standard deviation) of the α samples, for each simulation and estimation method. ese graphs are presented with the same scale on the vertical axe, in order to facilitate comparisons between methods. Globally, α appears less variable in series generated by the DH algorithm or with ARFIMA simulation, than in those obtained from SSM, which yielded the worst results. Concerning estimation methods, DFA results show a clear increase in α variability as true α increases. is effect is not present with PSD and ARFIMA modeling, in which variability remains stable over the whole true α range. PSD results reveal a high level of variability, and especially when series were generated by SSM. Conversely, ARFIMA modeling is characterized by a low variability.

Davies-Harte Algorithm.
e most important observation with the series generated with the DH algorithm is the strong bias for fBm series close to the 1/f boundary. is result, consistently obtained with the three analysis methods, should obviously be attributed to a bias in the generation technique. is result was already evidenced by Stadnitski [27], who suggested that "the observed discrepancies probably occurred due to the incapability of the Davies-Harte technique to provide fBm series with H close to 0. ese results underline the importance of a proper data generation in simulation studies and indicate a revision of results from studies that employed the Davies-Harte technique" (pp. 144-145).
Delignières [28] showed that this bias was not related to the Davies-Harte algorithm, but more fundamentally to the fGn/fBm model itself. Based on the premises of the model, and especially on the autocorrelation of fGn (equation (5)), the author derived an analytical expression for the autocorrelation of fBm: where N is the length of the series. We report in Figure 4 the theoretical lag-1 autocorrelations, computed for fGn series ranging from H � 0 to 1 according to equation (5)   breakdown of the correlation properties of series around the fGn/fBm boundary, and the limit behavior of fBm, H approaching 0, is uncorrelated white noise. is suggests that fBm series, obtained by the cumulative summation of the corresponding fGn series, should in fact be fGn for very low H values. On the basis of equation (33), we computed the expected ρ fBm (1) values, for H ranging from 0.01 to 0.1, and we compared these expected values to the mean lag-1 autocorrelations observed in the series simulated with the Davies-Harte algorithm. e results are illustrated in Figure 5 and show that the correlational structure of simulated series matched closely that expected from the fGn/fBm model.

Spectral Synthesis Method.
In contrast with the DH algorithm, SSM provides continuity around the 1/f boundary.
is continuity was expected, as SSM works indifferently over the whole range of β exponents. is result could be surprising, as fGn/fBm and spectral models are often considered equivalent, representing similar properties in the time and frequency domains, respectively.
Another important result is related to the variability of estimates, which appears systematically higher in series generated by SSM than that observed with other methods of generation. is represents an important problem with this method, which seems unable to provide series sufficiently close to the expected exponents. is result is interesting, as fBm series were in that case obtained by cumulative summation of their correspondent fGn, as for the DH algorithm. Additionally, the obtained sets of series presented the lowest variability in exponent estimations, whatever the used method.
ese results suggest that one could get a good confidence in series generated by ARFIMA, as compared with the two other methods.

Detrended Fluctuation
Analysis. DFA works quite well with fGn series but presents a systematic negative bias and a high level of variability for fBm series. For understanding these poor results with fBm series, as compared with other methods, it is important to keep in mind that DFA actually works on integrated series and in this case on integrated fBm.
is family of overdiffusive processes is not well known, and the diffusion property exploited by DFA seems moderately appropriate with such series [9].
DFA is a very popular method, which presents the advantage to be indifferently applicable to both fGn and fBm. Some other methods gave satisfying results but are only relevant for fGn (e.g., the dispersional analysis or the rescaled range analysis [33]) or for fBm (e.g., the scaled windowed variance analysis [34]). However, researchers are often unable to know if their signals refer to any of these families, especially if these signals are close to the 1/f boundary. DFA allows overcoming this problem. e poor performances of DFA with fBm signals could be neglected, considering that most physiological series fall into the fGn family. In some cases, however, for example, for postural sway or gaze fluctuations, series are clearly nonstationary and should be considered as fBm. In such cases, one could consider to apply DFA on differenced series or to omit the integration step in the algorithm (see, e.g., [35,36]).
It is worth noting that the present results were obtained with the evenly spaced version of DFA, which was proven to significantly improve the accuracy and to reduce the variability of the original method [29]. Results would have been even worse with the original algorithm.

Power Spectral Density.
Despite the application of the refinements proposed by Eke et al. [10], PSD provided the worst estimation results. In terms of accuracy, PSD tends to overestimate α in fBm series, and results are particularly deceptive in terms of variability. PSD remains a quite popular method, especially because it provides appealing graphical results. e bilogarithmic representation of the power spectrum, beyond the estimation of the 1/f slope, can give essential indications about the presence and the nature of short-term fluctuations in the series (see, e.g., [37,38]). However, for the specific purpose of exponent estimation, PSD seemed clearly outperformed by other methods.
5.6. ARFIMA Modeling. ARFIMA modeling has been neglected by most previous studies comparing fractal methods [9][10][11]. Rangarajan and Ding [39] developed an integrated approach of fractal analysis associating rescaled range and spectral analyses, but they never considered ARFIMA as a possible alternative or complement. However, this method provided the most accurate and less variable estimates. As indicated in the introduction, ARFIMA modeling was initially designed for accounting for stationary series. As proposed by Diebolt and Guiraud [32], we differentiated nonstationary series before the application of ARFIMA modeling, and this method gave satisfying results. We applied a very simple procedure for distinguishing stationary and nonstationary series, considering that series yielding d estimates equals to the upper limit obtainable with the algorithm (0.49999) should be considered nonstationary.
is simple procedure provided good results, as evidenced by   the nice continuity observed around the 1/f boundary. In the initial steps of this work, we tested the ARFIMA package [40,41] for the matrix computing language Ox [42]. is algorithm, however, gave worse results around the 1/f boundary. Liu et al. [43] recently proposed a first evaluation of diverse ARFIMA programs, available on various platforms (Matlab, R, SAS, and OxMetrics), providing solutions for simulation, parameter estimation, and forecasting. Further investigations should be necessary for providing effective guidelines for the selection of the most relevant solutions.
In the present paper, we limited ourselves to the estimation of ARFIMA (0, d, 0) models. ARFIMA modeling could also account for various autoregressive and/or moving average processes that could contaminate empirical series.
is method allows isolating long-range correlations in series and then to provide a better estimation of fractal exponents. Finally, a procedure using ARFIMA modeling has been proposed for gauging the effective presence of genuine longrange correlations in empirical series [26,38,44]. Indeed, short-term processes could sometimes mimic 1/f-like fluctuations, yielding false detections of long-range correlations with classical methods such as DFA [38]. e method proposed by Torre et al. [26] allows testing the relative likelihood of various ARMA and ARFIMA models and finally to conclude to the effective presence of long-range correlations.
In order to illustrate the experimental usefulness of these results, we applied the three estimation methods to a set of series of stride durations, collected during 15 min walking bouts with two groups of participants: the first one was composed of 22 young participants (28.07 yrs ± 8.88) and the second of 23 older persons (72.36 yrs ± 4.88). ese data were collected during recent experiments performed in our lab [45,46], and each series had a length of 512 data points. As indicated previously, Hausdorff et al. [47] showed that aging was characterized by a typical loss of complexity in walking dynamics, and one could expect to evidence significant differences in α estimates between the two groups. We present in Table 1 the results obtained with the three methods, converted in α metrics for comparison. In accordance with our previous results, PSD yielded a greater variability in estimates than DFA or ARFIMA. We compared the two samples with a oneway ANOVA, which did not evidence any difference between the two groups on the basis of PSD estimates. In contrast, a significant difference was found for DFA and ARFIMA estimates, with a stronger effect size for ARFIMA.
Beyond our observations about the superiority of ARFIMA modeling, in terms of accuracy and variability of exponent estimation on exact synthetic signals, this final result shows that ARFIMA modeling provides a better statistical power than DFA or PSD in the analysis of experimental data [48]. Further investigations remain necessary, however, for assessing the respective performances of these methods, especially with shorter time series, which represent an essential challenge in human behavioral experiments [9-12, 21, 29, 49].

Conclusion
is study provides a strong support for the use of ARFIMA methods, for simulation as well as for parameters estimation purposes. As previously indicated, ARFIMA methods were not considered in recent comparative studies and rarely exploited in empirical works, at least in the physiological and behavioral domains. e accuracy and the low variability of exponent estimation with ARFIMA modeling should convince researchers of the advantages of this method, especially for detecting mean differences between groups. is method should attract the attention of researchers, either for their future experiments or for revisiting their previous results.

Data Availability
e time series used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as potential conflicts of interest.