Smile, skew and term structure; volatility surface with SVI.

One volatility per strike and expiry

If Black-Scholes described the market well, all the options on an underlying would give the same implied volatility. They do not, and the way they do not has structure. For a fixed expiry, the implied volatility changes with the strike: that is the smile or, in equity indices, rather a smirk, the skew, because out-of-the-money puts trade much more expensive than equidistant calls. For a fixed relative strike, it changes with the tenor: that is the term structure. The two together form the volatility surface.

The cause of the skew is economic before it is mathematical. Indices fall faster than they rise and their volatility rises when they fall (the leverage effect), so the distribution of returns has a fat left tail; in addition, there is structural demand for puts as portfolio insurance. Both things make far puts more expensive than a lognormal would say. In the course data, the ORX index has precisely those two properties (−0.72 correlation between return and variance, and downward jumps), and its true surface has skew.

Three conventions for the whole lesson:

• Moneyness $k = \ln(K/F)$, with $F$ the forward of the expiry: $k < 0$ are strikes below the forward.
• Total variance $w(k, T) = \sigma^2_{\text{imp}}(k, T)\,T$, with $T$ in calendar days / 365. It is the natural magnitude for comparing tenors: the variance of a long tenor is the sum of those of its segments.
• Delta of 25 in the spot convention with continuous dividend, $e^{-qT}N(d_1)$ for the call; the “25-delta” put is the one with delta −0.25.

flowchart TB
  A[Weekly chain] --> B[OTM implied vol]
  B --> C[SVI per slice]
  C --> D{No arbitrage}
  D -- no --> C
  D -- yes --> E[Global SSVI]
  E --> F[Interpolation]
  F --> G[Usable surface]

The diagram is the order of the lesson. Each box is checked against the true surface, which in these data is known.

The surface on a Wednesday

opciones3_superficie.parquet is the Wednesday snapshot of the ORX chain with all the expiries up to two years: the six monthlies and the quarterlies. They are synthetic data; their true surface is a Gatheral and Jacquier (2014) SSVI whose parameters per slice are in opciones3_verdad.parquet. First, the implied volatility of each out-of-the-money option with a positive bid, with the safeguarded vectorized Newton of lesson 3.4.4. The forward and the rate of each slice are taken from the truth to separate the fitting problem from the parity one, which the previous lesson already solved:

import math
import numpy as np
import pandas as pd
from scipy.special import ndtr

sup = pd.read_parquet("opciones3_superficie.parquet")
ver = pd.read_parquet("opciones3_verdad.parquet")


def black(F, K, T, D, s, tipo):
    sd = s * np.sqrt(T)
    d1 = np.log(F / K) / sd + sd / 2
    c = D * (F * ndtr(d1) - K * ndtr(d1 - sd))
    return np.where(tipo == "C", c, c - D * (F - K))


def vega(F, K, T, D, s):
    sd = s * np.sqrt(T)
    d1 = np.log(F / K) / sd + sd / 2
    return D * F * np.exp(-d1 * d1 / 2) * np.sqrt(T / (2 * math.pi))


def iv_vec(p, F, K, T, D, tipo):
    """Newton with a bisection safeguard (lesson 3.4.4)."""
    lo, hi = np.full(len(p), 1e-4), np.full(len(p), 5.0)
    s = np.full(len(p), 0.3)
    for _ in range(60):
        d = black(F, K, T, D, s, tipo) - p
        lo, hi = np.where(d < 0, s, lo), np.where(d > 0, s, hi)
        paso = s - d / np.maximum(vega(F, K, T, D, s), 1e-12)
        mal = (paso <= lo) | (paso >= hi) | ~np.isfinite(paso)
        s = np.where(mal, (lo + hi) / 2, paso)
    return s


hoy = pd.Timestamp("2025-06-04")
x = sup[sup["fecha"] == hoy].merge(ver, on=["fecha", "vencimiento"])
x["T"] = x["dias"] / 365
x["k"] = np.log(x["strike"] / x["forward"])
otm = np.where(x["tipo"] == "C", x["k"] >= 0, x["k"] < 0)
x = x[otm & (x["bid"] > 0)].copy()
x["iv"] = iv_vec(((x["bid"] + x["ask"]) / 2).to_numpy(),
                 x["forward"].to_numpy(), x["strike"].to_numpy(),
                 x["T"].to_numpy(),
                 np.exp(-x["tipo_libre"] * x["T"]).to_numpy(),
                 x["tipo"].to_numpy())
