Implied volatility: Newton-Raphson, Brent, py_vollib and QuantLib.

From price to volatility

A pricing model such as Black-Scholes-Merton or Black-76 (lesson 3.4.2) takes five inputs and returns a price: forward (or spot and dividends), strike, time to expiry, interest rate and volatility. The first four are observed; the volatility is not. The market solves the problem the other way around: it quotes the price and lets everyone deduce the volatility that justifies it. That is the implied volatility: the σ that, plugged into the model, reproduces the market price.

For a European option on an index, with the forward $F$ at expiry and the discount factor $D = e^{-rT}$, Black-76 gives

and the implied volatility is the root of $f(\sigma) = C_{\text{model}}(\sigma) - C_{\text{market}} = 0$. There is no closed formula to solve for it: it has to be found with a numerical method. Two properties guarantee that the search makes sense:

• Monotonicity. The derivative of the price with respect to σ is the vega (lesson 3.4.3), $\mathcal{V} = D\,F\,\varphi(d_1)\sqrt{T} > 0$. The price increases strictly with σ, so there is at most one root.
• Bounds. As σ tends to 0, the price tends to the discounted intrinsic value, $D\max(F - K, 0)$ for the call; as σ grows without limit, to $D\,F$. A price outside that interval violates the arbitrage bounds of lesson 3.4.1 and has no implied volatility. Between the two, it exists and is unique.

Conventions matter as much as the method. In this lesson: prices in index points per option (the contract is worth 100 USD per point, but that does not change the volatility), time to expiry in calendar days / 365 because ORX options expire at the close, continuously compounded rate as a decimal and annual volatility as a decimal (0.25 = 25%). With the forward, the dividend is not needed: it is already inside $F$.

flowchart TB
  A[Chain bid and ask] --> B[Forward and rate]
  B --> C{Filters}
  C -- passes --> D[Mid price]
  C -- fails --> X[Discard]
  D --> E{Within bounds}
  E -- yes --> F[Root in sigma]
  E -- no --> X
  F --> G[Smile]

The order of the diagram is that of the lesson: first the forward, then the numerical methods and, at the end, the filters, which is where it is decided whether the number that comes out means anything.

One day's chain and its forward

The data are the daily option chain on ORX, a fictitious index that mimics the S&P 500, with European options that mimic SPX options (opciones3_cadena.parquet). They are synthetic and carry the known answer: each row has its iv_modelo, the true noise-free implied volatility, and the bids and asks were generated around the price for that volatility with a realistic spread and a bit of noise. It is therefore possible to check how wrong each step is. The chosen day is June 4, 2025, with the July 18 expiry:

import math
import numpy as np
import pandas as pd

cad = pd.read_parquet("opciones3_cadena.parquet")
ind = pd.read_csv("opciones3_indice.csv", parse_dates=["fecha"])
ver = pd.read_parquet("opciones3_verdad.parquet")

hoy, vto = pd.Timestamp("2025-06-04"), pd.Timestamp("2025-07-18")
dia = cad[(cad["fecha"] == hoy) & (cad["vencimiento"] == vto)]
S = ind.loc[ind["fecha"] == hoy, "cierre"].item()
dias = (vto - hoy).days
T = dias / 365            # calendar days: expires at the close
print(f"ORX {S:.2f} · {dias} days · T = {T:.4f}")
print(dia["tipo"].value_counts().to_dict())
cols = ["tipo", "strike", "bid", "ask", "iv_modelo"]
dia[dia["strike"].isin([4000, 5150, 6000])][cols]

ORX 5137.51 · 44 days · T = 0.1205
{'C': 45, 'P': 45}
       tipo  strike     bid     ask  iv_modelo
516554    C  4000.0  1153.9  1164.2     0.4130
516570    C  5150.0   160.2   161.7     0.2231
516583    C  6000.0     0.5     0.6     0.1709
516599    P  4000.0     9.8    10.1     0.4130
516615    P  5150.0   155.2   156.7     0.2231
516628    P  6000.0   838.2   844.7     0.1709

