Source code for neurokit2.complexity.fractal_hurst

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.special

from .fractal_dfa import _fractal_dfa_findscales

[docs] def fractal_hurst(signal, scale="default", corrected=True, show=False): """**Hurst Exponent (H)** This function estimates the Hurst exponent via the standard rescaled range (R/S) approach, but other methods exist, such as Detrended Fluctuation Analysis (DFA, see :func:`fractal_dfa`). The Hurst exponent is a measure for the "long-term memory" of a signal. It can be used to determine whether the time series is more, less, or equally likely to increase if it has increased in previous steps. This property makes the Hurst exponent especially interesting for the analysis of stock data. It typically ranges from 0 to 1, with 0.5 corresponding to a Brownian motion. If H < 0.5, the time-series covers less "distance" than a random walk (the memory of the signal decays faster than at random), and vice versa. The R/S approach first splits the time series into non-overlapping subseries of length n. R and S (sigma) are then calculated for each subseries and the mean is taken over all subseries yielding (R/S)_n. This process is repeated for several lengths *n*. The final exponent is then derived from fitting a straight line to the plot of :math:`log((R/S)_n)` vs :math:`log(n)`. Parameters ---------- signal : Union[list, np.array, pd.Series] The signal (i.e., a time series) in the form of a vector of values. or dataframe. scale : list A list containing the lengths of the windows (number of data points in each subseries) that the signal is divided into. See :func:`fractal_dfa` for more information. corrected : boolean if ``True``, the Anis-Lloyd-Peters correction factor will be applied to the output according to the expected value for the individual (R/S) values. show : bool If ``True``, returns a plot. See Also -------- fractal_dfa Examples ---------- .. ipython:: python import neurokit2 as nk # Simulate Signal with duration of 2s signal = nk.signal_simulate(duration=2, frequency=5) # Compute Hurst Exponent h, info = nk.fractal_hurst(signal, corrected=True, show=True) h References ---------- * Brandi, G., & Di Matteo, T. (2021). On the statistics of scaling exponents and the Multiscaling Value at Risk. The European Journal of Finance, 1-22. * Annis, A. A., & Lloyd, E. H. (1976). The expected value of the adjusted rescaled Hurst range of independent normal summands. Biometrika, 63(1), 111-116. * """ # Sanity checks if isinstance(signal, (np.ndarray, pd.DataFrame)) and signal.ndim > 1: raise ValueError( "Multidimensional inputs (e.g., matrices or multichannel data) are not supported yet." ) n = len(signal) scale = _fractal_dfa_findscales(n, scale) # get individual values for (R/S)_n rs_vals = np.array([_fractal_hurst_rs(signal, window) for window in scale]) # filter NaNs (zeros should not be possible, because if R is 0 then S is also zero) valid = ~np.isnan(rs_vals) rs_vals = rs_vals[valid] n_vals = np.asarray(scale)[valid] # it may happen that no rs vals are left (if all values of data are the same) if len(rs_vals) == 0: raise ValueError( "All (R/S) values are NaN. Check your data, or try a different window length." ) # Transform RS values rs_vals = np.log10(rs_vals) if corrected: rs_vals -= np.log10([expected_rs(n) for n in n_vals]) n_vals = np.log10(n_vals) # fit a line to the logarithm of the obtained (R/S) values poly = np.polyfit(n_vals, rs_vals, 1) h = poly[0] # Get the slope if corrected: h = poly[0] + 0.5 if show: _fractal_hurst_plot(poly, n_vals, rs_vals, corrected=corrected, ax=None) return h, {"Values": n_vals, "Scores": rs_vals, "Corrected": corrected, "Intercept": poly[1]}
# ============================================================================= # Utilities # ============================================================================= def expected_rs(n): """Calculates the expected (R/S)_n for white noise for a given n. This is used as a correction factor in the function hurst_rs. It uses the formula of Anis-Lloyd-Peters. """ # front = (n - 0.5) / n i = np.arange(1, n) back = np.sum(np.sqrt((n - i) / i)) if n <= 340: middle = scipy.special.gamma((n - 1) * 0.5) / np.sqrt(np.pi) / scipy.special.gamma(n * 0.5) else: middle = 1.0 / np.sqrt(n * np.pi * 0.5) return middle * back def _fractal_hurst_rs(signal, window): """Calculates an individual R/S value in the rescaled range approach for a given window size (the size of the subseries in which data should be split). """ n = len(signal) m = n // window # number of sequences # cut values at the end of data to make the array divisible by n signal = signal[: n - (n % window)] # split remaining data into subsequences of length n seqs = np.reshape(signal, (m, window)) # calculate means of subsequences means = np.mean(seqs, axis=1) # normalize subsequences by substracting mean y = seqs - means.reshape((m, 1)) # build cumulative sum of subsequences y = np.cumsum(y, axis=1) # find ranges r = np.max(y, axis=1) - np.min(y, axis=1) # find standard deviation # we should use the unbiased estimator, since we do not know the true mean s = np.std(seqs, ddof=1, axis=1) # some ranges may be zero and have to be excluded from the analysis idx = np.where(r != 0) r = r[idx] s = s[idx] # it may happen that all ranges are zero (if all values in data are equal) if len(r) == 0: return np.nan # return mean of r/s along subsequence index return np.mean(r / s) def _fractal_hurst_plot(poly, n_vals, rs_vals, corrected=False, ax=None): if ax is None: # ax option in case more plots need to be added later fig, ax = plt.subplots() fig.suptitle("Hurst Exponent via Rescaled Range (R/S) Analysis") else: fig = None ax.set_ylabel(r"$log$((R/S)_n)") ax.set_xlabel(r"$log$(n)") # Plot log((R/S)_n) vs log(n) ax.scatter( n_vals, rs_vals, marker="o", zorder=1, label="_no_legend_", ) fit_values = [poly[0] * i + poly[1] for i in n_vals] if corrected: label = "corrected h = {}".format(round(poly[0] + 0.5, 2)) else: label = "h = {}".format(round(poly[0], 2)) ax.plot(n_vals, fit_values, color="#E91E63", zorder=2, linewidth=3, label=label) ax.legend(loc="lower right") return fig # ============================================================================= # Generalized Hurst Exponent # ============================================================================= def _fractal_hurst_generalized(signal, q=2): """TO BE DONE. The Generalized Hurst exponent method is assesses directly the scaling properties of the time series via the qth-order moments of the distribution of the increments. Different exponents `q` are associated with different characterizations of the multi-scaling complexity of the signal. In contrast to the popular R/S statistics approach, it does not deal with max and min functions, and thus less sensitive to outliers. From""" n = len(signal) H = np.zeros((len(range(5, 20)), 1)) k = 0 for Tmax in range(5, 20): x = np.arange(1, Tmax + 1, 1) mcord = np.zeros((Tmax, 1)) for tt in range(1, Tmax + 1): dV = signal[np.arange(tt, n, tt)] - signal[np.arange(tt, n, tt) - tt] VV = signal[np.arange(tt, n + tt, tt) - tt] N = len(dV) + 1 X = np.arange(1, N + 1, dtype=np.float64) Y = VV mx = np.sum(X) / N SSxx = np.sum(X ** 2) - N * mx ** 2 my = np.sum(Y) / N SSxy = np.sum(np.multiply(X, Y)) - N * mx * my cc1 = SSxy / SSxx cc2 = my - cc1 * mx ddVd = dV - cc1 VVVd = VV - np.multiply(cc1, np.arange(1, N + 1, dtype=np.float64)) - cc2 mcord[tt - 1] = np.mean(np.abs(ddVd) ** q) / np.mean(np.abs(VVVd) ** q) mx = np.mean(np.log10(x)) SSxx = np.sum(np.log10(x) ** 2) - Tmax * mx ** 2 my = np.mean(np.log10(mcord)) SSxy = np.sum(np.multiply(np.log10(x), np.transpose(np.log10(mcord)))) - Tmax * mx * my H[k] = SSxy / SSxx k = k + 1 mH = np.mean(H) / q return mH