QuanTAlib/python/tests/test_pandas_ta_parity.py

"""test_pandas_ta_parity.py — Validate quantalib vs pandas-ta using the same
methodology as our C# ValidationHelper:

  1. Generate a LONG GBM series (5000 bars, seeded) so recursive indicators converge.
  2. Compare only the LAST 100 values (DefaultVerificationCount = 100).
  3. Skip lookback/warmup bars before comparison window.
  4. Tolerance: 1e-7 default, looser for known algorithmic differences (SPEC §9.3).

This mirrors lib/feeds/gbm/ValidationHelper.cs exactly.
"""
from __future__ import annotations

import numpy as np
import pandas as pd
import pandas_ta as ta
import pytest

from quantalib.indicators import (
    sma, ema, dema, tema, wma, hma, trima, alma, rsi, roc, mom,
    stddev, variance, zscore, bbands,
)

# ---------------------------------------------------------------------------
# Constants — match C# ValidationHelper
# ---------------------------------------------------------------------------
SEED = 42
N = 10000  # long series for convergence (C# uses 500-5000)
VERIFY_COUNT = 100  # DefaultVerificationCount in C#
DEFAULT_TOL = 1e-9  # ValidationHelper.DefaultTolerance


# ---------------------------------------------------------------------------
# GBM data generation — match C# GBM feed (Geometric Brownian Motion)
# ---------------------------------------------------------------------------
def _generate_gbm(n: int, seed: int = SEED, start_price: float = 100.0,
                   mu: float = 0.05, sigma: float = 0.2,
                   dt: float = 1 / 252) -> np.ndarray:
    """Generate GBM close prices matching C# GBM feed logic.

    S(t+1) = S(t) * exp((mu - sigma^2/2)*dt + sigma*sqrt(dt)*Z)
    """
    rng = np.random.default_rng(seed)
    z = rng.standard_normal(n - 1)
    drift = (mu - 0.5 * sigma ** 2) * dt
    diffusion = sigma * np.sqrt(dt) * z
    log_returns = drift + diffusion
    prices = np.empty(n, dtype=np.float64)
    prices[0] = start_price
    np.cumsum(log_returns, out=prices[1:])
    prices[1:] += np.log(start_price)
    np.exp(prices[1:], out=prices[1:])
    prices[0] = start_price
    return prices


# Module-level test data (generated once, reused across all tests)
CLOSE = _generate_gbm(N)
SERIES = pd.Series(CLOSE, name="close")


# ---------------------------------------------------------------------------
# Comparison helper — mirrors ValidationHelper.VerifyData logic
# ---------------------------------------------------------------------------
def _verify_last_n(
    qtl_arr: np.ndarray,
    pta_result: pd.Series | np.ndarray,
    *,
    verify_count: int = VERIFY_COUNT,
    tolerance: float = DEFAULT_TOL,
    label: str = "",
) -> None:
    """Compare only the last `verify_count` values where both are finite.

    This matches the C# pattern:
        int start = Math.Max(0, count - skip);
        for (int i = start; i < count; i++) { ... compare ... }
    """
    pta = pta_result.to_numpy() if isinstance(pta_result, pd.Series) else pta_result
    assert len(qtl_arr) == len(pta), (
        f"{label}: length mismatch qtl={len(qtl_arr)} vs pta={len(pta)}"
    )

    count = len(qtl_arr)
    start = max(0, count - verify_count)

    q_tail = qtl_arr[start:]
    p_tail = pta[start:]

    # Both must be finite in the tail (converged region)
    finite = np.isfinite(q_tail) & np.isfinite(p_tail)
    assert np.sum(finite) > 0, f"{label}: no finite values in last {verify_count}"

    q_vals = q_tail[finite]
    p_vals = p_tail[finite]

    max_diff = float(np.max(np.abs(q_vals - p_vals)))
    assert max_diff <= tolerance, (
        f"{label}: max_diff={max_diff:.2e} exceeds tolerance={tolerance:.0e} "
        f"(compared {len(q_vals)} values in last {verify_count})"
    )


# ===========================================================================
# FIR Trend indicators — exact match expected
# ===========================================================================

