MFDFA: Efficient Multifractal Detrended Fluctuation Analysis in Python

Multifractal detrended fluctuation analysis (MFDFA) has become a central method to characterise the variability and uncertainty in empiric time series. Extracting the fluctuations on different temporal scales allows quantifying the strength and correlations in the underlying stochastic properties, their scaling behaviour, as well as the level of fractality. Several extensions to the fundamental method have been developed over the years, vastly enhancing the applicability of MFDFA, e.g. empirical mode decomposition for the study of long-range correlations and persistence. In this article we introduce an efficient, easy-to-use python library for MFDFA, incorporating the most common extensions and harnessing the most of multi-threaded processing for very fast calculations.

Multifractal detrended fluctuation analysis (MFDFA) has become a central method to characterise the variability and uncertainty in empiric time series. Extracting the fluctuations on different temporal scales allows quantifying the strength and correlations in the underlying stochastic properties, their scaling behaviour, as well as the level of fractality. Several extensions to the fundamental method have been developed over the years, vastly enhancing the applicability of MFDFA, e.g. empirical mode decomposition for the study of long-range correlations and persistence. In this article we introduce an efficient, easy-to-use python library for MFDFA, incorporating the most common extensions and harnessing the most of multi-threaded processing for very fast calculations.
MFDFA is a numerical algorithm designed to determine the self-similarity of a stochastic process. Putting it simply, the algorithm examines the relation between the diffusion of the process and its propagation in time or space. Auto-regressive and stochastic processes with different power-law scaling will diffuse with different rates. Fluctuation Analysis (FA) provides a method to uncover these correlations, but fails in the presence of trends in the data, which, for example, are particularly present in weather and climate data. Detrending the data via polynomial fittings (DFA) allows one to uncover solely the relation between the inherent fluctuations and the time scaling of a process, thus circumventing the impact of non-stationarity in the data. Likewise, other methods-as empirical mode decomposition or moving average windows-are viable options to detrend the data. Another problem is that a process might be driven by more than one time scale, i.e., have more than one internal period, which can be removed either with local polynomial fittings or EMD. Moreover, a stochastic process might be of a monofractal or multifractal nature. By studying a continuum of power variations of DFA one extends into MFDFA, which permits the study of the fractality of the data by comparing power variations, i.e., a multifractal spectrum.
In this software we sought to design a computationally efficient code focused on computational speed and usability. There are currently no flexible and available implementations of MFDFA in python. Available are some MATLAB [55] as well as R packages [56,57]. There is a particularly thorough introductory guide to MFDFA in MATLAB with a source-code by Espen A. F. Ihlen [55], which is easy to implement but numerically inefficient. With this implementation efficiency was sought. This was achieved by making the most out of python, reshaping the code to allow for multi-threading, especially relying on numpy's polynomial, which scales easily with modern computers having more processor cores [58]. Moreover, this library contains the most commonly applied methods alongside with DFA and MFDFA: the added feature of empirical mode decomposition is implemented to sub-stitute the polynomial fittings; A moving window is included, especially valuable for shorter time series; The extended DFA (eDFA) method is also included, adding a second metric of fractal scaling, especially valuable for multifractal or aperiodic time series.
In the following sections we will introduce MFDFA alongside some of the aforementioned methods incorported into the MFDFA library. We will present two classical applications, one with a monofractal process and one with a multifractal noise, and show how to use MFDFA to extract their characteristics from a single one-dimensional time series. Python code is presented to explicate the use of the MFDFA library. Subsequently we study two real-world time series: the sunspot time series from 1818 to 2020 which accounts for the daily recorded sunspots and the quarter-hourly electricity trading market, which accounts for a small volume of electricity sell and purchase at 15 minute windows in Continental Europe. Lastly we address a few details of the library and contribute a few closing remarks.

II. THEORETICAL BASIS
In the following we briefly summarise the theoretical basis of Multifractal Detrended Fluctuation Analysis. Later we detail the different included extensions and which modifications these add to the original MFDFA algorithm.