There is one call and one put per strike, from 2,400 to 7,200 points. The 4,000 one is a put 22% below the index; the 5,150 one, almost at the money; the 6,000 one, a call 17% above. The iv_modelo column is the answer that has to be recovered.

The forward is missing. Instead of computing it with the dividend (which in reality is estimated), it is deduced from the prices themselves with the put-call parity of lesson 3.4.1, $C - P = D\,(F - K)$: a line in $K$ with slope $-D$ and intercept $D\,F$. With the strikes within 5% of the index:

medio = dia.assign(medio=(dia["bid"] + dia["ask"]) / 2)
tabla = medio.pivot(index="strike", columns="tipo", values="medio")
cerca = tabla[np.abs(tabla.index / S - 1) < 0.05]
pend, corte = np.polyfit(cerca.index, cerca["C"] - cerca["P"], 1)
D = -pend
F = corte / D
r = -math.log(D) / T
v = ver[(ver["fecha"] == hoy) & (ver["vencimiento"] == vto)].iloc[0]
print(f"{len(cerca)} strikes")
print(f"Parity: F = {F:.2f}, r = {r:.3%}")
print(f"Truth:  F = {v.forward:.2f}, r = {v.tipo_libre:.3%}")

10 strikes
Parity: F = 5155.15, r = 4.214%
Truth:  F = 5155.25, r = 4.203%

With ten strikes, parity gives a forward of 5,155.15 points, a tenth away from the true one (5,155.25), and a rate of 4.214% against 4.203%. The difference is the noise in the quotes. The forward is above the spot (5,137.51) because the rate exceeds the dividend yield. It is the same idea Cboe uses for the VIX (lesson 3.4.6): the options market says what its forward is, and that forward already carries the expected dividend, which thus does not need to be estimated.

Newton-Raphson: fast and fragile

Newton-Raphson takes advantage of the known derivative: at each step it replaces the curve by its tangent and jumps to the point where the tangent crosses the target price,

Near the root it converges quadratically: the number of correct decimals doubles at each iteration. First, the Black-76 price and vega vectorized with numpy:

from scipy.special import ndtr


def black(F, K, T, D, s, tipo):
    """Black-76 price (per index point)."""
    sd = s * np.sqrt(T)
    d1 = np.log(F / K) / sd + sd / 2
    call = D * (F * ndtr(d1) - K * ndtr(d1 - sd))
    return np.where(tipo == "C", call, call - D * (F - K))


def vega(F, K, T, D, s):
    """Derivative of the price with respect to σ (per unit of σ)."""
    sd = s * np.sqrt(T)
    d1 = np.log(F / K) / sd + sd / 2
    return D * F * np.exp(-d1 * d1 / 2) / math.sqrt(2 * math.pi) \
        * np.sqrt(T)


fila = dia[(dia["strike"] == 5150) & (dia["tipo"] == "C")].iloc[0]
p = black(F, 5150.0, T, D, fila.iv_modelo, "C")
print(f"Call 5150 with iv_modelo: {p:.2f}"
      f" (bid {fila.bid}, ask {fila.ask})")
print(f"Vega: {vega(F, 5150.0, T, D, fila.iv_modelo) / 100:.2f}"
      " price points per volatility point")

Call 5150 with iv_modelo: 160.96 (bid 160.2, ask 161.7)
Vega: 7.09 price points per volatility point

With the true volatility, the Black-76 price falls inside the spread. The vega says how much a volatility point is worth: 7.09 price points, 709 USD per contract. Read backward, it is the conversion that matters in this lesson: half a spread of 0.75 points is equivalent to about ±0.11 volatility points.

def newton(precio, K, tipo, s0, tol=1e-10, maxit=50):
    s, camino = s0, [s0]
    for _ in range(maxit):
        dif = float(black(F, K, T, D, s, tipo)) - precio
        if abs(dif) < tol:
            break
        s = s - dif / float(vega(F, K, T, D, s))
        camino.append(s)
        if not 0 < s < 10:          # outside reasonable values
            break
    return s, camino


