Fractal Methods in Nonequilibrium Systems

Abstract

The principles of complexity science can be applied to systems with natural hazards, and also human social and economic systems. Disruptive, extreme events result from emergent properties of complex, nonequilibrium systems that consist of interdependent components whose interactions result in a competition between organized, interaction-dominated behavior and irregular or stochastic, fluctuation-dominated behavior. The use of fractal analysis allows insights into the development of such extreme events, and provides input and constraints for models. The main goal of this paper is to share and expand the scope of some of the fractal methods commonly used in complex systems studies. A proper appreciation of the strengths and limitations of fractal methods can improve the assessment and analysis of risks associated with systems that exhibit extreme events.

1. Introduction

Many physical and social systems demonstrate the existence of properties that can be described as complex, featuring behaviors with both regular and irregular components. Complex behavior has been shown to arise when a system has many elements whose links, internal competition for resources, and interdependencies are inherently difficult to analyze. Because of these complex, often nonlinear, interconnections studying one part of the system may not necessarily lead to a better understanding or ability to accurately forecast the future development of any part of the system.

Historically science has progressed through a reductionist approach whereby the system is diminished to its simplest components. Then equations and models are sought which focus on adapting the values of these system parameters to match the observed data. Although this has been an extremely fruitful avenue of research, there is increasing recognition that many systems feature a wide variety of phenomena which do not act in isolation. Rather they should be understood in terms of multiple interactions between different internal system elements, as well as interactions from outside the system.

Whatever the discipline, there is an emerging perception that in our cosmos empirical data frequently demonstrate certain universal patterns of behavior [1, 2]. Part of the difficulty is captured in the probability distributions of system parameters. In equilibrium thermodynamics and statistical mechanics, the internal system dynamics constantly tend towards one or more equilibrium states with simple probability distributions. An example of an equilibrium probability distribution is the unimodal Gaussian distribution function in which most system events are within one standard deviation of the mean. This is the case for the classic Gibbs-Boltzmann statistical mechanics and the Gaussian distribution is one example in which the Gibbs-Boltzmann entropy is maximized to ensure the equilibrium state. Systems whose probability distribution assume the Gaussian form rarely have events that fall outside three standard deviations from the mean. Events like this, which deviate significantly from the mean, can be called ‘extreme’ events [3]. But observations in social and physical dynamics demonstrate that many systems have internal processes or external influences for which physical equilibrium descriptions of macroscopic behavior do not apply; extreme events are much more common than would be anticipated in equilibrium systems. This is most notably observed in heavy-tailed distribution functions which demonstrates that the mean is not an adequate measure of the future behavior [4].

It is not easy to settle on a definition of an extreme event, but there are certain common features [5]. Extreme events occur suddenly and are usually unexpected. They exceed a threshold value close to the upper or lower range of historical data, for instance being defined as those that occur in the top or bottom 5% or so of measurements. Other ways scientists describe extreme events include how different they are from the average, how often they occur, or their likelihood of happening again. As mentioned above, they produce heavy-tailed distributions. Because the probability of appearance of an extreme event is extremely small, the waiting time for its appearance is often too long for all practical purposes; the event occurs on a scale large compared to the system size. Thus, the sudden appearance of an extreme event can be considered exogenous to the system [6]. An open system is by nature a nonequilibrium system since there are exogenous inputs to, and losses from, the system. Thus, establishing the formal Gibbs-Boltzmann entropy can be difficult or impossible.

In nature, the most frequently observed, fascinating, and intricate phenomena are typically collective stochastic systems. These are systems where many individual components interact, with processes characterized by randomness and dissipative energetics. They show long-range correlations, intermiitency, and combinations of features which violate physical equilibrium descriptions of macroscopic behavior. This has led to concepts in terms of nonequilibrium phase transitions [7, 8, 9, 10], and nonextensive generalizations of Gibbs-Boltzmann entropy [11, 12]. The complexity of various physical and social systems is well recognized, arising mainly from the inherent nonlinearity and far from equilibrium nature, and generates fractal geometries in space and time.