A. Multifractal Detrended Fluctuation Analysis
Multifractal Detrended Fluctuation Analysis studies the variances of the fluctuations of a given process by considering increasing segments of a time series.
i) Take a time series X(t) (in time or space t) with N data points, discretised as X i , i = 1, 2, . . . , N . Find the "detrended" profile of the process by defining i.e., the cumulative sum of X i subtracting the mean µ X of the data. ii) Section the data into smaller non-overlapping segments of length s, obtaining therefore N s = int(N/s) segments. Given the total length of the data is not always a multiple of the segment's length s, discard the last points of the data.
iii) Consider the same data, apply the same procedure, but discard now instead the first points of the data. One has now 2N s segments of the time series. iv) To each of this segments fit a polynomial y v of order m and calculate the variance of the difference of the data to the polynomial fit for v = 1, 2, . . . , N s , where y (v−1)s+i is the polynomial fitting for the segment Y (v−1)s+i of length s, fitted via leastsquares. The order of the polynomial y v can be freely chosen, giving rise to the denotes (MF)DFA1, (MF)DFA2, . . . , (MF)DFAm, dependent on the chosen degree m of the polynomial. v) Notice now F (v, s) is a function of each variance of each v-segment of data and of the different s-length segments chosen. Define the q-th order fluctuation function by averaging over the N s variances of the segments of size s The fluctuation function F q (s) depends on two parameters: the segment size s and the q-th power. The fluctuation function F q (s) is the function we will focus on which the MFDFA algorithm developed extracts from the data. Two closely related algorithms are discussed and introduced here, DFA [1] and MFDFA [3]. DFA is a particular case of MFDFA for the choice of q = 2. What is presented above is the MFDFA algorithm as according to Kantelhardt et al. [3], for which a particular choice of q = 2 leads to the fluctuation function F 2 (s). The DFA fluctuation function F 2 (s) can unveil solely the monofractal spectrum of a time series. If the examined time series X i is monofractal, DFA is sufficient to describe and uncover the scaling relations in the data. If not, one must rely on MFDFA and the study of the spectrum unveiled by varying the q-th power.
We will later detail two changes: i) The first involving empirical mode decomposition (EMD) for detrending, where the local polynomial fittings are replaced and the trends of the data are subtracted by removing select Intrinsic Mode Functions (IMFs) obtained via empirical mode decomposition. ii) The second change involves substituting the non-overlapping segments with overlapping ones.
The inherent scaling properties of the data, if the data displays power-law correlations, can now be studied in a log-log plot of F q (s) versus s, where the scaling of the data obeys a power-law with exponent h(q) as where h(q) is the generalised Hurst exponent or selfsimilarity exponent, which will dependent on q if the data is multifractal, and relates directly to the Hurst index [59]. The generalised Hurst exponent h(q) is obtained by finding the slope of F q (s) curve in the log-log plots.
If the data is monofractal, the generalised Hurst expo- nent h(q) = H is independent of q and the generalised Hurst exponent is simply the Hurst index H. On the other hand, if the data is multifractal, the dependence on q can be understood by studying the multifractal scaling exponent τ (q), given by which depends on the generalised Hurst exponent h(q). Similarly, one can construct the singularity spectrum D(α) as the Legendre transform [60][61][62][63]. If τ (q) is sufficiently smooth, the singularity strength α is given by from which the singularity spectrum D(α) can be constructed as The singularity spectrum D(α) describes the dimension of the subset of the time series which is characterised by the singularity strength α [64]. The breadth of singularity strength α indicates the strength of the multifractality of the time series, centred around the most prominent scale of the time series, i.e., h. The singularity spectrum D(α) takes the shape of an inverted parabola with a maximum at D(α = 0) = D 0 , known as the boxcounting or Minkowski-Bouligand dimension, or sometimes simply fractal dimension [60]. D(α = 1) = D 1 is known as the information dimension and D(α = 2) = D 2 the correlation dimension [65]. For a clearer discussion of these properties, see Refs. [3,66]. An extensive and very illustrative representation of this can be found in Ref. [55]. For a careful analysis of the meaning and interpretation of the generalised Hurst coefficients extracted from (MF)DFA, see Ref. [67], where a description and clarification is given on what are persistent and antipersistent motions, stationary and non-stationarity time series, among other relevant details.

