Wavelet-based discrimination of isolated singularities masquerading as multifractals in detrended fluctuation analyses

The robustness of two widespread multifractal analysis methods, one based on detrended fluctuation analysis and one on wavelet leaders, is discussed in the context of time-series containing non-uniform structures with only isolated singularities. Signals generated by simulated and experimentally-realized chaos generators, together with synthetic data addressing particular aspects, are taken into consideration. The results reveal essential limitations affecting the ability of both methods to correctly infer the non-multifractal nature of signals devoid of a cascade-like hierarchy of singularities. Namely, signals harboring only isolated singularities are found to artefactually give rise to broad multifractal spectra, resembling those expected in the presence of a well-developed underlying multifractal structure. Hence, there is a real risk of incorrectly inferring multifractality due to isolated singularities. The careful consideration of local scaling properties and the distribution of H\"older exponent obtained, for example, through wavelet analysis, is indispensable for rigorously assessing the presence or absence of multifractality.


Introduction
The concept of multifractality, whereby not just one but an extended spectrum of exponents is required to account for the dynamics of a system, represents one of the pillars of complex signal analysis [1][2][3][4][5]. The term was coined in the context of fully-developed turbulence [6] and mathematically formalized, with the multifractal or singularity spectrum as principal characteristic, by Halsey et al. in the year 1986 [7]. Within the framework of multifractal formalism, a function is decomposed into subsets, which are characterized by a Hölder exponent α and a fractal dimension f (α) [7]. The identified set of Hölder exponents provides explicit information about the regularity of a time-series. When only one type of regularity is present in a signal, its statistical properties are quantified by a single Hölder exponent; this results in the convergence of the multifractal spectrum to a single point. A typical example of such a monofractal process is fractional Brownian motion (fBm) with a homogenous organization of the fluctuations. On the other hand, multifractal signals, related to intermittency phenomena and correlations heterogeneity, are characterized by an extended set of Hölder exponents and developed multifractal spectrum resembling an inverted parabola [8]. A representative example is the binomial cascade generated through the iterative and multiplicative procedure. Thus, the multifractal methodology offers an opportunity to distinguish between signals characterized by the same autocorrelation function or power spectrum, but having a different underlying organization. The pervasiveness of emergent fractal and multifractal structures has made it possible to apply this methodology fruitfully across diverse areas such as physics [9,10], biology [11][12][13], physiology and neuroscience [14,15], chemistry [16,17], economy [18][19][20][21][22][23][24], even linguistics [25,26] and music [27,28].
The study of the multifractal properties of time-series can be approached in two ways. The most common approach involves estimating the spectrum through global procedures that neglect the precise location of the singularities. These procedures, as we show in our study, work well for fractal structure with a well-developed hierarchy of singularities and densely roughness of the signals. An alternative one, based on assessing singularity locally, offers the possibility to determine also the location of the Hölder exponents. This opens up the possibility of quantifying the singular behavior of the signal whenever singularities appear isolated, e.g., when analyzed processes reveal local smoothness or outliers due to measurement error. Although this approach provides more information about the fluctuation organization, it is considerably more numerically unstable, and as such, rarely applied in practical time-series analysis. However, as we clearly demonstrate in our study, this is possibly the only way of distinguishing artefactual (or apparent) from genuine multifractality.
Several methods have been devised for numerically estimating the multifractal spectrum of a process given a simulated or recorded signal (in the present context, equivalent to time-series) [29][30][31]. The prevalent algorithms encompass the multifractal detrended fluctuation analysis (MFDFA) [32], the wavelet transform modulus maxima [33], and its development wavelet leaders (WL) [34]. Due to their numerical stability and usual accuracy, these two methods are commonly recognized as the most reliable means of estimating the multifractal spectrum [35][36][37]. Their accuracy has been repeatedly demonstrated via synthetic time-series having known multifractal properties [35]. Moreover, the results obtained with these two methods can be used as a cross-test of the validity of the assessed multifractality, since MFDFA and WL are based on different numerical methods of quantifying the multifractality. The former hinged around the scaling properties of the variance, whereas the latter is grounded on the wavelet transform for decomposing the signal and characterize its self-similar properties. Thus, in our study, we applied both approaches as complementary methodologies. Albeit distinct, these methods are all based on the well-known q-filtering technique, which decomposes a time-series into subsets predicated on the fluctuation amplitudes. As such, they similarly provide information about the average level of fractality across all time-series segments. While their practical usefulness is beyond question, the results should be interpreted cautiously due to the inherent complexity of both the signals under analysis and of the al-gorithms themselves. Here, several compelling examples are given of how even elementary systems can yield signals for which a naive interpretation of the multifractal spectrum, obtained via both the MFDFA and the WL, leads to completely flawed conclusions. To this end, consideration is given to simulated and experimentally-recorded time-series from the Saito chaos generator, a simple four-dimensional non-linear dynamical system with a strong hysteretic component, as well as to synthetic signals which by construction cannot possess any multifractal structure. The results show that relying only on the multifractal spectrum width as an indicator of multifractality can be profoundly misleading. At a minimum, other spectral attributes such as the asymmetry have to be taken into account, where the key is considering the local scaling properties and the distribution of Hölder exponents. This paper is organized as follows. In Sec. II, the notion of local and global Hölder exponents is introduced, together with the methods of MFDFA and WL. In Sec. III, a selection of examples of true multifractality is firstly introduced. Then, the case of a hysteretic oscillator is considered to exemplify the consequences of isolated singularities on both MFDFA and WL, followed by the presentation of paradigmatic synthetic signals leading to similar issues. Finally, in Sec. IV, practical recommendations for the proper application of these analyses are offered.