x["w"] = x["iv"] ** 2 * x["T"]
x.groupby("dias").agg(opciones=("iv", "size"), k_min=("k", "min"),
                      k_max=("k", "max")).round(3)

      opciones  k_min  k_max
dias
16          32 -0.414  0.103
44          40 -0.765  0.185
72          39 -0.767  0.245
107         37 -0.689  0.329
135         39 -0.692  0.354
170         38 -0.694  0.352
198         18 -0.622  0.322
289         19 -0.629  0.315
379         20 -0.636  0.363
471         20 -0.642  0.356
562         20 -0.649  0.350
653         19 -0.655  0.343

Twelve expiries, from 16 to 653 days, with between 18 and 40 usable options per slice. The short slices reach $k = -0.77$, puts 53% below the forward: after the April 2025 fall the low strikes are still listed (a strike does not disappear until its expiry). Those beyond six months have half the strikes because their spacing is double.

Smile, skew and term structure

Before fitting anything, it is worth knowing how to read a smile with the three figures the market uses:

• ATM volatility: that of the forward, $k = 0$. It measures the level.
• 25-delta risk reversal: $RR_{25} = \sigma_{25C} - \sigma_{25P}$. It measures the tilt; in an index, negative.
• 25-delta butterfly: $BF_{25} = (\sigma_{25C} + \sigma_{25P})/2 - \sigma_{ATM}$. It measures the curvature.

At the money, the true total variance is the theta parameter of the truth: $\sigma_{ATM} = \sqrt{\theta/T}$. The term structure is that volatility for each expiry, and it is best seen by comparing different days:

def theta_en(g, T):
    """θ linearly interpolated in T, with θ(0) = 0."""
    return np.interp(T, np.r_[0, g["dias"] / 365], np.r_[0, g["theta"]])


def atm(fecha, T):
    """True ATM vol at tenor T."""
    g = ver[ver["fecha"] == pd.Timestamp(fecha)]
    return np.sqrt(theta_en(g, T) / T)


vx = pd.read_csv("vol3_indices.csv", parse_dates=["fecha"])
plazos = np.array([30, 91, 182, 365, 548]) / 365
filas = {}
for f in ["2017-06-14", "2025-06-04", "2020-03-18"]:
    vorx = vx.loc[vx["fecha"] == f, "vorx"].item()
    filas[f"{f} (VORX {vorx:.1f})"] = 100 * atm(f, plazos)
pd.DataFrame(filas, index=["1m", "3m", "6m", "1y", "18m"]).round(1)

     2017-06-14 (VORX 11.1)  2025-06-04 (VORX 25.9)  2020-03-18 (VORX 72.8)
1m                     10.3                    22.1                    63.0
3m                     12.4                    22.4                    48.3
6m                     14.0                    22.4                    38.1
1y                     15.5                    22.0                    30.3
18m                    16.4                    21.6                    27.0

In calm, in June 2017 with the VORX (the VIX-like index of the data, lesson 3.4.6) at 11.1, the structure rises from 10.3% at one month to 16.4% at 18 months: the market charges for the calm not lasting. In March 2020 it inverts, from 63.0% at one month to 27.0% at 18 months: the crisis is expected to pass. In June 2025, two months after the April fall, it is almost flat, around 22%. Behind it is the mean reversion of variance: short tenors reflect today's variance and long ones, its long-term level. In the data model, the ATM total variance combines the instantaneous variance, a slow level and a long-term mean, multiplied by a variance risk premium.

The figure draws the three curves from 10 to 630 days; the March 2020 one, highlighted, falls to less than half between one month and one year.

mi = ver[ver["fecha"].dt.weekday == 2]
tabla = mi.groupby("fecha").apply(
    lambda g: pd.Series({"un_mes": np.sqrt(theta_en(g, 30 / 365)
                                           * 365 / 30),
                         "un_ano": np.sqrt(theta_en(g, 1.0))}))
invertida = tabla["un_mes"] > tabla["un_ano"]
print(f"{len(tabla)} Wednesdays; inverted: {invertida.mean():.1%}")
print(tabla[invertida].index.year.value_counts().sort_index()
      .to_dict())