Empirical mode decomposition
Empirical mode decomposition (EMD) is a method with a variety of applications in time series analysis [68]. It seeks to extract the modes of oscillation of a time series strictly from the data. One can harness the ability of the EMD, i.e., the Hilbert-Huang decomposition of a time series, to obtain the trend or trends of the time series and utilise those to transform non-stationary into stationary data. The central concept, developed by Qian, Gu, and Zhou [6], is to substitute the detrending method employed in the traditional MFDFA, i.e., polynomial fittings, by removing instead particular Intrinsic Mode functions extracted via EMD. A sketch of the method can be seen in Fig. 1.
EMD can be summarised in a few steps: a set of intrinsic mode functions (IMFs) are extracted from the time series, obeying: 1) the number of extrema and the number of zero crossings must maximally differ by one. 2) for any point, the mean value of the envelope defined by the local extrema is zero. Numerical methods-as cubic splinesare used to find the curve that best fits "between" the local extrema of the time series. The method is applied iteratively: i) Obtain an IMF by finding the "best" curve between the local extrema of the time series; ii) Subtract this IMF to the time series; iii) Repeat. Apply the process recursively to the time series until the final IMF contains solely a residual trend of the data.

Moving windows
The overlapping moving windows included in this library is not aimed at detrending, but instead for the analyses of rather short time series or very large scales in longer time series [8]. In the literature several applications of moving average windows have been proposed as methods to remove trends and ensure stationarity, by simply removing a windowed average to the time se-ries [8,9,42]. This is not what we do here. Here, we substitute the non-overlapping segmentation, as explain after Eq. (1), by a moving window, replacing the two separate segmentations by a moving window of each segment size s over the time series, as proposed by Zhou and Leung [8]. This is particularly relevant when examining short time series, where quickly the choice of larger lags s separates the data into a small number of segments, resulting in a poor statistics for the scaling at larger lags. A sketch of the methods can be seen in Fig. 1.

Extended Detrended Fluctuation Analysis
A new metric of similar nature as the fluctuation function F q (s), given in Eq. (3), has been proposed in Ref. [10] This measure supersedes the q-order powers and takes in solely the case of DFA where q = 2. Instead of finding the average of the variances over each choice of segments of size s, it considers the difference between the extrema of the fluctuation function at each segment s. Take F (v, s) as given in Eq. (2) and extract the maximum and minimum of the variances over all windows v for a certain window size s This new metric ∆F (s) is denoted Extended Fluctuation Analysis. In general, ∆F (s) can scale as a power law with a different exponent This metric takes into account aperiodicities in the data which, in some sense, would be accounted for as a multifractal behaviour. It can unravel a second scaling phenomenon due to local changes of a time series' period.

III. EXAMPLES
To exemplify the usage of MFDFA, we first take two common examples of stochastic processes, a fractional Ornstein-Uhlenbeck process and general process that has a symmetric Lévy α-stable distribution, with single parameter α. We will show how to extract the fluctuation function F q (s) and how to interpret the plots conventionally extracted to perform the analysis. Subsequently we test the algorithm with real-world data on sunspot time series, following Ref. [26], and later apply the algorithm to electricity price time series from the European Power Exchange.
A. Numerically generated data

Fractional Ornstein-Uhlenbeck process
To study the scaling effects in continuous stochastic processes, three exemplary fractional Ornstein-Uhlenbeck processes are taken, defined as [69] with a fractional Brownian motion B H t with the covariance function Eq. (10) fixed mean reverting strength θ = 1.0 and volatility σ = 0.5, with three distinct Hurst indices of H = 0.3, 0.5, and 0.7. Note here that the classic uncorrelated Brownian motion is given by H = 0.5. A fractional Brownian motion has a self-similarity exponent given by the Hurst index H, thus the three choices of fractional Ornstein-Uhlenbeck should result in a scaling of h(q) = H + 1. The +1 is due to the integration, which smooths the fluctuations and thus increases the regularity of the process. We will now numerically integrate these processes and utilise the MFDFA library to identify the Hurst coefficients and the presence of a monofractal vs multifractal spectrum in the time series. Let us exemplify how to numerically generate data and utilise the MFDFA library Load the MFDFA library alongside with the fractional Brownian noise generator fgn included in your python console or editor. To numerically integrate an Ornstein-Uhlenbeck process, given by Eq. (10), we utilise an Euler-Maruyama scheme with a stepsize ∆t = 0.001 for a total time of t = 10 4 (thus we have 10 7 data points). Here exemplified is the fractional Ornstein-Uhlenbeck process with H = 0.3.  To retrieve the MFDFA spectrum of the generated time series, define the set of q power variations and the lags s to examine, and call the MFDFA function. When not declaring the values of the q powers, q = 2 is assumed, thus resulting in the conventional DFA. Likewise, not declaring the order of the polynomial fitting, a first-order polynomial is assumed, i.e., order = 1.
In Fig. 2 the MFDFA of the three processes can be seen. In panel a) the fluctuation function F 2 (s), with q = 2, is shown for a polynomial detrending of first order. This is the conventional DFA. The slopes of each curve in the log-log plot reveal the Hurst indices of each process, i.e., the fractional Ornstein-Uhlenbeck with Hurst H = 0.3 scales with a slope of 1.3 = 0.3+1, the other two with H = 0.5 and H = 0.7 have a slope of 1.5 and 1.7, respectively. The inset in panel a) shows the fluctuation function F q (s), with q = −10, −2, 2, and 10, for the fractional Ornstein-Uhlenbeck process with H = 0.3. Note that the slope of all power variations is the same, i.e., the process is monofractal, as expected. The monofractality of the process is also evident in panel b). The multifractal scaling exponent τ (q) is shown and is purely a linear function. Likewise, in the inset, the generalised Hurst indices h(q) for the three processes for a set of power variations q ∈ [−10, 10] is displayed. The linear shape of τ (q) and constant value of h(q) in the inset indicates, as expected, that these three processes are monofractal. Small deviations are seen for very negative q powers (q 7), which arise due to the numeric (negative) powering operation, which highly depends on the numerical precision of the data.