def medio_de(K, tipo):
    f = dia[(dia["strike"] == K) & (dia["tipo"] == tipo)].iloc[0]
    return float(f.bid + f.ask) / 2


for K, tipo in [(5150.0, "C"), (4000.0, "P")]:
    s, camino = newton(medio_de(K, tipo), K, tipo, s0=0.20)
    print(K, tipo, [round(x, 4) for x in camino])

5150.0 C [0.2, 0.2231, 0.2231]
4000.0 P [0.2, 12.7831]

The near-the-money call converges in two steps. The 4,000 put jumps in the first one to σ = 12.78, a 1,278% annual volatility. The reason lies in the tangent: with σ = 20% and 44 days, the 4,000 strike is 3.65 standard deviations from the forward, the model price is 0.0099 points and the vega, 0.79 per unit of σ. The step is (0.0099 − 9.95) / 0.79 ≈ −12.6, and Newton adds 12.6 to σ. In the wings of the smile, a fixed starting point is a bad idea.

def arranque_mk(K):
    """Inflection point of the price in σ: sqrt(2·|ln(F/K)|/T)."""
    return math.sqrt(2 * abs(math.log(F / K)) / T)


for K, tipo in [(4000.0, "P"), (3400.0, "P"), (6000.0, "C")]:
    s0 = arranque_mk(K)
    s, camino = newton(medio_de(K, tipo), K, tipo, s0)
    print(f"{K:.0f} {tipo}: start {s0:.3f}, "
          f"{len(camino) - 1} iterations, σ = {s:.6f}")

4000 P: start 2.052, 7 iterations, σ = 0.412856
3400 P: start 2.628, 8 iterations, σ = 0.502269
6000 C: start 1.587, 9 iterations, σ = 0.170921

Manaster and Koehler (1982) proposed starting at $\sigma^* = \sqrt{2\,|\ln(F/K)|/T}$, the point where the vega is maximal (there $d_1 d_2 = 0$) and the price curve in σ changes from convex to concave. If the root is below, the function is convex between the two and the tangent always stays below the curve: each step approaches the root without overshooting. If it is above, it is concave and the same happens from the other side. With that starting point, Newton reaches $10^{-10}$ in the price in 7 to 9 iterations for the three options, the 4,000 put included.

The figure draws the two starting points on the 4,000 put, with the price on a log scale. From 2.05, the steps go down the curve without overshooting the root. The tangent at 0.20 is so flat that, over the whole plotted range, it does not reach 1.6 price points.

Brent: slow and safe

If the derivative is treacherous, it can be dispensed with. Bisection starts from an interval $[\sigma_a, \sigma_b]$ over which the function changes sign and halves it at each step: it never fails, but it gains only one bit per iteration. Brent's method combines bisection with interpolation (secant and inverse quadratic interpolation): when the interpolation looks promising, it uses it; when not, it bisects. It converges almost as fast as the secant and keeps the guarantee of bisection. In SciPy it is brentq:

from scipy.optimize import brentq


def iv_brent(precio, K, tipo, a=1e-4, b=5.0):
    f = lambda s: float(black(F, K, T, D, s, tipo)) - precio
    return brentq(f, a, b, xtol=1e-12, full_output=True)


for K, tipo in [(5150.0, "C"), (4000.0, "P"), (6000.0, "C")]:
    s, res = iv_brent(medio_de(K, tipo), K, tipo)
    print(f"{K:.0f} {tipo}: σ = {s:.6f}, {res.iterations} iterations,"
          f" {res.function_calls} evaluations")

5150 C: σ = 0.223091, 6 iterations, 7 evaluations
4000 P: σ = 0.412856, 15 iterations, 16 evaluations
6000 C: σ = 0.170921, 19 iterations, 20 evaluations

Brent reaches the same σ with 6, 15 and 19 iterations, one price evaluation at each (plus the initial one). Newton from Manaster-Koehler needed 7 to 9, with price and vega at each step. For one option the difference does not matter; for millions of rows, it does, and then a vectorized method (at the end of the lesson) or Jäckel's is used. The interval $[10^{-4}, 5]$ is the only assumption: if the price is outside the bounds, brentq finds no sign change and raises ValueError, which is exactly what is wanted.

