""" hurst-calculation.py ===================== Computes the Hurst exponent H (rescaled-range, R/S) and the DFA scaling exponent alpha on the daily Q drop count series. References ---------- - Hurst, H.E. (1951). "Long-term storage capacity of reservoirs." Transactions of the American Society of Civil Engineers, 116, 770-808. - Mandelbrot, B.B. & Wallis, J.R. (1968-1969). Series of papers on fractional Gaussian noise and the R/S statistic. - Peng, C.-K. et al. (1994). "Mosaic organization of DNA nucleotides." Physical Review E, 49(2), 1685-1689. (DFA introduction) - Kantelhardt, J.W. et al. (2002). "Multifractal detrended fluctuation analysis of nonstationary time series." Physica A, 316, 87-114. Interpretation -------------- - H = 0.5: random walk (white-noise increments, no memory) - H > 0.5: persistent / long-range positive correlation - H < 0.5: anti-persistent (mean-reverting) - alpha (DFA) maps 1:1 with H for stationary fGn under the same caveats. Methodological caveats are listed in README.md. This script computes the two estimators on whatever CSV is provided. It does not adjust for the selection bias of the corpus, the calibration window, or the social-media publication bias. Dependencies: numpy, pandas (Python 3.10+) """ from __future__ import annotations import math from pathlib import Path import numpy as np import pandas as pd CSV_PATH = Path(__file__).parent / "q-drops-dates.csv" N_BOOTSTRAP = 1000 def load_daily_series(csv_path: Path) -> np.ndarray: """Load drops, return the daily count series (one entry per calendar day).""" df = pd.read_csv(csv_path, parse_dates=["date"]) counts = df.groupby(df["date"].dt.date).size() full_index = pd.date_range(counts.index.min(), counts.index.max(), freq="D") series = counts.reindex(full_index.date, fill_value=0).to_numpy(dtype=float) return series, df["date"].min(), df["date"].max() def hurst_rs(x: np.ndarray, min_window: int = 8, max_window: int | None = None) -> float: """ Rescaled-range (R/S) Hurst exponent (Mandelbrot & Wallis 1968). For each window size n in a logarithmic grid: 1) Cut the series into k = floor(N/n) non-overlapping windows. 2) For each window: mean-adjust, take cumulative sum Y, compute R = max(Y) - min(Y), S = std(window). 3) Average R/S across the k windows -> (R/S)_n. The Hurst exponent is the slope of log((R/S)_n) vs log(n). """ x = np.asarray(x, dtype=float) N = len(x) if max_window is None: max_window = N // 4 # log-spaced window sizes ns = np.unique( np.floor( np.logspace(np.log10(min_window), np.log10(max_window), num=20) ).astype(int) ) ns = ns[ns >= min_window] rs_means = [] for n in ns: k = N // n if k < 1: continue rs_vals = [] for i in range(k): chunk = x[i * n : (i + 1) * n] mean = chunk.mean() y = np.cumsum(chunk - mean) r = y.max() - y.min() s = chunk.std(ddof=0) if s > 0: rs_vals.append(r / s) if rs_vals: rs_means.append((n, np.mean(rs_vals))) arr = np.array(rs_means) log_n = np.log(arr[:, 0]) log_rs = np.log(arr[:, 1]) slope, _intercept = np.polyfit(log_n, log_rs, 1) return float(slope) def dfa(x: np.ndarray, min_window: int = 8, max_window: int | None = None) -> float: """ Detrended Fluctuation Analysis exponent alpha (Peng et al. 1994). 1) Build the integrated profile Y = cumsum(x - mean(x)). 2) For each window size n: split Y into N/n non-overlapping segments, fit a linear trend on each segment, compute the RMS of the residuals across all segments -> F(n). 3) alpha is the slope of log(F(n)) vs log(n). """ x = np.asarray(x, dtype=float) N = len(x) if max_window is None: max_window = N // 4 profile = np.cumsum(x - x.mean()) ns = np.unique( np.floor( np.logspace(np.log10(min_window), np.log10(max_window), num=20) ).astype(int) ) ns = ns[ns >= min_window] f_vals = [] for n in ns: k = N // n if k < 1: continue residuals_squared = [] for i in range(k): seg = profile[i * n : (i + 1) * n] t = np.arange(n) slope, intercept = np.polyfit(t, seg, 1) trend = slope * t + intercept residuals_squared.append(np.mean((seg - trend) ** 2)) f_vals.append((n, math.sqrt(np.mean(residuals_squared)))) arr = np.array(f_vals) log_n = np.log(arr[:, 0]) log_f = np.log(arr[:, 1]) slope, _intercept = np.polyfit(log_n, log_f, 1) return float(slope) def bootstrap_pvalue(series: np.ndarray, observed_h: float, n_iter: int = N_BOOTSTRAP) -> float: """ Monte Carlo permutation test. Null hypothesis: the daily values are exchangeable (no temporal order). For each permutation, recompute H. p-value = fraction of permuted H that are >= observed_h. """ rng = np.random.default_rng(seed=42) count_ge = 0 for _ in range(n_iter): permuted = rng.permutation(series) h_perm = hurst_rs(permuted) if h_perm >= observed_h: count_ge += 1 return (count_ge + 1) / (n_iter + 1) # Laplace correction def main(): print(f"Loading {CSV_PATH.name} ...") series, dmin, dmax = load_daily_series(CSV_PATH) n_drops = int(series.sum()) n_days = len(series) print(f"Loaded {n_drops} drops between {dmin.date()} and {dmax.date()}") print(f"Daily series: {n_days} days") print() h = hurst_rs(series) a = dfa(series) print(f"Hurst R/S exponent H = {h:.4f}") print(f"DFA scaling exponent alpha = {a:.4f}") print() print(f"Bootstrap Monte Carlo (n={N_BOOTSTRAP} permutations of the daily series)") print(" H_0: daily values are exchangeable (random ordering).") p = bootstrap_pvalue(series, h, n_iter=N_BOOTSTRAP) if p <= 1.0 / N_BOOTSTRAP: print(f" p-value < {1.0 / N_BOOTSTRAP:.4f} (no permutation reached observed H)") else: print(f" p-value = {p:.4f}") print() if h > 0.7: print("Conclusion: H significantly different from 0.5 (random walk).") print(" Series shows long-range persistence.") elif h > 0.55: print("Conclusion: weak persistence (H above 0.5 but below threshold).") else: print("Conclusion: H consistent with random walk (no long-range memory).") if __name__ == "__main__": main()