Lévy-driven process
As a second example, take a collection of Lévy distributed random variables [70]. That is, each X t is drawn independently from a symmetric α-stable distribution, such that the probability density function of X(t) vanishes as a power-law P (x) ∼ |x| −(α−1) for large |x| [70]. These processes exhibit heavy tails, ill-defined variances, and multifractal scaling. In Fig. 3 three symmetric Lévy α-stable distributed processes with α = 1.75, 1.25, and 0.75 are studied with MFDFA.
The multifractal behaviour can be identified directly in panel a), where the fluctuation function F q (s) for α = 1.25 is shown. The lines of F q (s) are not parallel for different positive q powers, showing that the process is not mono-fractal. In fact, the process is bi-fractal, having a separate behaviour for q < 0 and q > α. For positive power variations q > α, the generalised Hurst exponent h(q) decays like 1/q. For values of q < 0, the generalised Hurst exponent h(q) = 1/α. This can be seen clearly in Fig. 3 b), where the generalised Hurst exponent h(q) is displayed. In the inset, one notices that the multifractal scaling exponent τ (q) = 0, for q > α (always zero for q > 2), once again showing that none of these processes are distinguishable for positive power variations.
In general, without the aid of the multifractal spectra, which we uncovered by studying MFDFA for a range of q values, it is not possible to distinguish between Lévy distributed processes. The particular choice of q = 2, i.e., conventional DFA, obscures the fractality of these processes, as they all show a similar scaling for q = 2, i.e., h(q = 2) = 1/2 for all Lévy motions, including (nonfractional) Brownian motions (where α = 2).

B. Real-world data, empirical mode decomposition, and extended DFA
In order to evaluate the efficiency of the algorithm, we test here two real-world data sets. Firstly, using MFDFA we will evaluate the multifractality of sunspots time series, a recurring phenomenon on the Sun's photosphere which can be observed with a telescope [26]. Secondly, we will analyse the German and Autrian spot market intraday quarter-hourly electricity price extracted from the European Power Electricity Exchange (EPEX SPOT) [71,72]. We illustrate the application of two advanced features of the developed python package, moving windows and the masking of missing data points.

