Generalization of Higuchi's fractal dimension for multifractal analysis of time series with limited length

We introduce a generalization of Higuchi's estimator of the fractal dimension as a new way to characterize the multifractal spectrum of univariate time series. The resulting multifractal Higuchi dimension analysis (MF-HDA) method considers the order-$q$ moments of the partition function provided by the length of the time series graph at different levels of subsampling. The results obtained for different types of stochastic processes as well as real-world examples of word length series from fictional texts demonstrate that MF-HDA provides a reliable estimate of the multifractal spectrum already for moderate time series lengths. Practical advantages as well as disadvantages of the new approach as compared to other state-of-the-art methods of multifractal analysis are discussed, highlighting the particular potentials of MF-HDA to distinguish mono- from multi-fractal dynamics based on relatively short time series.


Introduction
Since Benoit Mandelbrot applied the concept of fractality for the first time to real-world time series, different approaches have been proposed to detect fractal and multifractal scaling properties in diverse signals from complex systems [1,2,3,4,5,6] (for recent reviews, see [7,8]). In many real-world cases, this scaling behavior reflects the fact that embedded variability components cover a wide range of time scales despite the absence of a dominant periodicity [9,10]. For monofractal time series, a single scaling exponent is sufficient to completely characterize the fractal properties of the dynamics displayed by the analyzed sequence. By contrast, for multifractal signals a variety of scaling exponents are necessary to fully describe the dynamics.
In the context of complex time series, the two standard methods to detect multifractality are based on either the construction of a partition function or the concept of extended self-similarity associated with the notion of structure functions (for details, see [3,8]). Both approaches provide direct procedures to obtain the spectrum of multifractal scaling exponents involved in the multifractal measure. While the partition function method is mainly rooted in the box-counting approach [3] and may therefore require relatively long time series for proper evaluation, structure functions rely on the nonlinear scaling behaviors of increments bridging increasingly large time differences [21,22]. A more sophisticated method for multifractal analysis, the wavelet transform modulus maxima (WTMM) method, has been introduced by Arneodo et al. [23,24] and applied to a great variety of problems from various disciplines. Here, the wavelet transform is used to obtain the sequence of local maximum values over a range of scales, which also accounts for the wide-spread nonstationarity of real-world time series. In a similar spirit, a conceptually related approach based on time scale decomposition by means of a more data-adaptive technique (the empirical mode decomposition) has been developed more recently [25,26].
Another state of the art approach of multifractal analysis has been introduced by Kantelhardt et al., who proposed a generalization of detrended fluctuation analysis [27] to obtain an estimate of the multifractal spectrum of nonstationary time series. This multifractal detrended fluctuation analysis (MF-DFA), together with various algorithmic variants thereof differing in the way of how local time series detrending is performed, provides a robust way to determine the multifractal scaling characteristics in terms of generalized Hurst exponents h(q), where q represents the order of the statistical moment under study. However, both WTMM and MF-DFA exhibit certain practical challenges in determining the corresponding multifractal spectra from relatively short time series [28,29], which may lead to spurious identification of multifractality.
As a possible alternative, this paper introduces a generalization of Higuchi's fractal dimension (HFD) estimator for univariate time series [30,31], which has been widely used as a technique for quantifying monofractal scaling characteristics in nonstationary time series from different areas of science, ranging from geophysics [30,32,33,34,35] to biomedicine [36,37,38,39]. Specifically, we introduce the concept of multifractal Higuchi dimension analysis (MF-HDA), which is applicable to comparatively short sequences and provides very stable estimates of scaling exponents for the construction of the multifractal spectrum.
The remainder of this paper is organized as follows. In Section 2, we briefly review some essential details on the HFD method and introduce a thorough modification of the original approach to obtain the corresponding scaling exponent. Section 3 describes the multifractal generalization of the monofractal dimension estimator. Its application to simulated time series from two illustrative stochastic model systems and a subsequent comparison of the obtained performance with that of MF-DFA are presented in Sections 4 and 5. In Section 6, we then apply our method to some real-world examples of world length sequences in English fictional texts. Finally, a discussion and some conclusions and general remarks are provided in Sections 7 and 8, respectively.
2 The fractal dimension of irregular time series 2.1 Higuchi's estimator of the fractal dimension Consider a time-series X(t) with t = 1, 2, . . . , N . First, we construct a set of sub-sampled time series at different scales (determined by k) with m = 1, 2, . . . , k, where ⌊•⌋ denotes the (lower) integer part of a real number. Here, k and m are integers which describe a time lag and the initial time index, respectively. Next, the length of the curve associated with a given sequence X m k is defined as [30] where ∆X m k (j) = |X(m + jk) − X(m + (j − 1) · k)| are the absolute increments of the subsampled series, j max = N −m k is the total number of increments at scale k, and (N −1)/(j max ·k) represents a normalization factor.
We notice that the quantity ∆X m k represents the differences corresponding to each k-lag value, which when aggregated and normalized by j max gives the average of the increments for that scale k. For the purpose of the present study, consider for example X(t) to represent the position of a one-dimensional random walk at time t. Then, ∆X m k is the absolute displacement within a time step k with m = 1, 2, . . . , k, where the best possible temporal resolution corresponds to k = 1, the resolution decreases by one half for k = 2, and so on. Thus, the mean absolute displacement of X(t) starting at m at a resolution k is given by and the total mean displacement is calculated as ℓ m (k) = represents the typical number of displacements available at a resolution k. In order to make these total mean displacements obtained for different resolutions k mutually comparable, the total length of the curve describing the time series graph coarse-grained at scale k and starting at m is defined as L m (k) = ℓ m (k)/k, and the total mean curve length at scale k is given as L(k) = L m (k) [30,31].
In case of a self-affine (i.e., fractal) time series, it has been demonstrated that where D A ∈ (1, 2) is an estimate of the fractal dimension of the original time series graph and the subscript A stands for averaging [30,31].