619 Wednesdays; inverted: 23.4%
{2014: 2, 2015: 12, 2016: 17, 2017: 6, 2018: 19, 2019: 6, 2020: 20, 2021: 8, 2022: 28, 2023: 2, 2024: 7, 2025: 18}

23.4% of the 619 Wednesdays have the one-month ATM volatility above the one-year one (with θ linearly interpolated in tenor). They are concentrated in the stress years of the script: 28 in 2022, 20 in 2020, 19 in 2018 and 18 in 2025. A strategy that sells short-dated volatility and buys long-dated (a calendar, lesson 3.5.3) earns with contango and loses all at once when the curve inverts.

SVI: five parameters per expiry

The SVI (stochastic volatility inspired) parameterization, devised at Merrill Lynch in 1999 and presented publicly by Gatheral in 2004 (as Gatheral and Jacquier, 2014, recount), describes the total variance of a slice with five parameters:

Its wings are straight in total variance, as Lee's moment formula requires: far from the money, $w(k)$ cannot grow faster than $2|k|$. The fit is done by least squares on the volatility (not on the variance, which would weight the wings too much), with bounds so that the parameters do not leave their domain:

from scipy.optimize import least_squares


def svi(k, a, b, rho, m, sig):
    return a + b * (rho * (k - m) + np.sqrt((k - m) ** 2 + sig ** 2))


LIM = ([-1, 1e-6, -0.999, -1, 1e-4], [1, 5, 0.999, 1, 2])


def ajusta_svi(k, iv, T):
    def resid(p):
        return np.sqrt(np.maximum(svi(k, *p), 1e-10) / T) - iv
    p0 = [0.5 * (iv ** 2 * T).min(), 0.1, -0.5, 0.0, 0.1]
    return least_squares(resid, p0, bounds=LIM).x


cols = ["svi_a", "svi_b", "svi_rho", "svi_m", "svi_sigma"]
ajustes, filas = {}, []
for dias, g in x.groupby("dias"):
    k, T = g["k"].to_numpy(), g["T"].iloc[0]
    p = ajusta_svi(k, g["iv"].to_numpy(), T)
    ajustes[dias] = p
    verdad = g[cols].iloc[0].to_numpy()
    iv_v = np.sqrt(svi(k, *verdad) / T)
    filas.append({"dias": dias, "rho": p[2], "rho_v": verdad[2],
                  "b": p[1], "b_v": verdad[1],
                  "ruido_max": 100 * np.abs(g["iv"] - iv_v).max(),
                  "error_max": 100 * np.abs(
                      np.sqrt(svi(k, *p) / T) - iv_v).max()})
pd.DataFrame(filas).set_index("dias").round(3)

        rho  rho_v      b    b_v  ruido_max  error_max
dias
16   -0.673 -0.662  0.021  0.022      0.338      0.091
44   -0.665 -0.662  0.037  0.037      0.412      0.087
72   -0.678 -0.662  0.046  0.047      0.132      0.098
107  -0.656 -0.662  0.058  0.057      0.103      0.044
135  -0.660 -0.662  0.064  0.064      0.054      0.030
170  -0.661 -0.662  0.072  0.072      0.044      0.021
198  -0.657 -0.662  0.078  0.077      0.027      0.013
289  -0.661 -0.662  0.092  0.092      0.026      0.008
379  -0.657 -0.662  0.105  0.104      0.023      0.019
471  -0.688 -0.662  0.112  0.114      0.030      0.014
562  -0.655 -0.662  0.123  0.123      0.018      0.012
653  -0.634 -0.662  0.133  0.131      0.041      0.007

The _v columns are the truth. The fitted ρ ranges from −0.634 to −0.688, around the true −0.662, which is the same in all the slices because the surface is an SSVI; $b$ comes out almost exact. ruido_max is the largest distance between the volatility of an option and the truth: up to 0.41 points in the 44-day slice, the quote noise seen in lesson 3.4.4. error_max is that of the fitted SVI: 0.087 in that same slice and less than 0.1 in all of them. A five-parameter model cannot follow the noise of forty quotes: it acts as a filter.