Hölder exponents
Analyzing the regularity of a signal provides essential insight into its statistical properties and possible underlying geometrical structure. To characterize the local singular properties of a time-series, then, the point-wise Hölder exponent α can be considered. Given a function f , for each x 0 ∈ the same is defined as follows [38]: where f belongs to the Hölder space Values of the Hölder exponent approaching zero indicate increasing irregularity of the function f ; conversely, larger values of α denote more regular fluctuations. The multifractal formalism is a statistical description of functions through quantifying their distribution of the point-wise Hölder exponents, which is naturally extended to continuous signals and discretized time-series.

Multifractal detrended fluctuation analysis
Multifractal detrended fluctuation analysis (MFDFA) [32] is a method for detecting and quantifying the scaling properties of time-series which is widely applied across diverse areas of experimental and computational science [39][40][41][42][43][44]. It comprises multiple steps, which may be summarized as follows. Let us consider a time-series x i having length N , i = 1, 2...N and, as a first step, calculate its profile according to where x denotes the mean of time-series x i . Since fractality manifests as patterns which are self-similar across different temporal scales, the profile has to be analyzed over segments of different length. Thus, the time-series is next subdivided into N s non-overlapping segments ν of length s (N s = int(N/s)) starting from the beginning. However, since the length is not necessarily an integer multiple of the scale s, the procedure is also repeated starting from the end, yielding a total of 2N s segments. To remove possible trends in the timeseries, which can distort the results, in each segment ν a polynomial of order m (P (m) ν ) is fit and subsequently subtracted from the data. The effectiveness of this detrending step strongly depends on the polynomial order, and it is generally agreed that small values of m provide the most reliable results [45]. Here, no statistically discernible differences were found between results obtained for m = 2, 3 and 4, thus, for brevity, only results assuming m = 4 are given.
As a next step, the detrended variance is calculated within each segment according to Then, in order to quantify the fractal properties of the signal with respect to the amplitude, the q-order filtering technique is applied, obtaining the q-order fluctuations function where q operates as a filter which discriminates fluctuations based on their amplitude; more precisely, negative and positive settings of q respectively emphasize small and large changes. Fractality in a time-series manifests itself as power-law behavior of F q (s) over different scales, that is, where h(q) denotes the generalized Hurst exponent. Hence, h(q) represents the fractality of the fluctuations selected by a given setting of q. For monofractal time-series, h(q) is constant and equals the Hurst exponent h(q) = H [46,47]. This can be used to classify time-series with respect to linear correlations. Namely, H > 0.5 indicates persistent dynamics (i.e., positive long-range correlation), whereas for H < 0.5 a signal is anti-persistent (i.e., a tendency to reverse is observed); on the other hand, H = 0.5 denotes the absence of any linear correlation. For multifractal signals, h(q) is a decreasing function of q, and the Hurst exponent is retrieved at h(q = 2) = H. Thus, for better visualizing the results and interpreting the spectrum of the generalized Hurst exponents, the same can be converted into the multifractal spectrum via the Legendre transform of the scaling function τ (q) = qh(q)−1, or directly through where α is the Hölder exponent, and f (α) refers to the fractal dimension of the data supported by a particular α. The intensity of multifractality, and thus the degree of signal complexity, is often quantified through the width of the multifractal spectrum, that is, ∆α = α max − α min . The larger ∆α, the more developed a multifractal structure is deemed to be. Another important feature of the multifractal spectrum is its asymmetry. For the paradigmatic case of the binomial cascade, a mathematical multifractal, f (α) resembles a symmetric inverted parabola [48]. However, for real-world time-series, the spectrum is often asymmetric, having one side better developed than the other; this stems from a heterogeneous organization of the signal fluctuations across scales. Hence, through quantifying the spectral asymmetry, one can retrieve critical information about the temporal organization of a time-series. The asymmetry parameter is defined as [49] A where ∆α L and ∆α R stand, respectively, for the distances between the spectral maximum and the smallest and largest values of α. In turn, the degree of the asymmetry is quantified as |A α |, whereas the sign indicates the asymmetry direction. A positive value of A α hallmarks a leftwards-stretched spectrum and denotes a well-developed fractal organization of the large fluctuations, while smaller ones are governed by simpler dynamics. Contrariwise, for negative A α the spectrum is stretched towards the right, denoting more complex behavior of the small fluctuations compared to the larger ones.

