--- title: "3.4.5 Volatility smile, skew and surface with SVI" description: "Smile, skew and term structure; SVI and SSVI fitted against the true surface, arbitrage checks and interpolation between expiries." url: "https://pythonparatrading.com/en/python-trading-course/level-3/volatility-smile-skew-svi-surface/" lang: "en" image: "https://pythonparatrading.com/og/en/python-trading-course/level-3/volatility-smile-skew-svi-surface.png" --- # 3.4.5 · Smile, skew and term structure; volatility surface with SVI > An option chain does not give one volatility, but hundreds: one per strike and expiry. This lesson arranges them into a surface, fits it with SVI and SSVI on synthetic data whose true surface is known, checks that the fit leaves no arbitrage and measures how wrong each way of interpolating between expiries is. - Course: [Level 3: Derivatives: futures, spreads, forex and options](https://pythonparatrading.com/en/python-trading-course/level-3/) - Module: 3.4 · Option theory and pricing - Duration: 75 min - Prerequisites: [3.4.2 · Binomial trees, Black-Scholes-Merton, Black-76 and Monte Carlo](https://pythonparatrading.com/en/python-trading-course/level-3/binomial-black-scholes-black-76-monte-carlo/), [3.4.3 · Greeks: delta, gamma, vega, theta, rho, vanna, volga and charm](https://pythonparatrading.com/en/python-trading-course/level-3/first-and-second-order-option-greeks/), [3.4.4 · Implied volatility: Newton-Raphson, Brent, py_vollib and QuantLib](https://pythonparatrading.com/en/python-trading-course/level-3/implied-volatility-newton-brent-quantlib/) ## By the end of this lesson - Describe a smile with its market measures (ATM volatility, 25-delta risk reversal and butterfly) and the term structure in calm and in stress. - Fit the raw SVI parameterization to each expiry and SSVI to the whole surface, and compare the parameters with the true ones. - Check the absence of butterfly arbitrage (positive density) and calendar arbitrage (increasing total variance), and see how a badly posed fit violates them. - Interpolate the surface between expiries and choose the method by its measured error, not by habit. ## 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. ```mermaid 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: ```python 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) ``` ```output 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: ```python 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) ``` ```output 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. *Fig. — Term structure of the true ATM volatility of ORX on three Wednesdays: in calm markets it rises with maturity; under stress it inverts.* Three ATM volatility curves against maturity, from 10 to 730 days. June 2017: rises from 10% at one month to 17% at two years. June 2025: almost flat, around 22%. March 2020 (highlighted): falls from 63% at one month to 30% at one year and 25% at two years. 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. ```python 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()) ``` ```output 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: $$ w(k) = a + b\left[\rho\,(k - m) + \sqrt{(k - m)^2 + \sigma^2}\right]. $$ | Parameter | What it moves | Constraint | |---|---|---| | $a$ | The level of the variance | $a + b\,\sigma\sqrt{1-\rho^2} \ge 0$ (non-negative minimum variance) | | $b$ | The slope of the wings, $b(1 \pm \rho)$ | $b \ge 0$; by Lee (2004), $b(1 + \lvert\rho\rvert) \le 2$ | | $\rho$ | The asymmetry: negative, the put wing steeper | $-1 < \rho < 1$ | | $m$ | Shifts the smile in $k$ | — | | $\sigma$ | The curvature at the minimum (rounds the “V”) | $\sigma > 0$ | 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: ```python 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) ``` ```output 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: ```python 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}") ``` ```output 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: $$ g(k) = \left(1 - \frac{k\,w'}{2w}\right)^2 - \frac{w'^2}{4}\left(\frac{1}{w} + \frac{1}{4}\right) + \frac{w''}{2} \ \ge 0. $$ - **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. ```python 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) ``` ```output 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”: ```python 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])) ``` ```output 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. *Fig. — Risk-neutral density of ln(S_T/F) implied by the SVI of the 44-day slice: fitted with all strikes and only with the near-the-money ones, which turns it negative on the call wing.* Two density curves against k, from −0.8 to 0.35. With all strikes, a skewed bell with a long left tail, positive everywhere. With strikes within 6% of the forward (highlighted), almost the same in the centre, but it falls below zero between k = 0.17 and 0.35, with a minimum of −0.16: a negative density, a butterfly arbitrage. 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: $$ w(k, \theta_T) = \frac{\theta_T}{2}\left[1 + \rho\,\varphi(\theta_T)\,k + \sqrt{(\varphi(\theta_T)\,k + \rho)^2 + 1 - \rho^2}\right], \qquad \varphi(\theta) = \frac{\eta}{\theta^{\gamma}(1 + \theta)^{1-\gamma}}. $$ 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$: ```python 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") ``` ```output ρ = -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$. ```python 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) ``` ```output 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. *Fig. — Error when rebuilding the 44-day slice from the 16 and 72-day ones, against its truth, in volatility points (log scale): four ways to interpolate; z = k/√θ, moneyness in ATM standard deviations.* Horizontal bars of the maximum error. Volatility linear at fixed k: 4.91 points (mean 1.50). Total variance at fixed k: 1.79 (mean 0.55). SSVI with interpolated θ: 0.12 (mean 0.08). Total variance at fixed z = k/√θ (highlighted): 0.08 (mean 0.06). 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? **Solution.** ```python 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") ``` ```output 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 ``` Lee's condition grows with tenor, from 0.036 at 16 days to 0.217 at 653, always very far from 2, and the minimum variance is positive in all the slices. The reason lies in SSVI: each slice is an SVI with $b = \theta\,\varphi(\theta)/2 = \eta\sqrt{\theta}/(2\sqrt{1 + \theta})$, which grows with the square root of the total variance. Lee's bound only binds at very long tenors or with very high volatilities. 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. **Solution.** ```python 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}") ``` ```output unweighted: maximum error 0.087, mean 0.020 weighted: maximum error 0.026, mean 0.010 ``` With weights, the maximum error against the truth falls from 0.087 to 0.026 points and the mean, from 0.020 to 0.010. The options with a narrow spread in volatility, near the money, carry more information about the true volatility, and the weighted fit respects them more. 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. **Solution.** Two errors. The first is interpolating at fixed strike: the 4,500 has $k = -0.134$ with the 16-day forward and $k = -0.136$ with the 44-day one, almost the same, but in standard deviations they are −2.9 and −1.8: in the short slice it is in the wing and in the long one, near the center. The second is interpolating volatilities instead of total variances: total variance is what adds up over time, and a day-weighted average of volatilities does not preserve it and can create calendar arbitrage. The correct approach is to compute the $k$ of the strike with the 30-day forward, convert it to $z = k/\sqrt{\theta_{30}}$ (with $\theta_{30}$ interpolated), evaluate each SVI at its own $k_i = z\sqrt{\theta_i}$, interpolate the total variance linearly in $T$ and divide by $T$; or use the fitted SSVI directly. ## Self-assessment 1. In the data, the skew per unit of ln(K/F) is much steeper at 16 days than at 653. Does it mean that the market fears a short-term fall more? - A. Yes: the short-term skew is always the most informative - B. No: it is an error in the synthetic data - C. Not necessarily: a 10% fall is many more standard deviations at 16 days than at two years; measured per standard deviation (k/√θ) the skew is similar at all tenors - D. Yes, because short puts are cheaper *Answer C: Not necessarily: a 10% fall is many more standard deviations at 16 days than at two years; measured per standard deviation (k/√θ) the skew is similar at all tenors.* The same ln(K/F) is a much larger distance, in standard deviations, the shorter the tenor. That is why skew is compared in k/√θ or in delta. In SSVI, the ATM slope of the volatility in k is approximately ρ·η/(2√T): it grows like 1/√T, and measured in standard deviations it depends only on the ATM volatility. 2. What condition guarantees that a slice of the surface has no butterfly arbitrage? - A. That the implied volatility is positive at all strikes - B. That ρ is between −1 and 1 - C. That the ATM volatility is lower than that of the wings - D. That Gatheral's function g(k), proportional to the risk-neutral density, is non-negative for all k *Answer D: That Gatheral's function g(k), proportional to the risk-neutral density, is non-negative for all k.* The second derivative of the call price with respect to the strike is the discounted risk-neutral density. Written in terms of the total variance w(k), its sign is that of g(k). A negative g implies a butterfly with a negative price: an arbitrage. 3. An SVI fitted only with the strikes within 6% of the forward reproduces those points well. What is the risk? - A. That the wings are left unconstrained: in the data, such a slice gives negative density, another a negative total variance at the minimum and others cross the next expiry - B. None: if it fits well where there is data, it fits well everywhere - C. That the fit is slower - D. That the forward changes *Answer A: That the wings are left unconstrained: in the data, such a slice gives negative density, another a negative total variance at the minimum and others cross the next expiry.* With few points near the money, the five parameters are not identified: many combinations give the same curve in the center and very different wings. In the 44-day slice, g(k) reaches −0.23; in the 562-day one, the minimum of w is negative. 4. To interpolate the 44-day slice between the 16- and 72-day ones, which method is least wrong in the data? - A. Linearly interpolating the implied volatility at fixed strike - B. Interpolating the total variance at fixed standard deviations (k/√θ), or using SSVI with interpolated θ - C. Interpolating the total variance at fixed ln(K/F) - D. They all give the same result *Answer B: Interpolating the total variance at fixed standard deviations (k/√θ), or using SSVI with interpolated θ.* In volatility at fixed k the error reaches 4.9 points and in total variance at fixed k, 1.8. At fixed k/√θ it falls to 0.07 and with SSVI to 0.12: the shape of the smile scales with the standard deviation, not with the strike. ## Key points - The index smile is a skew: out-of-the-money puts have more volatility than calls. In the 44-day slice of June 4, 2025, the 25-delta risk reversal is −9.0 points. - The ATM term structure slopes upward in calm (10.3% at one month and 15.5% at one year in June 2017) and inverts in stress (63.0% and 30.3% in March 2020). - SVI fits each slice with five parameters and, with all the strikes, recovers the true surface with less than 0.1 points of error; SSVI fits the whole surface with two (ρ and η) and recovers ρ = −0.662 and η = 0.953 (truth: −0.662 and 0.95). - No arbitrage: g(k) ≥ 0 in each slice and total variance increasing with tenor. A fit with few strikes violates them in the wings even if it fits the center well. - To interpolate between expiries, total variance at fixed standard deviations or SSVI: maximum error of 0.07 to 0.12 points, against 4.9 in volatility at fixed strike. ## Resources ### Documentation - [SciPy: scipy.optimize.least_squares](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html) · docs.scipy.org ### Books - [Gatheral, J. and Jacquier, A. (2014): Arbitrage-free SVI volatility surfaces. Quantitative Finance 14(1), 59-71](https://doi.org/10.1080/14697688.2013.819986) · doi.org - [Gatheral, J. and Jacquier, A.: Arbitrage-free SVI volatility surfaces, arXiv version (1204.0646)](https://arxiv.org/abs/1204.0646) · arxiv.org - [Lee, R. W. (2004): The Moment Formula for Implied Volatility at Extreme Strikes. Mathematical Finance 14(3), 469-480](https://doi.org/10.1111/j.0960-1627.2004.00200.x) · doi.org ### Data - [Level 3 data: weekly option surface on ORX (Parquet)](https://pythonparatrading.com/curso/datos/nivel-3/opciones3_superficie.parquet) · Download · opciones3_superficie.parquet - [Level 3 data: true SVI parameters of each slice (Parquet)](https://pythonparatrading.com/curso/datos/nivel-3/opciones3_verdad.parquet) · Download · opciones3_verdad.parquet Previous lesson: [3.4.4 · Implied volatility: Newton-Raphson, Brent, py_vollib and QuantLib](https://pythonparatrading.com/en/python-trading-course/level-3/implied-volatility-newton-brent-quantlib/) Next lesson: [3.4.6 · Realized vs implied volatility; volatility risk premium; VIX and VSTOXX](https://pythonparatrading.com/en/python-trading-course/level-3/realized-vs-implied-volatility-vix/)