Greeks: delta, gamma, vega, theta, rho, vanna, volga and charm.

What each Greek measures

An option changes value because the underlying, the volatility, time and the interest rate move. The Greeks are the partial derivatives of the price with respect to each one: they say how much the price changes per unit change of each factor, with the others held still. They are the tool with which a portfolio is hedged (lesson 3.6.2), its risk measured and its P&L explained. The first-order ones are first derivatives; the second-order ones measure how much the first ones change.

flowchart TB
  V[Price V] --> D[Delta]
  V --> VE[Vega]
  V --> TH[Theta]
  V --> RH[Rho]
  D --> G[Gamma]
  D --> VA[Vanna]
  D --> CH[Charm]
  VE --> VO[Volga]

Vanna sits under delta and could sit under vega: it is the cross derivative, and the two readings are the same. The table fixes the units used throughout the course. Almost all errors with Greeks are about units, and each library has its own:

ORX prices are per index point with a multiplier of 100 USD: a delta of 0.52 on one contract is equivalent to 52 “indices” and, with the index at 5,148, to 268,000 USD of exposure; a theta of −2.32 is 232 USD per contract per day.

Closed-form formulas in Black-Scholes-Merton

With the same variables as in lesson 3.4.2 ($d_1$, $d_2$, continuous $q$) and $n(\cdot)$ the standard normal density, the first-order Greeks of a call are:

and the second-order ones, the same for call and put except for charm:

Gamma, vega, vanna and volga are the same for the call and the put of the same strike, because the difference between the two is a forward (parity, lesson 3.4.1), linear in $S$ and independent of $\sigma$. The function returns all of them already in the units of the table, for a call or a put:

import math
import numpy as np
import pandas as pd
from scipy.stats import norm

def griegas(S, K, T, r, q, sig, call=True):
    sq = np.sqrt(T)
    sd = sig * sq
    d1 = (np.log(S / K) + (r - q) * T) / sd + sd / 2
    d2 = d1 - sd
    eq, er, n1 = np.exp(-q * T), np.exp(-r * T), norm.pdf(d1)
    w = 1 if call else -1                 # call or put sign
    precio = w * (S * eq * norm.cdf(w * d1)
                  - K * er * norm.cdf(w * d2))
    delta = w * eq * norm.cdf(w * d1)
    theta = (-S * eq * n1 * sig / (2 * sq)
             - w * r * K * er * norm.cdf(w * d2)
             + w * q * S * eq * norm.cdf(w * d1))
    rho = w * K * T * er * norm.cdf(w * d2)
    charm = (w * q * eq * norm.cdf(w * d1)
             - eq * n1 * (2 * (r - q) * T - d2 * sd) / (2 * T * sd))
    vega = S * eq * n1 * sq
    return {"precio": precio, "delta": delta,
            "gamma": eq * n1 / (S * sd),
            "vega": vega / 100,               # per vol. point
            "theta": theta / 365,             # per calendar day
            "rho": rho / 100,                 # per rate point
            "vanna": -eq * n1 * d2 / sig / 100,
            "volga": vega * d1 * d2 / sig / 100**2,
            "charm": charm / 365}

# the ORX call from lesson 3.4.2
S, K, T, r, q, sig = 5148.06, 5150.0, 30 / 365, 0.04191, 0.013438, 0.217
call = griegas(S, K, T, r, q, sig)
put = griegas(S, K, T, r, q, sig, call=False)
pd.DataFrame({"call": call, "put": put}).map("{:.6g}".format)

                call           put
precio        132.59       122.503
delta       0.524404     -0.474492
gamma     0.00124182    0.00124182
vega         5.86996       5.86996
theta       -2.31833      -1.91836
rho          2.10992       -2.1084
vanna   -8.31677e-06  -8.31677e-06
volga    7.69206e-06   7.69206e-06
charm   -0.000476373  -0.000513149

The at-the-money call at 30 days is worth 132.59 points with a delta of 0.524 and the put, one of −0.474: the difference is the delta of the synthetic forward, $e^{-qT}$ = 0.9989. Gamma, vega, vanna and volga match for both. Vanna and volga are almost zero at the money ($d_2 \approx 0$); they show up in the wings, further down. The call's rho is positive (a higher rate makes the forward more expensive) and the put's, negative.

Three ways to compute them: formula, finite differences and libraries

Formulas can be copied wrong; finite differences cannot. A central difference approximates the derivative with two prices around the point:

The cell compares the formulas with finite differences on the pricing function itself, with py_vollib (vollib package, lesson 3.4.2) and with QuantLib, printing each library in its native units:

import QuantLib as ql
from vollib.black_scholes_merton.greeks import analytical as vl

def precio(S=S, T=T, r=r, sig=sig):
    return griegas(S, K, T, r, q, sig)["precio"]

def delta_de(S=S, T=T, sig=sig):
    return griegas(S, K, T, r, q, sig)["delta"]

h, hs, hr = 1.0, 1e-4, 1e-4
dif = {
    "delta": (precio(S + h) - precio(S - h)) / (2 * h),
    "gamma": (precio(S + h) - 2 * precio() + precio(S - h)) / h**2,
    "vega": (precio(sig=sig + hs) - precio(sig=sig - hs))
            / (2 * hs) / 100,
    "theta": (precio(T=T - 1 / 365) - precio()),      # one day
    "rho": (precio(r=r + hr) - precio(r=r - hr)) / (2 * hr) / 100,
    "vanna": (delta_de(sig=sig + hs) - delta_de(sig=sig - hs))
             / (2 * hs) / 100,
    "charm": delta_de(T=T - 1 / 365) - delta_de(),     # one day
}
F = S * math.exp((r - q) * T)
bc = ql.BlackCalculator(ql.PlainVanillaPayoff(ql.Option.Call, K),
                        F, sig * math.sqrt(T), math.exp(-r * T))
origen = {
    "delta": (vl.delta("c", S, K, T, r, sig, q), bc.delta(S)),
    "gamma": (vl.gamma("c", S, K, T, r, sig, q), bc.gamma(S)),
    "vega": (vl.vega("c", S, K, T, r, sig, q), bc.vega(T)),
    "theta": (vl.theta("c", S, K, T, r, sig, q), bc.theta(S, T)),
    "rho": (vl.rho("c", S, K, T, r, sig, q), bc.rho(T)),
    "vanna": (np.nan, bc.vanna(S, T)),
    "charm": (np.nan, np.nan),
}
pd.DataFrame({"formula": {k: call[k] for k in dif},
              "finite diff.": dif,
              "vollib": {k: v[0] for k, v in origen.items()},
              "QuantLib (native)": {k: v[1]
                                    for k, v in origen.items()}}
             ).map("{:.6g}".format)

            formula  finite diff.      vollib QuantLib (native)
delta      0.524404      0.524404    0.524404          0.524404
gamma    0.00124182    0.00124182  0.00124182        0.00124182
vega        5.86996       5.86996     5.86996           586.996
theta      -2.31833      -2.33649    -2.31833          -846.192
rho         2.10992       2.10992     2.10992           210.992
vanna  -8.31677e-06   -8.3169e-06         nan      -0.000831677
charm  -0.000476373  -0.000481304         nan               nan

Formula and finite differences match on all the instantaneous Greeks. Theta and charm by a one-day difference are somewhat larger in absolute value (−2.336 against −2.318 per day): the formula's theta is an instantaneous rate, and time value decays faster and faster as expiry approaches. py_vollib already gives the working units: vega and rho divided by 100 and theta by 365. QuantLib gives vega (587.0) and rho (210.99) per unit, that is, per jump of 100 volatility or rate points, and theta per year (−846.19); for daily theta it has thetaPerDay. None is wrong: they measure the same thing in other units, and a system that mixes two libraries without converting units can hedge a hundred times too much or too little.

Volga holds a surprise. QuantLib 1.44 has BlackCalculator.volga(T), but its value does not match the second derivative with respect to σ except with $T = 1$:

filas = []
for Tx in [0.25, 1.0, 4.0]:
    Fx = S * math.exp((r - q) * Tx)
    v = lambda s: ql.BlackCalculator(
        ql.PlainVanillaPayoff(ql.Option.Call, 5000.0), Fx,
        s * math.sqrt(Tx), math.exp(-r * Tx)).value()
    bcx = ql.BlackCalculator(ql.PlainVanillaPayoff(ql.Option.Call,
                                                   5000.0),
                             Fx, sig * math.sqrt(Tx),
                             math.exp(-r * Tx))
    fd = (v(sig + 1e-3) - 2 * v(sig) + v(sig - 1e-3)) / 1e-6
    filas.append({"T": Tx, "QuantLib volga": bcx.volga(Tx),
                  "finite diff.": fd,
                  "ratio": bcx.volga(Tx) / fd,
                  "1/sqrt(T)": 1 / math.sqrt(Tx)})
