Delta hedging, gamma scalping and dispersion.

What remains after hedging the delta

A long call wins if the underlying rises. If at the same time the amount of underlying given by its delta is sold, the first order of the move disappears: what the call gains on a small rise, the sale loses. What remains is the second order, the curvature. Lesson 3.4.3 does the full attribution of the P&L of an ATM put into delta, gamma, vega and theta; here its result for the hedged portfolio (the option and the $-\Delta$ position in the underlying) is enough. With the Taylor expansion of the option value $V(S, \sigma, t)$, the P&L of one day is

In Black-Scholes (without rates or dividends, to see the idea) theta and gamma are tied by the pricing equation: $\Theta = -\tfrac12 \Gamma S^2 \sigma_i^2$, with $\sigma_i$ the implied volatility at which the option is priced. Substituting and calling $r_t = dS/S$ the day's return:

It is the central formula of the lesson. The first term says that a long, hedged option wins on the day the squared move exceeds the one the implied volatility was discounting ($\sigma_i^2 \Delta t$) and loses on the other days, when it pays theta with nothing in return. Summed over the life of the option, the result is the difference between realized and implied variance weighted by each day's $\Gamma S^2$. The second term is the effect of the implied volatility itself changing (lesson 3.4.3). Bakshi and Kapadia (2003) use this same decomposition to measure the volatility risk premium in real S&P 500 options: the gains of long, hedged options are negative on average.

flowchart TB
  A[Buy straddle] --> B[Compute delta]
  B --> C[Sell delta]
  C --> D{S moves}
  D -- a lot --> E[Earns gamma]
  D -- little --> F[Pays theta]
  E --> G[Rebalance delta]
  F --> G
  G --> B

The loop is what desks call gamma scalping: with positive gamma, the delta rises when the price rises and falls when it falls, so rebalancing the hedge forces you to sell after rallies and buy after drops. Each turn of the loop captures a piece of the move. Theta is the price of that privilege.

Conventions in this lesson: prices per index unit (the real contract multiplies by 100 USD), $T$ in calendar days / 365 (ORX options expire at the close), continuous rate and dividend yield, annual volatility as a fraction, vega per unit of volatility (1.00 = 100 points; divided by 100, per point) and each trade's P&L in % of the index at opening.

The data: one ATM straddle every month

The data are those of the fictitious ORX index, which mimics the S&P 500 (lesson 3.4.6): opciones3_indice.csv with the close, the rate curve and the dividend yield, and opciones3_cadena.parquet, the daily chain with the three nearest monthly expirations. They are synthetic data with a planted variance risk premium: the model surface prices variance 1.30 times above the expected one in calm periods and 1.00 at peak stress (public/curso/datos/nivel-3/README.md). That is the known answer the lesson has to recover.

import numpy as np
import pandas as pd
from scipy.special import ndtr
from scipy.stats import norm

ind = pd.read_csv("opciones3_indice.csv", parse_dates=["fecha"])
ind = ind.set_index("fecha")
cad = pd.read_parquet("opciones3_cadena.parquet")
cad["medio"] = (cad["bid"] + cad["ask"]) / 2
S = ind["cierre"]
r = ind["tipo_1m"] / 100                 # continuous, as a fraction
q = ind["rentabilidad_dividendo"] / 100  # continuous
print(len(cad), "rows,", cad["fecha"].nunique(), "sessions")
dia = cad[(cad["fecha"] == "2020-02-21")
          & (cad["vencimiento"] == "2020-03-20")]
k = dia.loc[(dia["strike"] - S["2020-02-21"]).abs().idxmin(),
            "strike"]
print("S =", S["2020-02-21"], " ATM strike =", k)
dia[dia["strike"] == k][["tipo", "bid", "ask", "iv_modelo"]]

544548 rows, 3017 sessions
S = 2978.34  ATM strike = 3000.0
       tipo   bid   ask  iv_modelo
275893    C  28.1  28.5     0.1191
275911    P  52.4  52.9     0.1191

On 21 February 2020 the index closed at 2,978.34 and the nearest strike was 3,000. The call and the put of that strike, with 28 days, have the same iv_modelo, the model's true implied volatility, without the noise of the quotes: 11.9%. In real data that column does not exist and the implied volatility is computed from the mid with the methods of lesson 3.4.4; here it is used so that the error of the implied volatility does not mix with what we want to measure.

The Greeks, with the forward $F = S e^{(r-q)T}$, come from the Black-Scholes-Merton formulas of lesson 3.4.2. Before using them, a check against QuantLib with that day's call:

import QuantLib as ql