Libraries: vollib, py_vollib_vectorized and QuantLib

vollib (formerly py_vollib)

For years, py_vollib was the reference Python library for Black, Black-Scholes and Black-Scholes-Merton with Jäckel's (2015) Let's Be Rational algorithm, which obtains the volatility to machine precision in two iterations thanks to a rational approximation of the starting point. Its installed version, 1.0.12, no longer contains code: it is a transition package that depends on vollib (the new name) and warns when imported:

import warnings

with warnings.catch_warnings(record=True) as avisos:
    warnings.simplefilter("always")
    import py_vollib.black.implied_volatility as pv_iv
print(avisos[0].category.__name__, "·", avisos[0].message)

from vollib.black.implied_volatility import implied_volatility
print(pv_iv.implied_volatility is implied_volatility)

K, tipo = 4000.0, "P"
iv = implied_volatility(medio_de(K, tipo), F, K, r, T, tipo.lower())
print(f"vollib, put 4000: σ = {iv:.6f}")

DeprecationWarning · py_vollib is deprecated and will be removed in a future release; please import from vollib instead.
True
vollib, put 4000: σ = 0.412856

The warning appears once, on import. py_vollib.black.implied_volatility is the same object as in vollib: the old package only redirects. New code should import from vollib, as in the rest of the lesson. The function receives the discounted price, F, K, r, t and 'c' or 'p', and gives the same σ as Newton and Brent: 0.412856.

py_vollib_vectorized

py_vollib_vectorized wrapped the py_vollib functions with numba to compute whole chains at once. Its latest version, 0.1.1, is from February 2021, and it does not compile with the current versions of numba (0.68) and of py_lets_be_rational (1.1.2, from May 2026): the first call fails when typing the functions. This is checked in a separate process, so as not to leave numba loaded in this one:

import subprocess
import sys

prueba = """
import numpy as np
import py_vollib_vectorized as pvv
pvv.vectorized_implied_volatility_black(
    np.array([3.9]), np.array([100.0]), np.array([100.0]),
    np.array([0.03]), np.array([0.25]), np.array(["c"]))
"""
res = subprocess.run([sys.executable, "-c", prueba],
                     capture_output=True, text=True)
error = [l for l in res.stderr.splitlines() if "Error" in l][-1]
print("Exit code:", res.returncode)
print(error)

Exit code: 1
numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)

The process ends with code 1 and a numba TypingError. Reading the installed code shows why: the compiled function of py_vollib_vectorized calls normalised_black_call from py_lets_be_rational, which in version 1.1.2 is a normal Python function (there, numba compilation became optional, with the environment variable PY_LETS_BE_RATIONAL_ENABLE_NUMBA), and numba cannot type an uncompiled function inside a compiled one. Setting that variable is not enough either: compilation fails further on.

import os
import numba

os.environ["NUMBA_DISABLE_JIT"] = "1"   # numba as plain Python
numba.config.reload_config()            # it was already imported
import py_vollib_vectorized as pvv

otm = dia[np.where(dia["tipo"] == "P", dia["strike"] < F,
                   dia["strike"] >= F) & (dia["bid"] > 0)].copy()
otm["medio"] = (otm["bid"] + otm["ask"]) / 2
n = len(otm)
iv_vec = pvv.vectorized_implied_volatility_black(
    otm["medio"].to_numpy(), np.full(n, F), otm["strike"].to_numpy(),
    np.full(n, r), np.full(n, T), otm["tipo"].str.lower().to_numpy(),
    return_as="numpy")
print(n, "out-of-the-money options with bid > 0")
print(np.round(iv_vec[:4], 6))

40 out-of-the-money options with bid > 0
[0.214444 0.205999 0.19825  0.190714]

With NUMBA_DISABLE_JIT=1, numba leaves the decorated functions as plain Python and the library works, although without the speed that was its reason for being. Two details of the code. The variable is read when numba is imported, and here numba was already loaded (vollib imports it through py_lets_be_rational): hence reload_config(). And on import, py_vollib_vectorized replaces the vollib functions with its vectorized versions; implied_volatility was imported earlier and is still the original. An unmaintained dependency can break with any update of the environment; for large chains, your own vectorized Newton (below) or QuantLib are safer bets.

