Advanced option greeks in Python, on a live chain

Delta, gamma, theta and vega are the four everyone learns. This page is about what sits one and two derivatives past them, and about the part that usually gets skipped, which is checking that the formulas are right before trusting a single number they produce.

Everything here comes from a live NVDA option chain rather than a worked example. The engine prices under Black-Scholes-Merton with a continuous dividend yield, computes the closed forms, then differentiates the price surface numerically to confirm each one. That last step caught two real sign errors, which is why it lives in the script rather than in a comment.

1. The greeks, in one function

Order is how many times you differentiate, and with respect to what. Gamma is delta differentiated by spot. Vanna is delta differentiated by volatility, which is the same object as vega differentiated by spot. Zomma is gamma differentiated by volatility, color is gamma differentiated by time, speed is gamma differentiated by spot a second time.

One convention has to be stated before any of it means anything. Time derivatives here are with respect to calendar time, so they are reported per day and come out negative for a long option. Published tables disagree on this point, and that disagreement is exactly where the sign errors live.

greeks.py

import math
import numpy as np

SQ2PI = math.sqrt(2.0 * math.pi)

def _phi(x):
    return np.exp(-0.5 * np.asarray(x, dtype=float) ** 2) / SQ2PI

def _N(x):
    x = np.asarray(x, dtype=float)
    return 0.5 * (1.0 + np.vectorize(math.erf)(x / math.sqrt(2.0)))

def d1d2(S, K, tau, r, q, sig):
    v = sig * np.sqrt(tau)
    a = (np.log(S / K) + (r - q + 0.5 * sig ** 2) * tau) / v
    return a, a - v

def greeks(S, K, tau, r, q, sig, kind="call"):
    a, b = d1d2(S, K, tau, r, q, sig)
    rt, pa = np.sqrt(tau), _phi(a)
    eq, er = np.exp(-q * tau), np.exp(-r * tau)

    vega  = S * eq * pa * rt
    gamma = eq * pa / (S * sig * rt)

    # second order
    vanna = -eq * pa * b / sig                  # dDelta/dsig = dVega/dS
    volga = vega * a * b / sig                  # dVega/dsig
    # published veta and color tables are written in tau, charm is written in t.
    # negate both so every time derivative here is d/dt. Checked in section 2.
    veta  = S * eq * pa * rt * (q + (r - q) * a / (sig * rt) - (1 + a * b) / (2 * tau))

    # third order
    speed = -gamma / S * (a / (sig * rt) + 1.0)   # dGamma/dS
    zomma = gamma * (a * b - 1.0) / sig           # dGamma/dsig
    color = (eq * pa / (2 * S * tau * sig * rt)
             * (2 * q * tau + 1.0
                + (2 * (r - q) * tau - b * sig * rt) / (sig * rt) * a))
    ultima = -vega / sig ** 2 * (a * b * (1 - a * b) + a * a + b * b)

    return dict(vega=vega, gamma=gamma, vanna=vanna, volga=volga, veta=veta,
                speed=speed, zomma=zomma, color=color, ultima=ultima)

2. The check that matters

A closed form greek is only as good as its transcription, and the higher order ones are easy to get wrong. So before any number is used, the script differentiates the price surface numerically and compares. If the two disagree by more than a small tolerance it refuses to continue.

This is what caught the two sign errors. Veta and color both came out exactly negated against the numerical derivative, because the tables they were taken from are written in time to expiry while the charm formula beside them is written in calendar time. A relative error of 2.00 is the signature of a sign flip, and it is unmissable once you look for it.

validate.py

def _fd(f, x, h, order=1):
    if order == 1:
        return (f(x + h) - f(x - h)) / (2 * h)
    if order == 2:
        return (f(x + h) - 2 * f(x) + f(x - h)) / h ** 2

# differentiate the PRICE SURFACE numerically, compare to every closed form
S, K, tau, r, q, sig = 100.0, 100.0, 0.25, 0.042, 0.005, 0.28
g = greeks(S, K, tau, r, q, sig)
P = lambda s=S, t=tau, v=sig: float(bs_price(s, K, t, r, q, v, "call"))
hS, hV, hT = S * 2e-4, 1e-4, 1e-5

num = {
    "vanna": _fd(lambda v: _fd(lambda s: P(s=s, v=v), S, hS), sig, hV),
    "volga": _fd(lambda v: P(v=v), sig, hV, order=2),
    # note the MINUS: these are d/dt, and tau runs the other way
    "veta":  -_fd(lambda t: float(greeks(S, K, t, r, q, sig)["vega"]), tau, hT),
    "speed": _fd(lambda s: float(greeks(s, K, tau, r, q, sig)["gamma"]), S, hS),
    "zomma": _fd(lambda v: float(greeks(S, K, tau, r, q, v)["gamma"]), sig, hV),
    "color": -_fd(lambda t: float(greeks(S, K, t, r, q, sig)["gamma"]), tau, hT),
}
for name, nv in num.items():
    av = float(g[name])
    err = abs(av - nv) / max(abs(av), 1e-12)
    print(f"{name:<7}{av:>14.6g}{nv:>14.6g}{err:>11.1e}")

After negating veta and color the worst relative error across all 48 checks is 3.6e-05, which is finite difference noise rather than a mistake. Ten of the twelve formulas were right the first time. Two were not, and nothing except this check would have told me.