pd.DataFrame(filas).round(4)

      T  QuantLib volga  finite diff.  ratio  1/sqrt(T)
0  0.25        953.2139      476.6165    2.0        2.0
1  1.00        512.1042      512.1166    1.0        1.0
2  4.00        475.6863      951.4064    0.5        0.5

The ratio is exactly $1/\sqrt{T}$: the “volga” of BlackCalculator is the derivative of vega with respect to the total standard deviation $\sigma\sqrt{T}$, not with respect to $\sigma$. To use it as volga it has to be multiplied by $\sqrt{T}$. It is not an error of method but of convention (or of documentation), and it only shows up by checking it. That is the rule of the lesson: no Greek from a library is used without checking it once against finite differences.

The step of the finite difference

A large step $h$ gives a truncation error (the central difference is off by $h^2$ times the next derivative); a small one gives a rounding error: with double precision ($\varepsilon \approx 2.2\cdot 10^{-16}$), subtracting two nearly equal prices of 132 points and dividing by $h^2$ amplifies that epsilon. For gamma the total error is of the order of $c_1 h^2 + c_2\,\varepsilon V/h^2$, minimal with $h \propto \varepsilon^{1/4}$, of the order of one ten-thousandth of the underlying (about 0.6 points with ORX at 5,148):

pasos = 10.0 ** np.arange(-4, 2.75, 0.25)
gamma_fd = np.array([(precio(S + x) - 2 * precio() + precio(S - x))
                     / x**2 for x in pasos])
err_rel = np.abs(gamma_fd / call["gamma"] - 1)
mejor = pasos[err_rel.argmin()]
print(f"Step with the smallest error: {mejor:.2f} points "
      f"(relative error {err_rel.min():.1e})")
for x in [1e-4, 1e-3, 1e-2, 1.0, 100.0]:
    i = np.argmin(abs(pasos - x))
    print(f"h = {x:g}: relative error {err_rel[i]:.1e}")

Step with the smallest error: 0.06 points (relative error 1.1e-08)
h = 0.0001: relative error 1.1e-02
h = 0.001: relative error 2.9e-04
h = 0.01: relative error 3.8e-06
h = 1: relative error 7.9e-07
h = 100: relative error 7.9e-03

With steps of one ten-thousandth of a point the numerical gamma is garbage (an error of 1.1%); between 0.01 and a few points, the error is in the millionths or less (the minimum measured, at 0.06 points, is partly luck of rounding); with 100 points it grows again from truncation. The figure draws the whole V. For vega and rho the trade-off is the same on another scale: $10^{-4}$ in volatility or rate is close to the optimum.

Finite differences with QuantLib: moving the quotes

For what a library does not provide (vanna, volga, charm in QuantLib with an arbitrary engine), the input is moved and the option repriced. In QuantLib the inputs are observable quotes (SimpleQuote): when their value changes, everything that depends on them is recalculated automatically. For charm, the valuation date is moved one day with curves that move along with it (those created with settlementDays = 0):

hoy = ql.Date(18, 6, 2025)
ql.Settings.instance().evaluationDate = hoy
dc, cal = ql.Actual365Fixed(), ql.NullCalendar()
q_s, q_v = ql.SimpleQuote(S), ql.SimpleQuote(sig)

def plana(x):
    return ql.YieldTermStructureHandle(ql.FlatForward(
        0, cal, ql.QuoteHandle(ql.SimpleQuote(x)), dc))

proceso = ql.BlackScholesMertonProcess(
    ql.QuoteHandle(q_s), plana(q), plana(r),
    ql.BlackVolTermStructureHandle(
        ql.BlackConstantVol(0, cal, ql.QuoteHandle(q_v), dc)))
opcion = ql.VanillaOption(ql.PlainVanillaPayoff(ql.Option.Call, K),
                          ql.EuropeanExercise(hoy + 30))
opcion.setPricingEngine(ql.AnalyticEuropeanEngine(proceso))

def con(quote, valor, f):
    viejo = quote.value()
    quote.setValue(valor)
    x = f()
    quote.setValue(viejo)
    return x

vanna = (con(q_v, sig + hs, opcion.delta)
         - con(q_v, sig - hs, opcion.delta)) / (2 * hs) / 100