The term fractal was introduced by Mandelbrot and Van Ness [13, 14] to describe objects that display self-similarity. Since then data exploration in terms of fractal properties has found regular use in research areas as diverse as economics [15], engineering [16, 17], financial markets [18], geology [19, 20], hydrology [21, 22], medicine [23, 24, 25], physics [25], space science [26, 27, 28] and many more. Fractal properties in nature and human dynamics arguably serve to yield increased understanding and advancement in human society and stewardship.

Dynamic features in natural and social systems often show the scale independent behavior characteristic of fractals [19, 29, 30]. One can define a fractal set according to the following power law,

$$N\propto \frac{1}{r^D}$$E1

where N is the value of an element of the set with characteristic linear dimension r, and D is the associated fractal dimension. The fractal dimension is a real number. It is equivalent to an Euclidean dimension when it has integer values. Not all fractals are self-similar according to the definition of Eq. (1); when the scale invariance depends on the direction, one speaks of self-affinity. A self-affine fractal may be defined in terms of a two-dimensional space (x,y) where the function f(rx,r^Hy) is statistically similar to f(x,y). That is, self-affine fractals demonstrate self-similarity only if different coordinates are scaled by different magnification factors. The quantity r is a scaling factor and H is the dimensionless Hurst exponent, also known as the Hausdorff dimension, which is related to the fractal dimension (D) by H = 2 – D [31]. In the following sections we will provide an overview of various approaches to estimate fractal properties, including the fractal dimension, of various physical and social systems. These approaches use the tools of chaos, nonlinear dynamics, and statistical mechanics to demonstrate the wide applicability to extreme events and risk management.

2. Fractal methods

In the following subsections we will illustrate the use of various fractal methods in a range of nonequilibrium systems. Our main goal is to summarize how fractal methods have wide, interdisciplinary applications in systems that are open, dissipative, and whose source parameters and interactions are not well prescribed, and which have an unknown, but large, number of degrees of freedom.

2.1 Rescaled range (R/S) analysis

Range scaling analysis, also known as the R/S method, was developed by H.E. Hurst [32] to study time series in civil engineering whose underlying processes are independent, though not necessarily Gaussian.

Let us define a time series x(t) over the interval T with values x(t₀), x(t₁), x(t₂), ..., x(t_n) where t₀ = 0, t₁ = τ, t₂ = 2τ, … , t_N = Nτ; N ∈ ℕ. The time interval between measurements is τ, and t_n = T. The mean signal is

$$\overline{x}\left(T\right)=\frac{1}{T}{\int }_0^{T}x\left(t\right)\mathit{dt}$$E2

and signal variance is

$$V\left(T\right)=\frac{1}{T}{\int }_0^T{\left[x\left(t\right)-\overline{x}\left(T\right)\right]}^2\mathit{dt}$$E3

with standard deviation

$$\sigma \left(T\right)={\left[V\left(T\right)\right]}^{1/2}$$E4

A rescaled range is calculated by first normalizing the original data by subtracting the sample mean:

$$Z_r=x\left(t_r\right)-\overline{x};\mathit{rϵ}\mathbb{N}$$E5

Then an integrated, or cumulative, time series, y, is derived from Eq. (5),

$$y_l=\sum _{i=1}^l{Z}_i;l=2,3,\dots ,N,$$E6

and an adjusted range, R(T), is formed in terms of the maximum minus the minimum value of the cumulative series, Y,

$$R\left(T\right)=sup\left(y_1,y_2,\dots ,y_T\right)-inf\left(y_1,y_2,\dots ,y_T\right).$$E7

The rescaled range, R/S, is R(T)/σ(T). It scales with respect to T, by a power law.

$${\left(R/S\right)}_T\propto T^H$$E8

where H is the Hurst exponent [32]. Figure 1 shows a magnetogram (top) from the experiment described by [33]. The magnetic field was measured by the author, at a 10-second cadence, near the Hartebeesthoek Satellite Tracking Station and Radio Astronomy Observatory (HBS), near Pretoria, South Africa on March 1, 1994, using a Geometrics G-856 magnetometer. The bottom subplot shows the log-log plot of the rescaled range (R/S), the slope of which yields a scaling exponent H = 0.39 ± 0.01.

Figure 1.