Multifractal analysis based on wavelet leaders
Another class of techniques for estimating the multifractal characteristics of a non-stationary time-series is based on the wavelet transform [50,51]. According to these techniques, a signal is decomposed into the elementary space-scale wavelet coefficients by means of a family of functions stemming from a basic function, the so-called mother wavelet. By scaling and translation of the mother wavelet ψ a,s (x) = s −1/2 ψ( x−a s ), one can obtain a decomposition of the signal at each scale s corresponding to a frequency band, separately for all time-points a (a, s ∈ , s > 0). The wavelet transform of a function f (x) is defined as [52] W f (a, s) Importantly, visualization of the resulting wavelet spectrum W f (a, s) on the scale-time plane promptly reveals the skeleton of the hierarchical structure of the process being analyzed. The choice of the mother wavelet is dictated by it being well-localized in both the time and frequency domains (derivatives of a Gaussian function are often used as mother wavelets). A crucial property of the wavelet transform is its close relation with the Hölder exponent α [50], Hence, a local singularity α(x 0 ) can be characterized by the scaling behavior of the wavelet transform around the point x 0 . Moreover, the maxima of the wavelet transform produce maxima lines in space-scale half-plane, which converge towards loci of singularity. Thus, by retrieving the power-law behavior of the wavelet transform coefficients along these lines, one can estimate the Hölder exponents, and in turn, quantify the singularity strength [33]. Due to the instability of the canonical wavelet-based multifractal methods whenever a large number of coefficients are close to zero, and due to its insensitivity to oscillating singularities [53], the notion of wavelet leaders (WL) was introduced [34,51]. For a discrete scale parameter s j = 2 −j and time a j,k = 2 −jk (j, k ∈ Z), the signal can be recovered via the formula [54,55] f where the wavelet coefficient c j,k is given by In this study, the Daubechies wavelet with 4 vanishing moments was used [56]. The wavelet leader of x 0 at the level j denotes the largest wavelet coefficient among those existing in the spatial neighborhood of x 0 at finer scales [57]. Formally, for the dyadic interval where 3λ j,k (x 0 ) = λ j,k−1 ∪ λ j,k ∪ λ j,k+1 = [2 j (k − 1, 2 j (k + 2))) and contains x 0 . For a given scale 2 j , one can define structure functions S(q, j) based on the q-th order average of the leaders where q is a real number, and Λ j is a set of dyadic intervals at scale j.
Power-law behavior of the structure function in the limit of small scales S(q, j) ≈ C q 2 jζ(q) , (2 j → 0) is a manifestation of scale invariance. Thus, ζ(q) determines the scaling exponents, and can be numerically estimated by means of a log-log regression. Since the ζ(q) function is necessarily concave [58], the Legendre transform can be used to estimate the multifractal spectrum according to the formula