Modified HFD estimator
For the purpose of generalizing the HFD method to a multifractal formalism, it is necessary to first introduce an alternative formulation of an estimator of the fractal dimension that is based on the same rationale as Higuchi's original method. Here, our modified approach replaces the averaging procedure for calculating the mean absolute displacement by an expectation value in terms of probabilities of increment sizes ∆X m k (Eq. 3). We first recall that the number of increments ∆X m k at scale k is ⌊(N − m)/k⌋, and the total number is k m=1 ⌊(N − m)/k⌋, which equals N − k [40]. For numerically evaluating the corresponding expectation values, we follow here a simple histogram based approach (other more sophisticated density estimators may be used as well, yet we focus here on this simple strategy and outline corresponding follow-up studies on this aspect as parts of future work). Let us consider P n (•) representing the discrete probabilities of increments following into disjoint intervals ∆X n (k), where n represents the n−th interval of an N b equiprobable partition of the entire range of increment values in increasing order, that is, each interval (bin) contains approximately a fraction of 1/N b of the N − k increments, i.e., P n (∆X) ≈ 1/N b . Note that in practice, the total number of increments ∆X may not be an integer multiple of N b , which results in minor deviations among class frequencies that become gradually less important as N b /N → ∞. For the purpose of the present work, we employ numerical estimates of the respective n/N b quantiles (n = 1, . . . , N b − 1) to define mutually disjoint intervals with approximately the same population.
In order to obtain a stable estimate of the expectation value of the increments, we start with calculating the "local-mean" increment values, i.e., the mean values for each bin, as ∆X n (k) = 1 fn fn s=1 ∆X s (k), where f n represents the actual number of increments within the n-th bin and n f n = N −k. Then, the expectation value of the increments at scale k can be estimated as In analogy with Higuchi's original method, the total length of the curve at resolution k can then be written as: Again, we may expect that if the original time series is fractal, then the total length follows a power-law behavior L(k) ∼ k −DE , where D E is our modified Higuchi type estimator of the fractal dimension and the subscript E stands for expectation. In order to verify the consistency of our modified approach with the classical Higuchi estimator, we applied both our methodology and the original HFD method to some elementary stochastic time series with known properties. Specifically, we first consider realizations of three representative monofractal processes: uncorrelated noise (white noise), intermittent long-term correlated noise (1/f noise) and classical Brownian motion [41]. The respective results are summarized in Fig. 1. In general, we observe an excellent agreement between our proposed modification and Higuchi's original method for all three monofractal series.
To further evaluate the agreement between both methodologies, we systematically generated monofractal time series of fractional Gaussian noise with gradually varying spectral scaling exponent β [43] and then applied both approaches for estimating the corresponding fractal dimensions. Figure 2 shows the numerical results for different processes with β ∈ (0, 4). Note that we have to expect appropriate estimates to exist only Power spectral exponent β  Fig. 2 Estimates of the fractal dimension (D E and D A ) for monofractal series of fractional Gaussian noises with prescribed spectral exponent β. Both methodologies lead to very similar average values for the simulated stochastic processes. The dashed line indicates the theoretical relationship between the fractal dimension and the spectral exponent β for N → ∞, which is given by D = (5 − β)/2 and applies for 1 < β < 3 [31,42]. Error bars represent the numerical standard deviations obtained from 100 independent realizations of length N = 10, 000.
for β ∈ (1, 3) and the finite-sample estimates to become gradually biased as any of the boundaries of the latter interval are approached [44]. This expectation is met by our numerical results, which again demonstrate that both methodologies lead to very similar fractal dimension estimates.

