--- title: "3.4.4 Implied volatility with Newton, Brent and QuantLib" description: "Implied volatility with Newton, Brent, vollib and QuantLib on chains with known truth, and where it fails: zero bids, in the money, stale quotes." url: "https://pythonparatrading.com/en/python-trading-course/level-3/implied-volatility-newton-brent-quantlib/" lang: "en" image: "https://pythonparatrading.com/og/en/python-trading-course/level-3/implied-volatility-newton-brent-quantlib.png" --- # 3.4.4 · Implied volatility: Newton-Raphson, Brent, py_vollib and QuantLib > The market quotes prices; models and strategies speak in volatilities. This lesson converts one into the other with four methods, checks that they agree to the tenth decimal and measures, on twelve years of synthetic chains with the true volatility known, how much error the bid-ask spread introduces and where the calculation stops making sense. - Course: [Level 3: Derivatives: futures, spreads, forex and options](https://pythonparatrading.com/en/python-trading-course/level-3/) - Module: 3.4 · Option theory and pricing - Duration: 70 min - Prerequisites: [3.4.1 · Payoffs, put-call parity and arbitrage bounds](https://pythonparatrading.com/en/python-trading-course/level-3/option-payoffs-put-call-parity-arbitrage/), [3.4.2 · Binomial trees, Black-Scholes-Merton, Black-76 and Monte Carlo](https://pythonparatrading.com/en/python-trading-course/level-3/binomial-black-scholes-black-76-monte-carlo/), [3.4.3 · Greeks: delta, gamma, vega, theta, rho, vanna, volga and charm](https://pythonparatrading.com/en/python-trading-course/level-3/first-and-second-order-option-greeks/) ## By the end of this lesson - Pose implied volatility as the root of an equation with a unique solution and bound it with the arbitrage bounds. - Program Newton-Raphson with a safe starting point and Brent with a bracket, and explain when and why Newton fails. - Obtain the same volatility with `vollib` (formerly `py_vollib`), `py_vollib_vectorized` and QuantLib, and know the maintenance status of each one. - Measure the error of the volatility recovered from mid prices against the true one and relate it to the bid-ask spread and the vega. - Filter a chain before computing: out-of-the-money options, zero bids, stale and crossed quotes. ## 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 $$ C = D\,[F\,N(d_1) - K\,N(d_2)], \qquad P = D\,[K\,N(-d_2) - F\,N(-d_1)], \qquad d_{1,2} = \frac{\ln(F/K)}{\sigma\sqrt{T}} \pm \frac{\sigma\sqrt{T}}{2}, $$ 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$. ```mermaid 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: ```python 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] ``` ```output 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: ```python 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%}") ``` ```output 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, $$ \sigma_{n+1} = \sigma_n - \frac{C(\sigma_n) - C_{\text{market}}}{\mathcal{V}(\sigma_n)}. $$ 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: ```python 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") ``` ```output 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. ```python 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]) ``` ```output 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. ```python 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}") ``` ```output 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. *Fig. — Newton-Raphson on the price of a put 22% out of the money (strike 4,000, 44 days): from the inflection point it converges; from σ = 20% the almost flat tangent jumps to σ = 12.8.* Curve of the put price against volatility, from 0 to 2.2: almost zero up to 20% and rising afterwards. Log scale. A horizontal line marks the mid price, 9.95 points. Newton steps from σ = 2.05 move down the curve to the root, σ = 0.413 (highlighted). The tangent at σ = 20% barely rises: inside the figure it stays below 2 points, far from the target price. 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`: ```python 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") ``` ```output 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. | Method | Needs | Convergence | Fails if… | |---|---|---|---| | Newton-Raphson | Price and vega | Quadratic near the root | It starts far away with a small vega (wings) | | Newton with Manaster-Koehler start | Price and vega | Quadratic, monotone | Almost never in Black-76; watch very large σ | | Bisection | An interval with a sign change | Linear (one bit per step) | There is no interval: price outside the bounds | | Brent | An interval with a sign change | Superlinear | There is no interval | | Let's Be Rational (Jäckel) | Normalized price | Machine precision in two steps | Price outside the bounds | ## 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: ```python 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}") ``` ```output 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: ```python 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) ``` ```output 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. ```python 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)) ``` ```output 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 4[^ql]; here it is given $10^{-10}$, 200 evaluations and σ between $10^{-4}$ and 5: ```python 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%}") ``` ```output 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: ```python 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()) ``` ```output 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. ```python 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") ``` ```output 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: ```python 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%})") ``` ```output 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. ```python 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) ``` ```output 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). *Fig. — Median error of the volatility recovered from the mid price, in volatility points, by row type of the ORX daily chain, 2014-2025 (number of rows in brackets).* Horizontal bars on a log scale. Out-of-the-money options without incidents: 0.013 points (210,949 rows). In the money: 0.074, with 13.8% having no solution. Crossed quotes: 0.127. Parity breach: 0.494. Stale: 2.491. Zero bid (highlighted): 4.817 points over 59,818 rows. 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. ```python 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}") ``` ```output 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. ```mermaid 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. **Solution.** ```python 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") ``` ```output 4000 P: σ = 0.412856 in 41 iterations 6000 C: σ = 0.170921 in 39 iterations ``` Bisection needs 41 and 39 iterations, against 7 to 9 for Newton from Manaster-Koehler and 15 to 19 for Brent. Each bisection step gains one bit: to go from an interval of 5 to a precision of $10^{-10}$ in the price takes about forty halvings. 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? **Solution.** ```python 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)) ``` ```output 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 ``` The true volatility falls inside the band for 100% of the options. The band measures between 0.18 and 0.25 points in median across the three ranges of $k = \ln(K/F)$: with the zero-bid options out, the spread in volatility is narrow over the whole range of this expiry. For the zero-bid ones the band would go from zero up to the volatility of the ask. 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. **Solution.** It is a deep in-the-money call. Its discounted intrinsic value, $D\,(F - K) \approx 2{,}741.2$ points, is almost all of the price: the mid, 2,743.25, leaves about 2 points of time value, and the spread measures 25.7. With those numbers, the volatility of the mid comes out at 87.8%, against the true 65.3%, and it is enough to move the price 0.7 points (from 2,741.3 to 2,742.0) to go from 66.9% to 79.5%; a little lower, there is no solution. The volatility of that strike is obtained from the 2,400 put (bid 0.05, ask 0.10), which by parity has the same one and whose price is all time value: its mid gives 65.2%. [^ql]: Signature of `VanillaOption::impliedVolatility` in QuantLib 1.44, `ql/instruments/vanillaoption.hpp`: https://github.com/lballabio/QuantLib/blob/v1.44/ql/instruments/vanillaoption.hpp ## Self-assessment 1. Why is the implied volatility of a European option, if it exists, unique? - A. Because the market only quotes one volatility per expiry - B. Because Newton-Raphson always converges to the same root - C. It is not unique: there are two roots, one for the call and one for the put - D. Because the Black-Scholes price is strictly increasing in σ (the vega is positive), so each price within the arbitrage bounds corresponds to a single σ *Answer D: Because the Black-Scholes price is strictly increasing in σ (the vega is positive), so each price within the arbitrage bounds corresponds to a single σ.* The vega is positive for any σ > 0, so the price is a strictly increasing function of σ between the discounted intrinsic value (as σ tends to 0) and the upper bound (as σ grows without limit). The call and the put of the same strike have the same σ if parity holds. 2. Newton-Raphson starting from σ = 20% jumps to σ = 12.8 in the first iteration for a put 22% out of the money. What explains it? - A. A numpy rounding error - B. That the vega at σ = 20% is tiny for that put: the tangent is almost horizontal and crosses the target price very far away - C. That the put is mispriced - D. That Newton only works with calls *Answer B: That the vega at σ = 20% is tiny for that put: the tangent is almost horizontal and crosses the target price very far away.* With 44 days and σ = 20%, a put at 4,000 with the forward at 5,155 is almost five standard deviations away: its price and its vega are almost zero. The Newton step is (price − target)/vega, and with a tiny vega the step is huge. From the Manaster and Koehler inflection point it converges in a few iterations. 3. The bid-ask spread of an out-of-the-money put is 0.10 points and its vega is 3 price points per volatility point. What is, approximately, the uncertainty of its implied volatility when using the mid price? - A. 0.10 volatility points - B. None: the mid price is the true price - C. About ±0.017 volatility points: half the spread divided by the vega - D. 3 volatility points *Answer C: About ±0.017 volatility points: half the spread divided by the vega.* The volatility of the bid and that of the ask are ±(0.05 / 3) ≈ ±0.017 points from that of the mid. Far from the money, the vega falls much faster than the spread and that band widens to tens of points. 4. A far put is quoted bid 0 and ask 0.05. What do you do with it when building the smile? - A. Exclude it: with a zero bid nobody buys it and the mid is not a market price; in the data, its recovered volatility deviates a median of 4.8 points from the true one - B. Use the mid, 0.025: it is as valid as any other - C. Use the ask, which is the price at which you would buy - D. Set its volatility to zero *Answer A: Exclude it: with a zero bid nobody buys it and the mid is not a market price; in the data, its recovered volatility deviates a median of 4.8 points from the true one.* A zero bid says that the value is between 0 and the ask: the mid is a convention, not a price. Cboe excludes those options from the VIX calculation for the same reason (lesson 3.4.6). 5. Why are out-of-the-money options (puts below the forward, calls above) used for the smile? - A. Because they are cheaper - B. Because in-the-money options have no implied volatility in theory - C. Because put-call parity requires it - D. Because their price is all time value: in the money, the intrinsic value dominates and a spread of a few points can leave the mid below the discounted intrinsic value, with no solution *Answer D: Because their price is all time value: in the money, the intrinsic value dominates and a spread of a few points can leave the mid below the discounted intrinsic value, with no solution.* By parity, the call and the put of the same strike have the same volatility, so you can choose the more informative one. In the data, 13.8% of in-the-money options have a mid below the discounted intrinsic value; out of the money, none. ## Key points - Implied volatility is the σ that equates the Black-76 price to the market price; it exists and is unique within the arbitrage bounds because the vega is positive. - Newton converges in 7 to 9 iterations if it starts at the Manaster-Koehler inflection point; from a fixed value it can shoot off on far options. Brent never fails if there is a bracket, at the cost of more evaluations. - `vollib` (the current name of `py_vollib`) uses Jäckel's algorithm and agrees with QuantLib and with our own methods to the eleventh decimal; `py_vollib_vectorized` 0.1.1, from 2021, does not compile with current numba. - The error of the volatility recovered from the mid is of the order of half the spread divided by the vega: hundredths of a point near the money and tens of points in the wings. - Before computing, filter: only out of the money, no zero bids, no crossed or stale quotes. What passes the filter is recovered with a median error of 0.013 points. ## Resources ### Documentation - [vollib: documentation (Black, Black-Scholes and Black-Scholes-Merton with Let's Be Rational)](https://vollib.org/) · vollib.org - [PyPI: py-vollib-vectorized 0.1.1 (latest version, February 2021)](https://pypi.org/project/py-vollib-vectorized/) · pypi.org - [QuantLib 1.44: blackformula.hpp (blackFormulaImpliedStdDev and its default tolerance)](https://github.com/lballabio/QuantLib/blob/v1.44/ql/pricingengines/blackformula.hpp) · github.com - [SciPy: scipy.optimize.brentq](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html) · docs.scipy.org ### Books - [Jäckel, P. (2015): Let's Be Rational. Wilmott 2015(75), 40-53](https://doi.org/10.1002/wilm.10395) · doi.org - [Manaster, S. and Koehler, G. (1982): The Calculation of Implied Variances from the Black-Scholes Model: A Note. Journal of Finance 37(1), 227-230](https://doi.org/10.1111/j.1540-6261.1982.tb01105.x) · doi.org ### Data - [Level 3 data: daily option chain on ORX, 2014-2025 (Parquet)](https://pythonparatrading.com/curso/datos/nivel-3/opciones3_cadena.parquet) · Download · opciones3_cadena.parquet Previous lesson: [3.4.3 · Greeks: delta, gamma, vega, theta, rho, vanna, volga and charm](https://pythonparatrading.com/en/python-trading-course/level-3/first-and-second-order-option-greeks/) Next lesson: [3.4.5 · Smile, skew and term structure; volatility surface with SVI](https://pythonparatrading.com/en/python-trading-course/level-3/volatility-smile-skew-svi-surface/)