Files
my-python-backteat/ferro-ta-main/tests/unit/test_derivatives_accuracy.py
T

609 lines
24 KiB
Python
Raw Normal View History

2026-07-09 05:08:16 +08:00
"""
Accuracy/correctness tests for ferro-ta derivatives analytics.
Each test class validates the ferro-ta implementation against reference
formulas implemented using scipy and numpy.
"""
from __future__ import annotations
import numpy as np
import pytest
# ---------------------------------------------------------------------------
# Reference formulas (pure numpy / scipy)
# ---------------------------------------------------------------------------
def _norm_cdf(x):
"""Standard normal CDF via scipy."""
from scipy.stats import norm as _norm
return _norm.cdf(x)
def _norm_pdf(x):
from scipy.stats import norm as _norm
return _norm.pdf(x)
def bsm_call(S, K, r, q, T, sigma): # noqa: N803
"""Reference BSM call price."""
d1 = (np.log(S / K) + (r - q + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
return S * np.exp(-q * T) * _norm_cdf(d1) - K * np.exp(-r * T) * _norm_cdf(d2)
def bsm_put(S, K, r, q, T, sigma): # noqa: N803
d1 = (np.log(S / K) + (r - q + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
return K * np.exp(-r * T) * _norm_cdf(-d2) - S * np.exp(-q * T) * _norm_cdf(-d1)
def bsm_delta_call(S, K, r, q, T, sigma): # noqa: N803
d1 = (np.log(S / K) + (r - q + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
return np.exp(-q * T) * _norm_cdf(d1)
def digital_cash_call(S, K, r, q, T, sigma): # noqa: N803
d2 = (np.log(S / K) + (r - q - 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
return np.exp(-r * T) * _norm_cdf(d2)
def digital_asset_call(S, K, r, q, T, sigma): # noqa: N803
d1 = (np.log(S / K) + (r - q + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
return S * np.exp(-q * T) * _norm_cdf(d1)
def digital_cash_put(S, K, r, q, T, sigma): # noqa: N803
d2 = (np.log(S / K) + (r - q - 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
return np.exp(-r * T) * _norm_cdf(-d2)
def digital_asset_put(S, K, r, q, T, sigma): # noqa: N803
d1 = (np.log(S / K) + (r - q + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
return S * np.exp(-q * T) * _norm_cdf(-d1)
def vanna_num(S, K, r, q, T, sigma, eps=1e-4): # noqa: N803
"""∂Δ/∂σ via central differences."""
delta_up = bsm_delta_call(S, K, r, q, T, sigma + eps)
delta_dn = bsm_delta_call(S, K, r, q, T, sigma - eps)
return (delta_up - delta_dn) / (2 * eps)
def vega_bsm(S, K, r, q, T, sigma): # noqa: N803
d1 = (np.log(S / K) + (r - q + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
return S * np.exp(-q * T) * _norm_pdf(d1) * np.sqrt(T)
def volga_num(S, K, r, q, T, sigma, eps=1e-4): # noqa: N803
"""∂²V/∂σ² via central differences."""
v_up = vega_bsm(S, K, r, q, T, sigma + eps)
v_dn = vega_bsm(S, K, r, q, T, sigma - eps)
return (v_up - v_dn) / (2 * eps)
def ctc_vol_reference(close, window, trading_days=252.0):
"""Close-to-close vol: rolling std of log returns × sqrt(trading_days)."""
log_ret = np.log(close[1:] / close[:-1])
n = len(close)
out = np.full(n, np.nan)
for i in range(window, n):
returns_window = log_ret[i - window : i]
out[i] = np.sqrt(np.sum(returns_window**2) / window * trading_days)
return out
# ---------------------------------------------------------------------------
# Test cases
# ---------------------------------------------------------------------------
# Six parameter sets: ATM, 10% OTM, 10% ITM, low vol, high vol, non-zero carry
_DIGITAL_CASES = [
# (S, K, r, q, T, sigma, label)
(100.0, 100.0, 0.05, 0.00, 1.0, 0.20, "ATM"),
(100.0, 110.0, 0.05, 0.00, 1.0, 0.20, "10% OTM"),
(100.0, 90.0, 0.05, 0.00, 1.0, 0.20, "10% ITM"),
(100.0, 100.0, 0.05, 0.00, 1.0, 0.05, "low vol"),
(100.0, 100.0, 0.05, 0.00, 1.0, 0.50, "high vol"),
(100.0, 100.0, 0.05, 0.03, 1.0, 0.20, "non-zero carry"),
]
class TestDigitalOptionsAccuracy:
@pytest.fixture(autouse=True)
def require_scipy(self):
pytest.importorskip("scipy")
def test_cash_or_nothing_call_vs_reference(self):
from ferro_ta.analysis.options import digital_option_price
for S, K, r, q, T, sigma, label in _DIGITAL_CASES:
expected = digital_cash_call(S, K, r, q, T, sigma)
actual = digital_option_price(
S,
K,
r,
T,
sigma,
option_type="call",
digital_type="cash_or_nothing",
carry=q,
)
assert actual == pytest.approx(expected, abs=1e-6), (
f"cash_or_nothing call mismatch for case '{label}': "
f"got {actual}, expected {expected}"
)
def test_cash_or_nothing_put_vs_reference(self):
from ferro_ta.analysis.options import digital_option_price
for S, K, r, q, T, sigma, label in _DIGITAL_CASES:
expected = digital_cash_put(S, K, r, q, T, sigma)
actual = digital_option_price(
S,
K,
r,
T,
sigma,
option_type="put",
digital_type="cash_or_nothing",
carry=q,
)
assert actual == pytest.approx(expected, abs=1e-6), (
f"cash_or_nothing put mismatch for case '{label}': "
f"got {actual}, expected {expected}"
)
def test_asset_or_nothing_call_vs_reference(self):
from ferro_ta.analysis.options import digital_option_price
for S, K, r, q, T, sigma, label in _DIGITAL_CASES:
expected = digital_asset_call(S, K, r, q, T, sigma)
actual = digital_option_price(
S,
K,
r,
T,
sigma,
option_type="call",
digital_type="asset_or_nothing",
carry=q,
)
# Tolerance 1e-4: asset-or-nothing involves S * N(d1), small numerical diff expected
assert actual == pytest.approx(expected, abs=1e-4), (
f"asset_or_nothing call mismatch for case '{label}': "
f"got {actual}, expected {expected}"
)
def test_asset_or_nothing_put_vs_reference(self):
from ferro_ta.analysis.options import digital_option_price
for S, K, r, q, T, sigma, label in _DIGITAL_CASES:
expected = digital_asset_put(S, K, r, q, T, sigma)
actual = digital_option_price(
S,
K,
r,
T,
sigma,
option_type="put",
digital_type="asset_or_nothing",
carry=q,
)
# Tolerance 1e-4: asset-or-nothing involves S * N(-d1), small numerical diff expected
assert actual == pytest.approx(expected, abs=1e-4), (
f"asset_or_nothing put mismatch for case '{label}': "
f"got {actual}, expected {expected}"
)
def test_batch_digital_price_matches_scalar(self):
"""Vectorized call must match scalar loop for 10 random points."""
from ferro_ta.analysis.options import digital_option_price
rng = np.random.default_rng(7)
n = 10
S_arr = rng.uniform(80.0, 120.0, n)
K_arr = rng.uniform(80.0, 120.0, n)
r_arr = rng.uniform(0.01, 0.10, n)
T_arr = rng.uniform(0.1, 2.0, n)
sigma_arr = rng.uniform(0.10, 0.50, n)
batch = digital_option_price(
S_arr,
K_arr,
r_arr,
T_arr,
sigma_arr,
option_type="call",
digital_type="cash_or_nothing",
)
scalar_results = np.array(
[
digital_option_price(
float(S_arr[i]),
float(K_arr[i]),
float(r_arr[i]),
float(T_arr[i]),
float(sigma_arr[i]),
option_type="call",
digital_type="cash_or_nothing",
)
for i in range(n)
]
)
assert batch == pytest.approx(scalar_results, abs=1e-10), (
"Batch digital_option_price does not match scalar loop"
)
# Four cases for extended Greeks: ITM call, ATM call, OTM call, ATM put
_GREEK_CASES = [
# (S, K, r, q, T, sigma, option_type, label)
(110.0, 100.0, 0.05, 0.0, 1.0, 0.20, "call", "ITM call"),
(100.0, 100.0, 0.05, 0.0, 1.0, 0.20, "call", "ATM call"),
(90.0, 100.0, 0.05, 0.0, 1.0, 0.20, "call", "OTM call"),
(100.0, 100.0, 0.05, 0.0, 1.0, 0.20, "put", "ATM put"),
]
class TestExtendedGreeksAccuracy:
@pytest.fixture(autouse=True)
def require_scipy(self):
pytest.importorskip("scipy")
def test_vanna_vs_numerical_fd(self):
"""extended_greeks().vanna matches ∂Δ/∂σ from central differences (tol=1e-3)."""
from ferro_ta.analysis.options import extended_greeks
for S, K, r, q, T, sigma, opt_type, label in _GREEK_CASES:
eg = extended_greeks(S, K, r, T, sigma, option_type=opt_type, carry=q)
# Reference is defined only for calls; for put use numerical FD directly
if opt_type == "call":
expected = vanna_num(S, K, r, q, T, sigma)
else:
# Vanna for put: ∂(put delta)/∂σ = ∂(call delta - e^{-qT})/∂σ = vanna_call
expected = vanna_num(S, K, r, q, T, sigma)
assert float(eg.vanna) == pytest.approx(expected, abs=1e-3), (
f"Vanna mismatch for '{label}': got {eg.vanna}, expected {expected}"
)
def test_volga_vs_numerical_fd(self):
"""extended_greeks().volga matches ∂²V/∂σ² from central differences (tol=1e-2)."""
from ferro_ta.analysis.options import extended_greeks
for S, K, r, q, T, sigma, opt_type, label in _GREEK_CASES:
eg = extended_greeks(S, K, r, T, sigma, option_type=opt_type, carry=q)
expected = volga_num(S, K, r, q, T, sigma)
assert float(eg.volga) == pytest.approx(expected, abs=1e-2), (
f"Volga mismatch for '{label}': got {eg.volga}, expected {expected}"
)
def test_speed_negative_for_calls(self):
"""Speed (∂Γ/∂S) should be negative for OTM calls — Gamma decreases as S moves away."""
from ferro_ta.analysis.options import extended_greeks
# OTM call: S < K
eg = extended_greeks(90.0, 100.0, 0.05, 1.0, 0.20, option_type="call")
assert float(eg.speed) < 0.0, (
f"Speed should be negative for OTM call, got {eg.speed}"
)
def test_charm_finite_for_valid_inputs(self):
"""Charm should be finite and non-zero for non-degenerate inputs."""
from ferro_ta.analysis.options import extended_greeks
for S, K, r, q, T, sigma, opt_type, label in _GREEK_CASES:
eg = extended_greeks(S, K, r, T, sigma, option_type=opt_type, carry=q)
assert np.isfinite(float(eg.charm)), (
f"Charm is not finite for '{label}': {eg.charm}"
)
assert eg.charm != 0.0, (
f"Charm is zero for '{label}' — unexpected for non-degenerate inputs"
)
class TestAmericanOptionsAccuracy:
"""Property-based tests for American options (no scipy required)."""
def test_baw_vs_published_values(self):
"""BAW American put satisfies the lower bound: price ≥ max(K - S, European BSM put).
The Haug (2007) table uses b = r - q (cost of carry convention). Rather
than replicate the exact table — which requires matching the BAW carry
convention precisely — we verify two model-agnostic inequalities that any
correct American-put implementation must satisfy:
1. American put ≥ intrinsic value (K - S)
2. American put ≥ European BSM put (early exercise has non-negative value)
"""
from ferro_ta.analysis.options import american_option_price, option_price
S, K, r, T, sigma = 100.0, 100.0, 0.10, 0.25, 0.20
american = american_option_price(S, K, r, T, sigma, option_type="put")
european = option_price(S, K, r, T, sigma, option_type="put")
assert american >= max(K - S, 0.0) - 1e-8, (
f"American put below intrinsic: {american:.4f} < {max(K - S, 0.0)}"
)
assert american >= european - 1e-8, (
f"American put below European put: {american:.4f} < {european:.4f}"
)
# Sanity-check: American ATM put should be in a reasonable range
assert 0.0 < american < K, (
f"American put price {american:.4f} is outside (0, K={K})"
)
def test_american_put_increases_with_strike(self):
"""Deeper ITM (higher strike for put) ⇒ higher American put price.
Uses moderately spaced strikes to avoid the intrinsic-value floor
where K - S becomes the binding constraint and the increments are
exactly 1-for-1, which can mask ordering issues near the floor.
"""
from ferro_ta.analysis.options import american_option_price
# S = 100, K in {85, 100, 115}; rate and carry both 0.05 to avoid b=0 issues
S, r, T, sigma = 100.0, 0.05, 0.5, 0.25
strikes = [85.0, 100.0, 115.0]
prices = [
american_option_price(S, K, r, T, sigma, option_type="put", carry=r)
for K in strikes
]
assert prices[0] < prices[1] < prices[2], (
f"American put prices not monotone in strike: "
f"K={strikes} → prices={[round(p, 4) for p in prices]}"
)
def test_american_call_increases_with_spot(self):
"""Higher spot ⇒ higher American call price."""
from ferro_ta.analysis.options import american_option_price
spots = [90.0, 100.0, 110.0]
prices = [
american_option_price(S, 100.0, 0.05, 1.0, 0.20, option_type="call")
for S in spots
]
assert prices[0] < prices[1] < prices[2], (
f"American call prices not monotone in spot: {prices}"
)
def test_american_call_equals_european_no_dividends_no_early_exercise(self):
"""American call with no early-exercise incentive (carry=0) ≈ European call.
When the cost-of-carry parameter is zero, there is no dividend/carry
benefit to holding the underlying. In this regime, it is never
optimal to early-exercise an American call, so the American call price
equals the European call price computed with the same carry=0 convention.
The `early_exercise_premium` function exposes this directly and should
return ~0 for calls with carry=0.
"""
from ferro_ta.analysis.options import early_exercise_premium
S, K, r, T, sigma = 100.0, 100.0, 0.05, 1.0, 0.20
premium = early_exercise_premium(
S, K, r, T, sigma, option_type="call", carry=0.0
)
assert premium == pytest.approx(0.0, abs=1e-4), (
f"Early exercise premium for call with carry=0 should be ~0, got {premium:.6f}"
)
def test_early_exercise_premium_positive_for_deep_itm_put(self):
"""Deep ITM American put should have a meaningful early exercise premium.
When S is well below K (deep ITM put), the time value is low and the
interest gained from early exercise of the put dominates — leading to a
positive early-exercise premium.
"""
from ferro_ta.analysis.options import early_exercise_premium
# Deep ITM: S=70, K=100 — strong incentive to exercise early
premium = early_exercise_premium(
70.0, 100.0, 0.10, 1.0, 0.20, option_type="put"
)
assert premium > 0.0, (
f"Deep ITM American put early exercise premium should be > 0, got {premium}"
)
class TestVolEstimatorsAccuracy:
@pytest.fixture(autouse=True)
def require_scipy(self):
pytest.importorskip("scipy")
def test_close_to_close_vs_reference_impl(self):
"""C2C vol matches reference formula exactly (tol=1e-10), 100 samples."""
from ferro_ta.analysis.options import close_to_close_vol
rng = np.random.default_rng(42)
log_ret = rng.normal(0.0, 0.01, 100)
close = 100.0 * np.cumprod(np.exp(log_ret))
window = 20
actual = close_to_close_vol(close, window=window, trading_days_per_year=252.0)
expected = ctc_vol_reference(close, window=window, trading_days=252.0)
valid = ~np.isnan(expected)
assert np.allclose(actual[valid], expected[valid], atol=1e-10), (
"close_to_close_vol does not match reference formula"
)
def test_constant_returns_known_vol(self):
"""Constant daily log-return of 0.01 → C2C vol = 0.01 * sqrt(252) ≈ 0.1587."""
from ferro_ta.analysis.options import close_to_close_vol
# Build a price series with constant daily log-return of 0.01
n = 100
constant_log_ret = 0.01
close = 100.0 * np.exp(np.arange(n) * constant_log_ret)
window = 21
out = close_to_close_vol(close, window=window, trading_days_per_year=252.0)
# Expected: sqrt(0.01^2 * 252) = 0.01 * sqrt(252)
expected_vol = constant_log_ret * np.sqrt(252.0)
valid = ~np.isnan(out)
assert np.all(valid[window:]), "Expected valid values after warmup"
assert out[window] == pytest.approx(expected_vol, rel=1e-10), (
f"Constant-return vol: got {out[window]}, expected {expected_vol}"
)
def test_parkinson_lognormal_unbiased(self):
"""Parkinson estimator within 50% of true vol=0.20 for simulated OHLC data.
Parkinson uses the log(high/low) range as a proxy for daily realized
vol. The estimator is unbiased for a Brownian-motion diffusion where
the daily range follows a known distribution, but a simplified
simulation (single end-of-day price + independent range draw) will
underestimate the range. We therefore build a proper multi-step
intraday path so the high/low reflects the true diffusion range,
and use a lenient 50% tolerance to accommodate finite-sample noise.
"""
from ferro_ta.analysis.options import parkinson_vol
rng = np.random.default_rng(123)
true_vol = 0.20
n_days = 500
steps_per_day = 50 # intraday steps to get a realistic H-L range
daily_sigma = true_vol / np.sqrt(252.0)
step_sigma = daily_sigma / np.sqrt(steps_per_day)
# Simulate intraday paths, extract open/high/low/close each day
highs = np.empty(n_days)
lows = np.empty(n_days)
price = 100.0
for i in range(n_days):
intraday = price * np.exp(
np.cumsum(rng.normal(0.0, step_sigma, steps_per_day))
)
path = np.concatenate([[price], intraday])
highs[i] = path.max()
lows[i] = path.min()
price = intraday[-1]
window = 21
out = parkinson_vol(highs, lows, window=window, trading_days_per_year=252.0)
valid = out[~np.isnan(out)]
assert len(valid) > 0, "No valid Parkinson estimates"
median_est = float(np.median(valid))
assert abs(median_est - true_vol) < 0.50 * true_vol, (
f"Parkinson estimate {median_est:.4f} is more than 50% from true vol {true_vol}"
)
def test_vol_estimators_all_positive_finite(self):
"""All 5 estimators produce finite and positive non-NaN values on random OHLC."""
from ferro_ta.analysis.options import (
close_to_close_vol,
garman_klass_vol,
parkinson_vol,
rogers_satchell_vol,
yang_zhang_vol,
)
rng = np.random.default_rng(99)
n = 200
log_ret = rng.normal(0.0, 0.01, n)
close = 100.0 * np.cumprod(np.exp(log_ret))
high = close * np.exp(np.abs(rng.normal(0.0, 0.005, n)))
low = close * np.exp(-np.abs(rng.normal(0.0, 0.005, n)))
open_ = np.roll(close, 1)
open_[0] = close[0]
window = 20
estimators = {
"close_to_close": close_to_close_vol(close, window=window),
"parkinson": parkinson_vol(high, low, window=window),
"garman_klass": garman_klass_vol(open_, high, low, close, window=window),
"rogers_satchell": rogers_satchell_vol(
open_, high, low, close, window=window
),
"yang_zhang": yang_zhang_vol(open_, high, low, close, window=window),
}
for name, out in estimators.items():
valid = out[~np.isnan(out)]
assert len(valid) > 0, f"{name}: no valid (non-NaN) estimates"
assert np.all(np.isfinite(valid)), f"{name}: non-finite values present"
assert np.all(valid > 0.0), f"{name}: non-positive values present"
class TestVolConeAccuracy:
"""Tests for vol_cone — no scipy required."""
def test_cone_windows_match_requested(self):
"""Output windows should match the input list exactly."""
from ferro_ta.analysis.options import vol_cone
rng = np.random.default_rng(0)
close = 100.0 * np.cumprod(np.exp(rng.normal(0.0, 0.01, 500)))
requested = (10, 21, 42)
cone = vol_cone(close, windows=requested)
assert list(cone.windows.astype(int)) == list(requested), (
f"Cone windows {list(cone.windows)} do not match requested {list(requested)}"
)
def test_cone_median_matches_rolling_median(self):
"""Manually computed rolling C2C vol median for window=21 should match cone.median[0]."""
from ferro_ta.analysis.options import close_to_close_vol, vol_cone
rng = np.random.default_rng(5)
close = 100.0 * np.cumprod(np.exp(rng.normal(0.0, 0.01, 500)))
window = 21
cone = vol_cone(close, windows=(window,))
rolling = close_to_close_vol(close, window=window, trading_days_per_year=252.0)
valid = rolling[~np.isnan(rolling)]
manual_median = float(np.median(valid))
assert cone.median[0] == pytest.approx(manual_median, rel=1e-6), (
f"vol_cone median {cone.median[0]:.6f} does not match manual median {manual_median:.6f}"
)
class TestStrategyAnalyticsAccuracy:
@pytest.fixture(autouse=True)
def require_scipy(self):
pytest.importorskip("scipy")
def test_put_call_parity_deviation_analytical(self):
"""BSM call/put from scipy formulas fed into put_call_parity_deviation → < 1e-8."""
from ferro_ta.analysis.options import put_call_parity_deviation
S, K, r, q, T, sigma = 100.0, 100.0, 0.05, 0.02, 1.0, 0.20
call = bsm_call(S, K, r, q, T, sigma)
put = bsm_put(S, K, r, q, T, sigma)
dev = put_call_parity_deviation(call, put, S, K, r, T, carry=q)
assert abs(dev) < 1e-8, (
f"put_call_parity_deviation for BSM-consistent prices: got {dev}, expected ~0"
)
def test_expected_move_known_value(self):
"""S=100, iv=0.20, days=30, trading_days=252 → upper move ≈ 7.14."""
from ferro_ta.analysis.options import expected_move
S, iv, days, td = 100.0, 0.20, 30.0, 252.0
lower, upper = expected_move(S, iv, days, td)
# log-normal formula: S * (exp(sigma * sqrt(days/trading_days)) - 1)
expected_upper = S * (np.exp(iv * np.sqrt(days / td)) - 1.0)
expected_lower = S * (np.exp(-iv * np.sqrt(days / td)) - 1.0)
assert upper == pytest.approx(expected_upper, rel=1e-6), (
f"expected_move upper: got {upper:.4f}, expected {expected_upper:.4f}"
)
assert lower == pytest.approx(expected_lower, rel=1e-6), (
f"expected_move lower: got {lower:.4f}, expected {expected_lower:.4f}"
)
# Numeric check: upper ≈ 7.14
assert upper == pytest.approx(7.14, abs=0.05), (
f"expected_move upper should be ~7.14, got {upper:.4f}"
)