Theoretical Foundations

This document provides the theoretical foundations for the statistical tests, validation techniques, and mathematical principles underlying LRDBench.

Long-Range Dependence (LRD)

Definition and Properties

Long-range dependence (LRD), also known as long memory or persistence, is a statistical property of time series where correlations between observations decay slowly with increasing time lag.

For a stationary time series \(\{X_t\}\) with autocorrelation function \(\rho(k)\), LRD is characterized by:

\[\rho(k) \sim c|k|^{2H-2} \text{ as } k \to \infty\]

where: - \(H\) is the Hurst parameter (\(0.5 < H < 1\)) - \(c\) is a positive constant - \(\sim\) denotes asymptotic equivalence

The Hurst parameter \(H\) quantifies the degree of long-range dependence: - \(H = 0.5\): No long-range dependence (white noise) - \(0.5 < H < 1\): Positive long-range dependence (persistence) - \(0 < H < 0.5\): Negative long-range dependence (anti-persistence)

Power Spectral Density

For LRD processes, the power spectral density \(S(f)\) exhibits a power-law behavior at low frequencies:

\[S(f) \sim c|f|^{1-2H} \text{ as } f \to 0\]

This relationship forms the basis for spectral-based estimators like the Geweke-Porter-Hudak (GPH) estimator.

Fractional Brownian Motion (FBM)

Mathematical Definition

Fractional Brownian Motion \(B_H(t)\) is a continuous-time Gaussian process with stationary increments, defined by:

\[B_H(t) = \frac{1}{\Gamma(H + \frac{1}{2})} \int_{-\infty}^t (t-s)^{H-\frac{1}{2}} dB(s)\]

where: - \(H\) is the Hurst parameter (\(0 < H < 1\)) - \(\Gamma(\cdot)\) is the gamma function - \(B(s)\) is standard Brownian motion

Properties

1. Self-similarity: \(B_H(at) \stackrel{d}{=} a^H B_H(t)\) for all \(a > 0\)

2. Stationary increments: \(B_H(t+s) - B_H(s) \stackrel{d}{=} B_H(t) - B_H(0)\)

3. Variance scaling: \(\text{Var}[B_H(t)] = t^{2H}\)

4. Covariance function: \(\text{Cov}[B_H(s), B_H(t)] = \frac{1}{2}(|s|^{2H} + |t|^{2H} - |s-t|^{2H})\)

Generation Methods

Method 1: Cholesky Decomposition The covariance matrix \(\Sigma\) is constructed and decomposed as \(\Sigma = LL^T\). The FBM is then generated as \(B_H = LZ\) where \(Z\) is a vector of independent standard normal random variables.

Method 2: Circulant Embedding The covariance function is embedded in a circulant matrix, which can be diagonalized using the FFT, enabling efficient generation.

Method 3: Davies-Harte Method Uses the spectral representation of FBM to generate samples via FFT.

Fractional Gaussian Noise (FGN)

Definition

Fractional Gaussian Noise is the increment process of FBM:

\[X_t = B_H(t+1) - B_H(t)\]

Properties

1. Stationarity: FGN is a stationary process

2. Autocorrelation: \(\rho(k) = \frac{1}{2}(|k+1|^{2H} - 2|k|^{2H} + |k-1|^{2H})\)

3. Long-range dependence: For \(H > 0.5\), \(\rho(k) \sim H(2H-1)k^{2H-2}\) as \(k \to \infty\)

ARFIMA Models

Definition

An ARFIMA(p,d,q) process \(\{X_t\}\) satisfies:

\[\Phi(B)(1-B)^d X_t = \Theta(B)\epsilon_t\]

where: - \(\Phi(B) = 1 - \phi_1 B - \cdots - \phi_p B^p\) (AR polynomial) - \(\Theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q\) (MA polynomial) - \((1-B)^d\) is the fractional differencing operator - \(\epsilon_t\) is white noise

Fractional Differencing

The fractional differencing operator \((1-B)^d\) is defined by the binomial expansion:

\[(1-B)^d = \sum_{k=0}^{\infty} \binom{d}{k} (-B)^k\]

where \(\binom{d}{k} = \frac{d(d-1)\cdots(d-k+1)}{k!}\)

For \(|d| < 0.5\), the process is stationary and invertible. The relationship between \(d\) and the Hurst parameter is:

\[H = d + 0.5\]

Multifractal Random Walk (MRW)

Definition

The Multifractal Random Walk is defined as:

\[X(t) = \int_0^t e^{\omega(s)} dB(s)\]

where: - \(B(s)\) is standard Brownian motion - \(\omega(s)\) is a stationary Gaussian process with covariance:

\[\text{Cov}[\omega(s), \omega(t)] = \lambda^2 \log_+ \frac{T}{|s-t| + \ell}\]

