# =============================================================================
#  Project   : anomaly-atlas
#  File      : src/anomaly_atlas/stats/reversion.py
#  Purpose   : Mean-reversion tests: variance ratios, AC1, AR half-life
#  Author    : Simon-Pierre Boucher
#  Contact   : contact@spboucher.ai
#  Data src  : hfmarketdata.io (sole data source)
#  Created   : 2026-08-12
#  Modified  : 2026-08-12
#  Platform  : macOS / Apple Silicon (arm64)
#  License   : All rights reserved (research code)
# =============================================================================
"""Mean-reversion statistics on return series.

Validated on synthetic ground truth before touching real data
(benchmarks/synthetic/test_synthetic_gate.py — charter §8.1).
"""

from __future__ import annotations

import numpy as np


def ac1(returns: np.ndarray) -> float:
    """Lag-1 autocorrelation of a return series."""
    r = np.asarray(returns, dtype=float)
    if len(r) < 3:
        return float("nan")
    a, b = r[:-1], r[1:]
    sa, sb = a.std(), b.std()
    if sa == 0.0 or sb == 0.0:
        return float("nan")
    return float(((a - a.mean()) * (b - b.mean())).mean() / (sa * sb))


def autocov1(returns: np.ndarray) -> float:
    """Lag-1 autocovariance (input to the Roll spread estimator)."""
    r = np.asarray(returns, dtype=float)
    if len(r) < 3:
        return float("nan")
    a, b = r[:-1], r[1:]
    return float(((a - a.mean()) * (b - b.mean())).mean())


def variance_ratio(returns: np.ndarray, q: int) -> float:
    """Lo-MacKinlay variance ratio VR(q) with overlapping q-period sums.

    Ground truth: VR = 1 for a random walk, < 1 under mean reversion,
    > 1 under momentum. Unbiased variance estimators, demeaned.
    """
    r = np.asarray(returns, dtype=float)
    n = len(r)
    if n < q + 2 or q < 2:
        return float("nan")
    mu = r.mean()
    var1 = ((r - mu) ** 2).sum() / (n - 1)
    rq = np.convolve(r, np.ones(q), mode="valid")  # overlapping q-sums
    # Lo-MacKinlay bias-corrected PER-PERIOD variance of q-sums: the factor q
    # lives inside m, so the ratio below is varq/var1 (NOT varq/(q*var1)).
    m = q * (n - q + 1) * (1 - q / n)
    varq = ((rq - q * mu) ** 2).sum() / m
    if var1 == 0.0:
        return float("nan")
    return float(varq / var1)


def half_life(log_prices: np.ndarray) -> float:
    """Mean-reversion half-life from an AR(1) fit: dp_t = a + b*p_{t-1} + e.

    Returns ln(2)/-ln(1+b) in bars for b in (-1, 0); +inf if b >= 0
    (no reversion). Matches the OU generator's ln(2)/-ln(1-kappa).
    """
    p = np.asarray(log_prices, dtype=float)
    if len(p) < 10:
        return float("nan")
    x, y = p[:-1], np.diff(p)
    vx = x.var()
    if vx == 0.0:
        return float("nan")
    b = ((x - x.mean()) * (y - y.mean())).mean() / vx
    if b >= 0.0 or b <= -1.0:
        return float("inf")
    return float(np.log(2.0) / -np.log1p(b))