def bs(tipo, S, K, T, iv, r, q):
    """Price, delta, gamma and vega (per unit of volatility)."""
    F = S * np.exp((r - q) * T)
    sd = iv * np.sqrt(T)
    d1 = (np.log(F / K) + 0.5 * sd * sd) / sd
    df = np.exp(-r * T)
    call = df * (F * ndtr(d1) - K * ndtr(d1 - sd))
    precio = np.where(tipo == "C", call, call - df * (F - K))
    delta = np.exp(-q * T) * (ndtr(d1) - (tipo == "P"))
    gamma = np.exp(-q * T) * norm.pdf(d1) / (S * sd)
    vega = S * np.exp(-q * T) * norm.pdf(d1) * np.sqrt(T)
    return precio, delta, gamma, vega

f0 = pd.Timestamp("2020-02-21")
mio = bs("C", S[f0], k, 28 / 365, 0.1191, r[f0], q[f0])

hoy = ql.Date(21, 2, 2020)
ql.Settings.instance().evaluationDate = hoy
dc = ql.Actual365Fixed()
curva = lambda x: ql.YieldTermStructureHandle(
    ql.FlatForward(hoy, float(x), dc))
vol = ql.BlackVolTermStructureHandle(
    ql.BlackConstantVol(hoy, ql.NullCalendar(), 0.1191, dc))
proc = ql.BlackScholesMertonProcess(
    ql.QuoteHandle(ql.SimpleQuote(float(S[f0]))),
    curva(q[f0]), curva(r[f0]), vol)
op = ql.EuropeanOption(ql.PlainVanillaPayoff(ql.Option.Call, k),
                       ql.EuropeanExercise(hoy + 28))
op.setPricingEngine(ql.AnalyticEuropeanEngine(proc))
pd.DataFrame({"own": [float(x) for x in mio],
              "QuantLib": [op.NPV(), op.delta(), op.gamma(),
                           op.vega()]},
             index=["price", "delta", "gamma", "vega"]).round(6)

              own    QuantLib
price   28.323279   28.323279
delta    0.408275    0.408275
gamma    0.003949    0.003949
vega   320.040712  320.040712

They match to six digits. The vega of 320.0 is per unit of volatility: one point (0.01) moves the call 3.20 index points, 320 USD per contract.

The trade is the one in the README's realism table: on the day of each monthly expiration an ATM straddle of the next expiration with at least 20 days is opened, the delta is hedged at the close of each session and it is settled at expiration at its intrinsic value. The first trade enters on 2 January 2014 and each of the following ones on the day the previous one expires; there are 143:

vencs = sorted(cad["vencimiento"].unique())
fechas = S.index
con_cadena = set(cad["fecha"])
ops, j = [], 0
while True:
    f = fechas[j]
    v = next((x for x in vencs if (x - f).days >= 20), None)
    if v is None or v > fechas[-1]:
        break
    if f not in con_cadena:
        j += 1
        continue
    sub = cad[(cad["fecha"] == f) & (cad["vencimiento"] == v)]
    kk = sub.loc[(sub["strike"] - S[f]).abs().idxmin(), "strike"]
    for d in fechas[(fechas >= f) & (fechas <= v)]:
        ops.append((f, v, kk, d))
    j = fechas.get_loc(v)      # the next one enters when this one expires

pan = pd.DataFrame(ops, columns=["entrada", "vencimiento",
                                 "strike", "fecha"])
pan["op"] = pan.groupby("entrada").ngroup()
pc = cad.pivot_table(index=["fecha", "vencimiento", "strike"],
                     columns="tipo", values=["medio", "iv_modelo"])
pc.columns = [f"{a}_{b}" for a, b in pc.columns]
pan = pan.merge(pc.reset_index(), how="left",
                on=["fecha", "vencimiento", "strike"])
for col, serie in (("S", S), ("r", r), ("q", q)):
    pan[col] = serie.reindex(pan["fecha"]).to_numpy()
pan["T"] = (pan["vencimiento"] - pan["fecha"]).dt.days / 365
print(pan["op"].nunique(), "trades,", len(pan), "rows")
print("without a quote:", pan["medio_C"].isna().sum(),
      "rows (expiration days and 2019-07-23)")

143 trades, 3153 rows
without a quote: 144 rows (expiration days and 2019-07-23)

First, the known answer