Parameters: - \(\lambda^2\): Intermittency parameter - \(T\): Integral time scale - \(\ell\): Small-scale cutoff

Properties

1. Multifractality: The process exhibits different scaling exponents at different time scales

2. Log-normal multipliers: The process can be constructed using log-normal multipliers

3. Scaling: \(\langle |X(t+\tau) - X(t)|^q \rangle \sim \tau^{\zeta(q)}\) where \(\zeta(q)\) is the multifractal spectrum

Statistical Estimators

Temporal Domain Estimators

Detrended Fluctuation Analysis (DFA)

Algorithm:

1. Integration: \(y(i) = \sum_{k=1}^i (x_k - \bar{x})\)

2. Segmentation: Divide into \(N_s = N/s\) non-overlapping segments

3. Detrending: Fit polynomial \(y_n(i)\) to each segment

4. Fluctuation: Calculate \(F^2(s) = \frac{1}{s} \sum_{i=1}^s [y(i) - y_n(i)]^2\)

5. Scaling: \(F(s) \sim s^H\)

Theoretical Foundation: DFA measures the scaling of fluctuations around local trends, making it robust to non-stationarities.

Mathematical Formulation: For a time series of length \(N\), the DFA fluctuation function is:

\[F(s) = \sqrt{\frac{1}{N_s} \sum_{v=1}^{N_s} F^2(v,s)}\]

where \(F^2(v,s)\) is the mean squared fluctuation in segment \(v\) of size \(s\).

R/S Analysis

Algorithm:

1. Segmentation: Divide data into segments of length \(k\)

2. Rescaled Range: For each segment, calculate: - \(R = \max_{1 \leq i \leq k} S_i - \min_{1 \leq i \leq k} S_i\) (range) - \(S = \sqrt{\frac{1}{k} \sum_{i=1}^k (x_i - \bar{x})^2}\) (standard deviation)

3. Scaling: \(R/S \sim k^H\)

Theoretical Foundation: R/S analysis measures the scaling of the range of partial sums, normalized by the standard deviation.

Mathematical Formulation: For a segment of length \(k\), the rescaled range is:

\[R/S = \frac{\max_{1 \leq i \leq k} \sum_{j=1}^i (x_j - \bar{x}) - \min_{1 \leq i \leq k} \sum_{j=1}^i (x_j - \bar{x})}{\sqrt{\frac{1}{k} \sum_{i=1}^k (x_i - \bar{x})^2}}\]

Higuchi Method

Detrended Moving Average (DMA)

Algorithm:

1. Moving Average: Compute local moving averages at window length \(s\)

2. Residuals: :math:` ,epsilon(i,s) = x(i) - m_s(i) `

3. Fluctuation: :math:` F^2(s) = frac{1}{N-s+1}sum_{i=1}^{N-s+1} epsilon(i,s)^2 `

4. Scaling: \(F(s) \sim s^H\)

Theoretical Foundation: DMA measures scaling of residual variance after local averaging; for long‑memory processes, variance scales with window length.

Mathematical Formulation: For a centred series \(x(i)\), the moving average \(m_s(i)\) over window \(s\) defines residuals whose variance follows a power law in \(s\).

Algorithm:

1. Subseries Construction: For each \(k\), construct \(k\) subseries

2. Length Calculation: Calculate the length \(L_m(k)\) of each subseries

3. Average Length: \(L(k) = \frac{1}{k} \sum_{m=1}^k L_m(k)\)

4. Scaling: \(L(k) \sim k^{-D}\) where \(D = 2 - H\)

Theoretical Foundation: The Higuchi method estimates the fractal dimension by measuring how the length of the time series changes with different sampling intervals.

Mathematical Formulation: For a time series \(\{x_i\}\) and lag \(k\), the length is:

\[L_m(k) = \frac{1}{k} \left[ \frac{N-1}{k^2} \sum_{i=1}^{[(N-m)/k]} |x_{m+ik} - x_{m+(i-1)k}| \right]\]

Spectral Domain Estimators

Geweke-Porter-Hudak (GPH) Estimator

Algorithm:

1. Periodogram: Calculate \(I(f_j) = \frac{1}{2\pi N} |\sum_{t=1}^N x_t e^{-i2\pi f_j t}|^2\)

2. Log-Regression: Fit \(\log I(f_j) = c - d \log(4\sin^2(\pi f_j)) + \epsilon_j\)

3. Estimation: \(H = d + 0.5\)

Theoretical Foundation: The GPH estimator is based on the spectral representation of ARFIMA processes, where the log-periodogram follows a linear relationship with the log-frequency.

Mathematical Formulation: For frequencies \(f_j = j/N\), the regression model is:

\[\log I(f_j) = c - d \log(4\sin^2(\pi f_j)) + \epsilon_j\]

where \(d\) is the fractional differencing parameter and \(H = d + 0.5\).

Whittle Estimator

Algorithm:

1. Spectral Density: Assume parametric form \(S(f; \theta)\)

2. Whittle Likelihood: \(L(\theta) = \sum_{j=1}^{N/2} \left[ \log S(f_j; \theta) + \frac{I(f_j)}{S(f_j; \theta)} \right]\)

3. Optimization: Maximize \(L(\theta)\) to estimate parameters

Theoretical Foundation: The Whittle estimator maximizes an approximation to the likelihood function in the frequency domain, making it asymptotically efficient.

Mathematical Formulation: The Whittle likelihood function is:

\[L(\theta) = \sum_{j=1}^{N/2} \left[ \log S(f_j; \theta) + \frac{I(f_j)}{S(f_j; \theta)} \right]\]

where \(S(f; \theta)\) is the theoretical spectral density and \(I(f_j)\) is the periodogram.

Wavelet Domain Estimators

Wavelet Variance

Algorithm:

1. Wavelet Decomposition: Apply discrete wavelet transform

2. Variance Calculation: \(\sigma^2_j = \frac{1}{n_j} \sum_{k=1}^{n_j} d_{j,k}^2\)

3. Scaling: \(\sigma^2_j \sim 2^{j(2H-1)}\)

Theoretical Foundation: Wavelet variance measures the energy at different scales, providing a robust estimate of the scaling exponent.

Mathematical Formulation: For wavelet coefficients \(d_{j,k}\) at scale \(j\), the variance is:

\[\sigma^2_j = \frac{1}{n_j} \sum_{k=1}^{n_j} d_{j,k}^2\]

where \(n_j\) is the number of coefficients at scale \(j\).

Continuous Wavelet Transform (CWT)

Algorithm:

1. CWT Calculation: \(W_x(a,b) = \frac{1}{\sqrt{a}} \int_{-\infty}^{\infty} x(t) \psi^*\left(\frac{t-b}{a}\right) dt\)

2. Wavelet Spectrum: \(S_x(a) = \int_{-\infty}^{\infty} |W_x(a,b)|^2 db\)

3. Scaling: \(S_x(a) \sim a^{2H+1}\)

Theoretical Foundation: CWT provides a time-scale representation that preserves both temporal and frequency information.

Mathematical Formulation: The continuous wavelet transform is:

\[W_x(a,b) = \frac{1}{\sqrt{a}} \int_{-\infty}^{\infty} x(t) \psi^*\left(\frac{t-b}{a}\right) dt\]

where \(\psi(t)\) is the mother wavelet and \(a, b\) are scale and translation parameters.

Multifractal Estimators

Generalised Hurst Exponent (GHE)

Algorithm:

1. Increments: \(\Delta X_\tau(t) = X(t+\tau) - X(t)\)

2. q‑th Moments: \(K_q(\tau) = \langle |\Delta X_\tau|^q \rangle\)

3. Scaling: \(K_q(\tau) \sim \tau^{qH(q)}\); estimate \(H(q)\) by linear regression of \(\log K_q\) vs \(\log \tau\)

Notes: For monofractal processes, \(H(q)\) is constant and equals the standard Hurst parameter; variability in \(H(q)\) indicates multifractality.

Multifractal Detrended Fluctuation Analysis (MFDFA)

Algorithm:

1. Profile: \(Y(i) = \sum_{k=1}^i (x_k - \bar{x})\)

2. Segmentation: Divide into \(N_s = N/s\) segments

3. Detrending: Fit polynomial to each segment

4. Fluctuation: \(F_q(s) = \left[ \frac{1}{N_s} \sum_{v=1}^{N_s} F^2(v,s)^{q/2} \right]^{1/q}\)

5. Scaling: \(F_q(s) \sim s^{h(q)}\)

Theoretical Foundation: MFDFA extends DFA to capture multifractal scaling by considering different moments of the fluctuation function.

Mathematical Formulation: The qth order fluctuation function is:

\[F_q(s) = \left[ \frac{1}{N_s} \sum_{v=1}^{N_s} F^2(v,s)^{q/2} \right]^{1/q}\]

The multifractal spectrum \(f(\alpha)\) is obtained via Legendre transform:

\[\alpha = h(q) + qh'(q), \quad f(\alpha) = q[\alpha - h(q)] + 1\]

Theoretical Performance Considerations

Computational Complexity: Big-O notation for time and space complexity Convergence Rate: Rate at which estimator approaches true value Asymptotic Efficiency: Ratio of estimator variance to Cramér-Rao lower bound

For detailed performance metrics, statistical tests, and validation procedures, see the Validation Techniques and Statistical Tests section.

Practical Examples

Monte Carlo Simulation Example

For Monte Carlo validation examples and implementation details, see the dedicated Validation Guide.

Power Spectral Density Analysis

import numpy as np
from scipy import signal
from lrdbenchmark import FBMModel, FGNModel
import matplotlib.pyplot as plt

def power_spectral_density_example():
    """Demonstrate power spectral density analysis for LRD processes."""

    # Generate data with different Hurst parameters
    models = {
        'FBM (H=0.3)': FBMModel(H=0.3, sigma=1.0),
        'FBM (H=0.5)': FBMModel(H=0.5, sigma=1.0),
        'FBM (H=0.7)': FBMModel(H=0.7, sigma=1.0),
        'FBM (H=0.9)': FBMModel(H=0.9, sigma=1.0)
    }

    plt.figure(figsize=(12, 8))

    for model_name, model in models.items():
        # Generate data
        data = model.generate(2000, seed=42)

        # Compute power spectral density
        freqs, psd = signal.welch(data, fs=1.0, nperseg=256)

        # Plot PSD
        plt.loglog(freqs, psd, label=model_name, linewidth=2)

    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power Spectral Density')
    plt.title('Power Spectral Density of FBM Processes')
    plt.legend()
    plt.grid(True)
    plt.show()

    # Theoretical PSD for comparison
    plt.figure(figsize=(10, 6))
    freqs_theoretical = np.logspace(-3, 0, 100)

    for H in [0.3, 0.5, 0.7, 0.9]:
        # Theoretical PSD: S(f) ∝ f^(-2H+1)
        psd_theoretical = freqs_theoretical**(-2*H + 1)
        plt.loglog(freqs_theoretical, psd_theoretical,
                  label=f'Theoretical (H={H})', linestyle='--')

    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power Spectral Density')
    plt.title('Theoretical Power Spectral Density')
    plt.legend()
    plt.grid(True)
    plt.show()

# Run the example
if __name__ == "__main__":
    power_spectral_density_example()
    print("Power spectral density analysis completed!")

Autocorrelation Function Analysis

import numpy as np
from lrdbenchmark import FBMModel
import matplotlib.pyplot as plt

def autocorrelation_analysis_example():
    """Demonstrate autocorrelation function analysis for LRD processes."""

    # Generate FBM data with different H values
    H_values = [0.3, 0.5, 0.7, 0.9]
    max_lag = 100

    plt.figure(figsize=(12, 8))

    for H in H_values:
        # Generate data
        model = FBMModel(H=H, sigma=1.0)
        data = model.generate(2000, seed=42)

        # Compute autocorrelation function
        acf = np.correlate(data, data, mode='full')
        acf = acf[len(data)-1:len(data)-1+max_lag] / acf[len(data)-1]

        # Plot ACF
        lags = np.arange(max_lag)
        plt.plot(lags, acf, label=f'FBM (H={H})', linewidth=2)

    plt.xlabel('Lag')
    plt.ylabel('Autocorrelation')
    plt.title('Autocorrelation Function of FBM Processes')
    plt.legend()
    plt.grid(True)
    plt.show()

    # Theoretical ACF comparison
    plt.figure(figsize=(10, 6))
    lags_theoretical = np.arange(1, max_lag+1)

    for H in H_values:
        # Theoretical ACF: ρ(k) ∝ k^(2H-2)
        acf_theoretical = lags_theoretical**(2*H - 2)
        plt.loglog(lags_theoretical, acf_theoretical,
                  label=f'Theoretical (H={H})', linestyle='--')

    plt.xlabel('Lag')
    plt.ylabel('Autocorrelation')
    plt.title('Theoretical Autocorrelation Function')
    plt.legend()
    plt.grid(True)
    plt.show()

# Run the example
if __name__ == "__main__":
    autocorrelation_analysis_example()
    print("Autocorrelation analysis completed!")

Theoretical References

1. Beran, J. (1994). Statistics for Long-Memory Processes. Chapman & Hall.

2. Mandelbrot, B. B., & Van Ness, J. W. (1968). Fractional Brownian motions, fractional noises and applications. SIAM Review, 10(4), 422-437.

3. Peng, C. K., et al. (1994). Mosaic organization of DNA nucleotides. Physical Review E, 49(2), 1685.

4. Geweke, J., & Porter-Hudak, S. (1983). The estimation and application of long memory time series models. Journal of Time Series Analysis, 4(4), 221-238.

5. Kantelhardt, J. W., et al. (2002). Multifractal detrended fluctuation analysis of nonstationary time series. Physica A, 316(1-4), 87-114.

6. Abry, P., & Veitch, D. (1998). Wavelet analysis of long-range-dependent traffic. IEEE Transactions on Information Theory, 44(1), 2-15.