Sunspots
The sunspot numbers, also called Wolf numbers, are a rather simple measure of solar activity by counting in a weighted manner the number of groups of sunspots and single sunspots visible from the Earth in the solar photosphere, i.e. it is an integrated measure over space [74,75]. Hence, the sunspot numbers form a time series which has a mean period of about 11 years, but is far from being simply periodic [62]. Solar activity is the result of complex magneto-hydrodynamic processes in the Sun characterised by a highly complex spatio-temporal dynamics. It is of special interest to analyse the rather long series of sunspot numbers in order to explore some relations to the underlying spatio-temporal system.
The emergence of sunspots has a distinct statistics and a multifractal spectrum which has been examined in Ref. [26]. This publication has become a reference for multifractal studies as the data from the ILSO World Data Center, Royal Observatory of Belgium, Brussels is freely available [73]. Here, we will focus on numerical efficiency and how to deal with missing or corrupt data. We utilised another feature integrated in MFDFA that enables an efficient management of missing data points. In python's numpy arrays, missing or corrupt values in a time series can be handled with masked data, which logs the missing data points and takes these into account while performing averages, sums, and power operations. When calculating averages or the variance of a segment, or when taking powers, the masked entries are not taken into account. For the particular application with sunspot time series, which are recorded daily since 1818, there are 3247 missing values, over a total of 74145 entries, i.e., roughly 4.4% of the data is missing. To go around this, simply use The MFDFA will extract the variances as it is possible, taking into account the missing values in the time series. Here we highlight that MFDFA calculated 40 q-powers over 70 segments s in 426 ms ± 11.8 ms.
In Fig. 4 we display the fluctuation function F q (s) for q = −10, −5, −2, 2, 5, 10, in panel a), for s ∈ [70, 1000] These curves are not parallel, suggesting the time series is not monofractal. In panel b) the generalised Hurst exponent h(q) is shown as function of q, and similarly the multifractal scaling exponent τ (q) in inset. The generalised Hurst exponent h(q) is not constant over q and consequently the multifractal scaling exponent τ (q) is not linear, indicating clearly the time series in multifractal. In panel c) we display the singularity spectrum D(α) over the singularity strength α, as given by Eq. (7). The singularity strength α spans a wide range of values, over [1.25, 2.25], indicating how strongly multifractal the time series is. Here recall that a monofractal time series, as the fractional Ornstein-Uhlenbeck previously shown in Fig. 2, has a very narrow range of the singularity strength α, centred around H. For the case of sunspot time series we see a wide range of α, indicating the various scales of the phenomenon. Note as well that h(q) and α are always larger than 1, indicating that this is a non-stationary process.

German and Austrian spot market intraday quarter-hourly electricity price time series
We will examine now a 4-year long time series sampled at 15 minutes of the spot market intraday quarterhourly German and Austrian electricity price [71,72] of this particular data has been performed before, but other multifractal studies of price time series exist [44]. In Wang et al. [44], the authors examine different scaling properties for selected periods of low, regular, and high electricity price for some United States of America's electricity markets in the year 2000 and 2001. Here we propose a different analysis, studying the data and examining a short and long time scale of the data without separating different activity periods.
We know that the 15 minute trading electricity market amounts to a small volume of the overall exchange electricity sold, thus this market serves only electricity producers which can either extract or inject power from the power-grid system in a very fast manner (< 15 minutes) [76,77]. This will lead us to explore to separate scaling phenomena in the data: A short and a long timescale of market activities. The expectation is that the very short-time trading is highly volatile, given the necessity of the power grid in injecting or extracting power is a fairly speedy manner. In the long run, the quarter-hourly market is intrinsically linked to the larger hourly and daily electricity market, which has far less variability, as most of the power is sold in lengthier contracts, stabilising the value of the electricity price. Thus one expects a narrower multifractality at large temporal scales. Multifractality is nevertheless expected, as the system exhibits very large yet seldom negative prices, as well as an occasional four of five-fold increase of the (positive) prices, again occurring seldom and lasting very short periods.
In order to obtain a better statistics of the shorter time scales, we will employ MFDFA's moving windows method previously discussed. The moving window method requires the input of the number of steps used to "move" the windows. For the following example, the window parameter is set to 1, thus each overlapping window is displaced by solely 1 data point. This substantially increases the computational time as each averaging operation is repeated by the number of segments s − 1. For Fig. 5, the total calculation lasted 2 min 11 s ± 4.43 s for the windowed mode, compared with 764 ms ± 9.32 ms with the conventional non-overlapping windows, for 50 s segments. In Fig. 5 we display the MFDFA analysis of the price time series, as previously done for the sunpots in Fig. 4. We perform a similar analysis as above, thus we will condense the technical details and focus on the interpretation. Previous studies point to a clear separation of the scaling behaviour of price time series [43,78,79]. They separate two time scales for periods shorter and longer than 24 hours. These studies used pricing data from before 2004 for the Nordic grid (Nordpool). In our analysis, we similarly separate two scales, between 1 and 48 hours and between 48 hours and 10 days. These are indicated in purple (short timescale) and orange and green (long time scale). We first observe that negative q powers do not exist for the short time scale. This is not unusual, many processes do not show a multifractal spectrum with negative q values. Note that this involves taking negative powers of the average of the variances, which is not always well defined for short s segments. This also served as a threshold to assess the change in the fractal behaviour of the time series. For the large time scales (> 48 hours) the negative powers are well defined, and we can identify the full singularity spectrum D(α), as seen in panel c). We note that the short time scale (< 48 hours) has a very strong multifractality (in purple). The singularity strength α, which we can only extract for positive q values, has a very large breadth, especially with its equivalent for the large time scale (in orange). This is well grounded on the previous arguments of having a very volatile market at these short time scales, thus these results are in line with what is known about this market: The high volatility and occasional burst-into very large electricity prices or into negative prices-generate a wide range of the singularity strength. The long-term stability, connected with the larger intraday and day-ahead electricity markets, makes the process far less multifractal at large temporal scales.