The header figure draws five of the twelve fitted slices, with the volatility of each option of the 16-day slice as a dot. The 16-day one, highlighted, is the steepest in $k$: it falls from 62% to 16% in less than half a point of moneyness, while the 653-day one takes almost a whole point to fall from 35% to 16%. Measured in standard deviations, the difference almost disappears: $k = -0.41$ at 16 days is 9.0 standard deviations ($\sqrt{\theta} = 0.046$), and $k = -0.66$ at 653 days, only 2.3 ($\sqrt{\theta} = 0.286$). In SSVI, the ATM slope of the volatility is $\partial\sigma/\partial k = \rho\,\varphi(\theta)\sqrt{\theta}/(2\sqrt{T}) \approx \rho\,\eta/(2\sqrt{T})$: it grows like $1/\sqrt{T}$. At 16 days it is −1.5 volatility points for each 1% of $k$.

With one SVI per slice, the market measures can now be computed at any delta. The 25-delta strike is found with Brent on $k$, because the volatility of each strike depends on the strike itself:

from scipy.optimize import brentq

S = pd.read_csv("opciones3_indice.csv", parse_dates=["fecha"]) \
    .set_index("fecha").loc[hoy, "cierre"]


def strike_delta(p, F, T, r, objetivo):
    """k with spot delta e^(-qT)·N(d1) (call) or −e^(-qT)·N(−d1)."""
    eq = math.exp(-r * T) * F / S          # e^(-qT)
    def f(k):
        sd = math.sqrt(svi(k, *p))
        d1 = -k / sd + sd / 2
        delta = eq * (ndtr(d1) if objetivo > 0 else ndtr(d1) - 1)
        return delta - objetivo
    return brentq(f, -1.5, 1.0)


g = x[x["dias"] == 44].iloc[0]
T44 = g["T"]
for nombre, p in [("fit", ajustes[44]),
                  ("truth", g[cols].to_numpy(float))]:
    vol = lambda k: 100 * math.sqrt(svi(k, *p) / T44)
    k25p = strike_delta(p, g["forward"], T44, g["tipo_libre"], -0.25)
    k25c = strike_delta(p, g["forward"], T44, g["tipo_libre"], 0.25)
    print(f"{nombre}: ATM {vol(0):.2f}, 25P {vol(k25p):.2f} "
          f"(K = {g['forward'] * math.exp(k25p):.0f}), "
          f"25C {vol(k25c):.2f}, RR {vol(k25c) - vol(k25p):.2f}, "
          f"BF {(vol(k25c) + vol(k25p)) / 2 - vol(0):.2f}")

fit: ATM 22.21, 25P 27.52 (K = 4856), 25C 18.52, RR -9.01, BF 0.81
truth: ATM 22.22, 25P 27.52 (K = 4856), 25C 18.51, RR -9.01, BF 0.80

The ATM volatility comes out at 22.21%, against the true 22.22%. The 25-delta put, at the 4,856 strike (5.8% below the forward), is at 27.52%, and the 25-delta call, at 18.52%: a risk reversal of −9.01 points and a butterfly of 0.81, against −9.01 and 0.80 for the truth. The skew of an index is read in the risk reversal: the 25-delta put costs 9 volatility points more than the 25-delta call.

Static arbitrage: butterfly and calendar

A surface can fit the observed prices very well and still be impossible. Two conditions make it consistent with arbitrage-free prices:

• Butterfly (per slice). The second derivative of the call price with respect to the strike is the discounted risk-neutral density (Breeden and Litzenberger): a butterfly with strikes $K - h, K, K + h$ cannot be worth less than zero. Written in terms of $w(k)$, the condition is that Gatheral's function be non-negative:

• Calendar (between slices). With the forward as reference, the total variance has to increase with tenor at fixed $k$: $\partial w / \partial T \ge 0$. Otherwise, a long call would be worth less than an equivalent short one.

def g_svi(k, a, b, rho, m, sig):
    r = np.sqrt((k - m) ** 2 + sig ** 2)
    w = svi(k, a, b, rho, m, sig)
    w1 = b * (rho + (k - m) / r)
    w2 = b * sig ** 2 / r ** 3
    return (1 - k * w1 / (2 * w)) ** 2 \
        - w1 ** 2 / 4 * (1 / w + 0.25) + w2 / 2