volga = (con(q_v, sig + hs, opcion.vega)
         - con(q_v, sig - hs, opcion.vega)) / (2 * hs) / 100**2
d_hoy = opcion.delta()
ql.Settings.instance().evaluationDate = hoy + 1
charm = opcion.delta() - d_hoy
ql.Settings.instance().evaluationDate = hoy
print(f"vanna {vanna:.3e} (formula {call['vanna']:.3e})")
print(f"volga {volga:.3e} (formula {call['volga']:.3e})")
print(f"one-day charm {charm:.3e} (formula {call['charm']:.3e})")

vanna -8.317e-06 (formula -8.317e-06)
volga 7.692e-06 (formula 7.692e-06)
one-day charm -4.813e-04 (formula -4.764e-04)

All three match the formulas (charm, with the one-day difference versus the instantaneous rate already seen with theta). Moving a quote and repricing is also the way to compute Greeks in models without a formula: a tree for an American option (lesson 3.4.2) or Heston (lesson 3.4.7).

Profiles: where each Greek concentrates

The header figure draws gamma as a function of the index for three maturities, with the volatility of the ORX call: gamma is a bell centered near the strike that narrows and grows as expiry approaches. At the money it grows like $1/\sqrt{T}$: at 7 days it is twice that at 30 and almost four times that at 90. Theta follows gamma, because the Black-Scholes equation ties them together. For any European option in that model:

Whoever is long gamma (gains from large moves in either direction) pays theta every day; whoever sells it, collects it. It is the trade-off exploited by gamma scalping and premium selling (lessons 3.6.1 and 3.6.2). The cell computes the at-the-money Greeks for the three maturities and checks the equation with annual units:

filas = []
for dias in [7, 30, 90]:
    g = griegas(S, K, dias / 365, r, q, sig)
    pde = (365 * g["theta"] + 0.5 * sig**2 * S**2 * g["gamma"]
           + (r - q) * S * g["delta"] - r * g["precio"])
    filas.append({"days": dias, "precio": g["precio"],
                  "gamma": g["gamma"], "theta": g["theta"],
                  "vega": g["vega"],
                  "theta/price (%)": 100 * g["theta"] / g["precio"],
                  "BS equation": pde})
pd.DataFrame(filas).set_index("days").map("{:.4g}".format)

     precio     gamma   theta   vega theta/price (%) BS equation
days
7     62.13  0.002578  -4.603  2.843          -7.409  -2.078e-13
30    132.6  0.001242  -2.318   5.87          -1.749  -2.567e-13
90    237.1  0.000712  -1.409   10.1          -0.594  -1.492e-13

At 7 days the call is worth 62.13 points and loses 4.60 a day, 7.4% of its value; at 90 days it loses 1.41 a day, 0.59%. Vega behaves the opposite way to gamma: it grows with $\sqrt{T}$, so long-dated options are the ones with the most exposure to volatility and short-dated ones the most exposure to the move. The Black-Scholes equation gives zero at all three maturities up to floating-point rounding.

The Greeks of a chain: delta, probability and wings

In a real chain each strike has its own implied volatility (the smile, lesson 3.4.5), and the Greeks are computed with that of their strike. The cell does it with the 30-day ORX expiry of June 18, 2025 and the iv_modelo column:

cad = pd.read_parquet("opciones3_cadena.parquet")
c30 = cad[(cad["fecha"] == "2025-06-18")
          & (cad["vencimiento"] == "2025-07-18")
          & (cad["tipo"] == "C")].copy()
g = griegas(S, c30["strike"], T, r, q, c30["iv_modelo"])
for k in ["delta", "gamma", "vega", "vanna", "volga"]:
    c30[k] = g[k]
sd = c30["iv_modelo"] * math.sqrt(T)
d2 = (np.log(S / c30["strike"]) + (r - q) * T) / sd - sd / 2
c30["N(d2)"] = norm.cdf(d2)
c30["delta_usd"] = 100 * S * c30["delta"]   # per contract
sel = c30["strike"].isin([4600.0, 4900.0, 5150.0, 5400.0, 5700.0])
c30.loc[sel, ["strike", "iv_modelo", "delta", "N(d2)", "delta_usd",
              "gamma", "vega", "vanna", "volga"]].round(5)

        strike  iv_modelo    delta    N(d2)     delta_usd    gamma     vega    vanna    volga