Global vs. local Hölder exponents
The multifractality of a time-series manifests itself through sets of nontrivial Hölder exponents, which quantify the local variation in its irregularity [59]. These exponents may be collectively quantified by means of a "global" measure, obtained from the multifractal spectrum in (Eq. (7) or (Eq. 15) and denoted as α G , or directly through the analysis of local scaling properties of the signal by means of Eq. (10) and denoted as α L . Ideally, the two approaches should give consistent results. However, as demonstrated below, the multifractal analysis of complex time-series has limitations which, under certain circumstances, yield misleading signatures of multifractality.

Surrogates
To assess the statistical validity of the results, additional analyses were performed on surrogate time-series for all scenarios under consideration. Two commonly-used surrogate sets were generated. One relies on randomization of the Fourier phases and, as such, preserves only the spectral amplitudes while obliterating all non-linear inter-dependencies [60]. The other involves randomly permuting the time-points, destroying all temporal correlations while preserving the value distribution [21].

Examples of truly multifractal time-series
For comparison with the particular cases considered below, representative instances of real multifractals having diverse properties are firstly presented, based on analyses conducted over the range q ∈ [−4, 4] [61]. To this end, two mathematical multifractals are considered, namely the binomial cascade and the chaotic metronome derived from the Ikeda map [62], together with several real-world time-series: the inter-beat intervals extracted from electrocardiographic signals (103885 data points), the sentence length variability of the "Finnegans Wake" book by James Joyce, the logarithmic returns of the American stock market index S&P500 (7440 data points), and the sunspot number variability (43495 data points) [11,21,26,35,49]. In all these cases, the multifractal spectrum f (α G ) assumes the shape of a wide inverted parabola, spanning ∆α G > 0.2, indicating a multifractal organization of the data (Fig. 1, left). Yet, the spectra develop different degrees of asymmetry. For the binomial cascade, the inter-beat intervals, and the sentence length variability, the spectra appear almost symmetrical (A α ≈ 0), which suggests a homogeneous distribution of the correlations over small and large fluctuations. On the other hand, for stock market data and sunspot number variability, the asymmetry is, respectively, positively-and negatively-skewed.
Thus, multifractality of the S&P500 price variation is mainly the effect of a complex organization of the large fluctuations, whereas the arrangement of small fluctuations is primarily responsible for multifractality in the timeseries of sunspot numbers.
Importantly, the presence of true multifractality is confirmed, for all these cases, via analysis of the local scaling properties (Fig. 1, right). Therein, a continuous distribution of the estimated Hölder exponents spans a range of α L even broader than in the multifractal spectrum, incidentally revealing the higher sensitivity of the wavelet transform on the local scaling properties compared to the global methodology, which mainly reflects the prevalent singularities in the time-series. The following cases are presented: binomial cascade, chaotic metronome derived from the Ikeda map, inter-beat intervals, sentence length in "Finnegans Wake", logarithmic returns of the S&P500 index, and sunspot number variability.

Artefactual multifractality in the Saito chaos generator
To illustrate the potential pitfalls inherent in drawing hasty conclusions solely from global measures, let us now consider the case of the Saito chaos generator, which is a four-dimensional non-linear oscillator consisting of the following dimensionless state equations [63]: Despite its simple form and low dimensionality, this system readily generates rich dynamics spanning periodicity, quasi-periodicity, chaos, and eventually hyperchaos as a function of the parameters γ, δ, ε, η, and ρ. Here, it was initially deemed of interest from the perspective of its hypothetical ability to generate signals having a truly multifractal structure; however, in the course of numerical investigation, another feature was realized to be fundamentally important for the purposes of the present work, namely, the presence of the hysteresis function h(w), which only enters the equation of the state variable w. As a consequence of it, even though all state variables conjointly participate in the temporal dynamics, x, y, z have rather smooth an activity, whereas, in the limit of ε → 0, the temporal evolution w is characterized by sudden jumps. As shall become clear, it is these discontinuities, namely the combination of slow and fast motions corresponding to the continuous manifold and sudden jumps, which may lead to a mistaken inference of multifractality. Unless indicated otherwise, the parameters were set for operation in the hyperchaotic regime, that is, γ = 1, δ = 0.94, ε = 0.01, η = 1, and ρ = 14 [63,64].
Preliminary examination revealed differences between short-and longrange temporal correlations in the simulated time-series, giving rise to a cross-over in the fluctuation functions. In particular, the multiscale characteristics revealed a strong autocorrelation only over short time scales (i.e., s < 1000), occurring alongside a monofractal organization with weak linear correlations on the larger scales. Thus, in the analyses below, the focus is on the short-range correlations, which are relevant to the search for possible multifractality.

Numerical simulations
Time-series having a length of 10 6 points were simulated given equations (16)- (17), applying the adaptive step-size Runge-Kutta (4,5) method and returnig the results at a fixed step size of ∆t = 0.1 [65]. All simulations were repeated 10 times with randomized initial conditions. Representative segments for each variable in the hyperchaotic regime are depicted in Fig.  2a. Evidently, the dynamics of x, y, z are characterized by the markedly irregular behavior characteristic of chaotic systems. However, the dynamics of w are even more complex, featuring sharp upward and downward jumps. Even though the underlying system is the same, the multifractal properties of the signals, being influenced by the presence of singularities, could then be partially dependent on the variable under consideration. This observation is confirmed by the corresponding fluctuation functions F q (s) (Fig. 2b). Therein, it is clearly visible that the functions obey power-law behavior, which is a signature of fractal organization: however, while for x, y and z the scaling is rather homogeneous, for w a pronounced heterogeneity is apparent. Moreover, for the latter the majority of fluctuation functions have a slope close to those found close for the extreme values of q, i.e., q = ±4; only a minority assume intermediate levels, a fact that already points to a more bifractal-like organization of the data rather than to a well-developed multifractal structure.
Strikingly, rather similar characteristics of F q (s) are observed in the quasiperiodic regime, with δ = 0.65 ( Fig. 2c; for brevity, results are only shown for w). Though the dynamics are profoundly different compared to the hyperchaotic regime, the heterogeneity of the fluctuations functions remains most pronounced for w, with the distribution of slopes nearly unchanged and characteristic of a bifractal structure (Fig. 2d).
The multifractal analyses for the time-series of w generated as a function of the control parameter δ are depicted in Fig. 3. The parameter was swept in δ ∈ [0.6, 1], thus allowing the system to develop a wide range of dynamical behaviors comprising both chaotic motions and closed orbits [63]. The corresponding averaged Hurst exponent H and multifractal spectrum width ∆α G as estimated through the MFDFA and WL algorithms are depicted in Fig. 3a. It is evident that the multifractal characteristics are insensitive to the qualitative features of the system dynamics. The time-series remain strongly persistent, with H ≈ 1.5, and feature a wide spectrum with ∆α G ≈ 2.25: this could, at the surface, suggest a multifractal organization. In Fig. 3b, the multifractal spectra for the hyperchaotic and quasiperiodic regimes are compared. Their shape is almost identical, with a strong left-sided asymmetry A α ≈ 0.5: importantly, this coexists with an uneven distribution of the points along the spectrum, which concentrate mainly towards its ends. Here, analysis of the local scaling revealed fundamental subtleties of the data organization. The relative frequency histogram f r of the Hölder exponents α L forms two separable peaks, whose locations coincide with high-concentration points close to the minimal and maximal values of α G identified on the multifractal spectrum.
The locations of the singularities and their "strength" were recovered, as given by the Hölder exponents, through analysis of the local wavelet transform coefficients (Fig. 4). It is well-evident that for the time-series of the w variable, in both the quasiperiodic and hyperchaotic regimes (Fig. 4a), the maxima form separate lines on the space-scale half-plane (cf. Fig. 4b), which delineate isolated singularities (Fig. 4c). In the presence of a truly multifractal geometry, the maxima would follow a tree-like structure, stemming from the self-similar organization of the fluctuations. By contrast, consideration of the locations and strength of the singularities reveals that two discrete types are present in these time-series, and related to volatile portions of the signal: one reflects instants wherein the hysteretic behavior is apparent (α L ≈ 1.2), the other reflects the local extrema of the oscillatory component (α L ≈ 3). Thus, a faithful reconstruction of the multifractal spectrum would clustered around two separate points. However, as discussed below, the q-filtering method inherently yields a concave spectrum, and isolated peaks are impossible to obtain. The artifactual result, then, is purely the product of the averaging procedures inherent in the MFDFA methodology: together with the dense sampling of the q parameter, these generate a broad spectrum of Hölder exponents, even when only isolated singularities are present in the signal. Although it is markedly stretched towards the left-hand side, with a high concentration of points towards the two limit values of Hölder exponent, the estimated spectrum resembles an inverted parabola. Instead, the scaling properties of the time-series generated by this system should be represented by a single exponent for the x, y, z variables, and by a bifractal organization for the w variable.

Experimental confirmation
To independently confirm that the results presented above stem faithfully from the dynamics of this system, an experimental version of the same was conveniently constructed using two operational amplifiers (type TL082) and a non-linearity based on two anti-parallel series Zener diodes (type BZT52-C5V1). The corresponding circuit diagram is given in Fig. 5a, where r 1 = r 2 = R 1 = R 2 = R = 10 kΩ, r o = 820 Ω, C 1 = C 2 = 3.9 nF, L 0 = 3.3 mH, L = 32 mH (two inductors in series), and U Z = 5.1 V. These component values yield γ = C 1 /C 2 = 1, ε = L 0 /(r 2 1 C 1 ) = 0.0085, η = r 1 /r 2 = 1, and ρ = r 2 1 C 1 /L = 12.2. The signal corresponding to the variable w was digitized from the physical circuit board (Fig. 5b)    In agreement with the simulations, apparent multifractality only arises for the variable w. In Fig. 6, the time-series in the quasiperiodic and hyperchaos regimes are shown alongside the corresponding fluctuation functions. The heterogeneity of the latter is equally apparent in both cases, suggesting that the multifractal properties of the signals are only weakly dependent on the dynamics. The multifractal spectra are shown in Fig. 7: they are wide (∆α G > 2) but, in contrast to the simulations, more symmetric (A α ≈ 0.3) (cf. Fig. 7b). On the other hand, analysis of the local scaling properties reveals a singularity organization comparable to the simulations (cf. Fig. 7c). Thus, the variability of the Hölder exponents for the experimental signals is higher; however, the histogram still forms two separable clusters (albeit more dispersed than in the simulations), wherein high values reflect local signal maxima and small ones correspond to the sudden jumps. A possible explanation for the more diverse distribution of the singularities could be sought in the tolerances and non-ideal behaviors of the electronic components (e.g., finite quality factor and self-resonance of the inductors, smooth response of the diodes, etc.), which knowingly give rise to richer dynamics. Altogether, these results confirm the above conclusions, reassuring that they are not an artifact of the numerical integration.

Artefactual multifractality in the Rössler system
The next example of dynamical system that we consider is the Rössler system. Its dynamics are governed by the following system of three differential equations [67] Similarly to the Saito generator, the Rössler system reveals rich dynamics spanning periodic and chaotic behaviors. Notably, its dynamical properties depend on state equations without a hysteresis element. They are controlled by parameters that were set to a=0.3, b=0.2, c=5.7, knowingly realizing chaotic behavior with an intermediate level of folding. To ensure statistical reliability, we generated time-series having a length of 10 6 points, a representative fragment of which is visible in Fig. 8a. The distribution of the fluctuation functions already suggests a bi-fractal organization (Fig. 8b) [68]. However, the multifractal spectra estimated through the MFDFA and WL algorithms resemble typical multifractal characteristics, with a strong left-sided asymmetry (Fig. 8c). This is in contrast with the local scaling properties (Fig. 8d). The histogram of local Hölder exponents forms, similarly to the Saito generator case, two separate peaks that concentrate in the vicinity of the extreme α G values. Thus, the underlying structure reflects isolated singularities rather than a unified multifractal organization.

Further examples of artefactual multifractality in synthetic signals
The results presented above suggest that singular behavior in the Saito chaos generator can be quantified through just two scaling exponents. Moreover, the subsets corresponding to different singularities index separate components of the time-series, rather than constituting hierarchically-interwoven structures, which are the hallmark of true multifractality. Thus, a naive interpretation of the spectrum width ∆α G as a signature of multifractality can be faulty. To highlight this issue even more clearly, we finally consider processes that are, by construction, not multifractal. Yet, the methods based on q-filtering, namely the MFDFA and WL, yield misleadingly wide multifractal spectra in these cases.
As an instructive example, results from the multifractal analysis of the Lévy process, which possesses a well-recognized bifractal structure, are firstly presented. The multifractal spectrum of the Lévy time-series consists of two points, whose locations are directly related to the asymptotic behavior of the distribution tail P (x) ∼ x −(α Levy +1) and are given by [35,69]: where α Levy is the Lévy index and q is q-th moment of the fluctuation function F q (s). The multifractal analysis of the Lévy time-series having a length of 50000 points (Fig. 9a) with α Levy = 1.5 is reported in Fig. 10a,b. Therein, the bifractal nature of the data is clearly visible in the histogram of Hölder exponents: two peaks, corresponding to the theoretical values, can be readily identified. The dispersion of these peaks is artifact purely due to the finite time-series length. Yet, the multifractal spectrum estimated utilizing the MFDFA and WL methods is wide (i.e., ∆α G = 0.7) and strongly left-sided asymmetrical (A α ≈ 0.38): this could lead to the faulty conclusion that the process is multifractal when, in reality, it is not. Next, a minimal-complexity arrangement which can reproduce qualitative characteristics similar to those observed in the Saito chaos generator is considered: it simply consists of the linear superposition w(t) of two related signals. One is a pseudo-periodic signal given by, e.g., u(t) = i sin 2ωt/(p i / max p) where p i = {2, 3, 5, 7, 11}. The other is a sequence of binary fluctuations v(t) = W[u(t), ξ] generated by a hysteresis operator W acting on that signal, with v ∈ [−1, 1] and hysteresis parameter ξ = 0.2 max u. Their linear combination, e.g., w(t) = v(t) + u(t)/ max u is, by definition, not hierarchically interwoven and does not obey different scaling exponents [70,71], hence, the multifractal spectrum should consist of two separable points. To test this hypothesis, 10 time-series segments each having a lengths of 10 6 points were generated, and MFDFA was performed. The average of the multifractal spectra and the histogram of the Hölder exponents are depicted in Fig. 10c,d. In this case too, the multifractal spectrum appears well-developed (∆α G = 1.5), with a strong left-sided asymmetry (A α = 0.37), hence cursory interpretation of these results might suggest a complex multifractal structure. Again, the true nature of the process is revealed by the histogram of the Hölder exponents estimated through the wavelets, which shows that only two discrete types of singularities are present in the time-series. Thus, the analyzed structure is closer to a fractal structure than to multifractality. It is worth noting that the distribution of the Hölder exponents as well as the shape of the multifractal spectrum resemble closely the results for the variable w in the Saito chaos generator; here, however, there was no underlying non-linear dynamical system.
Finally, a pseudo-multifractal process, consisting of the superposition of a fractal time-series with periodic components, is considered. At a first glance, this process resembles the multifractal time-series of sunspot number variability (cf. Fig.1). However, as demonstrated below, careful inspection of the fractal characteristics illuminates its pseudo-multifractality. To this end, fractional Gaussian noise (fGn) [4] was generated: it represents a well-known example of a stochastic monofractal structure with possible long-range correlations quantified by the Hurst exponent, which has been applied to model phenomena across various fields of science. Namely, simulations produced a fractional Gaussian noise with an arbitrarily-chosen Hurst exponent of H = 0.8 (strongly persistent behavior), for which the amplitude of the pro-  Fig. 9c). Then, the periodic function F (i) was added to this amplitude modulated noise. In our simulations, N = 10 6 , A = 0.5 and T 0 = 4000 were set. The results of the local scaling analysis, as well as the multifractal spectrum, are shown in Fig. 10e,f. Analysis of the distribution of Hölder exponents confirms that the time-series is a composition of the two independent processes having different singular behaviors. The smaller values of α L concentrate around α L ≈ H = 0.8, corresponding to the fGn component, whereas the larger ones are related to the periodic trend. A cursory analysis of the multifractal spectrum indicates heterogeneity in the scaling properties. Particularly, the width of the spectrum, together with its strong right-sided asymmetry (A α = −0.85), could be taken as the hallmark of a multifractal time-series with a well-developed hierarchy of small fluctuations; however, this is purely an artifact. To avoid mistaking pseudo-multifractality as a valid multifractal structure, the analysis of the global fractal properties should always be corroborated by consideration of the local scaling properties as given by the wavelet transform. As a last example, we analyzed an artefactually-generated stochastic process with the singularity spectrum derived analytically. In this respect, we considered the square transform of fractional Brownian motion (Fig. 11a), which represents a bi-Hlder process whose spectrum is given by the following relation [72]: In our study, we generated fBm having a length of 10 6 points with a Hurst exponent H = 0.7. Fluctuations functions F q (s) (Fig. 11b) obtained through MFDFA show non-homogeneous scaling, which suggests a multiscaling behavior of the data. This is even more clearly visible in the singularity spectrum, which reproduces a concave hull supported by the interval in the range from H to 2H (Fig. 11c). Thus, based only on the Legendre-based methodology, a flawed conclusion on the multifractal structure of the data would be drawn. However, consideration of the histogram of Hölder exponents (Fig. 11d) estimated through wavelet analysis reveals the true bi-fractal na-ture of the process, with exponents corresponding to the theoretical expectations. This leads to the conclusion that the exponents identified through MFDFA and WL, except the two extreme values, are an artifact caused by a methodological limitation and do not contain any true information about the analysed process.

