## 3 MEDFA Algorithm  
In the early research studies, many researchers tried to remove the periodicity and trend in a time series to determine the true scale exponents. However, the removal of periodicity inevitably leads to unintended and subjective modification or a smoothing of the fluctuation. This deviation can be overcome by using MFDFA, since the user can choose fitting order \(m\), moment of variances \(q\) and scale \(s\). MEDFA is a well-established method to detect multifractality and the scaling behavior of noisy data in the presence of trends with unknown origin and shape. Detrended fluctuation analysis (DFA) was first proposed by [PERSON] et al. for detecting the long-range correlations of DNA sequences in 1995 [28]. Afterwards, an extension of DFA called multifractal detrending fluctuation analysis (MFDFA) was proposed by [PERSON] et al. for examining the multifractality of non-stationary time series in 2002 [29]. Currently, MFDFA has been successfully applied to analyze various data, such as hydrographic data [30], wind records [31], financial time series [32], traffic time series [33], control system assessment [34], mechanical vibration signals [35], etc. It has proven to be a powerful tool for uncovering the multifractality of non-stationary time series in the complex systems. The MFDFA method starts with a possibly non-stationary time series \(X_{i}\) for \(i=1,\ldots,N\), where \(N\) indicates its length. 1. Transform original data into mean-reduced cumulative sums, \[Y_{j}=\sum_{i=1}^{j}(X_{i}-\bar{X}),\ j=1,\ldots,N,\] (9) where \(\bar{X}\) is the mean of series, such that the aggregated time series are with zero mean. 2. Divide time series \(Y_{j}\) into \(N_{s}=\text{int}(N/s)\) non-overlapping segments of equal length \(s\), starting from the beginning. Since the length \(N\) of the series is often not a multiple of the considered time scale \(s\), to not miss any data, another set of segments starting from the end of data is made. As a result, \(2N_{s}\) segments are obtained covering the whole dataset. 3. Calculate the local trend for each of the segments \(k=1,\ldots,2N_{s}\) by a least-square fit of the series. Calculate the mean square error \(F^{2}(k,s)\) for the estimate of each segment \(k\) of length \(s\). \[F^{2}(k,s)=\frac{1}{s}\sum_{i=1}^{s}\left(E[(k-1)s+i]-p_{k}[i]\right)^{2},\] (10) for each segment \(k=1,\ldots,N_{s}\) and \[F^{2}(k,s)=\frac{1}{s}\sum_{i=1}^{s}\left(E[N-(k-N_{s})s+i]-p_{k}[i]\right)^{2},\] (11) for each segment \(k=N_{s}+1,\ldots,2N_{s}\). 5. Average all segments to obtain the \(q\)th order variance (or fluctuation) function \(F_{q}(s)\) for each size \(s\): \[F_{q}(s)=\left(\frac{1}{2N_{s}}\sum_{k=1}^{2N_{s}}[F^{2}(k,s)]^{q/2}\right)^{1 /q}.\] (12) For \(q=0\) use \[F_{0}(s)=\exp\left\{\frac{1}{4N_{s}}\sum_{k=1}^{2N_{s}}\ln(F^{2}(k,s))\right\},\] (13)
6. Repeat Steps (2)-(5) for different \(s\) evaluating new sets of variances \(F_{q}(s)\). 7. Plot \(F_{q}(s)\) for each \(q\) in log-log scale and estimate the linear fit with least squares. If slope \(h(q)\) varies with \(q\), multifractality is suspected. Single slope shows monofractal scaling. 8. Calculate multifractal exponent \(t(q)\) as \[t(q)=qh(q)-1.\] (14)
9. Use Legendre transform to evaluate the \(q\)-order singularity-Holder exponent \(h(q)\) and corresponding dimension \(D(q)\): \[ \begin{cases}h(q)=dt(q)/dq\\ D(q)=qh(q)-t(q)\end{cases} \] (15)  
For the above steps of calculation methodology, [PERSON] developed the MATLAB code in [36]. The approach is used to calculate multifractal spectrum, with graph \(\log(F_{q})\) versus \(\log(s)\) identifying crossovers with different \(q\) orders. The slope \(H_{q}\) of scaling function \(F_{q}\) with different \(q\) order demonstrates the LRD behavior of signals. The multifractal spectrum indicates how much dominant are the various fractal exponents present in the series.