def controles(ajustes, kk=np.linspace(-1.5, 0.6, 421)):
    dias = sorted(ajustes)
    w = np.array([svi(kk, *ajustes[d]) for d in dias])
    a, b, rho, m, sig = np.array([ajustes[d] for d in dias]).T
    return pd.DataFrame({
        "g_min": [g_svi(kk, *ajustes[d]).min() for d in dias],
        "w_min": a + b * sig * np.sqrt(1 - rho ** 2),
        "lee": b * (1 + np.abs(rho)),
        "dw_min_sig": np.r_[np.diff(w, axis=0).min(axis=1), np.nan]},
        index=pd.Index(dias, name="dias"))


controles(ajustes).round(4)

       g_min   w_min     lee  dw_min_sig
dias
16    0.2604  0.0012  0.0360      0.0019
44    0.2668  0.0033  0.0609      0.0019
72    0.2722  0.0055  0.0780      0.0025
107   0.2771  0.0083  0.0952      0.0019
135   0.2807  0.0105  0.1068      0.0023
170   0.2842  0.0132  0.1196      0.0018
198   0.2870  0.0152  0.1285      0.0058
289   0.2950  0.0219  0.1528      0.0055
379   0.3005  0.0282  0.1732      0.0048
471   0.3084  0.0343  0.1887      0.0058
562   0.3130  0.0403  0.2038      0.0056
653   0.3177  0.0463  0.2171         NaN

The twelve slices fitted with all the strikes pass the checks over the whole range from $k = -1.5$ to $0.6$, which goes beyond the listed strikes: minimum $g$ between 0.26 and 0.32, positive minimum variance, Lee's condition well below 2 and each slice below the next one (dw_min_sig, the smallest difference in total variance with the next slice, always positive).

What happens if the fit uses little data? A common case is keeping only the strikes close to the money, the most liquid ones, because the wings “are very noisy”:

cerca = {}
for dias, g in x[np.abs(x["k"]) < 0.06].groupby("dias"):
    cerca[dias] = ajusta_svi(g["k"].to_numpy(), g["iv"].to_numpy(),
                             g["T"].iloc[0])
c = controles(cerca)
print(c.round(4).to_string())
print("Slices with g < 0:", list(c.index[c["g_min"] < 0]))
print("With negative minimum variance:", list(c.index[c["w_min"] < 0]))
print("Crossing the next one:", list(c.index[c["dw_min_sig"] < 0]))

       g_min   w_min     lee  dw_min_sig
dias
16    0.2606  0.0012  0.0358      0.0025
44   -0.2289  0.0036  0.1778     -0.0798
72    0.2738  0.0044  0.0771      0.0031
107   0.2832  0.0075  0.0905      0.0021
135   0.2868  0.0098  0.1017      0.0029
170   0.2940  0.0127  0.1109      0.0035
198   0.3046  0.0220  0.1136      0.0100
289   0.1305  0.0320  0.2377     -0.1163
379   0.3270  0.0362  0.1485     -0.0152
471   0.2501  0.0106  0.2583     -0.0051
562   0.2320 -0.0132  0.2964     -0.0243
653   0.2779 -0.0058  0.2581         NaN
Slices with g < 0: [44]
With negative minimum variance: [562, 653]
Crossing the next one: [44, 289, 379, 471, 562]

With only 5 to 13 strikes per slice, the 44-day one has $g < 0$ (−0.23), the two longest have negative minimum variance (an imaginary volatility at some strike) and five cross the next one somewhere in the range. In the center they fit almost equally well: the problem is where there is no data, which is precisely where the surface is used to value far options or to compute a volatility index.

The figure translates $g$ into the risk-neutral density of $k$ in the 44-day slice. With all the strikes it is an asymmetric bell, positive over the whole range and with a long left tail (the skew). With the nearby strikes it is almost the same in the center, but it turns negative between $k = 0.17$ and $0.35$, with a minimum of −0.16: a negative probability that the index ends between 19% and 42% above the forward. A butterfly centered there would have a negative price.

SSVI: the whole surface with two parameters

SVI fits each slice separately, and nothing forces the slices to be consistent with each other. The SSVI (surface SVI) of Gatheral and Jacquier (2014) parameterizes the whole surface with the ATM total variance of each tenor, $\theta_T$, and a common shape function:

With $\theta_T$ increasing in $T$ there is no calendar arbitrage, and with $\eta\,(1 + |\rho|) \le 2$ (and $\gamma = 1/2$) no butterfly arbitrage either: it is a sufficient condition that can be checked in one line. Each slice of an SSVI is a raw SVI, which is why the parameters of the truth come in SVI form. The fit uses as $\theta_T$ the ATM variance of the SVI of each slice and looks for common $\rho$, $\eta$ and $\gamma$:

def w_ssvi(k, th, rho, eta, gam=0.5):
    phi = eta / (th ** gam * (1 + th) ** (1 - gam))
    return th / 2 * (1 + rho * phi * k
                     + np.sqrt((phi * k + rho) ** 2 + 1 - rho ** 2))


theta = {d: svi(0.0, *p) for d, p in ajustes.items()}
grupos = [(theta[d], g["k"].to_numpy(), g["iv"].to_numpy(),
           g["T"].iloc[0]) for d, g in x.groupby("dias")]


def resid(p):
    return np.concatenate([np.sqrt(w_ssvi(k, th, *p) / T) - iv
                           for th, k, iv, T in grupos])


fit = least_squares(resid, [-0.5, 0.5, 0.4],
                    bounds=([-0.999, 0.01, 0.01], [0.999, 5, 0.99]))
rho, eta, gam = fit.x
print(f"ρ = {rho:.4f} (truth {x['svi_rho'].iloc[0]:.4f}), "
      f"η = {eta:.3f} (0.95), γ = {gam:.3f} (0.5)")
print(f"η·(1 + |ρ|) = {eta * (1 + abs(rho)):.3f} <= 2")
ecm = 100 * np.sqrt(np.mean(fit.fun ** 2))
print(f"Root mean squared error: {ecm:.3f} points")

ρ = -0.6617 (truth -0.6616), η = 0.953 (0.95), γ = 0.499 (0.5)
η·(1 + |ρ|) = 1.584 <= 2
Root mean squared error: 0.037 points

SSVI recovers ρ = −0.6617 against −0.6616, η = 0.953 against 0.95 and γ = 0.499 against 0.5 (the values of the data model, in its README), with a root mean squared error of 0.037 points. Three shape parameters plus the twelve $\theta$ describe the whole surface. $\eta(1 + |\rho|) = 1.584 < 2$: no butterfly arbitrage by construction. And the fitted $\theta$ increase with tenor in the twelve slices (it is the ATM variance of each SVI, already checked when running the controls), so there is no calendar arbitrage at the money either.

Interpolating between expiries

The listed expiries are only a few; strategies and indices ask for constant tenors (30 days for the VIX, lesson 3.4.6) or the expiry of a custom option. You have to interpolate, and the way you do it matters. The cleanest test is to remove a known slice, the 44-day one, rebuild it from its 16- and 72-day neighbors and compare it with its truth. Four methods, with the weight $\alpha = (T - T_1)/(T_2 - T_1)$:

1. Implied volatility linear in $T$, at fixed $k$.
2. Total variance linear in $T$, at fixed $k$.
3. Total variance linear in $T$, at fixed standard deviations $z = k/\sqrt{\theta_T}$: each neighbor is evaluated at its own $k = z\sqrt{\theta_{T_i}}$.
4. SSVI with $\theta_T$ linearly interpolated and the fitted $\rho$, $\eta$, $\gamma$.

g = x[x["dias"] == 44]
k, T = g["k"].to_numpy(), g["T"].iloc[0]
verdad = np.sqrt(svi(k, *g[cols].iloc[0].to_numpy(float)) / T)
T1, T2 = 16 / 365, 72 / 365
p1, p2 = ajustes[16], ajustes[72]
al = (T - T1) / (T2 - T1)
th = (1 - al) * theta[16] + al * theta[72]
z = k / np.sqrt(th)

metodos = {
    "vol at fixed k": (1 - al) * np.sqrt(svi(k, *p1) / T1)
                      + al * np.sqrt(svi(k, *p2) / T2),
    "variance at fixed k": np.sqrt(((1 - al) * svi(k, *p1)
                                    + al * svi(k, *p2)) / T),
    "variance at fixed k/sqrt(θ)": np.sqrt((
        (1 - al) * svi(z * np.sqrt(theta[16]), *p1)
        + al * svi(z * np.sqrt(theta[72]), *p2)) / T),
    "SSVI": np.sqrt(w_ssvi(k, th, rho, eta, gam) / T)}
err = pd.DataFrame({m: 100 * np.abs(v - verdad)
                    for m, v in metodos.items()})