QuantLib

QuantLib offers two paths. The direct one is blackFormulaImpliedStdDev, which inverts the Black formula with the forward and the discount and returns the total standard deviation $\sigma\sqrt{T}$. The complete one builds the instrument (EuropeanOption) on a Black-Scholes-Merton process with rate and dividend curves, and calls impliedVolatility; it is the one used when the option is part of a portfolio valued with QuantLib. Its default values are a tolerance of $10^{-4}$ and σ between $10^{-7}$ and 411Note 1 · Signature of VanillaOption::impliedVolatility in QuantLib 1.44, ql/instruments/vanillaoption.hpp: https://github.com/lballabio/QuantLib/blob/v1.44/ql/instruments/vanillaoption.hpp; here it is given $10^{-10}$, 200 evaluations and σ between $10^{-4}$ and 5:

import QuantLib as ql

K, precio = 4000.0, medio_de(4000.0, "P")
sd = ql.blackFormulaImpliedStdDev(ql.Option.Put, K, F, precio, D)
print(f"blackFormulaImpliedStdDev: σ = {sd / math.sqrt(T):.6f}")
# displacement 0, start 0.3, tolerance 1e-12, 200 iterations
sd = ql.blackFormulaImpliedStdDev(ql.Option.Put, K, F, precio, D,
                                  0.0, 0.3, 1e-12, 200)
print(f"... with tolerance 1e-12: σ = {sd / math.sqrt(T):.6f}")

ql.Settings.instance().evaluationDate = ql.Date(4, 6, 2025)
dc = ql.Actual365Fixed()
hoy_ql = ql.Settings.instance().evaluationDate
q = r - math.log(F / S) / T     # dividend implied by the forward
curva = lambda x: ql.YieldTermStructureHandle(
    ql.FlatForward(hoy_ql, x, dc, ql.Continuous))
vol = ql.BlackVolTermStructureHandle(
    ql.BlackConstantVol(hoy_ql, ql.TARGET(), 0.20, dc))
proc = ql.BlackScholesMertonProcess(
    ql.QuoteHandle(ql.SimpleQuote(S)), curva(q), curva(r), vol)
put = ql.EuropeanOption(ql.PlainVanillaPayoff(ql.Option.Put, K),
                        ql.EuropeanExercise(ql.Date(18, 7, 2025)))
put.setPricingEngine(ql.AnalyticEuropeanEngine(proc))
iv_ql = put.impliedVolatility(precio, proc, 1e-10, 200, 1e-4, 5.0)
print(f"EuropeanOption.impliedVolatility: σ = {iv_ql:.6f}")
print(f"Implied dividend yield: {q:.3%}")

blackFormulaImpliedStdDev: σ = 0.412854
... with tolerance 1e-12: σ = 0.412856
EuropeanOption.impliedVolatility: σ = 0.412856
Implied dividend yield: 1.370%