Generalization of the modified Higuchi estimator
Based on our modification of Higuchi's fractal dimension estimator as presented in the previous section, our principal interest is the generalization of this monofractal analysis method to a multifractal framework. For this purpose, we consider the q-th order moment of the modified Higuchi's curve length as Here, q can take any real number (for a discussion of the special case of q = 0, see Appendix A). The behavior of this generalized Higuchi length L(k, q) for different orders of the moment q is indicative of a dependence of the characteristic (fractal) scaling exponents on q.
More precisely, if the time series under study displays power law correlations (i.e., long-term memory resulting in fractal scaling characteristics), we can expect to find the power law scaling behavior where d q is a generalized fractal dimension associated to the q-th moment. Notice that for q = 1, our modification of the standard Higuchi fractal dimension method is recovered. For Gaussian time series exhibiting an extended self-similarity (multifractal) property, d q is related to the generalized Hurst exponent h q (also known as the Hölder exponent) by the relationship h q = 2 − d q [10]. In our case, the singularities are identified by the increments ∆X s (k) q . For q > 0, the contributions of fluctuations corresponding to large increments are highlighted, while the opposite situation is obtained for q < 0, i.e., fluctuations associated to small increments are amplified.
We emphasize that the increments ∆X s (k) can be small enough to affect the numerical estimates of the moments of order q, especially for negative q, where reliable estimation becomes more and more difficult. To address this challenge, a suitable numerical regularization procedure may become necessary to allow a feasible calculation of the generalized lengths in Eq. (7). If the local behavior of the generalized length L(k, q) is found to be unstable (see Appendix B for details), we suggest here to replace the local mean ∆X n (k) = 1 fn fn s=1 ∆X s (k) in Eq. (7) by a regularized version trimmed at small increment values, ∆X n ′ (k) >pr , where n ′ represents the n ′ −th interval of a new equiprobable partition in which we have removed all increment values below their empirical r−th percentile p r . In this case, the number of increments at scale k is now given by denoting the number of elements below the r−th percentile. In our numerical experiments described in the remainder of this work, we typically set a percentile value within the range r ∈ [1, 15] (see Appendix B), while for the number of bins in the partition, we set N b = 16 (where for very long time series, finer partitions may also become numerically feasible).

Partition function and multifractal analysis
Considering the scaling behavior suggested in Eqs. (7) and (8), it follows that Using the fact that P n (∆X(k)) = f n /(N − k) and that for self-affine time series, h q = 2 − d q , Eq. (10) can be rewritten as which represents the well-known generalized structure function [3]. We notice that asymptotically (i.e., for large N ) 1/(N − k) ∼ k, and therefore, In the standard multifractal formalism, the partition function Z(q, k) is used to define the mass exponent τ (q), and the following scaling relation is observed: where we have assumed that the left-hand of Eq. (12) corresponds to the partition function Z(q, k). By comparing Eqs. (12) and (13), we obtain the relationship between the scaling exponents as Thus, we have related the exponent h q to the fractal scaling exponent τ (q), and the requirement of the fractal dimension value (when q = 0) for the geometric support of the multifractal measure −τ (0) = 1 = D f is fully satisfied [27]. Moreover, the canonical way to characterize a multifractal signal is the singularity spectrum f (α), which is obtained by means of τ (q) and its first derivative α as where f (α) represents the fractal dimension of the subset of the signal characterized by the local exponent α [3].
4 Numerical examples: stochastic time series