class TestTrendFIR:
    """FIR indicators: SMA, WMA, HMA — deterministic convolution, tight tolerance."""

    def test_sma(self) -> None:
        for length in (10, 20, 50):
            qtl = sma(CLOSE, length=length)
            pta = ta.sma(SERIES, length=length)
            _verify_last_n(qtl, pta, tolerance=1e-9,
                           label=f"SMA({length})")

    def test_wma(self) -> None:
        for length in (10, 14, 30):
            qtl = wma(CLOSE, length=length)
            pta = ta.wma(SERIES, length=length)
            _verify_last_n(qtl, pta, tolerance=1e-9,
                           label=f"WMA({length})")

    def test_hma(self) -> None:
        for length in (9, 14, 20):
            qtl = hma(CLOSE, length=length)
            pta = ta.hma(SERIES, length=length)
            _verify_last_n(qtl, pta, tolerance=1e-8,
                           label=f"HMA({length})")

    @pytest.mark.xfail(reason="TRIMA kernel differs: quantalib uses symmetric "
                               "triangular convolution, pandas-ta delegates to TA-Lib "
                               "which uses cascaded SMA (SPEC §9.3 known delta)")
    def test_trima(self) -> None:
        qtl = trima(CLOSE, length=14)
        pta = ta.trima(SERIES, length=14)
        _verify_last_n(qtl, pta, tolerance=1e-9, label="TRIMA(14)")

    @pytest.mark.xfail(reason="ALMA sigma/offset defaults differ between "
                               "quantalib and pandas-ta (SPEC §9.3 known delta)")
    def test_alma(self) -> None:
        qtl = alma(CLOSE, length=14)
        pta = ta.alma(SERIES, length=14)
        _verify_last_n(qtl, pta, tolerance=1e-6, label="ALMA(14)")


# ===========================================================================
# IIR Trend indicators — recursive, compare converged tail only
# ===========================================================================

class TestTrendIIR:
    """IIR indicators: EMA, DEMA, TEMA — recursive convergence.

    With 5000 bars the warmup difference is buried in the past.
    The last 100 bars should match tightly.
    """

    def test_ema(self) -> None:
        for length in (10, 20, 50):
            qtl = ema(CLOSE, length=length)
            pta = ta.ema(SERIES, length=length)
            _verify_last_n(qtl, pta, tolerance=1e-7,
                           label=f"EMA({length})")

    def test_dema(self) -> None:
        for length in (10, 20, 50):
            qtl = dema(CLOSE, length=length)
            pta = ta.dema(SERIES, length=length)
            _verify_last_n(qtl, pta, tolerance=1e-7,
                           label=f"DEMA({length})")

    def test_tema(self) -> None:
        for length in (10, 14, 30):
            qtl = tema(CLOSE, length=length)
            pta = ta.tema(SERIES, length=length)
            _verify_last_n(qtl, pta, tolerance=1e-7,
                           label=f"TEMA({length})")


# ===========================================================================
# Momentum indicators
# ===========================================================================

class TestMomentum:
    """Momentum: RSI (recursive), ROC, MOM."""

    def test_rsi(self) -> None:
        """RSI is recursive (Wilder smoothing). With 5000 bars, warmup
        convergence difference is negligible in the last 100."""
        for length in (7, 14, 21):
            qtl = rsi(CLOSE, length=length)
            pta = ta.rsi(SERIES, length=length)
            _verify_last_n(qtl, pta, tolerance=1e-6,
                           label=f"RSI({length})")

    @pytest.mark.xfail(reason="ROC formula differs: quantalib uses absolute "
                               "difference (close-prev), pandas-ta uses "
                               "percentage ((close/prev - 1)*100). "
                               "Known delta per SPEC §9.3.")
    def test_roc(self) -> None:
        """ROC: quantalib 'Roc' is Rate of Change (Absolute) = close - close[n].
        pandas-ta 'roc' is Rate of Change (Percentage) = ((c/c[n])-1)*100.
        These are fundamentally different indicators."""
        qtl = roc(CLOSE, length=10)
        pta = ta.roc(SERIES, length=10)
        _verify_last_n(qtl, pta, tolerance=1e-7, label="ROC(10)")

    def test_mom(self) -> None:
        for length in (5, 10, 20):
            qtl = mom(CLOSE, length=length)
            pta = ta.mom(SERIES, length=length)
            _verify_last_n(qtl, pta, tolerance=1e-9,
                           label=f"MOM({length})")


# ===========================================================================
# Statistics
# ===========================================================================

