"""Checks that a pricer should pass, grouped by the kind of knowledge they use.

Each check takes a pricer module and returns None when it passes, or a short
description of the first failing input. A crash counts as a failure.
"""
import random
from collections import namedtuple
from math import exp, sqrt

Check = namedtuple("Check", "id title family run")

TOLERANCE = 1e-7


def scenarios(count=300, seed=20261030, zero_dividend=False):
    """Random inputs that reach short maturities and far strikes."""
    rng = random.Random(seed)
    for _ in range(count):
        spot = 100.0 * exp(rng.uniform(-1.5, 1.5))
        strike = 100.0 * exp(rng.uniform(-1.5, 1.5))
        yield (spot, strike, rng.uniform(0.0, 0.15), 0.0 if zero_dividend else rng.uniform(0.0, 0.08),
               exp(rng.uniform(-3.0, 0.3)), exp(rng.uniform(-4.5, 1.0)))


def label(inputs):
    names = ("S", "K", "r", "q", "vol", "T")
    return ", ".join("{}={:.4g}".format(name, value) for name, value in zip(names, inputs))


def textbook(p):
    """One reference value, the way many test suites start."""
    value = p.call(100.0, 100.0, 0.05, 0.0, 0.2, 1.0)
    if abs(value - 10.4506) > 5e-5:
        return "call(100, 100, 5%, 0, 20%, 1y) gave {:.4f}, expected 10.4506".format(value)


def parity(p):
    for inputs in scenarios():
        spot, strike, rate, div, _, time = inputs
        difference = p.call(*inputs) - p.put(*inputs)
        expected = spot * exp(-div * time) - strike * exp(-rate * time)
        if abs(difference - expected) > TOLERANCE * max(1.0, spot, strike):
            return "{}: C-P={:.5f}, expected {:.5f}".format(label(inputs), difference, expected)


def bounds(p):
    for inputs in scenarios():
        spot, strike, rate, div, _, time = inputs
        value = p.call(*inputs)
        low = max(spot * exp(-div * time) - strike * exp(-rate * time), 0.0)
        high = spot * exp(-div * time)
        if value < low - 1e-8 or value > high + 1e-8:
            return "{}: call={:.5f} outside [{:.5f}, {:.5f}]".format(label(inputs), value, low, high)


def strike_shape(p):
    """Call prices fall with strike and are convex in it. Otherwise a density is negative."""
    for inputs in scenarios():
        spot, strike, rate, div, vol, time = inputs
        width = 0.02 * strike
        below = p.call(spot, strike - width, rate, div, vol, time)
        middle = p.call(*inputs)
        above = p.call(spot, strike + width, rate, div, vol, time)
        if above > below + 1e-8:
            return "{}: call rises with strike".format(label(inputs))
        if below - 2.0 * middle + above < -1e-8:
            return "{}: butterfly={:.6f}".format(label(inputs), below - 2.0 * middle + above)


def homogeneity(p):
    """Doubling spot and strike doubles the price."""
    for inputs in scenarios():
        spot, strike, rate, div, vol, time = inputs
        scaled = p.call(2.0 * spot, 2.0 * strike, rate, div, vol, time)
        if abs(scaled - 2.0 * p.call(*inputs)) > TOLERANCE * max(1.0, spot, strike):
            return "{}: not homogeneous".format(label(inputs))


def delta_bump(p):
    for inputs in scenarios():
        spot, strike, rate, div, vol, time = inputs
        bump = 1e-3 * spot * vol * sqrt(time)
        numeric = (p.call(spot + bump, strike, rate, div, vol, time) - p.call(spot - bump, strike, rate, div, vol, time)) / (2.0 * bump)
        analytic = p.delta(*inputs)
        if abs(numeric - analytic) > 1e-5:
            return "{}: delta {:.5f}, bump {:.5f}".format(label(inputs), analytic, numeric)