IV. THE MFDFA LIBRARY
The Multifractal Detrended Fluctuation Analysis library MFDFA in python presented is a standalone package based integrally on python's numpy [58]. It can be found in https://github.com/LRydin/MFDFA. It harnesses numpy's vectorised polynomial fittings, making it possible to utilise all computational cores in a computer's processor(s). Additionally, EMD is included as an extra feature, which is integrated into MFDFA by simply installing the python library PyEMD [80]. The conventional plots associated with multifractal analysis, i.e., F q (s) vs. s, h(q) and τ (q) vs. q, and D(α) vs. α, are available as well and require the plotting library matplotlib [81].
The MFDFA library accepts numpy's masked arrays, which is particularly convenient when dealing with time Ornstein-Uhlenbeck as given by Eq. (10). Included are firstorder and third-order polynomial fits, first-order fits with extended DFA, and first-order fits with a moving window with a step size of 5. A comparison with the distributed MATLAB code is included [55]. Tests ran on python 3.8.2 and MATLAB R2020b. MFDFA has a average speed-up compared with the MATLAB code, with a five-fold speed increase for first-order polynomial fits (m=1) and a ×27-fold increase for third-order polynomials fits (m=3). Both codes were tested for 100 segments s and 40 q powers. All tasks were performed on a laptop on two computer cores at 2.9 GHz each.
series with missing data, as exemplified in Sec. III B and Fig. 4.
The MFDFA library offers a considerable speed-up in comparison with the available MATLAB version. The library is fully developed to work with multi-threading, which shows an increase in the performance, while handling time series larger than 10 5 data points. In Fig. 6 we display the performance of the MFDFA library for time series of fractional Ornstein-Uhlenbeck processes given in Eq. (10) of increasing length. The MFDFA operation scales linearly with the number of points of the generated time series. The MFDFA algorithm runs in under 1 second for time series having up to 10 5 datapoints, with a first-order polynomial fittings, 100 segments s and 40 q powers, and outperforms the conventional library in MATLAB by up to a factor of ×10 3 in computational speed.
Estimation error and significance calculation have not been included in the library, as the focus lied on computational speed and the inclusion of several extra features, as discussed.

V. CONCLUSION
We have presented a numerically efficient python implementation of Multifractal Detrended Fluctuation Analysis called MFDFA. MFDFA has found extensive application in the past two decades, yet a reliable, allencompassing open-source software in python does not exist to this date. In this library we have harnessed the most of python's flexibility with handling matricial operations and multi-threaded polynomial fittings. In this implementation we have included some of the more common extensions of MFDFA, including a simple empirical mode decomposition as a mechanism to detrend the data, a moving window to handle very short time series, and the extended Detrended Fluctuation Analysis, which can track a different scaling mechanism for non-stationary time series. The MFDFA library can also handle missing values in the data with the aid of numpy's masked time series.
We have initially turned to two classic numerically generated stochastic processes, fractional Ornstein-Uhlenbeck processes and Lévy-distributed motions, and uncovered their monofractal and multifractal with MFDFA. Subsequently we have studied two real-world time series, the sunspot count from 1818 to 2020 and the quarterhourly electricity price time series from 2015 to 2019. For both we performed a multifractal analysis, unveiling their scaling properties and the strength of their multifractality. We focused on MFDFA's speed, the ability to handled missing data, and the integrated overlapping moving window. The analysis displayed here covered only part of MFDFA's integrated options, thus we leave the user to explore the other implemented methods, as the extended DFA and EMD detrending, as these are more specialised to particular research fields.
We hope with this contribution we open a door to fast MFDFA calculations that can the performed on a local machine without an extensive numerical effort and very long time runs, thus permitting in the future to analyse larger time series.