Monofractal fractional Gaussian noises
In order to demonstrate the performance of our multifractal Higuchi dimension analysis (MF-HDA) method when applied to simulated monofractal times series, we generated realizations of fractional Gaussian noises with prescribed Hurst exponents of H = 0.3, H = 0.5 and H = 0.75, respectively, by means of the Fourier transform method [43]. Figure 3 shows the results of the multifractal analysis for all three processes. As it can be seen in Fig. 3a, the generalized curve lengths for H = 0.5 exhibit a power law scaling behavior for which well defined exponents can be estimated (Eq. 8).
In Fig. 3b, the exponents h q = 2 − d q are plotted against the moment q. As expected, only a weak dependence of the numerical estimates on q is observed, indicating that the time series are, in fact, monofractal.
Moreover, Fig. 3c shows the behaviors of τ (q) versus q, where clearly linear dependencies on q are observed for all three processes. These results lead to very narrow multifractal spectra estimates (Fig. 3d), which are consistent with the expected monofractal properties of the original time series within the expected limitations of numerical finite sample estimates. The results show only a weak dependence on q, confirming the monofractality of the series. c) Behavior of τ (q) versus q, which indicates a linear dependence. d) Multifractal spectra f (α). Here, the singularity spectra f (α) were obtained from the Legendre transform of τ (q).

Binomial multifractal cascade
The binomial multifractal cascade (BMC) is one of the most representative models for which the multifractal features can be derived analytically, allowing a direct comparison with the numerical results of our method (for details on the BMC model, see [3]). In the following, we explain only briefly the essential details of the model [3,45,46]. In order to generate a series z j , with j = 1, . . . , 2 mmax , we consider the following expression: where a ∈ (0.5, 1) is a parameter, n(j − 1) is the number of digits equal to 1 in the binary representation of the index j − 1, and m max represents the number of iterations. One of the advantages of studying binomial cascades is that the key expressions of the multifractal properties can be obtained analytically in a straightforward manner as [27] τ and from which the multifractal spectrum can be obtained [3].
The numerical results of our MF-HDA method for BMC series with parameter a = 0.75 and N = 2 mmax = 2, 048 are shown in Fig. 4. The behavior of the generalized Higuchi curve length L(k, q) versus k for −5 ≤ q ≤ 5 is depicted in Fig. 4a, where evident changes in the slope (i.e., the generalized fractal dimensions d q ) are present as q varies. The corresponding h q exponents are shown in Fig. 4b. The fact that h q varies with q evidences the multifractal property of the BMC series. The mass function τ (q) (Fig. 4c) presents a clearly nonlinear behavior with q, as expected for multifractal series [3]. Finally, the resulting multifractal spectrum is shown in Fig. 4d, which was obtained from h q and τ (q) through the Legendre transform. We observe a very good agreement with the analytical predictions of the BMC model, except for an absence of the largest α values which may lead to an underestimation of the multifractal width most likely occurring due to the short length of the considered sequence.

Comparison of MF-HDA with MF-DFA
In the previous section, we have already presented results on the multifractal properties for selected time series with prescribed lengths and moments within a predefined interval of q. In the following, we will be interested in comparing our new MF-HDA method with the established MF-DFA approach, which presents one of the most widely recognized and applied methodologies to assess the presence of multifractality in nonstationary time series [8,27,28,47,48]. For a detailed study on the applicability of MF-DFA and the WTMM method as another benchmark approach in multifractal analysis, the reader is referred to [29].