3. Pulling a real chain

Toy examples hide the parts that bite. A live chain brings stale quotes, zero bids, missing implied vols and, in this case, a dividend yield that yfinance reported as 50 percent for NVDA. Left uncorrected that one field distorts every greek on the page, so the script rebuilds it from the actual trailing cash dividends instead.

chain.py

import yfinance as yf, datetime as dt

def fetch_chain(ticker="NVDA", max_expiries=14):
    tk = yf.Ticker(ticker)
    spot = float(tk.fast_info["last_price"])
    # info["dividendYield"] flips between fraction and percent and has returned
    # 50.0 for NVDA, so rebuild q from the actual trailing cash dividends
    div = tk.dividends
    cutoff = div.index.max() - dt.timedelta(days=365)
    q = float(div[div.index > cutoff].sum()) / spot
    if not (0.0 <= q < 0.25):
        q = 0.0
    r = float(yf.Ticker("^IRX").fast_info["last_price"]) / 100.0

    today, out = dt.date.today(), []
    for exp in tk.options[:max_expiries]:
        tau = (dt.date.fromisoformat(exp) - today).days / 365.0
        if tau <= 0.003:
            continue
        ch = tk.option_chain(exp)
        for kind, df in (("call", ch.calls), ("put", ch.puts)):
            for _, row in df.iterrows():
                iv = float(row.get("impliedVolatility") or 0.0)
                bid, ask = float(row.get("bid") or 0), float(row.get("ask") or 0)
                if not (0.01 < iv < 4.0) or bid <= 0 or ask <= 0:
                    continue
                out.append(dict(expiry=exp, tau=tau, kind=kind,
                                K=float(row["strike"]), iv=iv,
                                mid=0.5 * (bid + ask),
                                oi=int(row.get("openInterest") or 0)))
    return dict(ticker=ticker, spot=spot, r=r, q=q, rows=out)

The snapshot used throughout this page is NVDA at 221.88, a risk free rate of 3.73 percent from the three month bill, a dividend yield of 0.13 percent, and 1,546 contracts with a live two sided market. On the 30 day expiry the closest listed strike to spot is 220 at 43.8 percent implied, which puts one standard deviation at 27.85 dollars.

4. Where each greek lives on the chain

The most useful output is not a single number, it is the shape. Run the greeks across every listed strike and a geography appears that holds across names.

Vanna crosses zero at the money and peaks about one standard deviation out. Volga is effectively zero at the money and largest in the wings. Zomma runs negative through the middle and turns positive at both ends. There is a reason for the pattern. Any greek that is the spot derivative of something peaking at the money has to be zero there, because the derivative of a maximum is zero.

Screening eight names on the same day, SPY, QQQ, NVDA, TSLA, AAPL, MSFT, AMZN and META, the vanna peak sits within a rounding error of one standard deviation on every one of them, with implied vols running from 15 to 47 percent. Measure distance in dollars and the peaks scatter. Measure it in standard deviations and they line up.

5. Does any of it change a P&L

The honest test. Take NVDA worst single day of the last twelve months, minus 6.40 percent, add the five volatility points a day like that brings, reprice three live contracts exactly, then ask how much of the move each order of the Taylor expansion explains.

attribution.py

def taylor_attribution(S, K, tau, r, q, sig, kind, dS, dsig, days):
    g = greeks(S, K, tau, r, q, sig, kind)
    dt_yr = days / 365.0
    exact = (float(bs_price(S + dS, K, tau - dt_yr, r, q, sig + dsig, kind))
             - float(g["price"]))

    o1 = g["delta"] * dS + g["vega"] * dsig + g["theta"] * dt_yr
    o2 = (0.5 * g["gamma"] * dS ** 2
          + 0.5 * g["volga"] * dsig ** 2
          + g["vanna"] * dS * dsig
          + g["charm"] * dS * dt_yr
          + g["veta"] * dsig * dt_yr)
    o3 = ((1 / 6) * g["speed"] * dS ** 3
          + 0.5 * g["zomma"] * dS ** 2 * dsig
          + 0.5 * g["color"] * dS ** 2 * dt_yr
          + (1 / 6) * g["ultima"] * dsig ** 3)

    return dict(exact=exact, o1=o1, o2=o1 + o2, o3=o1 + o2 + o3)

Read those three rows and you have the answer. Delta and vega alone overshoot by a quarter at the money and by a third on the wing, because a straight line is the wrong shape for a day that actually moves. Second order does almost all of the repair. Third order takes the at the money contract from 97.6 to 99.9 percent, a difference nobody would ever notice in a P&L, and earns its place only out in the wings.

Which is the useful conclusion rather than the flattering one. Most of the time the third order greeks are a rounding error, and plenty of good traders never look at them. They start to matter in three specific places. Expiry week, when color goes vertical and a gamma number is stale within a day. A volatility regime shift, when zomma moves gamma out of the middle and into the wings. And large moves far from the strike. Outside those three, they are not where the money is.

Running it

The script needs numpy for the maths and yfinance for the chain. It validates first and refuses to print a single number if the check fails, which is the behaviour you want from anything that computes a hedge ratio.

terminal

pip install numpy yfinance
py -3 higher_order_greeks.py --validate      # formula checks only, no network
py -3 higher_order_greeks.py --ticker NVDA   # full run on a live chain

Quotes move, so your numbers will not match these exactly. The shapes will.