The README gives the result of the short version of this trade (collected at the mid, Black-Scholes delta with iv_modelo, cash at the 1-month rate and the hedge's dividends): 143 trades, mean P&L of 0.44% of the index, positive in 71% and worst trade of −11.8%, the one entering on 2020-02-21. Reproducing it is the proof that the panel and the Greeks are right:

def vendido_cubierto(g):
    g = g.reset_index(drop=True)
    kk, s0 = g["strike"][0], g["S"][0]
    caja, h = g["medio_C"][0] + g["medio_P"][0], 0.0
    for t in range(len(g) - 1):
        if np.isfinite(g["iv_modelo_C"][t]):  # there is a chain
            dl = sum(bs(tp, g["S"][t], kk, g["T"][t],
                        g[f"iv_modelo_{tp}"][t], g["r"][t],
                        g["q"][t])[1] for tp in ("C", "P"))
            caja -= (dl - h) * g["S"][t]   # buys the delta
            h = dl
        caja += h * g["S"][t] * g["q"][t] / 252   # dividends
        caja *= np.exp(g["r"][t] / 252)           # cash at 1 m
    s_fin = g["S"].iloc[-1]
    caja += h * s_fin - abs(s_fin - kk)           # settlement
    return 100 * caja / s0

corto = pan.groupby("op").apply(vendido_cubierto)
entradas = pan.groupby("op")["entrada"].first()
print(f"{len(corto)} trades; mean {corto.mean():.2f}%; "
      f"positive {(corto > 0).mean():.0%}")
print(f"worst {corto.min():.1f}% (entry "
      f"{entradas[corto.idxmin()].date()})")

143 trades; mean 0.44%; positive 71%
worst -11.8% (entry 2020-02-21)

The four figures match the README. Selling the straddle and hedging it wins on average, and wins almost three out of every four months. The worst trade wipes out in one blow what 27 average trades earn: it is the shape of premium selling from lesson 3.6.1, now with direction neutralized. What remains in the result is the difference between the volatility that was collected and the one that occurred.

The P&L decomposition, trade by trade

From here on the lesson looks at the long side, that of the gamma scalper, and without cash or dividends (they add less than 0.01 points per trade) so that each term is seen cleanly. For each day the actual P&L is computed (change in the straddle mid minus the delta times the change in the index; on the last day, the intrinsic value) together with the two terms of the formula: the gamma one, $\tfrac12 \Gamma (\Delta S^2 - S^2 \sigma_i^2 \Delta t)$, with $\Delta t$ in calendar days between sessions (the weekend counts three days of theta), and the vega one, $\mathcal{V} \Delta \sigma_i$ with the implied volatility of the same strike.

p = pan.dropna(subset=["medio_C"]).copy()
p["V"] = p["medio_C"] + p["medio_P"]
p["iv"] = p["iv_modelo_C"]          # same for the call and the put
args = [p[c].to_numpy() for c in ("S", "strike", "T", "iv",
                                  "r", "q")]
gc, gp = bs("C", *args), bs("P", *args)
p["delta"], p["gamma"] = gc[1] + gp[1], gc[2] + gp[2]
p["vega"] = gc[3] + gp[3]
fin = pan[pan["fecha"] == pan["vencimiento"]].copy()
fin["V"] = (fin["S"] - fin["strike"]).abs()   # intrinsic
a = (pd.concat([p, fin]).sort_values(["op", "fecha"])
     .reset_index(drop=True))
sig = a.groupby("op").shift(-1)       # the next session
dS = sig["S"] - a["S"]
dt = (sig["fecha"] - a["fecha"]).dt.days / 365
a["real"] = (sig["V"] - a["V"]) - a["delta"] * dS
a["gamma_pnl"] = 0.5 * a["gamma"] * (dS**2
                                     - a["S"]**2 * a["iv"]**2 * dt)
a["vega_pnl"] = (a["vega"] * (sig["iv"] - a["iv"])).fillna(0)
dias = a.dropna(subset=["real"])    # without the expiration row

s0 = pan.groupby("op")["S"].first()
por = (dias.groupby("op")[["real", "gamma_pnl", "vega_pnl"]].sum()
       .div(s0, axis=0) * 100)
por["resto"] = por["real"] - por["gamma_pnl"] - por["vega_pnl"]
print(por.describe().loc[["mean", "std", "min", "max"]].round(2))
c1 = por["real"].corr(por["gamma_pnl"])
c2 = por["real"].corr(por["gamma_pnl"] + por["vega_pnl"])
print(f"correlation with the actual: gamma {c1:.2f}; "
      f"gamma + vega {c2:.2f}")

       real  gamma_pnl  vega_pnl  resto
mean  -0.45      -0.64      0.35  -0.15
std    1.66       1.43      0.98   0.45
min   -3.46      -3.32     -1.81  -2.55
max   11.72       9.62      5.34   2.09
correlation with the actual: gamma 0.89; gamma + vega 0.97

The header figure is this table accumulated trade by trade: the actual P&L and the gamma term fall almost in parallel, month after month, and jump up in crises. The biggest jump is that of the straddle bought on 2020-02-21: 11.7 points of P&L, of which gamma explains 9.6. At the end of the sample the actual P&L accumulates −64.3 points (simple sum of the 143 trades) and the gamma term −92.1; the distance between the two lines is, above all, the vega term.

Three readings of the table:

• The mean of the gamma term is negative (−0.64% per trade). It is the planted variance risk premium: implied volatility tends to sit above realized, and whoever buys the option pays more theta than they recover with the moves. The mean actual P&L, −0.45%, mirrors the +0.44% of the short side (without cash and dividends).
• The vega term is not noise. Its mean is positive (+0.35%) because of a smile effect: the implied volatility of a fixed strike rises when the index moves away from it upward (the strike ends up below the money, in the expensive zone of the skew), and in a sample with the index rising that happens more often. With vega, the correlation with the actual P&L rises to 0.97.
• The residual picks up the higher-order terms (vanna and volga, lesson 3.4.3), the noise of quotes around the mid and the error of hedging only at the close. Its mean is −0.15 points and its standard deviation 0.45; in some trades it reaches 2.5 points.

Realized versus implied: the path matters

The pocket version of the formula says that the P&L is proportional to $\sigma_R^2 - \sigma_i^2$, the difference between the variance realized over the life of the option and the implied one at purchase. It is a good first approximation, but it assumes a constant gamma, and gamma is not constant: it is highest at the money and grows as expiration nears. The same 2% move is worth much more in the last week and at the strike than on the first day and with the index far away.

lr = np.log(S).diff()
filas = []
for o, g in pan.groupby("op"):
    f, v = g["entrada"].iloc[0], g["vencimiento"].iloc[0]
    rr = lr.loc[f:v].iloc[1:]
    rv = np.sqrt((rr**2).sum() / ((v - f).days / 365))
    filas.append((f.date(), g["iv_modelo_C"].iloc[0], rv))
e = pd.DataFrame(filas, columns=["entrada", "iv0", "rv"])
e["pnl"] = por["real"].to_numpy()
e["dif_var"] = e["rv"]**2 - e["iv0"]**2
print(f"realized > implied in {(e.rv > e.iv0).mean():.0%}"
      f" of trades")
print(f"P&L correlation with rv² − iv²: "
      f"{e['pnl'].corr(e['dif_var']):.2f}")
e.sort_values("pnl").iloc[[0, 1, 2, -3, -2, -1]].round(3)

realized > implied in 23% of trades
P&L correlation with rv² − iv²: 0.81
        entrada    iv0     rv     pnl  dif_var
81   2020-10-16  0.242  0.123  -3.460   -0.044
84   2021-01-15  0.222  0.116  -3.229   -0.036
83   2020-12-18  0.253  0.125  -2.707   -0.049
29   2016-06-17  0.132  0.271   3.290    0.056
134  2025-03-21  0.228  0.518   7.617    0.216
73   2020-02-21  0.119  0.857  11.719    0.720

Realized volatility exceeds implied in only 23% of trades, and the correlation of the P&L with the month's variance difference is 0.81, less than the 0.89 of the daily gamma term. The table shows the two extremes. The buyer's three worst trades are from late 2020 and early 2021: they were bought with implied volatilities of 22 to 25% after the crisis and the market moved 12%. The best ones are the crises: the 2020-02-21 entry was bought at 11.9% and the index moved at 86% annualized (+11.7%); the 2025-03-21 one was bought at 22.8% and realized volatility was 51.8% (+7.6%). What the buyer earns depends on the distance between the two, and in most months that distance is negative.

Hedging frequency: what is known and what the data measure

How often should the delta be rebalanced? The formula comes from hedging in continuous time; in practice you hedge every so often and each hedge costs bid-ask spread. To separate the two things, first a simulation with a known answer: a lognormal walk with realized volatility equal to implied (20%), a 28-day ATM call sold at its Black-Scholes price and hedged $N$ times over its life. With volatility correctly estimated, the expected value of the result is 0 whatever $N$ is; what changes is its dispersion:

rng = np.random.default_rng(362)
S0, K0, T0, sg, n_tray, pasos = 100.0, 100.0, 28 / 365, 0.20, 4000, 160

def call_bs(s, t):
    t = np.maximum(t, 1e-12)
    sd = sg * np.sqrt(t)
    d1 = (np.log(s / K0) + 0.5 * sd * sd) / sd
    return s * ndtr(d1) - K0 * ndtr(d1 - sd), ndtr(d1)

prima = call_bs(S0, T0)[0]
dtp = T0 / pasos
z = rng.standard_normal((n_tray, pasos))
cam = S0 * np.exp(np.cumsum(-0.5 * sg**2 * dtp
                            + sg * np.sqrt(dtp) * z, axis=1))
cam = np.c_[np.full(n_tray, S0), cam]
sim = []
for N in [1, 5, 20, 80, 160]:
    i = np.arange(0, pasos + 1, pasos // N)
    s_, t_ = cam[:, i], T0 - i * dtp
    d_ = call_bs(s_[:, :-1], t_[:-1])[1]
    pnl = (prima + (d_ * np.diff(s_, axis=1)).sum(axis=1)
           - np.maximum(cam[:, -1] - K0, 0))
    sim.append((N, pnl.mean(), pnl.std(), pnl.std() * np.sqrt(N)))
sim = pd.DataFrame(sim, columns=["coberturas", "media", "desv",
                                "desv_x_raizN"])
sim = sim.set_index("coberturas")
print(f"call premium: {prima:.3f}")
sim.round(3)

call premium: 2.210
            media   desv  desv_x_raizN
coberturas
1           0.020  1.657         1.657
5           0.016  0.814         1.820
20          0.011  0.428         1.916
80          0.004  0.212         1.899
160         0.001  0.149         1.890

The mean stays around zero in every row (between +0.001 and +0.02 points on a premium of 2.21: sampling error with 4,000 paths), and the standard deviation falls with $1/\sqrt{N}$: the last column, the deviation times $\sqrt N$, stays near 1.7-1.9. Going from 5 to 80 hedges (16 times more) divides the dispersion by 4. Each hedging error over an interval is $\tfrac12 \Gamma S^2 \sigma^2 \Delta t\,(Z^2 - 1)$, with variance proportional to $\Delta t^2$; summing $N$ independent intervals of length $T/N$, the total variance ends up at $N (T/N)^2 \propto 1/N$.

The data only have daily closes, so the maximum frequency is once a day: what happens within the session cannot be measured here. What can be compared is less frequent hedging and delta bands (rebalancing the hedge only when the portfolio delta drifts beyond a threshold):

def comprado(g, cada=1, banda=None, coste=0.0):
    """Long straddle; sells the delta every `cada` sessions or
    when it drifts beyond `banda`. `coste`: fraction of the
    notional traded in each hedge."""
    s, val, dl = (g[c].to_numpy() for c in ("S", "V", "delta"))
    h, pnl, n, cst = 0.0, -val[0], 0, 0.0
    for t in range(len(s) - 1):
        if np.isfinite(dl[t]):
            obj = -dl[t]
            toca = (t % cada == 0) if banda is None \
                else abs(obj - h) > banda
            if toca:
                cst += abs(obj - h) * s[t] * coste
                h, n = obj, n + 1
        pnl += h * (s[t + 1] - s[t])
    return 100 * (pnl + val[-1] - cst) / s[0], n

def resumen(**kw):
    res = [comprado(g, **kw) for _, g in a.groupby("op")]
    x = np.array([v for v, _ in res])
    return {"media": x.mean(), "desv": x.std(),
            "coberturas": np.median([n for _, n in res])}

reglas = {"every session": {}, "every 2": {"cada": 2},
          "every 5": {"cada": 5}, "every 10": {"cada": 10},
          "at open only": {"cada": 10_000},
          "band 0.05": {"banda": 0.05}, "band 0.10": {"banda": 0.1},
          "band 0.20": {"banda": 0.2}}
pd.DataFrame({k: resumen(**v) for k, v in reglas.items()}).T.round(2)

               media  desv  coberturas
every session  -0.45  1.66        20.0
every 2        -0.38  1.91        10.0
every 5        -0.17  3.06         4.0
every 10       -0.17  3.50         2.0
at open only   -0.03  3.63         1.0
band 0.05      -0.44  1.67        14.0
band 0.10      -0.43  1.75         9.0
band 0.20      -0.41  1.83         5.0

The standard deviation rises when hedging less: 1.66 points with the daily hedge, 3.06 every five sessions and 3.63 if the delta is hedged only at the open. The mean also changes (−0.45% daily, −0.17% every five, −0.03% without rebalancing), and it is tempting to read it as hedging less “being better”. It is not: hedging less leaves exposure to the direction of the index, and in a sample in which the index rises 12.5% a year that exposure has added up. It is a hidden directional bet, not an improvement of the hedge. Bands are an improvement: with a delta threshold of 0.20 the standard deviation is 1.83, almost that of the daily hedge, with 5 hedges per trade instead of 20.

The cost completes the picture. With the hedge done in an index future or in a basket, each rebalance pays bid-ask spread and commission on the notional traded. The figures that follow are assumptions (1 and 5 basis points of the notional traded), not prices from any market:

cost = {f"{k} · {c * 1e4:.0f} bp": resumen(coste=c, **v)
        for c in (0.0001, 0.0005)
        for k, v in (("every session", {}),
                     ("band 0.20", {"banda": 0.2}))}
pd.DataFrame(cost).T.round(2)

                      media  desv  coberturas
every session · 1 bp  -0.47  1.65        20.0
band 0.20 · 1 bp      -0.42  1.83         5.0
every session · 5 bp  -0.57  1.65        20.0
band 0.20 · 5 bp      -0.49  1.82         5.0

With a cost of 1 bp, the daily hedge subtracts about 2-3 bp per trade; with 5 bp, more than a tenth. Leland (1985) studied replication with proportional costs and revision at fixed intervals: shortening the interval reduces the replication error but makes the cost soar, and he proposed pricing the option with a volatility adjusted for those costs. In a one-month index straddle, with futures bid-ask spreads of a fraction of a basis point, the cost matters little; in illiquid single-stock options or with intraday hedging, a lot.

Gamma scalping: when the buyer wins

The gamma scalper buys the option and lives off rebalancing the hedge. The decomposition table already gave the average answer: they lose, because they pay the variance risk premium. They win when realized volatility exceeds implied, and that happens almost only at the start of crises. By terciles of VORX (the ORX volatility index, which mimics the VIX) on the entry day:

vx = pd.read_csv("vol3_indices.csv", parse_dates=["fecha"])
vx = vx.set_index("fecha")["vorx"]
e["vorx"] = vx.reindex(pd.to_datetime(e["entrada"])).to_numpy()
e["tercil"] = pd.qcut(e["vorx"], 3, labels=["low", "medium", "high"])
e.groupby("tercil", observed=True)["pnl"].agg(
    ["count", "mean", "median",
     ("positive", lambda s: (s > 0).mean())]).round(2)

        count  mean  median  positive
tercil
low        48  0.11   -0.30      0.38
medium     47 -0.83   -0.88      0.17
high       48 -0.64   -0.96      0.27

With VORX low at entry, the mean is positive (+0.11%) but the median is negative (−0.30%) and it wins only 38% of the time: the mean is lifted by the 2020-02-21 trade, bought with VORX at 15. With VORX medium or high, the buyer loses on average and wins in one out of every four or five. It is the profile of insurance: small losses almost every month and a large gain very occasionally, which arrives when volatility was cheap and the market breaks.

On 2020-02-21, gamma scalping is plain to see. The strike (3,000) was somewhat above the index, so the straddle started with a negative delta and the hedge was a small purchase of the index. As the index fell in March, the straddle's delta reached −0.89 and the hedge became a purchase of almost 0.9 units of index, rebalanced at each close: on each day of strong movement the position gains from the curvature and, when rebalancing the hedge, buys after drops and sells after rebounds. The four best days (each row's P&L is that of the interval starting that day):

d0 = dias[dias["op"] == e["pnl"].idxmax()]
dd = d0.set_index("fecha")[["S", "delta", "gamma_pnl", "real"]]
dd[["gamma_pnl", "real"]] *= 100 / d0["S"].iloc[0]
print(f"{len(dd)} days; actual P&L {dd['real'].sum():.2f}%; "
      f"from gamma {dd['gamma_pnl'].sum():.2f}%")
dd.sort_values("real").tail(4).round(3)

20 days; actual P&L 11.72%; from gamma 9.62%
                  S  delta  gamma_pnl   real
fecha
2020-03-11  2713.79 -0.891      1.404  0.753
2020-03-09  2803.60 -0.669      0.598  1.158
2020-03-10  2933.60 -0.223      2.697  2.286
2020-03-06  3034.07  0.205      4.636  5.556

The four intervals add up to 9.8 of the 11.7 points; the one of 2020-03-06 captures the 7.6% drop to Monday the 9th, with the delta still at +0.21 after the rebound on the 6th. Without those days the month would have been the buyer's usual one: theta paid day after day. That is why a long gamma portfolio is presented with its distribution, not with its mean, and why whoever sells it (lesson 3.6.1) must size the position to survive that month.

Term structure and skew

The two dimensions of the surface other than the level can also be traded with the delta hedged. The term structure compares implied volatilities of different expirations: a calendar sells the short expiration and buys the long one (lesson 3.5.3) and bets that the slope returns to its usual shape. The skew compares strikes: a delta-hedged risk reversal (buying the call and selling the put of the same expiration, or the other way around) bets that the difference in implied volatility between puts and calls is expensive or cheap relative to what the index will do.

opciones3_verdad.parquet has the true surface: the ATM total variance (theta) of each expiration. With it the ATM implied volatility at constant 30 and 60 days is interpolated (linear interpolation of total variance, lesson 3.4.5) and the slope is measured:

ver = pd.read_parquet("opciones3_verdad.parquet")
ver["T"] = ver["dias"] / 365
diaria = (ver.sort_values(["fecha", "vencimiento"])
          .groupby("fecha").head(3))

def atm_cte(g, dias):
    t, w = g["T"].to_numpy(), g["theta"].to_numpy()
    tt = dias / 365
    if t[0] >= tt:
        return np.sqrt(w[0] / t[0])
    b = min(np.searchsorted(t, tt), len(t) - 1)
    wi = w[b - 1] + (w[b] - w[b - 1]) * (tt - t[b - 1]) \
        / (t[b] - t[b - 1])
    return np.sqrt(wi / tt)

ts = diaria.groupby("fecha").apply(
    lambda g: pd.Series({"atm30": atm_cte(g, 30),
                         "atm60": atm_cte(g, 60)})) * 100
ts["pendiente"] = ts["atm60"] - ts["atm30"]
ts["vorx"] = vx.reindex(ts.index)
print(f"inverted (30 d > 60 d): {(ts.pendiente < 0).mean():.0%}"
      f" of sessions")
print(f"correlation of the slope with VORX: "
      f"{ts['pendiente'].corr(ts['vorx']):.2f}")
ts["pendiente"].describe().round(2)

inverted (30 d > 60 d): 28% of sessions
correlation of the slope with VORX: -0.67
count    3018.00
mean        0.38
std         1.21
min        -8.26
25%        -0.09
50%         0.56
75%         1.09
max         3.55
Name: pendiente, dtype: float64

The slope is positive almost whenever volatility is low (60-day implied volatility above 30-day: the structure in contango) and inverts in crises, with a correlation of −0.67 with VORX. It is inverted in 28% of sessions, in line with the 27% of Wednesdays the README measures between 1 month and 1 year. It is the model's known answer: short-term variance returns quickly to a slow mean, so after a shock the short implied volatility ends up well above the long one.

For the skew, the surface gives the implied volatility at any $k = \ln(K/F)$ with the SVI parameters. At about 30 days, the difference between the implied volatility 5% below and 5% above the forward, and the correlation between the daily index return and the change in the 30-day ATM implied volatility (the spot-vol correlation, which the skew discounts):

def iv_svi(g, k):
    w = g["svi_a"] + g["svi_b"] * (
        g["svi_rho"] * (k - g["svi_m"])
        + np.sqrt((k - g["svi_m"])**2 + g["svi_sigma"]**2))
    return np.sqrt(w / g["T"])

c30 = (ver[ver["dias"].between(20, 40)].drop_duplicates("fecha")
       .set_index("fecha"))
skew = 100 * (iv_svi(c30, -0.05) - iv_svi(c30, 0.05))
print(f"skew −5% / +5%: mean {skew.mean():.1f} points, "
      f"from {skew.min():.1f} to {skew.max():.1f}")
print(f"correlation of the skew with VORX: "
      f"{skew.corr(vx.reindex(skew.index)):.2f}")
rend = np.log(S).diff()
print(f"correlation return / change in the 30-day ATM: "
      f"{rend.corr(ts['atm30'].diff()):.2f}")

skew −5% / +5%: mean 9.1 points, from 6.8 to 11.9
correlation of the skew with VORX: 0.54
correlation return / change in the 30-day ATM: -0.62

The put 5% below the forward trades on average 9.1 volatility points above the call 5% above, and the difference grows with the volatility level (correlation of 0.54 with VORX). The measured spot-vol correlation is −0.62: when the index falls, implied volatility rises. The model plants it with a correlation of −0.72 between variance shocks and returns; ATM implied volatility is a smoothed average of future variance and picks up a smaller part.

With these data both trades can be tested, always delta-hedged: a short calendar when the structure is very inverted (a bet that it normalizes) or a risk reversal when the skew departs from its usual relationship with the level. Two warnings before doing it. First, the data have only three daily expirations and the full surface on Wednesdays: a calendar of more than three months can only be priced once a week. Second, both trades carry vega of opposite sign in each leg and move with the volatility level: a calendar sold “because the structure is inverted” is opened right in the crisis, when the level may keep rising, and loses on the short leg what the long leg does not earn.

Dispersion: what these data do not allow

The variance of an index depends on the variances of its components and on their correlations. With weights $w_i$ and implied volatilities $\sigma_i$ of the components, and a single average correlation $\overline{\rho}$:

It is the implied correlation. Cboe publishes S&P 500 implied correlation indices with this formula, with the ATM implied volatilities of the index and of its 50 largest components by float-adjusted capitalization (methodology white paper in the resources). A dispersion trade sells index volatility and buys that of the components, all delta-hedged: it wins if realized correlation ends up below implied. Driessen, Maenhout and Vilkov (2009) study options on the S&P 100 and on all its components and find that correlation risk is priced: it explains well the returns of index options versus those of stock options. According to their abstract, with realistic trading frictions that premium cannot be exploited.

flowchart TB
  I[Sell index vol] --> C[Hedged portfolio]
  A[Buy stock vol] --> C
  C --> R{Realized correlation}
  R -- low --> G[Wins]
  R -- high --> P[Loses: crisis]

The course data do not allow testing it: there are options on the ORX index, but not on its components (the KRN stock, the only one with options, is not part of the index and was not generated from it). What can be done is to compute the implied correlation as a formula. The numbers in the next cell are made up, just to see how it moves:

def correlacion_implicita(sig_indice, w, sig):
    w, sig = np.asarray(w), np.asarray(sig)
    sin_corr = np.sum(w**2 * sig**2)
    cruzados = np.sum(np.outer(w * sig, w * sig)) - sin_corr
    return (sig_indice**2 - sin_corr) / cruzados

w = [0.30, 0.25, 0.20, 0.15, 0.10]    # weights (made up)
sig = [0.28, 0.32, 0.25, 0.35, 0.40]  # implied volatilities of the stocks
for sig_i in (0.15, 0.20, 0.25):
    rho = correlacion_implicita(sig_i, w, sig)
    print(f"index at {sig_i:.0%}: implied correlation {rho:.2f}")

index at 15%: implied correlation 0.03
index at 20%: implied correlation 0.27
index at 25%: implied correlation 0.57

With the stocks at 25-40%, an index at 15% implies an average correlation of 0.03; at 20%, of 0.27, and at 25%, of 0.57: raising the index's implied volatility without the stocks' moving is raising the correlation. In a crisis, realized correlations soar toward 1, which is exactly when the seller of index volatility loses: dispersion is another way of selling premium with tail risk, and the costs (dozens of stock options, each with its bid-ask spread) eat up a good part of the premium.

Exercises

1. In the simulation, make realized volatility 30% with the call sold at 20%. Compare hedging with the delta of the implied volatility (20%) and with that of the realized one (30%) with 160 hedges: what happens to the mean and to the dispersion?

z2 = np.random.default_rng(7).standard_normal((n_tray, pasos))
s_r = 0.30
cam2 = S0 * np.exp(np.cumsum(-0.5 * s_r**2 * dtp
                             + s_r * np.sqrt(dtp) * z2, axis=1))
cam2 = np.c_[np.full(n_tray, S0), cam2]
t_ = T0 - np.arange(pasos + 1) * dtp

def delta_con(s, t, sv):
    sd = sv * np.sqrt(np.maximum(t, 1e-12))
    return ndtr((np.log(s / K0) + 0.5 * sd * sd) / sd)

for nombre, sv in (("implied 20%", 0.20), ("realized 30%", 0.30)):
    d_ = delta_con(cam2[:, :-1], t_[:-1], sv)
    pnl = (prima + (d_ * np.diff(cam2, axis=1)).sum(axis=1)
           - np.maximum(cam2[:, -1] - K0, 0))
    print(f"{nombre}: mean {pnl.mean():.3f}, std {pnl.std():.3f}")

implied 20%: mean -1.081, std 0.546
realized 30%: mean -1.099, std 0.224

2. Repeat the P&L decomposition only with the trades that enter with the term structure inverted (negative 30-60 day slope on the entry day). Does the buyer earn more or less than in the whole set?

e["pend"] = ts["pendiente"].reindex(
    pd.to_datetime(e["entrada"])).to_numpy()
e.groupby(e["pend"] < 0)["pnl"].agg(
    ["count", "mean", "median",
     ("positive", lambda s: (s > 0).mean())]).round(2)

       count  mean  median  positive
pend
False    102 -0.44   -0.75      0.25
True      41 -0.48   -0.84      0.32

3. A desk proposes selling index volatility and buying that of its ten largest components because the implied correlation is at 0.60 and the realized one over the last year was 0.35. List what can go wrong.