Multifractal detrended fluctuation analysis
Given a time series X(t), we compute the associated profile [27] with t = 1, . . . , N , where N is the number of values in the time series andX denotes the corresponding mean. Next, the integrated series Y (t) is divided into N s boxes of equal size ⌊N/s⌋. The local trend Y r (t) is calculated separately for data from each of the segments r = 1, . . . , N s by least square fitting of a polynomial of prescribed degree (i.e., the order of detrending) and removed from the profile. Then, the mean square fluctuation of the detrended profile in each segment is given as The resulting (squared) detrended fluctuation function F 2 (s) is defined as the average mean square fluctuation taken over all segments r. For monofractal series, its square-root behaves like F (s) ∼ s H with H being an estimate of the process' characteristic Hurst exponent.
To generalize this detrended fluctuation analysis to a multifractal analysis framework, we consider the orderq moments of the mean square fluctuations at scale s by setting [F 2 (s, r)] q/2 1/q and evaluating this property for several values of q (retaining the standard detrended fluctuation analysis for q = 2) [27]. Finally, the scaling behavior is described by F (q, s) ∼ s hq , where h q is the generalized Hurst exponent, which is related to the mass exponent τ (q) by means of Eq. (14).
In the following, we will compare the performance of the resulting MF-DFA method with our new MF-HDA approach by considering (i) the effect of different time series lengths and (ii) different ranges of q values.

Fractional Gaussian noise
As in the previous section, we first consider monofractal time series with prescribed Hurst exponents of H = 0.3, H = 0.5 and H = 0.75 with different lengths N = 500, 1, 000, 5, 000, 15, 000, and 65, 000. We apply both MF-HDA and MF-DFA to the generated series. Figure 5 shows the corresponding estimates of the associated multifractal spectra. Especially for short time series lengths, MF-DFA displays wider spectra of values as compared to those obtained by means of MF-HDA. This behavior is especially noticeable for N ≤ 5, 000, where the MF-DFA spectrum is markedly extended towards larger α. As the time series length increases, both methods reveal gradually more similar multifractal features, except for the fact that the spectra estimated by means of MF-HDA are less symmetric, due to an absence of the largest α values (corresponding to very negative moment orders q), which are most affected by the employed regularization procedure removing the smallest increments.
Concerning the effect of the selected range of q values, we also evaluated the changes in f (α) in terms of different intervals of q for the same monofractal series, here with a fixed length of N = 130, 000. The corresponding results are shown in Fig. 6. We observe that, for all three monofractal sequences, both methods lead to similar spectra for the interval |q| ≤ 3, while for intervals |q| ≤ 5 and |q| ≤ 10, MF-HDA provides narrower spectra compared to those obtained by means of MF-DFA for the reasons already mentioned above related to the regularization affecting the smallest fluctuations.

Binomial multifractal cascade
In addition to simple monofractal processes, we also compared the MF-HDA and MF-DFA methods for the BMC model for several time series lengths. Figure 7 shows the spectra estimated with both methods along with the analytical results. We find that the obtained estimates are generally in good agreement with the theoretical spectra. Some minor deviations at small α values (i.e., for large positive q values) are observed with both methods (with MF-DFA generally showing slightly larger deviations from the theoretical values), while MF-HDA does not manage to capture the rightmost part of the spectra (large α associated with large negative q values) especially for very short time series lengths.
We also tested the dependence of the obtained estimates on the range of q values considered. Figure 8 shows that, as expected, the effect of widening the q interval is to increase the width of the spectrum, with fairly similar behaviors of both methods.
For a more systematic comparison between both approaches, we finally computed the absolute differences between the theoretical values of the singularity strength exponents α T (obtained from Eqs. (17) and (18)) and the corresponding numerical estimates α E obtained using either MF-HDA or MF-DFA. Figure 9 shows these absolute differences (i.e., the estimation bias) as a function of the moments q for several representative time series lengths. For very short series and for q ≥ −1, the MF-HDA estimates are closer to the theoretical values, while differences between the theoretical and estimated values are smaller for q ≤ −2 when using MF-DFA (lacking the shortcomings due to the regularization procedure within our numerical implementation of MF-HDA). For long series, both methods lead to similar results with slightly smaller deviations when using the new MF-HDA method. From these results, we conclude that in the BMC model MF-HDA exhibits certain advantages (lower bias) over MF-DFA especially for short time series and when only a small interval of q values is considered.

