lets-be-precise 0.1.3

Numerically precise Black-Scholes pricing and implied-volatility primitives.

lets-be-precise

Numerically precise Black-Scholes pricing and implied-volatility primitives.

lbp.lbp is a compact compiled kernel whose forward (price) and inverse (implied volatility) paths target the lower edge of double-precision error across the whole (sigma, log(K/F)) plane -- including the deep-OTM and small-t corners where the standard Jaeckel and Cody-erfc paths shed several digits. The wheel also ships lbp.lbr, a thin Cython binding for Peter Jaeckel's "Let's Be Rational" 2013-2023 reference, exposed for honest side-by-side speed and accuracy comparison.

Install

pip install lets-be-precise

Optional extras:

• lets-be-precise[demo] -- pulls in mpmath so python -m lbp.demos.* comparison modules can run their high-precision truth side.
• lets-be-precise[dev] -- adds pytest and scipy for the (separately hosted) test suite.

Quick start

import math
import lbp                                # compiled BSM kernel
import lbp.lbr                            # Jaeckel LBR reference (same flags)

sigma, logk = 0.20, 0.15                  # logk = ln(K/F)
mlogk = -logk                             # lbp's k_or_mlogk on the mlogk branch
                                          # (matches LBR's x = ln(F/K))

C    = lbp.black(sigma, mlogk)            # call price under our BSM kernel
sig1 = lbp.iv(C, mlogk)                   # implied vol, lbp's Householder solver

# LBR uses Jaeckel's normalised price b = C / sqrt(F K) and the reflected
# moneyness x = -logk = ln(F/K); see lbp.lbr's docstrings for full detail.
beta = C * math.exp(-logk / 2.0)          # C -> Jaeckel's normalised price
sig2 = lbp.lbr.iv(beta, logk)             # same answer to ~eps

assert abs(sig1 - sig2) < 1e-13

The two implied-vol entry points are deliberately the same shape so the algorithmic differences are timing-clean for benchmarks; see src/lbp/lbr.pyx for the conversion table.

Numerical reliability

The compiled BSM kernel is built on top of mathpf's Mills-ratio primitives, including the cancellation-free symmetric divided difference R_DD(x, dx, +1) = (R(x - dx) - R(x + dx)) / (2 dx). This is the algorithmic counterpart to Jaeckel's two single-purpose expansions inside normalised_black_call_over_vega. They satisfy

b/vega = R(-h-t) - R(-h+t) = 2 t * R_DD(-h, t, +1)

so each Jaeckel expansion matches one R_DD branch. The headline difference:

region	what Jaeckel uses	R_DD branch	accuracy
h < -10, h+t < -9.79 (deep OTM)	17th-order asymptotic series in q	CF asymp ladder (n=2/4/6/8) or Taylor	both ~eps; AEXP a touch cleaner at the Taylor corner
-10 <= h <= 0, t < 0.21 (small t)	12th-order Taylor in t	5-term Taylor seeded by R'''(x)	mathpf wins by up to ~7 bits near h = -10

The small-t gap comes from Jaeckel's coefficient a := 1 + h * Y(h) (with Y(h) := Phi(h)/phi(h)), which suffers cancellation because Y(h) -> -1/h as h -> -infinity. At h = -10, |a| ~ 1/h^2 ~ 0.01 -- two decimal digits are killed before the series even runs. mathpf's R_DD Taylor branch is seeded from R'''(x) and descends to R'(x) via the cancellation-free (r_d3 + 1) / (x^2 + 3) relation, so the offending coefficient never appears.

Worst single cell, inside Jaeckel's actual STEXP dispatch region:

(h, t) = (-9.4, 0.17)            mpmath truth (80 dps)  =  3.72524275689953624e-03
LBR small_t_expansion_of_normalised_black_call_over_vega = 3.72524275689937664e-03
mathpf 2 t * R_DD(-h, t, +1)                             = 3.72524275689953624e-03

LBR    relerr : 4.28e-14   (~192 eps, ~7 bits lost)
mathpf relerr : 0          (bit-exact)

Reproduce on your machine:

pip install lets-be-precise[demo]
python -m lbp.demos.lbr_vs_mathpf

The demo prints both the AEXP-in-gate and STEXP-in-gate tables and explains the dispatch geometry. lbp.lbr also exposes Jaeckel's two file-scope expansions directly as lbp.lbr.aexp_over_vega(h, t) and lbp.lbr.stexp_over_vega(h, t) so users can probe further cells of their own choosing.

License

Proprietary -- see LICENSE. The package is distributed as binary wheels only via PyPI; the underlying Cython / C++ source is not included in the wheel.