Discussion
The present study was initially conceived as the search for a deterministic non-linear dynamical system, relatively low-dimensional, which could genuinely generate multifractal structures. The Saito chaos generator, representing a hysteretic oscillator which can be easily simulated numerically as well as realized experimentally in an electronic circuit, was identified as a promising candidate. Indeed, an initial investigation of its dynamics based on the commonly-accepted MFDFA technique yielded a spectrum suggestive of fully-developed multifractality; similar conclusions were drawn from the WL analysis. However, concerningly, these results could only be obtained for the dynamical variable w, which features discontinuities in the form of jumps, and not for the others. A clear indication of a flaw came from the observation of largely overlapping multifractal spectra between the hyperchaotic and quasiperiodic regimes: since in the latter there is no turbulence, this is unexpected.
Closer visual inspection of the time-series was instrumental in resolving this issue, because it demonstrated the absence of a cascaded hierarchy of singularities, fractally nested over the consecutive scales of magnification; instead, it revealed merely a sequence of isolated singularities, associated with a restricted number of distinct Hölder exponents. In the Saito chaos generator, these isolated singularities are produced by twofold dynamics of the system, namely a continuous manifold and sudden jumps. The dynamics of these components, whose organization is quantified by single scaling exponents, are not strongly interrelated: this results in a fluctuation structure devoid of hierarchical organization. An explicit step-by-step wavelet-based estimation of the same indicated values limited to two-well separated, narrow intervals ( Fig. 3): this is in stark contrast with the results given by both the MFDFA and WL algorithms, which generated a broad distribution of f (α), comprising the two Hölder exponents seen locally, but mistakenly suggesting coverage of the entire interval between these extremes. Tentatively, a combination of processes without intrinsic convolution can be indicated as a necessary condition for the emergence of artefactual multifractality related to isolated singularities. Thus, our methodology can be applied to dynamical systems revealing different forms of dynamics (cf. Fig. 7); for instance, the Saito system in the (quasi-)periodic regime is characterized by two types of behavior. One, related to the periodic component, and the other, produced by hysteresis. However, these dynamics are not hierarchically nested, which leads to a set of isolated singularities, which in turn masquerade as a uniform multifractal structure. Further synthetic signals revealed that the danger of over-evaluating the multifractal content is not confined to the dynamics of this particular oscillator: quite on the contrary, the diversity of these signals points to a potentially pervasive problem.
At the same time, for the genuine multifractals considered, such as the binomial cascade, the MFDFA and WL methods reproduced remarkably closely the explicitly-determined distribution of underlying Hölder exponents. Similarly, a good correspondence was found for the other cases of fully-developed multifractality, such as the series of heart inter-beat intervals or sentence length variability in "Finnegans Wake". These methods, then, are obviously not per-se flawed, but their results need to be interpreted much more cautiously than is generally done.
However, precise numerical determination of the singularity spectrum for experimental signals knowingly remains an open problem [73]. Due to the intricate nature of the multifractal formalism, there are no theoretical proofs of the mathematical validity of the algorithms used to estimate the singularity spectrum: these algorithms are treated simply as reasonable numerical approximations of the underlying ground truth [74]. In particular, the commonly used formulas for the multifractal spectrum have been proposed heuristically, through analogy with thermodynamic formalism, and verified numerically for specific cases based on multifractal measures [75]. Generally, the multifractal formalism yields the upper bound of the singularity spectrum as was shown in [76], but not the exact spectrum. Therefore, it does not lend itself to formal proof. The present demonstration of what cases under which the MFDFA and WL algorithms produce a flawed indication of multifractality provides essential inspiration and input for further, more mathematically rigorous inquiry into the validity of these methods.
It is important to point out that the observed disagreement between the multifractal spectra and the true distribution of the underlying Hölder exponents is not merely down to avoidable algorithmic implementation choices, but stems straight from the foundational assumptions upon which the multifractal analyses under consideration are built. The MFDFA and WL approaches represent two related advances of the Parisi-Frisch procedure [77], aiming for a compromise between reflecting the mathematical definition of the Hölder exponent and obtaining an estimate which is practically viable in terms of both computational load and stability. An averaging procedure related to the partition function (or its equivalent, cf. Eqs. (5) and (14)) and Legendre transform (cf. Eq. (15)) is central to their framework. Although numerically stable, this methodology has one serious drawback. Namely, the Legendre transform by construction imposes the concavity of the multifractal spectrum [73] a priori even when this is not true: if non-concave Hölder characteristics are considered, it provides only the upper bound of the Hölder spectrum [78]. Therefore, it is patently impossible to say without additional tests whether an observed concave spectrum is a valid representation of the data, or simply reflects the limitations of the methodology. To overcome this problem, multifractal methods not based on the Legendre transform have to be applied [37,79]. However, in many cases, these methods suffer from issues of practical implementation and computational stability. Therefore, in this work, the distribution of the Hölder exponents was instead directly analyzed. In the presence of non-concave multifractal characteristics, it formed two clearly-separable clusters: this inexorably illuminates the true organization of the data and thus demonstrates the risk of faulty interpretation of MFDFA-and WL-based results.
To conclude, while these techniques offer very efficient and practical tools for quantifying the principal characteristics of multifractal patterns both in time and in space, extreme care needs to be taken when a suspicion arises that the pattern does not stem from true, uniform multifractality. The reason is that the methods are, by design, forced to assume a multifractal form for the singularity spectrum. This aspect points to a definite necessity to be addressed in future developments, namely, mitigating this strong prior assumption. For now, inspection of the distribution of the local Hölder exponents as given by wavelet-based techniques is crucial in identifying such instances and classifying them correctly.