Estimators API

lrdbenchmark provides a comprehensive suite of 18 estimators for detecting and quantifying long-range dependence in time series data.

Base Estimator

class lrdbenchmark.analysis.base_estimator.BaseEstimator(**kwargs: Any)[source]

Bases: object

Minimal common interface for estimators.

Conventions - parameters: dict of user-specified configuration - results: dict of computed outputs. At minimum for H estimators:

  • ‘hurst_parameter’: float

  • ‘method’: str

Optional keys: ‘r_squared’, ‘confidence_interval’, ‘diagnostics’

__init__(**kwargs: Any) → None[source]
set_params(**kwargs: Any) → BaseEstimator[source]
get_params() → Dict[str, Any][source]
get_results() → Dict[str, Any][source]
get_parameters() → Dict[str, Any][source]

Alias for get_params() for backward compatibility.

set_parameters(**kwargs: Any) → None[source]

Alias for set_params() for backward compatibility.

_validate_parameters() → None[source]

Validate estimator parameters. Override in subclasses if needed.

Temporal Estimators

Detrended Fluctuation Analysis (DFA)

class lrdbenchmark.analysis.temporal.dfa_estimator.DFAEstimator(min_scale: int = 10, max_scale: int | None = None, num_scales: int = 10, order: int = 1, use_optimization: str = 'auto')[source]

Bases: BaseEstimator

Unified DFA Estimator with modular backend support.

Input convention. Pass a stationary increment series (e.g. fractional Gaussian noise or first differences of fractional Brownian motion). For such processes the integrated profile scales as \(F(s) \sim s^{H}\) in the mid-scale range, and the hurst_parameter field is the slope of \(\log F(s)\) versus \(\log s\) (i.e. the DFA scaling exponent, identified with \(H\) for fGn-like inputs). Feeding fBm levels without differencing does not satisfy the same scaling law and can mislead the estimate.

Backends: - ‘jax’: GPU/TPU accelerated (if available) - ‘numba’: CPU JIT compiled - ‘numpy’: Reference implementation