def gamma_bump(p):
    for inputs in scenarios():
        spot, strike, rate, div, vol, time = inputs
        bump = 1e-3 * spot * vol * sqrt(time)
        numeric = (p.delta(spot + bump, strike, rate, div, vol, time) - p.delta(spot - bump, strike, rate, div, vol, time)) / (2.0 * bump)
        analytic = p.gamma(*inputs)
        if abs(numeric - analytic) > 1e-4 * abs(analytic) + 1e-7:
            return "{}: gamma {:.6f}, bump of delta {:.6f}".format(label(inputs), analytic, numeric)


def vega_bump(p):
    for inputs in scenarios():
        spot, strike, rate, div, vol, time = inputs
        bump = 1e-4
        numeric = (p.call(spot, strike, rate, div, vol + bump, time) - p.call(spot, strike, rate, div, vol - bump, time)) / (2.0 * bump)
        analytic = p.vega(*inputs)
        if abs(numeric - analytic) > 1e-4 * max(1.0, abs(analytic)):
            return "{}: vega {:.5f}, bump {:.5f}".format(label(inputs), analytic, numeric)


def implied_round_trip(p):
    for inputs in scenarios():
        spot, strike, rate, div, vol, time = inputs
        vega_value = p.vega(*inputs)
        if vega_value < 1e-3 or not 1e-3 < vol < 4.0:
            continue
        target = p.call(*inputs)
        recovered = p.implied_vol(target, spot, strike, rate, div, time)
        if recovered is None or abs(recovered - vol) > 1e-5:
            return "{}: recovered vol {}".format(label(inputs), recovered)


def monte_carlo(p):
    """An independent method. The simulation does not share code with the formula."""
    rng = random.Random(7)
    paths = 40000
    for inputs in list(scenarios(count=6, seed=11)):
        spot, strike, rate, div, vol, time = inputs
        total, squares = 0.0, 0.0
        for _ in range(paths // 2):
            draw = rng.gauss(0.0, 1.0)
            payoffs = 0.0
            for sign in (1.0, -1.0):
                final = spot * exp((rate - div - 0.5 * vol * vol) * time + vol * sqrt(time) * sign * draw)
                payoffs += max(final - strike, 0.0)
            value = 0.5 * exp(-rate * time) * payoffs
            total += value
            squares += value * value
        count = paths // 2
        mean = total / count
        error = sqrt(max(squares / count - mean * mean, 0.0) / count)
        price = p.call(*inputs)
        if abs(price - mean) > 4.0 * error + 1e-3:
            return "{}: formula {:.4f}, simulation {:.4f} ± {:.4f}".format(label(inputs), price, mean, error)


def limits(p):
    """Edges of the input space, where a formula divides by something small."""
    spot, strike, rate, div = 100.0, 90.0, 0.05, 0.02
    for vol, time in ((0.0, 1.0), (1e-14, 1.0), (0.2, 0.0)):
        value = p.call(spot, strike, rate, div, vol, time)
        expected = max(spot * exp(-div * time) - strike * exp(-rate * time), 0.0)
        if abs(value - expected) > 1e-8:
            return "vol={}, T={}: call={:.5f}, expected {:.5f}".format(vol, time, value, expected)
    far = p.call(spot, strike, rate, div, 0.2, 1e-10)
    if abs(far - (spot - strike)) > 1e-3:
        return "near expiry: call={:.5f}, expected payoff {:.5f}".format(far, spot - strike)


CHECKS = [
    Check("textbook", "One reference price", "Example", textbook),
    Check("parity", "Put-call parity", "No-arbitrage", parity),
    Check("bounds", "Price bounds", "No-arbitrage", bounds),
    Check("strike-shape", "Falling and convex in strike", "No-arbitrage", strike_shape),
    Check("homogeneity", "Doubling spot and strike", "No-arbitrage", homogeneity),
    Check("delta-bump", "Delta against a bump", "Self-consistency", delta_bump),
    Check("gamma-bump", "Gamma against a bump of delta", "Self-consistency", gamma_bump),
    Check("vega-bump", "Vega against a bump", "Self-consistency", vega_bump),
    Check("implied-round-trip", "Implied volatility round trip", "Self-consistency", implied_round_trip),
    Check("monte-carlo", "Monte Carlo", "Independent method", monte_carlo),
    Check("limits", "Zero volatility and expiry", "Limits", limits),
]