Real-world application: Word length sequences
As an application to empirical data sets with intrinsic (time) order, we consider the sequences of subsequent word lengths (i.e., the numbers of letters forming each word [49]) obtained from a selection of fictional texts in English language, which include mainly novels from different authors (see Tab correlations and multifractality in such word length sequences, which have been attributed to the fact that written language is the conformation of grammatical properties and semantic connotations with the idea of expressing information [50]. In order to allow for a feasible application of the MF-HDA method to the different word length series, the original length sequence has first been integrated to obtain a profile compatible with a fractional Brownian motion regime. Then, we applied both MF-HDA and MF-DFA to word length sequences with different segment lengths (N = 2, 500, 5, 000, 10, 000). In order to improve the reliability of the obtained statistics, several non-overlapping segments of the same length N have been considered for each book. Then, we have used the following procedure to perform MF-HDA: (i) For each q and the m-th segment, the generalized curve length L m (k, q) is calculated.
is evaluated, where M denotes the number of segments for each book. (iii) The value d q,av is estimated from the fit of the scaling behavior L av (k, q) ∼ k dq,av . Figure 10 shows the results obtained with both MF-HDA and MF-DFA. We observe that unlike for the simple stochastic model cases studied in the previous section, the spectra f (α) obtained with MF-HDA are broader than those provided by MF-DFA, especially for short sequences. Notoriously, the expected convex shape of multifractal spectra is not well expressed in the MF-DFA results. Moreover, as the length of the indicates that the accumulated word length series exhibit correlated behavior, while for shuffled and integrated increments very narrow spectra centred at h = 0.5 are observed. (Right) Multifractal spectra estimated using MF-DFA. A similar procedure as for MF-HDA was used to estimate the scaling exponents hq,av [48]. The fitting region was 50 ≤ s ≤ T /10. considered segments is increased, the width of the spectra decreases but still remains relatively broad for MF-HDA as compared to MF-DFA. This result is remarkable since, as we have noticed in our controlled numerical experiments above, we may have expected MF-DFA tending to rather overestimate the multifractal spectral width for short time series. As a consistency check, we finally generated shuffled versions of the word length series by applying a random permutation of the individual word lengths, thereby destroying eventual correlations in the sequence. The results of applying both multifractal analysis methods to the accordingly randomized sequences are also presented in Fig. 10. As expected, we obtain narrow spectra centered at α = 0.5 for all books and different segment lengths, which do not differ much between the two methods.

Discussion
The identification of multifractal properties in time series is a challenge especially when dealing with short sequences, as are present in many real-world cases. The proposed MF-HDA method offers an alternative to assess the multifractal properties of relatively short series Table 1 Books used for the multifractal analysis of word length sequences. The table lists the title and author, code name and the total number of words, along with estimates of the global Hurst exponent h q=1 , the spectral mode α * and the spectral width ∆α as obtained using MF-HDA averaged over different segments with length N = 2, 500 (see Fig. 10). Although we have not followed a particular selection strategy for the studied titles beyond just considering well-known books, all estimated values for the global Hurst exponent (h q=1 ) are clearly larger than 0.5, confirming the presence of long range correlations in the word length series. These values are also in good agreement with previous monofractal analyses reported in [49,51,52]. The α * values at which the multifractal spectra take their maxima are also consistently larger than 0.5 for all books. Moreover, most of the books are characterized by relatively broad multifractal spectra with ∆α ≥ 0.1, except for AAW, ANC and TTM (highligthed in bold face), for which the spectra exhibit a discontinuous behavior. Our results are consistent with previous analyses based on traditional multifractal procedures [53,54] and provide a better characterization of the complexity displayed in written texts. Interestingly, the book with the largest multifractal spectral width (∆α = 0.36) corresponds to Joyce's Ulysses (ULY, underlined), which has been described as a text with particularly great diversity in language [55]. by constructing generalized curve lengths, the scaling of which encodes the multifractal characteristics of the signal. The implementation of our methodology is relatively simple, as it is based on the expected q-moment values of absolute increment distributions. We have described several cases where the MF-HDA method identifies a multifractal scaling behavior of the time series with very similar results as the celebrated MF-DFA approach, which provides a well-established methodology for the analysis of non-stationary signals. For very short series, the results for the MF-HDA indicate a certain gain in accuracy as compared to those obtained using the MF-DFA method, except for very negative moment orders q.
As with the majority of studies focusing on the detection of multifractality in non-stationary time series, the proposed methodology is subject to certain limitations: (i) dominant scaling exponents of the time series can be located outside of the region where the MF-HDA allows us to detect multifractal scaling in a reliable and robust manner; (ii) the limited sampling rate and numerical precision of the time series makes the expected moments of the generalized curve length quite difficult to estimate (especially for negative moments); (iii) the removal of the p r percentile of the absolute increments has a regularizing effect but can potentially lead to another bias in the estimated scaling exponents, especially for very negative q values. Further refining the presented numerical procedures to address the aforementioned challenges should be addressed in follow-up studies.

Conclusions
We have presented a generalization of the Higuchi fractal dimension analysis to a multifractal framework. The main advantage of the new MF-HDA method in comparison with other established multifractal analysi techniques lies in the property of obtaining stable estimates of multifractal scaling characteristics even when the analyzed sequences are relatively short. However, certain precautions must be considered when analyzing realworld time series, for which the origin of multifractality can be associated with different factors such as power law correlations or a heavy-tailed probability density function of the data [9,27]. For a more reliable estimation of multifractal characteristics, an integrated approach involving complementary methods such as MF-HDA, MF-DFA or WTMM is recommended. We envision to further explore the corresponding potentials and more systematically apply our methodology to a range of real-world and simulated time series in our future work.

Code availablity
Numerical implementations and examples for the application of MF-HDA can be found at https://github.com/carrizales90/MF-HDA.
A The case q = 0 In multifractal analysis, it is common that the scaling exponent dq is not well defined when q → 0. In our case, this value cannot be directly determined by means of the generalized curve lengths (Eq. (7)) due to the presence of a divergence in the exponent. More formally, we have Using some algebraic operations applied to the latter equation, which are omitted here for brevity, and applying L'Hôpital's rule, we find that a logarithmic transformation is required in order to determine the scaling exponent d 0 as where E[ ln{∆X(k)} ] = N b n=1 ln{∆X(k)} Pn(∆X(k)).

B Regularization effect of removing the lower percentile of absolute increments
As discussed in Section 3, numerical instabilities can appear in the evaluation of the moments of the generalized curve length L(k, q) (Eq. 7), especially for q < −1. In this case, we have suggested that the local mean could be replaced by a one-sided trimmed version ∆X n ′ (k) >pr in Eq. (7), where n ′ represents the n ′ -th interval of a new equiprobable partition in which we have removed the r-th percentile pr of the empirical distribution of the absolute increments. We note that numerical experiments with both, one-sided (asymmetrically) and two-sided (symmetrically) trimmed means revealed no qualitative differences in the resulting estimates (not shown), while removing the uppermost percentiles (i.e., very large increments) appears unnecessary since those values have no negative effects on the stability of the numerical estimates of the generalized curve lengths.
To address the problem of selecting a specific percentile to be removed, we focus here on just one statistical property, the confidence interval (CI) of the estimated slope (i.e., the scaling exponent dq in Eq. (8)) in the linear regression of the double-logarithmic generalized curve length versus scale relationship, at a certain confidence level (here, α = 0.05), and for the most negative value of q. For two-sided confidence intervals, the CI width (CIW) (measured in units of the associated standard error) is given byÎq,p r ≡ Iq,p r /S dq , with S dq being the standard error of the estimated slope dq [56]. Figure 11 shows the behavior of the rescaled CIW (in units ofÎ pr =0 ) as a function of the removed percentile pr, for some of the simulated stochastic processes and real-world data sets discussed in Sections 4 and 6, respectively, for q = −5. The results show that, as the removed percentile is increased, the CIW decreases in such a way that, for fractional Gaussian noises with H = 0.3, H = 0.5 and H = 0.75, the rescaled CIW has decayed by more than one half of its initial value when pr = 1, while for the world length data (exemplified here by the ULY book) the observed decay is slower. For practical purposes, we suggest that a criterion for selecting the value of the percentile to be removed should consider empirically a value of the percentile for which the rescaled CIW has stabilized, that is, even if higher percentiles are removed, there are no substantial further changes. We observe that pr ≈ 5 and pr ≈ 12 would be desirable in the cases of the simulated fractional Gaussian noises and word length data, respectively.
While the suggested strategy presents just a first attempt to improving the practical estimation of the generalized Higuchi fractal dimensions, we emphasize that there may be cases in which the rescaled CIW may behave in a more unstable way with increasing percentile pr. In such cases, additional numerical tests with larger pr values may become necessary to determine a reliable value leading to sufficiently stable estimates.