519219  4600.0     0.3300  0.89551  0.87844  461011.41198  0.00037  2.65304 -0.00636  0.11841
519225  4900.0     0.2705  0.75899  0.73507  390734.43659  0.00078  4.58487 -0.00721  0.07515
519230  5150.0     0.2170  0.52440  0.50018  269966.39702  0.00124  5.86996 -0.00001  0.00001
519235  5400.0     0.1722  0.18500  0.17231   95237.01786  0.00105  3.93792  0.01464  0.19359
519240  5700.0     0.1615  0.01673  0.01492    8614.73892  0.00017  0.61380  0.00559  0.17552

Two readings of the table:

• Delta is not the probability of exercise. The risk-neutral probability of ending in the money is $N(d_2)$, not $N(d_1)$: for the 5,400 call the delta is 0.185 and $N(d_2)$, 0.172. The difference grows with volatility and time to expiry. And both are risk-neutral probabilities, not forecasts (lesson 3.4.2).
• Vanna and volga live in the wings. At the money they are almost zero; for the 4,600 call (deep in the money) vanna is −0.0064 and for the 5,400 one (out of the money), +0.0146: raising volatility lowers the delta of in-the-money calls and raises that of out-of-the-money ones. Volga is positive in the wings: far options gain from volatility convexly, and that is why the market pays more for them (part of the explanation of the smile, lesson 3.4.5).

Explaining a P&L with the Greeks

The daily use of the Greeks is to break down what a position has gained or lost. With a Taylor expansion of the price up to second order:

On April 3, 2025 ORX fell 4.9% in one day (one of the jumps planted in the data, like that of the S&P 500 on that day). We take the at-the-money put of the May expiry, priced at the close of April 2 and of April 3 with Black-Scholes-Merton, the index, the curve, the dividend and the iv_modelo of its strike on each day:

ind = pd.read_csv("opciones3_indice.csv", parse_dates=["fecha"],
                  index_col="fecha")
plazos = np.array([1 / 12, 0.25, 0.5, 1.0, 2.0])
cols = ["tipo_1m", "tipo_3m", "tipo_6m", "tipo_1a", "tipo_2a"]

def estado(fecha, venc, strike, tipo):
    f, v = pd.Timestamp(fecha), pd.Timestamp(venc)
    Tf = (v - f).days / 365
    fila = cad[(cad["fecha"] == f) & (cad["vencimiento"] == v)
               & (cad["strike"] == strike)
               & (cad["tipo"] == tipo)].iloc[0]
    return {"S": float(ind.loc[f, "cierre"]), "T": Tf,
            "r": float(np.interp(Tf, plazos,
                                 ind.loc[f, cols].to_numpy(float)
                                 / 100)),
            "q": float(ind.loc[f, "rentabilidad_dividendo"]) / 100,
            "sig": float(fila["iv_modelo"]),
            "medio": (fila["bid"] + fila["ask"]) / 2}

def atribucion(f0, f1, venc, strike, tipo):
    a, b = estado(f0, venc, strike, tipo), estado(f1, venc, strike,
                                                  tipo)
    es_call = tipo == "C"
    g = griegas(a["S"], strike, a["T"], a["r"], a["q"], a["sig"],
                es_call)
    v1 = griegas(b["S"], strike, b["T"], b["r"], b["q"], b["sig"],
                 es_call)["precio"]
    dS, dsig = b["S"] - a["S"], 100 * (b["sig"] - a["sig"])
    dt = (pd.Timestamp(f1) - pd.Timestamp(f0)).days
    dr = 100 * (b["r"] - a["r"])
    partes = pd.Series({
        "delta": g["delta"] * dS, "gamma": 0.5 * g["gamma"] * dS**2,
        "vega": g["vega"] * dsig, "theta": g["theta"] * dt,
        "rho": g["rho"] * dr, "vanna": g["vanna"] * dS * dsig,
        "volga": 0.5 * g["volga"] * dsig**2,
        "charm": g["charm"] * dt * dS})
    partes["residual"] = (v1 - g["precio"]) - partes.sum()
    partes["total (model)"] = v1 - g["precio"]
    partes["total (mids)"] = b["medio"] - a["medio"]
    print(f"S {a['S']} to {b['S']} ({dS / a['S']:+.2%}); vol. "
          f"{a['sig']} to {b['sig']}; {dt} day(s)")
    return partes.round(2)

atribucion("2025-04-02", "2025-04-03", "2025-05-16", 5200.0, "P")