err.agg(["mean", "max"]).T.round(3)

                              mean    max
vol at fixed k               1.499  4.914
variance at fixed k          0.546  1.786
variance at fixed k/sqrt(θ)  0.062  0.075
SSVI                         0.076  0.116

Interpolating the volatility at fixed $k$ is off by up to 4.9 points, and the total variance at fixed $k$, by up to 1.8. The reason is the one in the header figure: the shape of the smile scales with the standard deviation, and between 16 and 72 days the scale changes a lot. At fixed $k$, points of the two smiles that do not correspond are mixed: the same $k$ is the far wing at 16 days and almost the center at 72. At fixed $z = k/\sqrt{\theta}$, the maximum error falls to 0.075 points, and with SSVI to 0.12; the error left in SSVI comes from interpolating $\theta$ linearly in $T$ (0.005990 against the true 0.005950), because the true ATM variance is not linear in tenor.

The figure ranks the four methods by their maximum error, on a log scale: between the worst and the best there is a factor of 65.

What changes in a real system

• The forward rules. The whole surface is built in $k = \ln(K/F)$: with a badly estimated forward (dividends, rates) the smile shifts and false parity violations appear. Parity on the mids (lesson 3.4.4) is the standard way to obtain it.
• Fit weights. Fitting in volatility with equal weights gives as much importance to a far put with a spread of 2 volatility points as to an ATM one with 0.1. It is common to weight by the inverse of the spread in volatility or by the vega.
• Calibrating with constraints. In production, checking arbitrage after the fact is not enough: it is imposed in the fit (bounds like those of Lee and Gatheral-Jacquier, or an SSVI fit followed by a per-slice SVI that does not stray from it).
• The surface moves. Between two days, the smile shifts with the index differently depending on the regime (sticky strike or sticky delta); a skew strategy (lesson 3.6.2) depends on which one holds.

Exercises

1. Check Lee's condition, $b(1 + |\rho|) \le 2$, and that of non-negative minimum variance in the twelve slices fitted with all the strikes. Which slice comes closest to Lee's bound and why?

c = controles(ajustes)
print(c[["lee", "w_min"]].round(4).T.to_string())
print("Lee maximum:", round(c["lee"].max(), 4), "at",
      c["lee"].idxmax(), "days")

dias      16      44      72      107     135     170     198     289     379     471     562     653
lee    0.0360  0.0609  0.0780  0.0952  0.1068  0.1196  0.1285  0.1528  0.1732  0.1887  0.2038  0.2171
w_min  0.0012  0.0033  0.0055  0.0083  0.0105  0.0132  0.0152  0.0219  0.0282  0.0343  0.0403  0.0463
Lee maximum: 0.2171 at 653 days

2. Repeat the SVI fit of the 44-day slice weighting each residual by the inverse of the half-spread in volatility (compute it as in lesson 3.4.4) and compare the maximum error against the truth with that of the unweighted fit.

g = x[x["dias"] == 44]
k, T, iv = g["k"].to_numpy(), g["T"].iloc[0], g["iv"].to_numpy()
D = np.exp(-g["tipo_libre"].iloc[0] * T)
vg = vega(g["forward"].to_numpy(), g["strike"].to_numpy(), T, D, iv)
semi = ((g["ask"] - g["bid"]) / 2).to_numpy() / vg
peso = 1 / np.maximum(semi, 1e-4)


def resid_p(p):
    return peso * (np.sqrt(np.maximum(svi(k, *p), 1e-10) / T) - iv)


pp = least_squares(resid_p, [0.001, 0.1, -0.5, 0.0, 0.1],
                   bounds=LIM).x
verdad = np.sqrt(svi(k, *g[cols].iloc[0].to_numpy(float)) / T)
for nombre, p in [("unweighted", ajustes[44]), ("weighted", pp)]:
    e = 100 * np.abs(np.sqrt(svi(k, *p) / T) - verdad)
    print(f"{nombre}: maximum error {e.max():.3f}, mean {e.mean():.3f}")

unweighted: maximum error 0.087, mean 0.020
weighted: maximum error 0.026, mean 0.010

3. A colleague interpolates the 30-day volatility of the 4,500-strike put by taking the volatility of that same put at the 16- and 44-day expiries and averaging them weighted by days. Explain the two errors of the method and how you would do it.