(Top) Total magnetic field measured at Hartebeesthoek (HBS) in March 1994. (Bottom) Log-log plot of the corresponding rescaled range (R/S). The blue crosses show the data and the red line is the best-fit, 95% confidence intervals.

2.2 Spectral analysis

Fourier transforms are an efficient technique for decomposing a signal from the temporal or spatial domain into the frequency domain. Fractal scaling can be examined, under certain conditions, using a power spectral density (PSD) analysis. If the PSD, given as a function of frequency, follows a power-law, in other words P(ν) ∝ ν^−β. The regularity of the scaling exponent is connected with the physical nature of the signal. That is, the slope β obtained from the plot of power spectral density against frequency provides a measure of the correlation level of the signal. When −1 < β < 1, the signal is termed a fractional Gaussian noise (fGn). fGn can be obtained as successive increments of a signal that displays fractional Brownian motion (fBm), which is a signal where both the real and imaginary components of the Fourier amplitudes are Gaussian-distributed random variables [34, 35]. The mean of the signal Fourier amplitudes $\overline{\phi \left(\nu \right)}=0$ and $\overline{\phi \left(\nu \right)\phi {\left({\nu }^{\prime }\right)}^{\ast }}=P\left(\nu \right)\delta \left(\nu -{\nu }^{\prime }\right).$ This means that for β = 2, fBm is the same as a random walk and the power law spectrum varies as an inverse square. When the scaling exponent is above β = 2, the signal is called persistent, because if the data at some point have x(t_i + 1) > x(t_i), for instance, then the probability is greater than 0.5 that x(t_i + 2) > x(t_i + 1). Signals with exponents below 2 are called antipersistent because if x(t_i + 1) > x(t_i), the probability is greater than 0.5 that x(t_i + 2) < x(t_i + 1). The special case where β = 2 (H = 0.5) indicates Brownian motion. When the signal is fBm, it exhibits power-law scaling with slope 1 < β < 3. In this case the signal is nonstationary but has stationary increments over a range of scales. For fBm one can generally relate the power-spectral scaling index to the Hurst index, introduced in the previous section:

$$\beta =2H+1.$$E9

The blue crosses in Figure 2 show the PSD for the HBS data shown before. This yields a linear best-fit fractal scaling exponent for all the scales of β = 2.21 ± 0.10 (H = 0.61 ± 0.05). This best-fit is not shown in the figure because there is a crossover region in the figure differentiating two different scaling regimes; the best-fit straight line for all frequencies, corresponding to a single monofractal scaling exponent, does not capture the statistics evident in the PSD, particularly at the lower frequencies. Instead, we plot the two red curves which show the best-fits, with the 95% confidence intervals, for low frequency and higher frequency ranges on each side of the crossover. For low frequencies, β = 2.20 ± 0.10 (H = 0.60 ± 0.05) and for the higher frequencies the best-fit yields β = 1.43 ± 0.02 (H = 0.22 ± 0.01). These results show what is often observed in dynamical data from different research areas, viz. that the correlations of recorded data often do not follow the same scaling law for all time scales [36, 37]. This can be because the statistical nature of the data is more complex than a simple fBm [38]. However, in the case illustrated in Figure 2 it may be the result of nonstationarity of the data. The weakness inherent in spectral analysis is related to the problem of data nonstationarity [39]. Nonstationarity means that statistical properties are not constant through the signal, and traditional analysis methods that assume stationarity cannot be used, except with extreme care.

Figure 2.

The power spectral density versus frequency for the Hartebeesthoek (HBS) data (blue crosses). The red lines show the best-fit, 95% confidence intervals, for low frequency and higher frequency ranges. For low frequencies, β = 2.20 (H = 0.60) and for the higher frequencies the best-fit is β = 1.43 (H = 0.22).

Although ubiquitous in science and engineering applications, Fourier analysis is often not suitable for the scrutiny of nonlinear and nonstationary data, which is precisely what is so common in the real world [40, 41]. Fourier techniques, without appropriate preconditioning of data, can yield spurious results for fractal properties. Thus, it behooves the researcher, when studying complex phenomena, whether that be turbulence or any other such systems with intermittent and extreme events, to carefully consider the nature of the data before choosing a method to study fractal properties, since data may feature complex statistics, including multifractional, multifractal and other models [27, 42, 43, 44].

