""" Compare ferro_ta derivatives analytics against analytical references and optional third-party libraries available in the current environment. The suite focuses on selected core workflows: - Black-Scholes-Merton call pricing - implied-volatility recovery - first-order Greeks - Black-76 call pricing Outputs include: - speed timings with per-run samples and variance stats - analytical accuracy metrics - Python-tracked peak allocation snapshots - machine, runtime, build, and package metadata """ from __future__ import annotations import argparse import importlib.util import json import math import sys import time import tracemalloc from collections.abc import Callable from dataclasses import dataclass from pathlib import Path from typing import Any import numpy as np from ferro_ta.analysis.options import ( black_76_price as ft_black_76_price, ) from ferro_ta.analysis.options import ( greeks as ft_greeks, ) from ferro_ta.analysis.options import ( implied_volatility as ft_implied_volatility, ) from ferro_ta.analysis.options import ( option_price as ft_option_price, ) try: from benchmarks.metadata import benchmark_metadata, package_versions except ModuleNotFoundError: # pragma: no cover - script execution fallback from metadata import benchmark_metadata, package_versions N_WARMUP = 1 N_RUNS = 7 DEFAULT_SIZES = [1_000, 10_000] DEFAULT_ACCURACY_SIZE = 512 DEFAULT_SEED = 42 INV_SQRT_2PI = 1.0 / math.sqrt(2.0 * math.pi) @dataclass(frozen=True) class Case: name: str label: str expected_key: str component_names: tuple[str, ...] | None = None accuracy_target: str = "expected_output" @dataclass(frozen=True) class Provider: name: str kind: str note: str functions: dict[str, Callable[[dict[str, np.ndarray]], np.ndarray]] max_speed_size: int | None = None def supports(self, case_name: str) -> bool: return case_name in self.functions CASES = [ Case("bsm_call_price", "BSM Call Price", "call_price"), Case( "bsm_call_iv", "BSM Call IV Recovery", "volatility", accuracy_target="reconstructed_price", ), Case( "bsm_call_greeks", "BSM Call Greeks", "call_greeks", component_names=("delta", "gamma", "vega", "theta", "rho"), ), Case("black76_call_price", "Black-76 Call Price", "black76_call_price"), ] def _median(values: list[float]) -> float: ordered = sorted(values) mid = len(ordered) // 2 if len(ordered) % 2: return ordered[mid] return (ordered[mid - 1] + ordered[mid]) / 2.0 def _summary_stats(samples_ms: list[float]) -> dict[str, float]: if not samples_ms: return { "median_ms": 0.0, "mean_ms": 0.0, "min_ms": 0.0, "max_ms": 0.0, "stddev_ms": 0.0, "cv_pct": 0.0, } mean_ms = sum(samples_ms) / len(samples_ms) variance = ( sum((sample - mean_ms) ** 2 for sample in samples_ms) / (len(samples_ms) - 1) if len(samples_ms) > 1 else 0.0 ) stddev_ms = math.sqrt(variance) cv_pct = (stddev_ms / mean_ms * 100.0) if mean_ms else 0.0 return { "median_ms": round(_median(samples_ms), 4), "mean_ms": round(mean_ms, 4), "min_ms": round(min(samples_ms), 4), "max_ms": round(max(samples_ms), 4), "stddev_ms": round(stddev_ms, 4), "cv_pct": round(cv_pct, 3), } def _timed_runs_ms( fn: Callable[[dict[str, np.ndarray]], np.ndarray], chain: dict[str, np.ndarray], ) -> list[float]: for _ in range(N_WARMUP): fn(chain) samples_ms: list[float] = [] for _ in range(N_RUNS): t0 = time.perf_counter() fn(chain) samples_ms.append((time.perf_counter() - t0) * 1000.0) return samples_ms def _python_peak_bytes( fn: Callable[[dict[str, np.ndarray]], np.ndarray], chain: dict[str, np.ndarray], ) -> int | None: try: tracemalloc.start() tracemalloc.reset_peak() fn(chain) _, peak = tracemalloc.get_traced_memory() return int(peak) except Exception: return None finally: tracemalloc.stop() def _throughput_contracts_s(size: int, median_ms: float) -> float: if median_ms <= 0: return 0.0 return size / (median_ms / 1000.0) def _normal_pdf_numpy(x: np.ndarray) -> np.ndarray: return INV_SQRT_2PI * np.exp(-0.5 * x * x) def _normal_cdf_numpy(x: np.ndarray) -> np.ndarray: abs_x = np.abs(x) t = 1.0 / (1.0 + 0.2316419 * abs_x) poly = ( (((((1.330274429 * t) - 1.821255978) * t) + 1.781477937) * t - 0.356563782) * t + 0.319381530 ) * t cdf = 1.0 - _normal_pdf_numpy(abs_x) * poly return np.where(x >= 0.0, cdf, 1.0 - cdf) def _normal_cdf_scalar(x: float) -> float: return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0))) def _normal_pdf_scalar(x: float) -> float: return INV_SQRT_2PI * math.exp(-0.5 * x * x) def _bsm_price_numpy( spot: np.ndarray, strike: np.ndarray, rate: np.ndarray, time_to_expiry: np.ndarray, volatility: np.ndarray, *, option_type: str, carry: np.ndarray, ) -> np.ndarray: sqrt_t = np.sqrt(time_to_expiry) sigma_sqrt_t = volatility * sqrt_t d1 = ( np.log(spot / strike) + (rate - carry + 0.5 * volatility * volatility) * time_to_expiry ) / sigma_sqrt_t d2 = d1 - sigma_sqrt_t spot_df = np.exp(-carry * time_to_expiry) strike_df = np.exp(-rate * time_to_expiry) if option_type == "call": out = spot * spot_df * _normal_cdf_numpy( d1 ) - strike * strike_df * _normal_cdf_numpy(d2) else: out = strike * strike_df * _normal_cdf_numpy( -d2 ) - spot * spot_df * _normal_cdf_numpy(-d1) return np.ascontiguousarray(out, dtype=np.float64) def _black76_price_numpy( forward: np.ndarray, strike: np.ndarray, rate: np.ndarray, time_to_expiry: np.ndarray, volatility: np.ndarray, *, option_type: str, ) -> np.ndarray: sqrt_t = np.sqrt(time_to_expiry) sigma_sqrt_t = volatility * sqrt_t d1 = ( np.log(forward / strike) + 0.5 * volatility * volatility * time_to_expiry ) / sigma_sqrt_t d2 = d1 - sigma_sqrt_t discount = np.exp(-rate * time_to_expiry) if option_type == "call": out = discount * ( forward * _normal_cdf_numpy(d1) - strike * _normal_cdf_numpy(d2) ) else: out = discount * ( strike * _normal_cdf_numpy(-d2) - forward * _normal_cdf_numpy(-d1) ) return np.ascontiguousarray(out, dtype=np.float64) def _bsm_greeks_numpy( spot: np.ndarray, strike: np.ndarray, rate: np.ndarray, time_to_expiry: np.ndarray, volatility: np.ndarray, *, option_type: str, carry: np.ndarray, ) -> np.ndarray: sqrt_t = np.sqrt(time_to_expiry) sigma_sqrt_t = volatility * sqrt_t d1 = ( np.log(spot / strike) + (rate - carry + 0.5 * volatility * volatility) * time_to_expiry ) / sigma_sqrt_t d2 = d1 - sigma_sqrt_t pdf = _normal_pdf_numpy(d1) carry_df = np.exp(-carry * time_to_expiry) strike_df = np.exp(-rate * time_to_expiry) if option_type == "call": delta = carry_df * _normal_cdf_numpy(d1) theta = ( -(spot * carry_df * pdf * volatility) / (2.0 * sqrt_t) - rate * strike * strike_df * _normal_cdf_numpy(d2) + carry * spot * carry_df * _normal_cdf_numpy(d1) ) rho = strike * time_to_expiry * strike_df * _normal_cdf_numpy(d2) else: delta = carry_df * (_normal_cdf_numpy(d1) - 1.0) theta = ( -(spot * carry_df * pdf * volatility) / (2.0 * sqrt_t) + rate * strike * strike_df * _normal_cdf_numpy(-d2) - carry * spot * carry_df * _normal_cdf_numpy(-d1) ) rho = -strike * time_to_expiry * strike_df * _normal_cdf_numpy(-d2) gamma = carry_df * pdf / (spot * sigma_sqrt_t) vega = spot * carry_df * pdf * sqrt_t return np.ascontiguousarray( np.column_stack([delta, gamma, vega, theta, rho]), dtype=np.float64, ) def _bsm_price_scalar( spot: float, strike: float, rate: float, time_to_expiry: float, volatility: float, *, option_type: str, carry: float, ) -> float: sqrt_t = math.sqrt(time_to_expiry) sigma_sqrt_t = volatility * sqrt_t d1 = ( math.log(spot / strike) + (rate - carry + 0.5 * volatility * volatility) * time_to_expiry ) / sigma_sqrt_t d2 = d1 - sigma_sqrt_t spot_df = math.exp(-carry * time_to_expiry) strike_df = math.exp(-rate * time_to_expiry) if option_type == "call": return spot * spot_df * _normal_cdf_scalar( d1 ) - strike * strike_df * _normal_cdf_scalar(d2) return strike * strike_df * _normal_cdf_scalar( -d2 ) - spot * spot_df * _normal_cdf_scalar(-d1) def _black76_price_scalar( forward: float, strike: float, rate: float, time_to_expiry: float, volatility: float, *, option_type: str, ) -> float: sqrt_t = math.sqrt(time_to_expiry) sigma_sqrt_t = volatility * sqrt_t d1 = ( math.log(forward / strike) + 0.5 * volatility * volatility * time_to_expiry ) / sigma_sqrt_t d2 = d1 - sigma_sqrt_t discount = math.exp(-rate * time_to_expiry) if option_type == "call": return discount * ( forward * _normal_cdf_scalar(d1) - strike * _normal_cdf_scalar(d2) ) return discount * ( strike * _normal_cdf_scalar(-d2) - forward * _normal_cdf_scalar(-d1) ) def _bsm_greeks_scalar( spot: float, strike: float, rate: float, time_to_expiry: float, volatility: float, *, option_type: str, carry: float, ) -> tuple[float, float, float, float, float]: sqrt_t = math.sqrt(time_to_expiry) sigma_sqrt_t = volatility * sqrt_t d1 = ( math.log(spot / strike) + (rate - carry + 0.5 * volatility * volatility) * time_to_expiry ) / sigma_sqrt_t d2 = d1 - sigma_sqrt_t pdf = _normal_pdf_scalar(d1) carry_df = math.exp(-carry * time_to_expiry) strike_df = math.exp(-rate * time_to_expiry) if option_type == "call": delta = carry_df * _normal_cdf_scalar(d1) theta = ( -(spot * carry_df * pdf * volatility) / (2.0 * sqrt_t) - rate * strike * strike_df * _normal_cdf_scalar(d2) + carry * spot * carry_df * _normal_cdf_scalar(d1) ) rho = strike * time_to_expiry * strike_df * _normal_cdf_scalar(d2) else: delta = carry_df * (_normal_cdf_scalar(d1) - 1.0) theta = ( -(spot * carry_df * pdf * volatility) / (2.0 * sqrt_t) + rate * strike * strike_df * _normal_cdf_scalar(-d2) - carry * spot * carry_df * _normal_cdf_scalar(-d1) ) rho = -strike * time_to_expiry * strike_df * _normal_cdf_scalar(-d2) gamma = carry_df * pdf / (spot * sigma_sqrt_t) vega = spot * carry_df * pdf * sqrt_t return delta, gamma, vega, theta, rho def _implied_vol_bisection_numpy( price: np.ndarray, spot: np.ndarray, strike: np.ndarray, rate: np.ndarray, time_to_expiry: np.ndarray, *, option_type: str, carry: np.ndarray, lower: float = 1e-6, upper: float = 5.0, tolerance: float = 1e-8, max_iterations: int = 100, ) -> np.ndarray: lo = np.full_like(price, lower, dtype=np.float64) hi = np.full_like(price, upper, dtype=np.float64) mid = np.full_like(price, 0.2, dtype=np.float64) for _ in range(max_iterations): mid = (lo + hi) / 2.0 estimate = _bsm_price_numpy( spot, strike, rate, time_to_expiry, mid, option_type=option_type, carry=carry, ) too_low = estimate < price lo = np.where(too_low, mid, lo) hi = np.where(too_low, hi, mid) if float(np.max(np.abs(estimate - price))) < tolerance: break return np.ascontiguousarray(mid, dtype=np.float64) def _implied_vol_bisection_scalar( price: float, spot: float, strike: float, rate: float, time_to_expiry: float, *, option_type: str, carry: float, lower: float = 1e-6, upper: float = 5.0, tolerance: float = 1e-10, max_iterations: int = 100, ) -> float: lo = lower hi = upper mid = 0.2 for _ in range(max_iterations): mid = (lo + hi) / 2.0 estimate = _bsm_price_scalar( spot, strike, rate, time_to_expiry, mid, option_type=option_type, carry=carry, ) if abs(estimate - price) < tolerance: return mid if estimate < price: lo = mid else: hi = mid return mid def _reference_python_loop( chain: dict[str, np.ndarray], fn: Callable[..., float | tuple[float, ...]], *, include_forward: bool = False, include_price: bool = False, ) -> np.ndarray: rows: list[Any] = [] for idx in range(len(chain["strike"])): kwargs: dict[str, float | str] = { "strike": float(chain["strike"][idx]), "rate": float(chain["rate"][idx]), "time_to_expiry": float(chain["time_to_expiry"][idx]), "volatility": float(chain["volatility"][idx]), "carry": float(chain["carry"][idx]), } if include_forward: kwargs["forward"] = float(chain["forward"][idx]) else: kwargs["spot"] = float(chain["spot"][idx]) if include_price: kwargs["price"] = float(chain["call_price"][idx]) rows.append(fn(**kwargs)) return np.asarray(rows, dtype=np.float64) def _build_chain(n: int, *, seed: int) -> dict[str, np.ndarray]: rng = np.random.default_rng(seed) spot = np.ascontiguousarray(rng.uniform(80.0, 120.0, size=n), dtype=np.float64) strike = np.ascontiguousarray(rng.uniform(70.0, 130.0, size=n), dtype=np.float64) rate = np.ascontiguousarray(rng.uniform(0.0, 0.07, size=n), dtype=np.float64) carry = np.ascontiguousarray(rng.uniform(0.0, 0.03, size=n), dtype=np.float64) time_to_expiry = np.ascontiguousarray( rng.uniform(7.0 / 365.0, 2.0, size=n), dtype=np.float64 ) volatility = np.ascontiguousarray(rng.uniform(0.08, 0.65, size=n), dtype=np.float64) forward = np.ascontiguousarray( spot * np.exp((rate - carry) * time_to_expiry), dtype=np.float64, ) chain = { "spot": spot, "strike": strike, "rate": rate, "carry": carry, "time_to_expiry": time_to_expiry, "volatility": volatility, "forward": forward, } call_price = _reference_python_loop( chain, lambda *, spot, strike, rate, time_to_expiry, volatility, carry: ( _bsm_price_scalar( spot, strike, rate, time_to_expiry, volatility, option_type="call", carry=carry, ) ), ) call_greeks = _reference_python_loop( chain, lambda *, spot, strike, rate, time_to_expiry, volatility, carry: ( _bsm_greeks_scalar( spot, strike, rate, time_to_expiry, volatility, option_type="call", carry=carry, ) ), ) black76_call_price = _reference_python_loop( chain, lambda *, forward, strike, rate, time_to_expiry, volatility, carry: ( _black76_price_scalar( forward, strike, rate, time_to_expiry, volatility, option_type="call", ) ), include_forward=True, ) chain["call_price"] = np.ascontiguousarray(call_price, dtype=np.float64) chain["call_greeks"] = np.ascontiguousarray(call_greeks, dtype=np.float64) chain["black76_call_price"] = np.ascontiguousarray( black76_call_price, dtype=np.float64 ) return chain def _ferro_ta_call_price(chain: dict[str, np.ndarray]) -> np.ndarray: return np.asarray( ft_option_price( chain["spot"], chain["strike"], chain["rate"], chain["time_to_expiry"], chain["volatility"], option_type="call", model="bsm", carry=chain["carry"], ), dtype=np.float64, ) def _ferro_ta_call_iv(chain: dict[str, np.ndarray]) -> np.ndarray: return np.asarray( ft_implied_volatility( chain["call_price"], chain["spot"], chain["strike"], chain["rate"], chain["time_to_expiry"], option_type="call", model="bsm", carry=chain["carry"], ), dtype=np.float64, ) def _ferro_ta_call_greeks(chain: dict[str, np.ndarray]) -> np.ndarray: result = ft_greeks( chain["spot"], chain["strike"], chain["rate"], chain["time_to_expiry"], chain["volatility"], option_type="call", model="bsm", carry=chain["carry"], ) return np.ascontiguousarray( np.column_stack( [ np.asarray(result.delta, dtype=np.float64), np.asarray(result.gamma, dtype=np.float64), np.asarray(result.vega, dtype=np.float64), np.asarray(result.theta, dtype=np.float64), np.asarray(result.rho, dtype=np.float64), ] ), dtype=np.float64, ) def _ferro_ta_black76_call_price(chain: dict[str, np.ndarray]) -> np.ndarray: return np.asarray( ft_black_76_price( chain["forward"], chain["strike"], chain["rate"], chain["time_to_expiry"], chain["volatility"], option_type="call", ), dtype=np.float64, ) def _reference_numpy_call_price(chain: dict[str, np.ndarray]) -> np.ndarray: return _bsm_price_numpy( chain["spot"], chain["strike"], chain["rate"], chain["time_to_expiry"], chain["volatility"], option_type="call", carry=chain["carry"], ) def _reference_numpy_call_iv(chain: dict[str, np.ndarray]) -> np.ndarray: return _implied_vol_bisection_numpy( chain["call_price"], chain["spot"], chain["strike"], chain["rate"], chain["time_to_expiry"], option_type="call", carry=chain["carry"], ) def _reference_numpy_call_greeks(chain: dict[str, np.ndarray]) -> np.ndarray: return _bsm_greeks_numpy( chain["spot"], chain["strike"], chain["rate"], chain["time_to_expiry"], chain["volatility"], option_type="call", carry=chain["carry"], ) def _reference_numpy_black76_call_price(chain: dict[str, np.ndarray]) -> np.ndarray: return _black76_price_numpy( chain["forward"], chain["strike"], chain["rate"], chain["time_to_expiry"], chain["volatility"], option_type="call", ) def _reference_python_call_price(chain: dict[str, np.ndarray]) -> np.ndarray: return _reference_python_loop( chain, lambda *, spot, strike, rate, time_to_expiry, volatility, carry: ( _bsm_price_scalar( spot, strike, rate, time_to_expiry, volatility, option_type="call", carry=carry, ) ), ) def _reference_python_call_iv(chain: dict[str, np.ndarray]) -> np.ndarray: return _reference_python_loop( chain, lambda *, price, spot, strike, rate, time_to_expiry, volatility, carry: ( _implied_vol_bisection_scalar( price, spot, strike, rate, time_to_expiry, option_type="call", carry=carry, ) ), include_price=True, ) def _reference_python_call_greeks(chain: dict[str, np.ndarray]) -> np.ndarray: return _reference_python_loop( chain, lambda *, spot, strike, rate, time_to_expiry, volatility, carry: ( _bsm_greeks_scalar( spot, strike, rate, time_to_expiry, volatility, option_type="call", carry=carry, ) ), ) def _reference_python_black76_call_price(chain: dict[str, np.ndarray]) -> np.ndarray: return _reference_python_loop( chain, lambda *, forward, strike, rate, time_to_expiry, volatility, carry: ( _black76_price_scalar( forward, strike, rate, time_to_expiry, volatility, option_type="call", ) ), include_forward=True, ) def _reprice_bsm_call_from_iv( chain: dict[str, np.ndarray], implied_vols: np.ndarray, ) -> np.ndarray: rows = [ _bsm_price_scalar( float(spot), float(strike), float(rate), float(time_to_expiry), max(float(iv), 1e-12), option_type="call", carry=float(carry), ) for spot, strike, rate, time_to_expiry, iv, carry in zip( chain["spot"], chain["strike"], chain["rate"], chain["time_to_expiry"], np.asarray(implied_vols, dtype=np.float64), chain["carry"], ) ] return np.asarray(rows, dtype=np.float64) def _py_vollib_provider() -> Provider | None: if importlib.util.find_spec("py_vollib") is None: return None from py_vollib.black_scholes_merton import black_scholes_merton as py_vollib_bsm from py_vollib.black_scholes_merton.implied_volatility import ( implied_volatility as py_vollib_iv, ) def _price(chain: dict[str, np.ndarray]) -> np.ndarray: return np.asarray( [ py_vollib_bsm( "c", float(s), float(k), float(t), float(r), float(vol), float(q), ) for s, k, r, t, vol, q in zip( chain["spot"], chain["strike"], chain["rate"], chain["time_to_expiry"], chain["volatility"], chain["carry"], ) ], dtype=np.float64, ) def _iv(chain: dict[str, np.ndarray]) -> np.ndarray: return np.asarray( [ py_vollib_iv( float(price), "c", float(s), float(k), float(t), float(r), float(q), ) for price, s, k, r, t, q in zip( chain["call_price"], chain["spot"], chain["strike"], chain["rate"], chain["time_to_expiry"], chain["carry"], ) ], dtype=np.float64, ) return Provider( name="py_vollib", kind="third_party", note="Scalar Black-Scholes-Merton baseline from py_vollib.", functions={ "bsm_call_price": _price, "bsm_call_iv": _iv, }, max_speed_size=1_000, ) def available_providers() -> list[Provider]: providers = [ Provider( name="ferro_ta", kind="project", note="Rust-backed vectorized implementation.", functions={ "bsm_call_price": _ferro_ta_call_price, "bsm_call_iv": _ferro_ta_call_iv, "bsm_call_greeks": _ferro_ta_call_greeks, "black76_call_price": _ferro_ta_black76_call_price, }, ), Provider( name="reference_numpy", kind="reference", note="Pure NumPy analytical formulas with vectorized IV bisection.", functions={ "bsm_call_price": _reference_numpy_call_price, "bsm_call_iv": _reference_numpy_call_iv, "bsm_call_greeks": _reference_numpy_call_greeks, "black76_call_price": _reference_numpy_black76_call_price, }, ), Provider( name="reference_python_loop", kind="reference", note="Scalar math-loop analytical baseline; useful for accuracy sanity checks.", functions={ "bsm_call_price": _reference_python_call_price, "bsm_call_iv": _reference_python_call_iv, "bsm_call_greeks": _reference_python_call_greeks, "black76_call_price": _reference_python_black76_call_price, }, max_speed_size=1_000, ), ] optional = _py_vollib_provider() if optional is not None: providers.append(optional) return providers def _accuracy_metrics(actual: np.ndarray, expected: np.ndarray) -> dict[str, Any]: actual_arr = np.asarray(actual, dtype=np.float64) expected_arr = np.asarray(expected, dtype=np.float64) abs_error = np.abs(actual_arr - expected_arr) rel_error = abs_error / np.maximum(np.abs(expected_arr), 1e-12) return { "output_shape": list(actual_arr.shape), "max_abs_error": round(float(np.max(abs_error)), 12), "mean_abs_error": round(float(np.mean(abs_error)), 12), "rmse": round( float(np.sqrt(np.mean(np.square(actual_arr - expected_arr)))), 12 ), "max_rel_error": round(float(np.max(rel_error)), 12), } def _accuracy_summary(rows: list[dict[str, Any]]) -> list[dict[str, Any]]: summary: list[dict[str, Any]] = [] for case in CASES: case_rows = [ row for row in rows if row["case"] == case.name and "max_abs_error" in row ] if not case_rows: continue best = min(case_rows, key=lambda row: float(row["max_abs_error"])) worst = max(case_rows, key=lambda row: float(row["max_abs_error"])) summary.append( { "case": case.name, "best_provider": best["provider"], "best_max_abs_error": best["max_abs_error"], "worst_provider": worst["provider"], "worst_max_abs_error": worst["max_abs_error"], } ) return summary def _speed_summary(rows: list[dict[str, Any]]) -> list[dict[str, Any]]: summary: list[dict[str, Any]] = [] for case in CASES: for size in sorted( { int(row["size"]) for row in rows if row["case"] == case.name and "median_ms" in row } ): case_rows = [ row for row in rows if row["case"] == case.name and row.get("size") == size and "median_ms" in row ] if not case_rows: continue fastest = min(case_rows, key=lambda row: float(row["median_ms"])) ranking = [ { "provider": row["provider"], "median_ms": row["median_ms"], "contracts_per_s": row["contracts_per_s"], } for row in sorted(case_rows, key=lambda row: float(row["median_ms"])) ] summary.append( { "case": case.name, "size": size, "fastest_provider": fastest["provider"], "fastest_median_ms": fastest["median_ms"], "ranking": ranking, } ) return summary def _print_provider_inventory(providers: list[Provider]) -> None: print("Providers:") for provider in providers: supported = ", ".join( case.name for case in CASES if provider.supports(case.name) ) cap = ( f" (speed cap {provider.max_speed_size})" if provider.max_speed_size else "" ) print(f" - {provider.name} [{provider.kind}] {provider.note}{cap}") print(f" supported: {supported}") print() def _print_accuracy_table(rows: list[dict[str, Any]], accuracy_size: int) -> None: print(f"Accuracy ({accuracy_size} contracts)") print( "IV accuracy is measured as price reconstruction error from the recovered IV." ) header = f"{'Case':<22} {'Provider':<24} {'Max abs err':<14} {'RMSE':<14} {'Max rel err':<14}" print(header) print("-" * len(header)) for row in rows: if "max_abs_error" not in row: continue print( f"{row['label']:<22} {row['provider']:<24} " f"{row['max_abs_error']:<14.6g} {row['rmse']:<14.6g} {row['max_rel_error']:<14.6g}" ) print() def _print_speed_table(rows: list[dict[str, Any]]) -> None: print(f"Speed (median of {N_RUNS} measured runs after {N_WARMUP} warmup)") header = f"{'Case':<22} {'Size':<8} {'Provider':<24} {'Median ms':<12} {'Contracts/s':<14} {'Peak alloc':<12}" print(header) print("-" * len(header)) for row in rows: if "median_ms" not in row: continue peak = row.get("python_peak_allocation_bytes") peak_label = "n/a" if peak is None else str(peak) print( f"{row['label']:<22} {row['size']:<8} {row['provider']:<24} " f"{row['median_ms']:<12.4f} {row['contracts_per_s']:<14.2f} {peak_label:<12}" ) print() def run_benchmark( *, sizes: list[int], accuracy_size: int, json_path: str | None, ) -> dict[str, Any]: providers = available_providers() speed_chains = { size: _build_chain(size, seed=DEFAULT_SEED + size) for size in sizes } accuracy_chain = _build_chain(accuracy_size, seed=DEFAULT_SEED) accuracy_rows: list[dict[str, Any]] = [] speed_rows: list[dict[str, Any]] = [] _print_provider_inventory(providers) for case in CASES: for provider in providers: if not provider.supports(case.name): continue fn = provider.functions[case.name] actual = fn(accuracy_chain) if case.name == "bsm_call_iv": compared_actual = _reprice_bsm_call_from_iv(accuracy_chain, actual) expected = np.asarray(accuracy_chain["call_price"], dtype=np.float64) else: compared_actual = actual expected = np.asarray( accuracy_chain[case.expected_key], dtype=np.float64 ) row = { "case": case.name, "label": case.label, "provider": provider.name, "provider_kind": provider.kind, "sample_size": accuracy_size, "accuracy_target": case.accuracy_target, } row.update(_accuracy_metrics(compared_actual, expected)) if case.component_names is not None: row["component_names"] = list(case.component_names) accuracy_rows.append(row) _print_accuracy_table(accuracy_rows, accuracy_size) for case in CASES: for size in sizes: chain = speed_chains[size] for provider in providers: if not provider.supports(case.name): continue if ( provider.max_speed_size is not None and size > provider.max_speed_size ): continue fn = provider.functions[case.name] samples_ms = _timed_runs_ms(fn, chain) stats = _summary_stats(samples_ms) median_ms = float(stats["median_ms"]) contracts_per_s = _throughput_contracts_s(size, median_ms) peak_bytes = _python_peak_bytes(fn, chain) speed_rows.append( { "case": case.name, "label": case.label, "size": size, "provider": provider.name, "provider_kind": provider.kind, "median_ms": round(median_ms, 4), "contracts_per_s": round(contracts_per_s, 2), "runs_ms": [round(sample, 4) for sample in samples_ms], "stats": stats, "python_peak_allocation_bytes": peak_bytes, "input_layout": { "dtype": "float64", "contiguous": True, }, } ) _print_speed_table(speed_rows) package_names = ["numpy", "ferro-ta", "py_vollib"] metadata = benchmark_metadata( "benchmark_derivatives_compare", extra={ "dataset": { "generator": "synthetic_option_chain", "speed_sizes": sizes, "accuracy_size": accuracy_size, "dtype": "float64", "array_layout": "C-contiguous", "seed": DEFAULT_SEED, "ranges": { "spot": [80.0, 120.0], "strike": [70.0, 130.0], "rate": [0.0, 0.07], "carry": [0.0, 0.03], "time_to_expiry_years": [7.0 / 365.0, 2.0], "volatility": [0.08, 0.65], }, }, "methodology": { "warmup_runs": N_WARMUP, "measured_runs": N_RUNS, "reported_metric": "median_ms", "speed_metric": "contracts_per_second", "accuracy_reference": ( "Scalar analytical Black-Scholes-Merton and Black-76 formulas " "using math.erf; IV accuracy is measured as repriced error " "from the recovered volatility because direct volatility " "differences can be unstable on low-vega contracts." ), "input_layout_notes": ( "Benchmarks use contiguous float64 arrays. If your workload " "passes non-contiguous arrays or mixed dtypes, benchmark that " "path separately." ), "allocation_notes": ( "python_peak_allocation_bytes is a tracemalloc snapshot of " "Python-tracked allocations only; it does not measure native RSS." ), "provider_notes": ( "reference_python_loop and py_vollib are scalar baselines and " "are size-capped in the speed table to keep runtime reasonable." ), }, "providers": [ { "name": provider.name, "kind": provider.kind, "note": provider.note, "max_speed_size": provider.max_speed_size, "supported_cases": [ case.name for case in CASES if provider.supports(case.name) ], } for provider in providers ], "packages": package_versions(*package_names), }, ) result = { "schema_version": 1, "command": " ".join(["python", *sys.argv]), "n_warmup": N_WARMUP, "n_runs": N_RUNS, "accuracy_size": accuracy_size, "sizes": sizes, "metadata": metadata, "accuracy": { "summary": _accuracy_summary(accuracy_rows), "results": accuracy_rows, }, "speed": { "summary": _speed_summary(speed_rows), "results": speed_rows, }, } if json_path: output_path = Path(json_path) output_path.parent.mkdir(parents=True, exist_ok=True) output_path.write_text(json.dumps(result, indent=2), encoding="utf-8") print(f"Results written to {output_path}") return result def main() -> int: parser = argparse.ArgumentParser( description="Compare ferro_ta derivatives analytics against reference implementations" ) parser.add_argument( "--json", default=None, help="Write the benchmark artifact to JSON", ) parser.add_argument( "--sizes", type=int, nargs="+", default=DEFAULT_SIZES, help="Contract counts to benchmark (default: 1000 10000)", ) parser.add_argument( "--accuracy-size", type=int, default=DEFAULT_ACCURACY_SIZE, help="Contract count used for the accuracy pass (default: 512)", ) args = parser.parse_args() run_benchmark( sizes=args.sizes, accuracy_size=args.accuracy_size, json_path=args.json, ) return 0 if __name__ == "__main__": raise SystemExit(main())