__init__(min_scale: int = 10, max_scale: int | None = None, num_scales: int = 10, order: int = 1, use_optimization: str = 'auto')[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter using DFA. Delegates calculation to the selected backend.

The returned hurst_parameter is the log-log slope of the fluctuation function vs. scale; see the class docstring for the recommended input (increments / fGn vs. fBm levels).

__init__(min_scale: int = 10, max_scale: int | None = None, num_scales: int = 10, order: int = 1, use_optimization: str = 'auto')[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter using DFA. Delegates calculation to the selected backend.

The returned hurst_parameter is the log-log slope of the fluctuation function vs. scale; see the class docstring for the recommended input (increments / fGn vs. fBm levels).

_get_compute_function(backend: str)[source]

Factory method for backend strategy.

get_optimization_info() → Dict[str, Any][source]
plot_analysis(figsize: Tuple[int, int] = (12, 8), save_path: str | None = None) → None[source]

Plot the DFA analysis results.

get_method_recommendation(n: int) → Dict[str, Any][source]

Get method recommendation for a given data size.

Detrended Moving Average (DMA)

class lrdbenchmark.analysis.temporal.dma_estimator.DMAEstimator(min_scale: int | None = None, max_scale: int | None = None, num_scales: int | None = None, use_optimization: str = 'auto', *, min_window_size: int | None = None, max_window_size: int | None = None, num_windows: int | None = None, window_sizes: Sequence[int] | None = None, overlap: bool = True)[source]

Bases: BaseEstimator

Unified DMA (Detrended Moving Average) Estimator for Long-Range Dependence Analysis.

DMA analyzes the root-mean-square fluctuation of detrended time series data using moving average detrending to estimate the Hurst parameter.

Features: - Automatic optimization framework selection (JAX, Numba, NumPy) - GPU acceleration with JAX when available - JIT compilation with Numba for CPU optimization - Graceful fallbacks when optimization frameworks fail

Parameters:

  • min_scale (int, optional (default=10)) – Minimum scale for analysis

  • max_scale (int, optional (default=None)) – Maximum scale for analysis. If None, uses data length / 4

  • num_scales (int, optional (default=10)) – Number of scales to test

  • use_optimization (str, optional (default='auto')) – Optimization framework to use: ‘auto’, ‘jax’, ‘numba’, ‘numpy’

__init__(min_scale: int | None = None, max_scale: int | None = None, num_scales: int | None = None, use_optimization: str = 'auto', *, min_window_size: int | None = None, max_window_size: int | None = None, num_windows: int | None = None, window_sizes: Sequence[int] | None = None, overlap: bool = True)[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate the Hurst parameter using DMA with automatic optimization.

Parameters:

data (array-like) – Input time series data

Returns:

Dictionary containing estimation results including: - hurst_parameter: Estimated Hurst parameter - r_squared: R-squared value of the fit - scales: Scales used in the analysis - fluctuation_values: Fluctuation values for each scale - log_scales: Log of scales - log_fluctuations: Log of fluctuation values

Return type:

dict

__init__(min_scale: int | None = None, max_scale: int | None = None, num_scales: int | None = None, use_optimization: str = 'auto', *, min_window_size: int | None = None, max_window_size: int | None = None, num_windows: int | None = None, window_sizes: Sequence[int] | None = None, overlap: bool = True)[source]
_select_optimization_framework(use_optimization: str) → str[source]

Select the optimal optimization framework.

_validate_parameters() → None[source]

Validate estimator parameters.

estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate the Hurst parameter using DMA with automatic optimization.

Parameters:

data (array-like) – Input time series data

Returns:

Dictionary containing estimation results including: - hurst_parameter: Estimated Hurst parameter - r_squared: R-squared value of the fit - scales: Scales used in the analysis - fluctuation_values: Fluctuation values for each scale - log_scales: Log of scales - log_fluctuations: Log of fluctuation values

Return type:

dict

_estimate_numpy(data: ndarray) → Dict[str, Any][source]

NumPy implementation of DMA estimation.

_estimate_numba(data: ndarray) → Dict[str, Any][source]

Numba-optimized implementation of DMA estimation.

_estimate_jax(data: ndarray) → Dict[str, Any][source]

JAX-optimized implementation of DMA estimation.

_calculate_fluctuation_numpy(data: ndarray, scale: int, overlap: bool) → float[source]

Calculate fluctuation value for a given scale using NumPy.

get_optimization_info() → Dict[str, Any][source]

Get information about available optimizations and current selection.

get_confidence_intervals(confidence_level: float = 0.95) → Dict[str, Tuple[float, float]][source]
get_estimation_quality() → Dict[str, Any][source]
plot_scaling(**kwargs) → None[source]
_calculate_fluctuation(data: ndarray | list, window_size: int) → float[source]

Backward-compatible helper for direct fluctuation calculation.

Get the recommended optimization framework.

static _should_suppress_fallback_warning(error: Exception) → bool[source]

Return True when a fallback is expected and shouldn’t raise a warning.

plot_analysis(figsize: Tuple[int, int] = (12, 8), save_path: str | None = None) → None[source]

Plot the DMA analysis results.

get_method_recommendation(n: int) → Dict[str, Any][source]

Get method recommendation for a given data size.

Higuchi Method

class lrdbenchmark.analysis.temporal.higuchi_estimator.HiguchiEstimator(min_k: int = 2, max_k: int | None = None, k_values: List[int] | None = None, use_optimization: str = 'auto')[source]

Bases: BaseEstimator

Estimate fractal dimension and Hurst exponent using the Higuchi method.

__init__(min_k: int = 2, max_k: int | None = None, k_values: List[int] | None = None, use_optimization: str = 'auto') → None[source]
estimate(data: ndarray | List[float], copy: bool = True) → Dict[str, Any][source]
__init__(min_k: int = 2, max_k: int | None = None, k_values: List[int] | None = None, use_optimization: str = 'auto') → None[source]
estimate(data: ndarray | List[float], copy: bool = True) → Dict[str, Any][source]
get_confidence_intervals(confidence_level: float = 0.95) → Dict[str, Tuple[float, float]][source]
get_estimation_quality() → Dict[str, Any][source]
get_optimization_info() → Dict[str, Any][source]
plot_scaling(**kwargs) → None[source]
plot_results(save_path: str | None = None) → None[source]

Generalised Hurst Exponent (GHE)

class lrdbenchmark.analysis.temporal.ghe_estimator.GHEEstimator(**kwargs)[source]

Bases: BaseEstimator

Unified GHE (Generalized Hurst Exponent) Estimator for Long-Range Dependence Analysis.

The GHE estimator analyzes the scaling properties of time series data by computing q-th order moments of increments and estimating the generalized Hurst exponent.

Features: - Automatic framework selection (JAX, Numba, NumPy) - Multiple q values for multifractal analysis - Robust error handling and fallback mechanisms - Comprehensive result reporting - Visualization capabilities

Parameters:

  • q_values (array-like, optional) – Array of q values for multifractal analysis. Default: [1, 2, 3, 4, 5]

  • tau_min (int, optional) – Minimum time lag. Default: 2

  • tau_max (int, optional) – Maximum time lag. Default: min(N//4, 50) where N is data length

  • tau_step (int, optional) – Step size for time lags. Default: 1

  • use_jax (bool, optional) – Force JAX backend if available. Default: None (auto-select)

  • use_numba (bool, optional) – Force Numba backend if available. Default: None (auto-select)

__init__(**kwargs)[source]

Initialize the GHE estimator.

Parameters:

**kwargs (dict) – Estimator parameters

estimate(data: ndarray) → Dict[str, Any][source]

Estimate the generalized Hurst exponent using the GHE method.

Parameters:

data (np.ndarray) – Time series data to analyze

Returns:

Dictionary containing estimation results

Return type:

Dict[str, Any]

__init__(**kwargs)[source]

Initialize the GHE estimator.

Parameters:

**kwargs (dict) – Estimator parameters

_validate_parameters() → None[source]

Validate estimator parameters.

_compute_qth_moments_numpy(data: ndarray, q_values: ndarray, tau_values: ndarray) → ndarray[source]

Compute q-th order moments using NumPy.

Parameters:

  • data (np.ndarray) – Time series data

  • q_values (np.ndarray) – Array of q values

  • tau_values (np.ndarray) – Array of time lags

Returns:

Array of shape (len(q_values), len(tau_values)) containing q-th moments

Return type:

np.ndarray

_compute_qth_moments_numba(data: ndarray, q_values: ndarray, tau_values: ndarray) → ndarray[source]

Compute q-th order moments using Numba.

Parameters:

  • data (np.ndarray) – Time series data

  • q_values (np.ndarray) – Array of q values

  • tau_values (np.ndarray) – Array of time lags

Returns:

Array of shape (len(q_values), len(tau_values)) containing q-th moments

Return type:

np.ndarray

_compute_qth_moments_jax(data: ndarray, q_values: ndarray, tau_values: ndarray) → ndarray[source]

Compute q-th order moments using JAX.

Parameters:

  • data (np.ndarray) – Time series data

  • q_values (np.ndarray) – Array of q values

  • tau_values (np.ndarray) – Array of time lags

Returns:

Array of shape (len(q_values), len(tau_values)) containing q-th moments

Return type:

np.ndarray

_estimate_hurst_exponents(tau_values: ndarray, moments: ndarray, q_values: ndarray) → Tuple[ndarray, ndarray, ndarray][source]

Estimate Hurst exponents from q-th moments.

Parameters:

  • tau_values (np.ndarray) – Array of time lags

  • moments (np.ndarray) – Array of q-th moments

  • q_values (np.ndarray) – Array of q values

Returns:

Hurst exponents, R-squared values, and standard errors

Return type:

Tuple[np.ndarray, np.ndarray, np.ndarray]

_select_backend(data_length: int) → str[source]

Select the best available backend for computation.

Parameters:

data_length (int) – Length of the input data

Returns:

Selected backend (‘jax’, ‘numba’, or ‘numpy’)

Return type:

str

estimate(data: ndarray) → Dict[str, Any][source]

Estimate the generalized Hurst exponent using the GHE method.

Parameters:

data (np.ndarray) – Time series data to analyze

Returns:

Dictionary containing estimation results

Return type:

Dict[str, Any]

plot_scaling_behavior(figsize: Tuple[int, int] = (12, 8)) → Figure[source]

Plot the scaling behavior of q-th moments.

Parameters:

figsize (Tuple[int, int], optional) – Figure size. Default: (12, 8)

Returns:

The created figure

Return type:

matplotlib.figure.Figure

get_multifractal_spectrum() → Dict[str, ndarray][source]

Compute the multifractal spectrum from generalized Hurst exponents.

Returns:

Dictionary containing multifractal spectrum parameters

Return type:

Dict[str, np.ndarray]

R/S Analysis

class lrdbenchmark.analysis.temporal.rs_estimator.RSEstimator(min_block_size: int | None = None, max_block_size: int | None = None, num_blocks: int | None = None, use_optimization: str = 'auto', *, min_window_size: int | None = None, max_window_size: int | None = None, num_windows: int | None = None, window_sizes: Sequence[int] | None = None, overlap: bool = False)[source]

Bases: BaseEstimator

Unified R/S (Rescaled Range) Estimator for Long-Range Dependence Analysis.

Backends: - ‘jax’: GPU/TPU accelerated (if available) - ‘numba’: CPU JIT compiled - ‘numpy’: Reference implementation

__init__(min_block_size: int | None = None, max_block_size: int | None = None, num_blocks: int | None = None, use_optimization: str = 'auto', *, min_window_size: int | None = None, max_window_size: int | None = None, num_windows: int | None = None, window_sizes: Sequence[int] | None = None, overlap: bool = False) → None[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter using R/S analysis.

__init__(min_block_size: int | None = None, max_block_size: int | None = None, num_blocks: int | None = None, use_optimization: str = 'auto', *, min_window_size: int | None = None, max_window_size: int | None = None, num_windows: int | None = None, window_sizes: Sequence[int] | None = None, overlap: bool = False) → None[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter using R/S analysis.

_resolve_block_sizes(n: int) → ndarray[source]

Construct block sizes.

get_optimization_info() → Dict[str, Any][source]
get_confidence_intervals(confidence_level=0.95)[source]

Spectral Estimators

Geweke-Porter-Hudak (GPH)

class lrdbenchmark.analysis.spectral.gph_estimator.GPHEstimator(min_freq_ratio: float = 0.01, max_freq_ratio: float = 0.1, apply_bias_correction: bool = True, use_optimization: str = 'auto')[source]

Bases: BaseEstimator

Unified Geweke-Porter-Hudak (GPH) Hurst parameter estimator.

This estimator uses log-periodogram regression with the regressor log(4*sin^2(ω/2)) to estimate the fractional differencing parameter d, then converts to Hurst parameter as H = d + 0.5.

Features: - Automatic optimization framework selection (JAX, Numba, NumPy) - GPU acceleration with JAX when available - JIT compilation with Numba for CPU optimization - Graceful fallbacks when optimization frameworks fail

Parameters:

  • min_freq_ratio (float, optional (default=0.01)) – Minimum frequency ratio (relative to Nyquist) for fitting.

  • max_freq_ratio (float, optional (default=0.1)) – Maximum frequency ratio (relative to Nyquist) for fitting.

  • apply_bias_correction (bool, optional (default=True)) – Whether to apply bias correction for finite sample effects.

  • use_optimization (str, optional (default='auto')) – Optimization framework to use: ‘auto’, ‘jax’, ‘numba’, ‘numpy’

__init__(min_freq_ratio: float = 0.01, max_freq_ratio: float = 0.1, apply_bias_correction: bool = True, use_optimization: str = 'auto')[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter using GPH method with automatic optimization.

Parameters:

data (array-like) – Time series data.

Returns:

Dictionary containing estimation results including: - hurst_parameter: Estimated Hurst parameter - d_parameter: Estimated fractional differencing parameter - intercept: Intercept of the linear fit - r_squared: R-squared value of the fit - m: Number of frequency points used in fitting - log_regressor: Log regressor values - log_periodogram: Log periodogram values

Return type:

dict

__init__(min_freq_ratio: float = 0.01, max_freq_ratio: float = 0.1, apply_bias_correction: bool = True, use_optimization: str = 'auto')[source]
_validate_parameters() → None[source]

Validate estimator parameters.

estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter using GPH method with automatic optimization.

Parameters:

data (array-like) – Time series data.

Returns:

Dictionary containing estimation results including: - hurst_parameter: Estimated Hurst parameter - d_parameter: Estimated fractional differencing parameter - intercept: Intercept of the linear fit - r_squared: R-squared value of the fit - m: Number of frequency points used in fitting - log_regressor: Log regressor values - log_periodogram: Log periodogram values

Return type:

dict

_estimate_numpy(data: ndarray) → Dict[str, Any][source]

NumPy implementation of GPH estimation.

_estimate_jax(data: ndarray) → Dict[str, Any][source]

JAX-optimized implementation of GPH estimation.

_estimate_numba(data: ndarray) → Dict[str, Any][source]

Numba-optimized implementation of GPH estimation.

get_optimization_info() → Dict[str, Any][source]

Get information about available optimizations and current selection.

Get the recommended optimization framework.

plot_scaling(save_path: str | None = None) → None[source]

Plot the scaling relationship and PSD.

get_method_recommendation(n: int) → Dict[str, Any][source]

Get method recommendation for a given data size.

Periodogram

class lrdbenchmark.analysis.spectral.periodogram_estimator.PeriodogramEstimator(min_freq_ratio: float = 0.01, max_freq_ratio: float = 0.1, use_welch: bool = True, window: str = 'hann', nperseg: int | None = None, use_multitaper: bool = False, n_tapers: int = 3, use_optimization: str = 'auto')[source]

Bases: BaseEstimator

Unified Periodogram-based Hurst parameter estimator.

This estimator computes the power spectral density (PSD) of the time series and fits a power law to the low-frequency portion to estimate the Hurst parameter. The relationship is: PSD(f) ~ f^(-beta) where beta = 2H - 1.

Backends: - ‘jax’: GPU-accelerated FFT (Welch/Periodogram). - ‘numpy’: SciPy-based (Periodogram/Welch/Multitaper). - ‘numba’: Defaults to NumPy backend.

__init__(min_freq_ratio: float = 0.01, max_freq_ratio: float = 0.1, use_welch: bool = True, window: str = 'hann', nperseg: int | None = None, use_multitaper: bool = False, n_tapers: int = 3, use_optimization: str = 'auto')[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter using periodogram analysis.

__init__(min_freq_ratio: float = 0.01, max_freq_ratio: float = 0.1, use_welch: bool = True, window: str = 'hann', nperseg: int | None = None, use_multitaper: bool = False, n_tapers: int = 3, use_optimization: str = 'auto')[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter using periodogram analysis.

_fit_power_law(freqs, psd, backend_name)[source]

Fit power law to the PSD.

get_optimization_info() → Dict[str, Any][source]
plot_scaling(save_path: str | None = None) → None[source]

Plot the scaling relationship and PSD.

Whittle Estimator

class lrdbenchmark.analysis.spectral.whittle_estimator.WhittleEstimator(min_freq_ratio: float = 0.01, max_freq_ratio: float = 0.5, use_optimization: str = 'auto', **kwargs)[source]

Bases: BaseEstimator

Whittle estimator for Hurst parameter estimation.

Minimizes the Whittle likelihood function based on the fGn spectral density. Supports JAX acceleration for likelihood evaluation.

__init__(min_freq_ratio: float = 0.01, max_freq_ratio: float = 0.5, use_optimization: str = 'auto', **kwargs)[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter using Whittle method.

__init__(min_freq_ratio: float = 0.01, max_freq_ratio: float = 0.5, use_optimization: str = 'auto', **kwargs)[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter using Whittle method.

_estimate_jax(data: ndarray) → Dict[str, Any][source]

JAX-accelerated Whittle estimation.

get_optimization_info() → Dict[str, Any][source]

Wavelet Estimators

Continuous Wavelet Transform (CWT)

class lrdbenchmark.analysis.wavelet.cwt_estimator.CWTEstimator(wavelet: str = 'morl', scales: ndarray | None = None, confidence: float = 0.95, use_optimization: str = 'auto', robust: bool = False, scale_range: Tuple[float, float] | None = None, trim_ends: int = 0)[source]

Bases: BaseEstimator

Unified Continuous Wavelet Transform (CWT) Estimator for Long-Range Dependence Analysis.

This estimator uses continuous wavelet transforms to analyze the scaling behavior of time series data and estimate the Hurst parameter for fractional processes.

Features: - Automatic optimization framework selection (JAX, Numba, NumPy) - GPU acceleration with JAX when available - JIT compilation with Numba for CPU optimization - Graceful fallbacks when optimization frameworks fail

Parameters:

  • wavelet (str, optional (default='cmor1.5-1.0')) – Wavelet type for continuous transform

  • scales (np.ndarray, optional (default=None)) – Array of scales for analysis. If None, uses automatic scale selection

  • confidence (float, optional (default=0.95)) – Confidence level for confidence intervals

  • use_optimization (str, optional (default='auto')) – Optimization framework to use: ‘auto’, ‘jax’, ‘numba’, ‘numpy’

__init__(wavelet: str = 'morl', scales: ndarray | None = None, confidence: float = 0.95, use_optimization: str = 'auto', robust: bool = False, scale_range: Tuple[float, float] | None = None, trim_ends: int = 0)[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate the Hurst parameter using Continuous Wavelet Transform analysis with automatic optimization.

Parameters:

data (array-like) – Input time series data

Returns:

Dictionary containing estimation results including: - hurst_parameter: Estimated Hurst parameter - confidence_interval: Confidence interval for the estimate - r_squared: R-squared value of the fit - scales: Scales used in the analysis - wavelet_type: Wavelet type used - slope: Slope of the log-log regression - intercept: Intercept of the log-log regression - scale_powers: Power at each scale

Return type:

dict

__init__(wavelet: str = 'morl', scales: ndarray | None = None, confidence: float = 0.95, use_optimization: str = 'auto', robust: bool = False, scale_range: Tuple[float, float] | None = None, trim_ends: int = 0)[source]
_select_optimization_framework(use_optimization: str) → str[source]

Select the optimal optimization framework.

_validate_parameters() → None[source]

Validate estimator parameters.

estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate the Hurst parameter using Continuous Wavelet Transform analysis with automatic optimization.

Parameters:

data (array-like) – Input time series data

Returns:

Dictionary containing estimation results including: - hurst_parameter: Estimated Hurst parameter - confidence_interval: Confidence interval for the estimate - r_squared: R-squared value of the fit - scales: Scales used in the analysis - wavelet_type: Wavelet type used - slope: Slope of the log-log regression - intercept: Intercept of the log-log regression - scale_powers: Power at each scale

Return type:

dict

_estimate_numpy(data: ndarray) → Dict[str, Any][source]

NumPy implementation of CWT estimation.

_estimate_numba(data: ndarray) → Dict[str, Any][source]

Numba-optimized implementation of CWT estimation.

_estimate_jax(data: ndarray) → Dict[str, Any][source]

JAX-optimized implementation of Wavelet Log Variance estimation.

_regression_statistics(x: ndarray, y: ndarray, weights: ndarray, slope: float, intercept: float, XtWX: ndarray) → Tuple[float, float][source]

Compute weighted regression diagnostics for slope.

_get_confidence_interval(estimated_hurst: float, slope_se: float, n_points: int) → Tuple[float, float][source]

Calculate confidence interval for the Hurst parameter estimate.

get_optimization_info() → Dict[str, Any][source]

Get information about available optimizations and current selection.

Get the recommended optimization framework.

plot_analysis(data: ndarray, figsize: Tuple[int, int] = (15, 10), save_path: str | None = None) → None[source]

Plot the CWT analysis results.

get_method_recommendation(n: int) → Dict[str, Any][source]

Get method recommendation for a given data size.

Wavelet Variance

class lrdbenchmark.analysis.wavelet.variance_estimator.WaveletVarianceEstimator(wavelet='db4', scales=None, confidence=0.95, use_optimization='auto', robust=False, j_min=2, j_max=None)[source]

Bases: BaseEstimator

Unified Wavelet Variance Estimator. Identical logic to LogVarianceEstimator, kept for compatibility/completeness.

__init__(wavelet='db4', scales=None, confidence=0.95, use_optimization='auto', robust=False, j_min=2, j_max=None)[source]
estimate(data)[source]
__init__(wavelet='db4', scales=None, confidence=0.95, use_optimization='auto', robust=False, j_min=2, j_max=None)[source]
estimate(data)[source]
get_optimization_info()[source]

Wavelet Log-Variance

class lrdbenchmark.analysis.wavelet.log_variance_estimator.WaveletLogVarianceEstimator(wavelet: str = 'db4', scales: List[int] | None = None, confidence: float = 0.95, use_optimization: str = 'auto', robust: bool = False, j_min: int = 2, j_max: int | None = None)[source]

Bases: BaseEstimator

Unified Wavelet Log Variance Estimator.

Uses wavelet decomposition to estimate Hurst parameter via log-variance scaling. Supports JAX acceleration for DWT.

__init__(wavelet: str = 'db4', scales: List[int] | None = None, confidence: float = 0.95, use_optimization: str = 'auto', robust: bool = False, j_min: int = 2, j_max: int | None = None)[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter.

__init__(wavelet: str = 'db4', scales: List[int] | None = None, confidence: float = 0.95, use_optimization: str = 'auto', robust: bool = False, j_min: int = 2, j_max: int | None = None)[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate Hurst parameter.

get_optimization_info() → Dict[str, Any][source]
plot_analysis(figsize=(12, 8), save_path=None)[source]

Wavelet Whittle

class lrdbenchmark.analysis.wavelet.whittle_estimator.WaveletWhittleEstimator(wavelet: str = 'db4', scales: List[int] | None = None, confidence: float = 0.95, use_optimization: str = 'auto', bootstrap_samples: int = 64, bootstrap_block_size: int | None = None)[source]

Bases: BaseEstimator

Unified Wavelet Whittle Estimator for Long-Range Dependence Analysis.

This estimator combines wavelet decomposition with Whittle likelihood estimation to provide robust estimation of the Hurst parameter for fractional processes.

Features: - Automatic optimization framework selection (JAX, Numba, NumPy) - GPU acceleration with JAX when available - JIT compilation with Numba for CPU optimization - Graceful fallbacks when optimization frameworks fail

Parameters:

  • wavelet (str, optional (default='db4')) – Wavelet type to use for decomposition

  • scales (List[int], optional (default=None)) – List of scales for wavelet analysis. If None, uses automatic scale selection

  • confidence (float, optional (default=0.95)) – Confidence level for confidence intervals

  • use_optimization (str, optional (default='auto')) – Optimization framework to use: ‘auto’, ‘jax’, ‘numba’, ‘numpy’

__init__(wavelet: str = 'db4', scales: List[int] | None = None, confidence: float = 0.95, use_optimization: str = 'auto', bootstrap_samples: int = 64, bootstrap_block_size: int | None = None)[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate the Hurst parameter using wavelet Whittle analysis with automatic optimization.

Parameters:

data (array-like) – Input time series data

Returns:

Dictionary containing estimation results including: - hurst_parameter: Estimated Hurst parameter - confidence_interval: Confidence interval for the estimate - whittle_likelihood: Whittle likelihood value - scales: Scales used in the analysis - wavelet_type: Wavelet type used - optimization_success: Whether optimization succeeded

Return type:

dict

__init__(wavelet: str = 'db4', scales: List[int] | None = None, confidence: float = 0.95, use_optimization: str = 'auto', bootstrap_samples: int = 64, bootstrap_block_size: int | None = None)[source]
_select_optimization_framework(use_optimization: str) → str[source]

Select the optimal optimization framework.

_validate_parameters() → None[source]

Validate estimator parameters.

estimate(data: ndarray | list) → Dict[str, Any][source]

Estimate the Hurst parameter using wavelet Whittle analysis with automatic optimization.

Parameters:

data (array-like) – Input time series data

Returns:

Dictionary containing estimation results including: - hurst_parameter: Estimated Hurst parameter - confidence_interval: Confidence interval for the estimate - whittle_likelihood: Whittle likelihood value - scales: Scales used in the analysis - wavelet_type: Wavelet type used - optimization_success: Whether optimization succeeded

Return type:

dict

_estimate_numpy(data: ndarray, compute_ci: bool = True) → Dict[str, Any][source]

NumPy implementation of Wavelet Whittle estimation.

_estimate_numba(data: ndarray) → Dict[str, Any][source]

Numba-optimized implementation of Wavelet Whittle estimation.

_bootstrap_confidence_interval(data: ndarray, point_estimate: float) → Tuple[float, float][source]

Approximate confidence interval using circular block bootstrap.

_circular_block_resample(data: ndarray, block_size: int) → ndarray[source]

Generate circular block bootstrap resample.

_estimate_jax(data: ndarray) → Dict[str, Any][source]

JAX-optimized implementation of Wavelet Whittle estimation.

_theoretical_spectrum_fgn(frequencies: ndarray, H: float, sigma: float = 1.0) → ndarray[source]

Calculate theoretical spectrum for fractional Gaussian noise.

_get_confidence_interval(estimated_hurst: float, whittle_likelihood: float, js: List[int]) → Tuple[float, float][source]

Calculate confidence interval for the Hurst parameter estimate.

get_optimization_info() → Dict[str, Any][source]

Get information about available optimizations and current selection.

Get the recommended optimization framework.

plot_analysis(figsize: Tuple[int, int] = (12, 8), save_path: str | None = None) → None[source]

Plot the wavelet Whittle analysis results.

get_method_recommendation(n: int) → Dict[str, Any][source]

Get method recommendation for a given data size.

Multifractal Estimators

Multifractal Detrended Fluctuation Analysis (MFDFA)

class lrdbenchmark.analysis.multifractal.mfdfa_estimator.MFDFAEstimator(q_orders: List[float] | ndarray | None = None, scales: List[int] | None = None, order: int = 1, use_optimization: str = 'auto', **kwargs)[source]

Bases: BaseEstimator

Multifractal Detrended Fluctuation Analysis (MFDFA) Estimator.

Estimates the generalized Hurst exponent H(q) and the multifractal spectrum D(h). Supports JAX acceleration.

__init__(q_orders: List[float] | ndarray | None = None, scales: List[int] | None = None, order: int = 1, use_optimization: str = 'auto', **kwargs)[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Perform MFDFA.

__init__(q_orders: List[float] | ndarray | None = None, scales: List[int] | None = None, order: int = 1, use_optimization: str = 'auto', **kwargs)[source]
estimate(data: ndarray | list) → Dict[str, Any][source]

Perform MFDFA.

get_optimization_info()[source]

Multifractal Wavelet Leaders

class lrdbenchmark.analysis.multifractal.wavelet_leaders_estimator.MultifractalWaveletLeadersEstimator(wavelet: str = 'db3', scales: List[int] | None = None, min_scale: int = 2, max_scale: int = 32, num_scales: int = 10, q_values: List[float] | None = None, use_optimization: str = 'numpy')[source]

Bases: BaseEstimator

Unified Multifractal Wavelet Leaders Estimator.

__init__(wavelet: str = 'db3', scales: List[int] | None = None, min_scale: int = 2, max_scale: int = 32, num_scales: int = 10, q_values: List[float] | None = None, use_optimization: str = 'numpy')[source]
estimate(data)[source]
__init__(wavelet: str = 'db3', scales: List[int] | None = None, min_scale: int = 2, max_scale: int = 32, num_scales: int = 10, q_values: List[float] | None = None, use_optimization: str = 'numpy')[source]
estimate(data)[source]
get_optimization_info()[source]

Machine Learning Estimators

For detailed documentation of machine learning estimators, see Machine Learning Estimators.

The following ML estimators are available:

  • Random Forest: Ensemble tree-based estimation with feature importance

  • Support Vector Regression: SVM-based estimation with RBF kernel

  • Gradient Boosting: Boosted tree estimation with comprehensive feature engineering

See the dedicated Machine Learning Estimators page for complete API documentation, performance metrics, and usage examples.

Neural Network Estimators

For detailed documentation of neural network estimators, see Neural Network Factory API.

The following neural network architectures are available:

  • CNN: Convolutional Neural Networks for spatial pattern recognition

  • LSTM: Long Short-Term Memory networks for temporal sequences

  • GRU: Gated Recurrent Units for efficient temporal modeling

  • Transformer: Attention-based architectures for complex patterns

See the dedicated Neural Network Factory API page for complete API documentation, architecture details, and usage examples.

High Performance Estimators

Note

The unified estimators automatically select optimal computation frameworks (JAX, Numba, or NumPy) based on data characteristics and hardware availability. For advanced users, the high-performance modules in lrdbenchmark.analysis.high_performance.jax and lrdbenchmark.analysis.high_performance.numba provide direct access to optimized implementations (e.g., DFAEstimatorJAX, DFAEstimatorNumba), but the unified estimators are recommended for most use cases as they automatically select the best framework based on data characteristics and hardware availability.

Backend Modules (Strategy Pattern)

The following backend modules provide modular implementations for each estimator family, enabling JAX GPU acceleration or Numba JIT compilation:

Temporal Backends

  • lrdbenchmark.analysis.temporal.dfa_backends: NumPy, JAX, Numba implementations for DFA.

  • lrdbenchmark.analysis.temporal.rs_backends: NumPy, JAX, Numba implementations for R/S.

Spectral Backends

  • lrdbenchmark.analysis.spectral.spectral_backends: NumPy, JAX implementations for Periodogram, Welch, Whittle.

Wavelet Backends

  • lrdbenchmark.analysis.wavelet.wavelet_backends: NumPy, JAX implementations for DWT-based variance estimators.

Multifractal Backends

  • lrdbenchmark.analysis.multifractal.mfdfa_backends: NumPy, JAX implementations for MFDFA.

  • lrdbenchmark.analysis.multifractal.wavelet_leaders_backends: NumPy, JAX implementations for Wavelet Leaders.

Usage Example (Backend Selection)

from lrdbenchmark import DFAEstimator

# Auto-select best backend (default)
estimator = DFAEstimator(use_optimization='auto')

# Force JAX for GPU acceleration
estimator_jax = DFAEstimator(use_optimization='jax')

# Force NumPy for maximum compatibility
estimator_numpy = DFAEstimator(use_optimization='numpy')

Usage Examples

Basic Estimator Usage

from lrdbenchmark import DFAEstimator
from lrdbenchmark import GPHEstimator
from lrdbenchmark import FBMModel

# Generate test data with known Hurst parameter
model = FBMModel(H=0.7, sigma=1.0)
data = model.generate(1000, seed=42)

print(f"Generated FBM data with true H = 0.7")
print(f"Data length: {len(data)}")
print(f"Data mean: {data.mean():.3f}, std: {data.std():.3f}")

# Use DFA estimator
dfa = DFAEstimator()
H_dfa = dfa.estimate(data)
print(f"DFA H estimate: {H_dfa:.3f}")
print(f"DFA error: {abs(H_dfa - 0.7):.3f}")

# Use GPH estimator
gph = GPHEstimator()
H_gph = gph.estimate(data)
print(f"GPH H estimate: {H_gph:.3f}")
print(f"GPH error: {abs(H_gph - 0.7):.3f}")

# Compare estimates
print(f"\nEstimate comparison:")
print(f"True H: 0.700")
print(f"DFA:    {H_dfa:.3f} (error: {abs(H_dfa - 0.7):.3f})")
print(f"GPH:    {H_gph:.3f} (error: {abs(H_gph - 0.7):.3f})")

Multiple Estimators Comparison

from lrdbenchmark import DFAEstimator
from lrdbenchmark import RSEstimator
from lrdbenchmark import GPHEstimator
from lrdbenchmark import WaveletVarianceEstimator
from lrdbenchmark import HiguchiEstimator
from lrdbenchmark import FBMModel, FGNModel
import pandas as pd

# Define estimators to test
estimators = {
    'DFA': DFAEstimator(),
    'R/S': RSEstimator(),
    'GPH': GPHEstimator(),
    'Wavelet Variance': WaveletVarianceEstimator(),
    'Higuchi': HiguchiEstimator()
}

# Test on different data models
test_cases = {
    'FBM (H=0.7)': FBMModel(H=0.7, sigma=1.0),
    'FBM (H=0.3)': FBMModel(H=0.3, sigma=1.0),
    'FGN (H=0.8)': FGNModel(H=0.8, sigma=1.0)
}

# Store results
all_results = []

for case_name, model in test_cases.items():
    print(f"\n=== Testing {case_name} ===")
    data = model.generate(1000, seed=42)
    true_H = model.H

    case_results = {'Case': case_name, 'True_H': true_H}

    for name, estimator in estimators.items():
        try:
            H_est = estimator.estimate(data)
            error = abs(H_est - true_H)
            case_results[name] = H_est
            case_results[f'{name}_error'] = error
            print(f"  {name}: H = {H_est:.3f} (error: {error:.3f})")
        except Exception as e:
            print(f"  {name}: Error - {e}")
            case_results[name] = None
            case_results[f'{name}_error'] = None

    all_results.append(case_results)

# Create summary DataFrame
df = pd.DataFrame(all_results)
print(f"\n=== SUMMARY ===")
print(df.round(3))

# Calculate average errors
error_columns = [col for col in df.columns if col.endswith('_error')]
avg_errors = df[error_columns].mean()
print(f"\n=== AVERAGE ERRORS ===")
for col in error_columns:
    estimator = col.replace('_error', '')
    print(f"{estimator}: {avg_errors[col]:.3f}")

Machine Learning Estimators

from lrdbenchmark import RandomForestEstimator
from lrdbenchmark import GradientBoostingEstimator
from lrdbenchmark import CNNEstimator
from lrdbenchmark import FBMModel, FGNModel, ARFIMAModel
import numpy as np
from sklearn.model_selection import train_test_split

# Generate comprehensive training dataset
print("Generating training dataset...")
training_data = []
training_labels = []

# Create diverse training data
H_values = np.linspace(0.3, 0.9, 15)  # 15 different H values
models = {
    'FBM': FBMModel,
    'FGN': FGNModel,
    'ARFIMA': lambda H: ARFIMAModel(d=H-0.5, p=1, q=1)
}

for H in H_values:
    for model_name, model_class in models.items():
        if model_name == 'ARFIMA':
            model = model_class(H)
        else:
            model = model_class(H=H, sigma=1.0)

        # Generate multiple realizations
        for i in range(20):
            data = model.generate(1000, seed=int(H*1000 + i))
            training_data.append(data)
            training_labels.append(H)

print(f"Generated {len(training_data)} training samples")
print(f"H range: {min(training_labels):.1f} to {max(training_labels):.1f}")

# Split into training and validation sets
X_train, X_val, y_train, y_val = train_test_split(
    training_data, training_labels, test_size=0.2, random_state=42
)

# Train Random Forest estimator
print("\nTraining Random Forest estimator...")
rf_estimator = RandomForestEstimator(
    n_estimators=100,
    max_depth=10,
    random_state=42
)
rf_estimator.fit(X_train, y_train)

# Train Gradient Boosting estimator
print("Training Gradient Boosting estimator...")
gb_estimator = GradientBoostingEstimator(
    n_estimators=100,
    learning_rate=0.1,
    max_depth=5,
    random_state=42
)
gb_estimator.fit(X_train, y_train)

# Evaluate on validation set
print("\n=== Validation Results ===")
rf_val_pred = rf_estimator.estimate(X_val)
gb_val_pred = gb_estimator.estimate(X_val)

rf_mae = np.mean(np.abs(np.array(rf_val_pred) - np.array(y_val)))
gb_mae = np.mean(np.abs(np.array(gb_val_pred) - np.array(y_val)))

print(f"Random Forest MAE: {rf_mae:.3f}")
print(f"Gradient Boosting MAE: {gb_mae:.3f}")

# Test on new data
print("\n=== Test on New Data ===")
test_cases = [
    ('FBM (H=0.6)', FBMModel(H=0.6, sigma=1.0)),
    ('FGN (H=0.4)', FGNModel(H=0.4, sigma=1.0)),
    ('ARFIMA (H=0.75)', ARFIMAModel(d=0.25, p=1, q=1))
]

for test_name, test_model in test_cases:
    test_data = test_model.generate(1000, seed=999)

    H_rf = rf_estimator.estimate([test_data])[0]
    H_gb = gb_estimator.estimate([test_data])[0]

    if 'H=' in test_name:
        true_H = float(test_name.split('H=')[1].split(')')[0])
    else:
        true_H = 0.75  # For ARFIMA

    print(f"{test_name}:")
    print(f"  True H: {true_H:.3f}")
    print(f"  RF estimate: {H_rf:.3f} (error: {abs(H_rf - true_H):.3f})")
    print(f"  GB estimate: {H_gb:.3f} (error: {abs(H_gb - true_H):.3f})")

High Performance Estimators

from lrdbenchmark import DFAEstimator, RSEstimator, FBMModel
import numpy as np

# Generate large dataset
model = FBMModel(H=0.7, sigma=1.0)
data = model.generate(10000, seed=42)

# Unified estimators automatically use optimal backend (JAX/Numba/NumPy)
dfa = DFAEstimator(use_optimization='auto')  # Auto-selects best framework
result = dfa.estimate(data)
print(f"DFA H estimate: {result['hurst_parameter']:.3f}")

# R/S estimator with automatic optimization
rs = RSEstimator(use_optimization='auto')
result = rs.estimate(data)
print(f"R/S H estimate: {result['hurst_parameter']:.3f}")

Parameter Tuning

from lrdbenchmark import DFAEstimator
from lrdbenchmark import GPHEstimator
from lrdbenchmark import FBMModel

# Generate test data
model = FBMModel(H=0.7, sigma=1.0)
data = model.generate(1000, seed=42)

# DFA with custom parameters
dfa = DFAEstimator(
    min_scale=4,
    max_scale=100,
    num_scales=20,
    polynomial_order=2
)
H_dfa = dfa.estimate(data)

# GPH with custom parameters
gph = GPHEstimator(
    num_frequencies=50,
    min_frequency=0.01,
    max_frequency=0.5
)
H_gph = gph.estimate(data)

print(f"DFA (custom): H = {H_dfa:.3f}")
print(f"GPH (custom): H = {H_gph:.3f}")

# Note: All estimators are documented above in their respective sections # No duplicate documentation needed

Error Handling

from lrdbenchmark import DFAEstimator
from lrdbenchmark import GPHEstimator

# Test with insufficient data
short_data = [1, 2, 3, 4, 5]  # Too short for most estimators

dfa = DFAEstimator()
try:
    H_dfa = dfa.estimate(short_data)
    print(f"DFA H estimate: {H_dfa:.3f}")
except ValueError as e:
    print(f"DFA error: {e}")

gph = GPHEstimator()
try:
    H_gph = gph.estimate(short_data)
    print(f"GPH H estimate: {H_gph:.3f}")
except ValueError as e:
    print(f"GPH error: {e}")

Performance Comparison

import time
from lrdbenchmark import DFAEstimator, FBMModel

# Generate test data
model = FBMModel(H=0.7, sigma=1.0)
data = model.generate(5000, seed=42)

# DFA with automatic optimization (will use JAX/Numba if available)
dfa = DFAEstimator(use_optimization='auto')
start_time = time.time()
result = dfa.estimate(data)
dfa_time = time.time() - start_time

# DFA with explicit NumPy backend (no optimization)
dfa_numpy = DFAEstimator(use_optimization='numpy')
start_time = time.time()
result_numpy = dfa_numpy.estimate(data)
dfa_numpy_time = time.time() - start_time

print(f"Optimized DFA: H = {result['hurst_parameter']:.3f}, Time = {dfa_time:.4f}s")
print(f"NumPy DFA: H = {result_numpy['hurst_parameter']:.3f}, Time = {dfa_numpy_time:.4f}s")
if dfa_numpy_time > 0:
    print(f"Speedup: {dfa_numpy_time/dfa_time:.2f}x")

Best Practices

  1. Data Length: Use at least 1000 samples for reliable estimates

  2. Parameter Selection: Choose appropriate scale ranges for your data

  3. Multiple Estimators: Compare results from different estimator types

  4. Error Handling: Always handle potential estimation errors

  5. Performance: Use high-performance estimators for large datasets

  6. Validation: Test on synthetic data with known Hurst parameters

Note

Different estimators may give slightly different results due to their underlying assumptions and methodologies. It’s recommended to use multiple estimators and compare their results.

Warning

Some estimators require specific data characteristics (e.g., stationarity for spectral methods). Always check the estimator documentation for requirements and limitations.