2.3 Higuchi method

Higuchi [45] developed a widely applied time-domain technique to determine fractal properties of nonstationary physical data that is complex and non-periodic [46, 47, 48]. It is deemed a good tool to accurately estimate the fractal dimension and is linked to the fractal dimension derived from Fourier spectra when the signal is a fBm [49]. The Higuchi method was first applied to study large-scale turbulent fluctuations of magnetic fields in space. In that context it allowed modification of existing methods to study the fluctuation properties of turbulent plasma properties beyond the inertial range [50]. It is simple to implement, efficient, and can rapidly achieve accurate and stable values of fractal dimension. Importantly, this is true even in noisy, nonstationary data [51].

The Higuchi method takes a time series, with properties already defined in Eqs. (1)-(4) and, from this series, obtains ${\mathit{X}}_{\mathit{k}}^{\mathit{m}}$, defined as:

$${\mathit{X}}_{\mathit{k}}^{\mathit{m}}:\mathit{x}\left(\mathit{m}\right),\mathit{x}\left(\mathit{m}+\mathit{k}\right),\mathit{x}\left(\mathit{m}+\mathbf{2}\mathit{k}\right),\dots ,\mathit{x}\left(\mathit{m}+\left[\frac{\mathit{N}-\mathit{k}}{\mathit{k}}\right]·\mathit{k})\right),$$E10

Here [] represents the integer part of the enclosed value. The integer $\mathit{m}\mathbf{=}\mathbf{1}\mathbf{,}\mathbf{2}\mathbf{,}\dots \mathbf{,}\mathit{k}$is the zero-point time and $\mathit{k}$ is the time interval, with $\mathit{k}\mathbf{=}\mathbf{1}\mathbf{,}\dots \mathbf{,}{\mathit{k}}_{\mathit{max}}\mathbf{;}{\mathit{k}}_{\mathit{max}}$ is a free tuning parameter. This means that given time interval equal to $\mathit{k}$, spawns $\mathit{k}$-sets of new time series with lengths defined by:

$${\mathit{L}}_{\mathit{m}}\left(\mathit{k}\right)=\frac{\left\{\left({\sum }_{\mathit{i}=\mathbf{1}}^{\left[\frac{\mathit{N}-\mathit{m}}{\mathit{k}}\right]}\left|\mathit{x}\left(\mathit{m}+\mathit{ik}\right)-\mathit{x}(\mathit{m}+\left(\mathit{i}-\mathbf{1}\right)·\mathit{k}\right|\right)\frac{\mathit{N}-\mathbf{1}}{\left[\frac{\mathit{N}-\mathit{m}}{\mathit{k}}\right]·\mathit{k}}\right\}}{\mathit{k}}.$$E11

The length of the curve for the time interval $\mathit{k}$ is then defined as the average over the $\mathit{k}$ sets of ${\mathit{L}}_{\mathit{m}}\left(\mathit{k}\right)\mathbf{,}$ i.e. $\mathit{L}\left(\mathit{k}\right)\mathbf{=}⟨{\mathit{L}}_{\mathit{m}}\left(\mathit{k}\right)⟩\mathbf{.}$ In cases when this equation scales according to the rule $\mathit{L}\left(\mathit{k}\right)\mathbf{\propto }{\mathit{k}}^{\mathbf{-}\mathbf{HF}\mathit{D}}\mathbf{,}$ the time series behaves as a fractal with dimension HFD. A problem with the method, not apparently well-appreciated in the research community, is that a poor choice of the tuning parameter (k_max); can generate spurious results. Recent research has suggested various methodologies to choose an appropriate tuning parameter based on the input data [52, 53, 54, 55]. We follow the methodology of [52], with k_max = 37 to derive the results in Figure 3, which yield HFD = 1.75 ± 0.02, and since H = 2-HFD this gives H = 0.25 ± 0.02.

Figure 3.

Average curve length versus scale size, k, for the HBS time series yielding HFD = 1.75, and corresponding H = 0.25. The blue crosses show the data and the red line is the best-fit, with 95% confidence intervals.

2.4 Structure function analysis

Structure functions have traditionally been used in studies of plasma turbulence [56]. Compared to techniques such as Fourier analysis, structure functions can derive results with irregularly sampled data, and makes no assumptions about the missing data, such as stationarity. In terms of spatial fields, the structure function S_p(x, r) of order p of a field τ(x) is given by:

$$S_p\left(\stackrel{⃑}{x},\stackrel{⃑}{r}\right)=⟨{\left|\tau \left(\stackrel{⃑}{x}+\stackrel{⃑}{r}\right)-\tau \left(\stackrel{⃑}{x}\right)\right|}^p⟩$$E12

Here x is the position and r is a lag with respect to x, and the angled braces, ⟨·⟩, indicate an ensemble average. This can easily be generalized to a time-series. Structure functions of high order (large values of p) yield information concerning the distribution of the large fluctuations of the spatial or temporal field. Smaller fluctuations generally dominate the lower-order structure functions.

If increments of the field τ(x + r) − τ(x) are assumed to be statistically homogeneous, S_p(x, r) is only a function of the lag r and may be written as S_p(r). In this situation, the ergodic theorem implies that a spatial average can be used to estimate the ensemble average in the equation.

For white noise [57], the second-order structure function scales trivially in time as $S_2\sim$c, a constant. In the case of diffusive processes, such as the well-known stochastic process of Brownian motion, the scaling obeys as $S_2\sim \tau$. For fractional Brownian motion (having white noise and Brownian motion as subclasses) the scaling obeys the fractal scaling law $S_2\sim {\tau }^{2H},$ where H is the Hurst exponent. Figure 4 shows the second-order structure function results, which yield H = 0.20 ± 0.01. The curve deviates significantly from a power-law at the smallest and largest scales. The Hurst exponent estimates can clearly be strongly influenced by the choice of the minimum and maximum fitting scale [58]. In addition, the length of the time-series can influence scaling estimates, thus finite sample properties of estimators should not be ignored [55, 58].

Figure 4.

Plot of second-order structure function (blue crosses) and best-fit red line with 95% confidence intervals for the Hartebeesthoek (HBS) data, yielding H = 0.20 ± 0.01.

2.5 Detrended fluctuation analysis (DFA)

The final method we consider is detrended fluctuation analysis (DFA) which was initially introduced to study DNA nucleotides and variations in the human heart rate, which display the kind of long-range correlations typical of nonequilibrium dynamical systems [59, 60]. It has since then been employed in an enormous range of fields including economics, engineering, space science, and many others [61, 62, 63, 64]. The method is a modified root mean squared analysis of a random walk designed specifically to be able to deal with nonstationarities in nonlinear data, and is considered to be among the most robust of statistical techniques designed to detect long-range correlations in time series [65, 66, 67]. Because DFA focuses on fluctuations around trend, rather than the signal range, it has been shown to be robust to the presence of trends [68] and nonstationary time series [37, 69].

Briefly, the methodology begins by removing the mean,$\bar{x}$, from the time series, x(t), and then integrating as in Eqs. (1)-(6). At this point the methodology deviates from the R/S analysis. Now the new time-series is divided into boxes of equal length, n. A trend, represented by a least-squares fit to the data, is removed from each box; the trend is typically a linear, quadratic, or cubic function. Box-size n has its abscissa denoted by $y_n\left(k\right)$. Next the trend is removed from the integrated time series, y(k), by subtracting the local trend, $y_n\left(k\right)$, in each box.

For a given box size n, the characteristic size of the fluctuations, designated by the fluctuation function F(n), is calculated as the RMS deviation between y(k) and the trend in each box.

$$F\left(n\right)=\sqrt{\frac{1}{N}{\sum }_{k=1}^N{\left[y\left(k\right)-y_n\left(k\right)\right]}^2}$$E13

The calculation in Eq. (13) is finally performed over all time scales (box sizes), to potentially establish a power-law scaling between F(n) and n, viz.

$$F\left(n\right)\propto n^{\alpha }$$E14

Here the $\alpha$ parameter is a scaling exponent. Given a fractional Brownian motion series, the scaling exponent $\alpha$ = H – the familiar Hurst exponent. Figure 5 shows the fluctuation function results, which yield H = 0.19 ± 0.02. At the largest scales there is considerable larger standard deviations for F(n), than, for instance, compared to the R/S results in Figure 1. However, because of difference in methodology and finite sample effects the R/S estimate yields much larger Hurst exponent.

Figure 5.

Blue crosses show the fluctuation versus box size for the HBS data obtained by detrended fluctuation analysis. The best-fit red line with 95% confidence intervals has a slope corresponding to H = 0.19 ± 0.01.

3. Summary and conclusion

In this paper we have looked at the use cases, in multiple fields of research, for various estimators of fractal long-range dependence. Analysis of data in terms of complexity science can help understand systems with natural hazards, and also human social and economic systems, where disruptive, extreme events result from emergent properties. This is because the value of system parameters at a certain time is related not just to the immediate past, but also to fluctuations in the remote past. Fractal analysis can thus allow insights into the origin and development of extreme events, and provide limits for theoretical models. The main goal of this paper has been to share and expand the scope of some of the fractal methods commonly used in complex systems studies. A proper appreciation of the strengths and limitations of fractal methods can improve the assessment and analysis of risks associated with systems that exhibit extreme events.

The different classical estimators of fractal scaling behavior we considered were the R/S, spectral, Higuchi, structure function, and DFA methods. Table 1 gives a comparison of the fractal scaling exponent results from the five different methods that we considered for the analysis of a short time series (8500 data points) of a ground magnetic field measured by the author near the Hartebeesthoek Satellite Tracking Station and Radio Astronomy Observatory (HBS), near Pretoria, South Africa on March 1, 1994.

Method	Hurst exponent, H
R/S	0.39 ± 0.01
Spectral	0.60 ± 0.05/0.22 ± 0.01
Higuchi	0.25 ± 0.02
Structure Function	0.20 ± 0.01
DFA	0.19 ± 0.02

Table 1.

Comparison of Hurst exponent from different fractal methods.

The results we have obtained show how different Hurst scaling values may be obtained by employing different fractal methods on the same input data. The PSD analysis appears to show two different scaling regimes and a crossover (Figure 2), which can arise due to nonstationarity in the data. This illustrates the importance of determining first whether the data to be studied are stationary or nonstationary. Techniques, such as R/S and spectral analysis assume stationary data and can produce spurious and misleading results if applied without care on nonstationary data. This is evident in the large difference between the estimate for the Hurst exponent for the R/S method, compared to the other methods.

Another issue, demonstrated in Figure 4, relates to the length of the data series under investigation and other finite sample properties of estimators. Figure 4 shows the second-order structure function results can deviate significantly from a power-law at the smallest and largest scales. Thus, Hurst exponent estimates can be influenced by the choice of the minimum and maximum fitting scale.

Next, although a single scaling law is possible, over a wide or infinite range of scales, due to the system dynamics there can be a multitude of actual scaling exponents and scaling regimes. In this case the system is not simply fractal, but multifractal [70, 71]. In such a complex system phase transitions in the regulation behavior can be connected to changes in internal fractal dynamics. This supports various hypotheses that connect extreme events to dynamical phase transitions in out of equilibrium systems [9, 61, 72].

Although we have not undertaken to determine which methods are most accurate, the cautions we have highlighted are important for researchers in the field to consider. Various studies have done relevant comparative studies of the different methods. For instance, some authors have used Monte Carlo methods of very specific synthetically generated data to test method accuracies at capturing scaling exponents [55, 58]. Others have considered the effects of trends, the length of the series, periodicities, and other considerations that influence determination of scaling behaviors [64, 66, 67, 68, 73, 74, 75]. Finally, there are various other methods of fractal dimension calculation that we have not discussed in this paper. These include the basic box-counting, dispersional, Hall-Wood, semi-periodogram, discrete cosine transform, and wavelet estimators [44, 76, 77, 78].

Acknowledgments

Research reported in this publication was supported by the National Science Foundation under Award Number AGS-2053689, and the National Institute of General Medical Sciences of the National Institutes of Health under Award Number P20GM103499. The content is solely the responsibility of the authors and does not necessarily represent the official views of any funding agency.