class TestStatistics:
    """STDDEV, VARIANCE, ZSCORE.

    Note: pandas-ta uses sample stddev (ddof=1), quantalib may use population.
    With 5000 bars and period=20, the difference is ~5% for ddof effect.
    We use relative tolerance where needed.
    """

    def test_stddev(self) -> None:
        for length in (10, 20, 50):
            qtl = stddev(CLOSE, length=length)
            pta = ta.stdev(SERIES, length=length)
            # pandas-ta uses ddof=1 (sample), quantalib may use ddof=0 (population)
            # With long lookback, ratio = sqrt((n-1)/n) ≈ 1 - 1/(2n)
            # For n=20: ratio ≈ 0.975, so try both
            pta_np = pta.to_numpy()
            count = len(qtl)
            start = max(0, count - VERIFY_COUNT)
            q_tail = qtl[start:]
            p_tail = pta_np[start:]
            finite = np.isfinite(q_tail) & np.isfinite(p_tail)

            if np.sum(finite) == 0:
                pytest.fail(f"STDDEV({length}): no finite in tail")

            q_f = q_tail[finite]
            p_f = p_tail[finite]

            # Try direct
            max_diff = float(np.max(np.abs(q_f - p_f)))
            if max_diff <= 1e-8:
                return

            # Try population-to-sample adjustment
            # sample_std = pop_std * sqrt(n/(n-1))
            adjusted = q_f * np.sqrt(length / (length - 1))
            adj_diff = float(np.max(np.abs(adjusted - p_f)))
            if adj_diff <= 1e-8:
                return

            # Try inverse adjustment
            adjusted_inv = q_f * np.sqrt((length - 1) / length)
            inv_diff = float(np.max(np.abs(adjusted_inv - p_f)))
            if inv_diff <= 1e-8:
                return

            pytest.fail(
                f"STDDEV({length}): direct={max_diff:.2e}, "
                f"pop→sample={adj_diff:.2e}, sample→pop={inv_diff:.2e}"
            )

    def test_variance(self) -> None:
        for length in (10, 20, 50):
            qtl = variance(CLOSE, length=length)
            pta = ta.variance(SERIES, length=length)
            pta_np = pta.to_numpy()
            count = len(qtl)
            start = max(0, count - VERIFY_COUNT)
            q_tail = qtl[start:]
            p_tail = pta_np[start:]
            finite = np.isfinite(q_tail) & np.isfinite(p_tail)

            if np.sum(finite) == 0:
                pytest.fail(f"VARIANCE({length}): no finite in tail")

            q_f = q_tail[finite]
            p_f = p_tail[finite]

            # Try direct
            max_diff = float(np.max(np.abs(q_f - p_f)))
            if max_diff <= 1e-8:
                return

            # Try ddof adjustment: var_sample = var_pop * n/(n-1)
            adjusted = q_f * (length / (length - 1))
            adj_diff = float(np.max(np.abs(adjusted - p_f)))
            if adj_diff <= 1e-8:
                return

            adjusted_inv = q_f * ((length - 1) / length)
            inv_diff = float(np.max(np.abs(adjusted_inv - p_f)))
            if inv_diff <= 1e-8:
                return

            pytest.fail(
                f"VARIANCE({length}): direct={max_diff:.2e}, "
                f"pop→sample={adj_diff:.2e}, sample→pop={inv_diff:.2e}"
            )

    def test_zscore(self) -> None:
        """ZSCORE = (x - mean) / stddev.

        quantalib uses population stddev (ddof=0), pandas-ta uses sample (ddof=1).
        The ratio is sqrt(n/(n-1)). We verify after applying the correction factor.
        """
        for length in (10, 20, 50):
            qtl = zscore(CLOSE, length=length)
            pta = ta.zscore(SERIES, length=length)
            pta_np = pta.to_numpy()
            count = len(qtl)
            start = max(0, count - VERIFY_COUNT)
            q_tail = qtl[start:]
            p_tail = pta_np[start:]
            finite = np.isfinite(q_tail) & np.isfinite(p_tail)

            if np.sum(finite) == 0:
                pytest.fail(f"ZSCORE({length}): no finite in tail")

            q_f = q_tail[finite]
            p_f = p_tail[finite]

            # Correct for ddof difference:
            # z_pop = (x - mean) / std_pop
            # z_sample = (x - mean) / std_sample
            # std_sample = std_pop * sqrt(n/(n-1))
            # so z_pop = z_sample * sqrt(n/(n-1))
            ddof_ratio = np.sqrt(length / (length - 1))

            # Try both correction directions
            diff_direct = float(np.max(np.abs(q_f - p_f)))
            diff_corrected = float(np.max(np.abs(q_f - p_f * ddof_ratio)))
            diff_inv = float(np.max(np.abs(q_f / ddof_ratio - p_f)))

            best = min(diff_direct, diff_corrected, diff_inv)
            assert best <= 1e-7, (
                f"ZSCORE({length}): best_diff={best:.2e} "
                f"(direct={diff_direct:.2e}, corrected={diff_corrected:.2e}, "
                f"inv={diff_inv:.2e})"
            )


# ===========================================================================
# Multi-output: Bollinger Bands
# ===========================================================================