Three details. blackFormulaImpliedStdDev returns σ√T, not σ. Its default tolerance, $10^{-6}$ (visible in QuantLib's blackformula.hpp header), gets the sixth decimal wrong: 0.412854; with $10^{-12}$ it agrees with the others. And the complete instrument needs a process with spot, dividend and rate; the constant volatility of 0.20 it is built with has no influence, because impliedVolatility replaces it while searching. The dividend implied by the parity forward, 1.370%, is close to the 1.3435% of the index file: the difference is the noise with which forward and rate were estimated.

The five methods, option by option

On the out-of-the-money options of the day, the five methods should give the same result. The comparison serves as proof that our own code is correct:

filas = []
for f in otm.itertuples():
    K, tipo, pm = f.strike, f.tipo, f.medio
    filas.append({
        "newton": newton(pm, K, tipo, arranque_mk(K))[0],
        "brent": iv_brent(pm, K, tipo)[0],
        "vollib": implied_volatility(pm, F, K, r, T, tipo.lower()),
        "quantlib": ql.blackFormulaImpliedStdDev(
            ql.Option.Call if tipo == "C" else ql.Option.Put,
            K, F, pm, D, 0.0, 0.3, 1e-12, 200) / math.sqrt(T)})
ivs = pd.DataFrame(filas, index=otm.index)
ivs["vectorized"] = iv_vec
dif = ivs.sub(ivs["brent"], axis=0).abs().max()
print("Maximum difference from Brent:")
print(dif.map(lambda x: f"{x:.1e}").to_string())

Maximum difference from Brent:
newton        4.6e-12
brent         0.0e+00
vollib        3.2e-13
quantlib      2.4e-12
vectorized    3.2e-13

Newton, vollib and the vectorized version agree with Brent below $10^{-11}$; QuantLib with a strict tolerance, at $2.4 \cdot 10^{-12}$. These are rounding differences. The numerical method is not what separates a good volatility from a bad one: what separates them is the price it is given.

What the spread costs

With the method settled, the question that matters is a different one: how close is the volatility that comes out of the mid price to the true one? The mid is not the model price: it is inside a spread, with noise, and rounded to the tick (0.05 points below 3 and 0.10 above). A price error $\Delta C$ becomes a volatility error of approximately $\Delta C / \mathcal{V}$, and the vega falls very fast away from the money. That is why it is worth measuring also the half-spread in volatility, $(a - b) / (2\mathcal{V})$: the band within which any volatility is compatible with the quotes.

otm["iv"] = ivs["brent"]
otm["k"] = np.log(otm["strike"] / F)
vg = vega(F, otm["strike"].to_numpy(), T, D,
          otm["iv_modelo"].to_numpy())
otm["error_pts"] = 100 * (otm["iv"] - otm["iv_modelo"])
otm["semi_pts"] = 100 * (otm["ask"] - otm["bid"]) / 2 / vg
sel = otm["strike"].isin([3400, 4000, 4600, 5150, 5600, 6000, 6200])
print(otm.loc[sel, ["tipo", "strike", "bid", "ask", "iv_modelo",
                    "iv", "error_pts", "semi_pts"]]
      .round(4).to_string(index=False))
print(f"Mean absolute error: {otm['error_pts'].abs().mean():.3f}"
      " points")

tipo  strike   bid    ask  iv_modelo     iv  error_pts  semi_pts
   C  5600.0  11.1  11.50     0.1690 0.1692     0.0247    0.0731
   C  6000.0   0.5   0.60     0.1709 0.1709     0.0021    0.1719
   C  6200.0   0.1   0.15     0.1772 0.1763    -0.0867    0.2887
   P  3400.0   2.0   2.10     0.5018 0.5023     0.0469    0.1509
   P  4000.0   9.8  10.10     0.4130 0.4129    -0.0144    0.1150
   P  4600.0  41.9  42.70     0.3188 0.3191     0.0269    0.1014
   P  5150.0 155.2 156.70     0.2231 0.2233     0.0170    0.1057
Mean absolute error: 0.032 points

Near the money the error is hundredths of a point: the 5,150 put recovers 22.33% against the true 22.31%, with a half-spread of 0.11 points. In the 6,200 call, with bid 0.10 and ask 0.15, the half-spread rises to 0.29 points and the error to −0.09. Across the 40 options, the mean absolute error is 0.032 points, and in all of them it is smaller than the half-spread: the mid is within the band of volatilities compatible with the bid and the ask (exercise 2 checks it).

The header figure draws the full smile of the day. The mids fall on the truth line except for the five far calls with a zero bid, from 6,400 to 7,200, which the previous cell had left out. Their mid is half the ask (0.025 points) and their volatility deviates up to 5.8 points from the true one.

Where it fails: zero bids, in the money and stale quotes

One day is not enough to see the rare failures. With the whole chain, 544,548 rows over twelve years, the scalar methods are slow, but Newton with a bisection safeguard vectorizes in a few lines: each option keeps its interval $[\sigma_a, \sigma_b]$ and, if the Newton step falls outside it, it bisects. The forward and the rate of each slice are taken from opciones3_verdad.parquet to isolate the error of the quotes from that of the parity:

x = cad.merge(ver[["fecha", "vencimiento", "dias", "tipo_libre",
                   "forward"]], on=["fecha", "vencimiento"])
Tx = x["dias"].to_numpy() / 365
Dx = np.exp(-x["tipo_libre"].to_numpy() * Tx)
Fx, Kx = x["forward"].to_numpy(), x["strike"].to_numpy()
tx = x["tipo"].to_numpy()
px = ((x["bid"] + x["ask"]) / 2).to_numpy()
intr = Dx * np.maximum(np.where(tx == "C", Fx - Kx, Kx - Fx), 0)
techo = Dx * np.where(tx == "C", Fx, Kx)
valida = (px > intr) & (px < techo)

lo, hi = np.full(len(x), 1e-4), np.full(len(x), 5.0)
s = np.full(len(x), 0.3)
for _ in range(60):
    d = black(Fx, Kx, Tx, Dx, s, tx) - px
    lo, hi = np.where(d < 0, s, lo), np.where(d > 0, s, hi)
    paso = s - d / np.maximum(vega(Fx, Kx, Tx, Dx, s), 1e-12)
    fuera = (paso <= lo) | (paso >= hi) | ~np.isfinite(paso)
    s = np.where(fuera, (lo + hi) / 2, paso)
x["iv"] = np.where(valida, s, np.nan)
x["error_pts"] = 100 * (x["iv"] - x["iv_modelo"])
x["otm"] = np.where(tx == "C", Kx >= Fx, Kx < Fx)
print(f"{len(x)} rows; no solution: {(~valida).sum()}"
      f" ({(~valida).mean():.1%})")

544548 rows; no solution: 37662 (6.9%)

The 37,662 rows with no solution (6.9% of the chain) are in-the-money options whose mid falls below the discounted intrinsic value. It is not a failure of the method: with those prices no implied volatility exists.

inc = pd.read_csv("opciones3_incidencias.csv",
                  parse_dates=["fecha", "vencimiento"])
x = x.merge(inc, on=["fecha", "vencimiento", "tipo", "strike"],
            how="left")
x["clase"] = x["incidencia"]
x.loc[x["clase"].isna() & (x["bid"] <= 0), "clase"] = "bid_nulo"
x.loc[x["clase"].isna() & ~x["otm"], "clase"] = "dentro_del_dinero"
x["clase"] = x["clase"].fillna("fuera_del_dinero")
res = x.groupby("clase").agg(
    filas=("iv", "size"),
    sin_solucion=("iv", lambda z: z.isna().mean()),
    mediana_abs=("error_pts", lambda z: z.abs().median()),
    p95_abs=("error_pts", lambda z: z.abs().quantile(0.95)))
res.sort_values("mediana_abs").round(3)

                       filas  sin_solucion  mediana_abs  p95_abs
clase
fuera_del_dinero      210949         0.000        0.013    0.170
dentro_del_dinero     272261         0.138        0.074   16.265
cotizacion_cruzada        12         0.000        0.127   15.770
paridad_violada            8         0.000        0.494    2.834
cotizacion_congelada    1500         0.000        2.491   27.635
bid_nulo               59818         0.000        4.817   59.862

Out of the money and with no incidents, the median error is 0.013 points and the 95th percentile, 0.17: the volatility is recovered almost exactly. In the money, 13.8% of the rows have no solution and the 95th percentile reaches 16 points, because the time value is small compared with the spread. Zero bids are 59,818 rows, 11% of the chain, with a median of 4.8 points and a 95th percentile of 60. Stale quotes, with a median of 2.5 points: the index moved more than 1% and the quote did not. Crossed quotes and parity violations are few; a parity violation is discovered by comparing the put with the call of the same strike (lesson 3.4.1).

The figure ranks the classes by their median error on a log scale: between the clean row and the zero-bid one there is a factor of 370.

orden = ["vencimiento", "tipo", "strike", "fecha"]
y = x.sort_values(orden).copy()
grupo = y.groupby(["vencimiento", "tipo", "strike"])
igual = ((grupo["bid"].shift() == y["bid"])
         & (grupo["ask"].shift() == y["ask"]))
rend = ind.set_index("fecha")["cierre"].pct_change()
movido = y["fecha"].map(rend).abs() > 0.01
y["sospechosa"] = igual & movido & (y["bid"] > 0)
real = y["clase"] == "cotizacion_congelada"
marca = y["sospechosa"]
print(f"Stale: {real.sum()}, with bid > 0: "
      f"{(real & (y['bid'] > 0)).sum()}")
print(f"Flagged: {marca.sum()}; hits: {(marca & real).sum()}")
err = y["error_pts"].abs()
print(f"Median error, hits: {err[marca & real].median():.2f}"
      f" points; false alarms: {err[marca & ~real].median():.2f}")

Stale: 1500, with bid > 0: 1090
Flagged: 1923; hits: 1081
Median error, hits: 2.07 points; false alarms: 0.16

The rule (same quote as the previous session, index moved more than 1% and positive bid) flags 1,923 rows and catches 1,081 of the 1,090 stale quotes with a positive bid; the other 410 stale quotes have a zero bid and are already removed by the previous filter. The 842 false alarms are cheap options whose price, rounded to the tick, had no reason to change: their median error is 0.16 points, so discarding them costs little. That of the hits is 2.07.

flowchart TB
  A[Chain row] --> B{Bid greater than 0}
  B -- no --> X[Discard]
  B -- yes --> C{Out of the money}
  C -- no --> P[Use the other leg]
  C -- yes --> D{Same quote}
  D -- yes, index moved --> X
  D -- no --> E[Compute sigma]

The diagram summarizes the filters in the order in which they should be applied. “Use the other leg” means taking the put of the same strike instead of the in-the-money call (or vice versa): by parity they have the same volatility.

Exercises

1. Program pure bisection on $[10^{-4}, 5]$ with a tolerance of $10^{-10}$ on the price and count the iterations it needs for the 4,000 put and the 6,000 call of the day. Compare them with Newton from Manaster-Koehler and with Brent.

def biseccion(precio, K, tipo, a=1e-4, b=5.0, tol=1e-10):
    n = 0
    while True:
        m = (a + b) / 2
        dif = float(black(F, K, T, D, m, tipo)) - precio
        n += 1
        if abs(dif) < tol or n == 200:
            return m, n
        a, b = (m, b) if dif < 0 else (a, m)


for K, tipo in [(4000.0, "P"), (6000.0, "C")]:
    s, n = biseccion(medio_de(K, tipo), K, tipo)
    print(f"{K:.0f} {tipo}: σ = {s:.6f} in {n} iterations")

4000 P: σ = 0.412856 in 41 iterations
6000 C: σ = 0.170921 in 39 iterations

2. Compute the implied volatility of the bid and the ask of each out-of-the-money option of the day (with a positive bid) and measure in what fraction iv_modelo falls inside that band. Where does the band widen the most?

banda = []
for f in otm.itertuples():
    lado = [iv_brent(p, f.strike, f.tipo)[0] for p in (f.bid, f.ask)]
    banda.append(lado)
otm[["iv_bid", "iv_ask"]] = banda
dentro = otm["iv_modelo"].between(otm["iv_bid"], otm["iv_ask"])
otm["ancho_pts"] = 100 * (otm["iv_ask"] - otm["iv_bid"])
print(f"iv_modelo inside the band: {dentro.mean():.1%}")
print(otm.groupby(pd.cut(otm["k"], [-0.5, -0.1, 0.1, 0.25]),
                  observed=True)["ancho_pts"].median().round(2))

iv_modelo inside the band: 100.0%
k
(-0.5, -0.1]    0.21
(-0.1, 0.1]     0.18
(0.1, 0.25]     0.25
Name: ancho_pts, dtype: float64

3. A 2,400-strike call of the day is quoted bid 2,730.40 and ask 2,756.10. Explain why its implied volatility should not be computed directly and how you would obtain the volatility of that strike.