S 5211.08 to 4956.94 (-4.88%); vol. 0.2252 to 0.2314; 1 day(s)
delta            115.59
gamma             31.37
vega               4.44
theta             -1.62
rho               -0.01
vanna              0.09
volga              0.00
charm              0.07
residual          -0.06
total (model)    149.87
total (mids)     150.05
dtype: float64

The put earned 149.87 points in the model (150.05 with the mid screen prices: 15,005 USD per contract). Delta explains 115.59; gamma, 31.37 more, because with a 254-point drop the linear approximation falls short; vega, 4.44 (the volatility of that strike rose 0.62 points); theta subtracts 1.62. Rho, vanna, volga and charm do not reach 0.1. What is left unexplained is −0.06: the third-order terms. The figure shows it in bars:

That day the breakdown is the textbook one: a short position in that put with a delta hedge would have lost the gamma (31 points) plus the vega (4.4) and would only have collected 1.6 of theta. It is the risk of selling options in one sentence, and the reason why lesson 3.6.1 presents premium selling with its worst episode and not with its Sharpe.

Common errors

• Library units. QuantLib's vega and rho per unit, theta per year; py_vollib per point and per day; the volga of BlackCalculator with respect to $\sigma\sqrt{T}$. They are converted in a single place in the code and tested.
• Theta per business or calendar day. With $T$ in calendar days / 365, the weekend consumes three days of theta. If the system uses business days for $T$, the volatility has to be on the same basis.
• Finite differences with absurd steps. $h = 10^{-8}$ is not more precise: it is noise. For gamma, of the order of one ten-thousandth of the underlying; for vega, $10^{-4}$ in volatility.
• Greeks with the wrong volatility. A Greek is computed with the volatility of its strike. And when the index moves, the volatility of each strike moves too (in ORX it rose as the index fell): the “formula” delta assumes the strike stuck to its volatility (sticky strike) and does not include that effect. It is discussed with the surface in lesson 3.4.5.
• Dollar Greeks without the multiplier. A gamma of 0.00124 per point looks like nothing, but on a 100-USD contract on an index at 5,148 the dollar delta changes by about 640 USD per point (100 × 5,148 × 0.00124) and about 33,000 USD with a 1% move: the scale of the contract decides whether a Greek is small.

Exercises

1. Repeat the attribution for the at-the-money call and put of the May expiry between April 8 and 9, 2025, when ORX rose 9.5% and the volatility of those strikes rose almost 9 points. How much is left unexplained and why?

pd.DataFrame({
    "call": atribucion("2025-04-08", "2025-04-09", "2025-05-16",
                       4800.0, "C"),
    "put": atribucion("2025-04-08", "2025-04-09", "2025-05-16",
                      4800.0, "P")})

S 4820.64 to 5279.9 (+9.53%); vol. 0.3532 to 0.4403; 1 day(s)
S 4820.64 to 5279.9 (+9.53%); vol. 0.3532 to 0.4403; 1 day(s)
                 call     put
delta          251.29 -207.29
gamma           75.92   75.92
vega            53.58   53.58
theta           -3.04   -2.67
rho              0.01   -0.01
vanna           -0.29   -0.29
volga            0.01    0.01
charm           -0.10   -0.12
residual       -19.26  -19.92
total (model)  358.11 -100.80
total (mids)   358.80 -101.15

2. Check the Black-Scholes equation on all the 30-day calls and puts of the lesson's chain, each with its iv_modelo. Does it hold with different volatilities at each strike?

c30_all = cad[(cad["fecha"] == "2025-06-18")
              & (cad["vencimiento"] == "2025-07-18")]
max_err = 0.0
for tipo in ["C", "P"]:
    x = c30_all[c30_all["tipo"] == tipo]
    sx = x["iv_modelo"]
    gx = griegas(S, x["strike"], T, r, q, sx, call=tipo == "C")
    pde = (365 * gx["theta"] + 0.5 * sx**2 * S**2 * gx["gamma"]
           + (r - q) * S * gx["delta"] - r * gx["precio"])
    max_err = max(max_err, float(np.abs(pde).max()))
print(f"{len(c30_all)} options; maximum residual {max_err:.1e}")

90 options; maximum residual 3.1e-13

3. You have a portfolio of ORX options with delta 0, gamma of −0.05 per point, vega of −800 per volatility point and theta of +90 per day (all in index points, without the multiplier). If tomorrow the index falls 100 points and volatility rises 3 points, roughly how much do you lose? Which Greek dominates?