class TestMultiOutput:
    """Multi-output indicators: BBands returns (upper, middle, lower) tuple."""

    def _get_bbands(self, length: int = 20, std: float = 2.0):
        """Get both quantalib and pandas-ta BBands results."""
        qtl = bbands(CLOSE, length=length, std=std)
        # quantalib uses population stddev (ddof=0); tell pandas-ta to match
        pta = ta.bbands(SERIES, length=length, std=std, ddof=0)

        # quantalib returns tuple of 3 numpy arrays: (upper, middle, lower)
        if isinstance(qtl, tuple):
            qtl_upper, qtl_mid, qtl_lower = qtl[0], qtl[1], qtl[2]
        elif hasattr(qtl, 'ndim') and qtl.ndim == 2:
            qtl_upper, qtl_mid, qtl_lower = qtl[:, 0], qtl[:, 1], qtl[:, 2]
        else:
            pytest.fail(f"Unexpected bbands return type: {type(qtl)}")

        # pandas-ta column names vary by version:
        # v0.4+: "BBL_20_2.0_2.0", "BBM_20_2.0_2.0", "BBU_20_2.0_2.0"
        # older:  "BBL_20_2.0",     "BBM_20_2.0",     "BBU_20_2.0"
        cols = list(pta.columns)
        bbu = [c for c in cols if c.startswith("BBU")]
        bbm = [c for c in cols if c.startswith("BBM")]
        bbl = [c for c in cols if c.startswith("BBL")]
        assert bbu and bbm and bbl, f"BBands columns not found: {cols}"

        pta_upper = pta[bbu[0]].to_numpy()
        pta_mid = pta[bbm[0]].to_numpy()
        pta_lower = pta[bbl[0]].to_numpy()

        return (qtl_upper, qtl_mid, qtl_lower), (pta_upper, pta_mid, pta_lower)

    def test_bbands_middle(self) -> None:
        """Middle band = SMA, should match exactly."""
        (_, q_mid, _), (_, p_mid, _) = self._get_bbands()
        _verify_last_n(q_mid, p_mid, tolerance=1e-9,
                       label="BBands middle")

    def test_bbands_upper(self) -> None:
        (q_upper, _, _), (p_upper, _, _) = self._get_bbands()
        # Tolerance depends on stddev ddof agreement
        _verify_last_n(q_upper, p_upper, tolerance=1e-6,
                       label="BBands upper")

    def test_bbands_lower(self) -> None:
        (_, _, q_lower), (_, _, p_lower) = self._get_bbands()
        _verify_last_n(q_lower, p_lower, tolerance=1e-6,
                       label="BBands lower")


# ===========================================================================
# Shape contract tests — output length must match input length
# ===========================================================================

class TestShape:
    """Verify output shapes match input for single-output indicators."""

    @pytest.mark.parametrize("indicator,length", [
        ("sma", 20), ("ema", 14), ("wma", 10), ("rsi", 14),
        ("mom", 10), ("roc", 10), ("stddev", 20), ("hma", 14),
    ])
    def test_output_length(self, indicator: str, length: int) -> None:
        fn = globals().get(indicator) or locals().get(indicator)
        if fn is None:
            fn = eval(indicator)  # noqa: S307
        result = fn(CLOSE, length=length)
        assert len(result) == N, (
            f"{indicator}({length}) output={len(result)} != input={N}"
        )


# ===========================================================================
# Performance comparison (informational, no assertions)
# ===========================================================================

class TestPerformance:
    """Throughput comparison. Uses 10K bars, 100 iterations.
    Results are printed, not asserted — documenting speedup only."""

    N_PERF = 10_000
    PERF_CLOSE = _generate_gbm(N_PERF, seed=99)
    PERF_SERIES = pd.Series(PERF_CLOSE, name="close")
    N_ITER = 100

    @pytest.mark.parametrize("name,qtl_fn,pta_fn,kwargs", [
        ("SMA(20)", sma, lambda s: ta.sma(s, length=20), {"length": 20}),
        ("EMA(20)", ema, lambda s: ta.ema(s, length=20), {"length": 20}),
        ("RSI(14)", rsi, lambda s: ta.rsi(s, length=14), {"length": 14}),
        ("WMA(14)", wma, lambda s: ta.wma(s, length=14), {"length": 14}),
        ("MOM(10)", mom, lambda s: ta.mom(s, length=10), {"length": 10}),
    ])
    def test_throughput(self, name: str, qtl_fn, pta_fn, kwargs: dict) -> None:
        import time

        data = self.PERF_CLOSE
        series = self.PERF_SERIES

        # quantalib
        t0 = time.perf_counter()
        for _ in range(self.N_ITER):
            _ = qtl_fn(data, **kwargs)
        qtl_us = (time.perf_counter() - t0) / self.N_ITER * 1e6

        # pandas-ta
        t0 = time.perf_counter()
        for _ in range(self.N_ITER):
            _ = pta_fn(series)
        pta_us = (time.perf_counter() - t0) / self.N_ITER * 1e6

        ratio = pta_us / qtl_us if qtl_us > 0 else float("inf")

        print(f"\n  {name} on {self.N_PERF:,} bars:")
        print(f"    quantalib : {qtl_us:8.1f} µs/call")
        print(f"    pandas-ta : {pta_us:8.1f} µs/call")
        print(f"    speedup   : {ratio:.1f}x")