Deep learning the SPY volatility surface

A volatility surface is a function of two things. How far an option's strike sits from the spot price, and how long until it expires. Feed those in, get an implied volatility back. The market redraws that surface for SPY every day the exchange is open, and it has a shape everyone in options knows by feel: a steep wall of fear on the downside puts, a gentle slope out to the long-dated calls.

The question here is whether a small neural network, handed nothing but raw quotes and no Black-Scholes structure, can learn the whole thing. And once one network can, which kind learns it best. The live dashboard at the foot of this page holds dozens of them, different depths, different activation functions, different feature sets, each trained on real chains and each watchable epoch by epoch as it converges. This page is the script underneath. The data, the network built by hand, the split that keeps the numbers honest, and the finding that made me rewrite my priors.

The walkthrough below trains a single network so the code stays readable. The full script that sweeps every architecture and activation, the version that feeds the dashboard, is on GitHub.

The data

The chains come from a public Kaggle set of SPY end-of-day options, spy-options-eod-volatility-surface-2010-2023, one parquet file per year. Fourteen years, close to six million quotes once cleaned, with implied vol and Greeks already computed. Each raw row is wide, carrying the call and the put side of one strike together under bracketed names like [C_IV] and [P_IV], so the first job is to melt it into one tidy row per quote and drop the junk.

Three filters do the cleaning. Keep the out-of-the-money side at each strike, puts below spot and calls above, because that is where implied vol reads cleanest, away from the intrinsic value of deep in-the-money contracts. Throw out anything with a zero bid, a spread wider than half its mid price, or an implied vol outside a sane band. And keep maturities between one week and one year, where the quotes are dense and liquid.

import numpy as np, pandas as pd, math
from pathlib import Path

DTE_MIN, DTE_MAX = 7, 365        # one week to one year, the liquid part
IV_MIN,  IV_MAX  = 0.03, 2.00    # kills negatives and 40.0 blow-ups
MAX_REL_SPREAD   = 0.50          # drop quotes wider than 50% of mid

def load_year(path):
    df = pd.read_parquet(path)
    df.columns = [c.strip().strip("[]") for c in df.columns]     # [C_IV] -> C_IV
    spot   = pd.to_numeric(df["UNDERLYING_LAST"], errors="coerce")
    strike = pd.to_numeric(df["STRIKE"],          errors="coerce")

    def side(prefix, otm, kind):
        bid = pd.to_numeric(df[f"{prefix}_BID"], errors="coerce")
        ask = pd.to_numeric(df[f"{prefix}_ASK"], errors="coerce")
        iv  = pd.to_numeric(df[f"{prefix}_IV"],  errors="coerce")
        mid = (bid + ask) / 2
        out = pd.DataFrame({"dte": pd.to_numeric(df["DTE"], errors="coerce"),
                            "spot": spot, "strike": strike, "type": kind, "iv": iv})
        keep = (otm & bid.gt(0) & mid.gt(0) & iv.between(IV_MIN, IV_MAX)
                & ((ask - bid) / mid).le(MAX_REL_SPREAD)
                & out["dte"].between(DTE_MIN, DTE_MAX))
        return out.loc[keep]

    puts  = side("P", strike < spot, "put")     # OTM puts sit below spot
    calls = side("C", strike > spot, "call")    # OTM calls sit above spot
    q = pd.concat([puts, calls], ignore_index=True)
    q["tau"]   = q["dte"] / 365.0
    q["log_m"] = np.log(q["strike"] / q["spot"])         # log-moneyness
    return q

DATA = Path("data/historical")
quotes = pd.concat([load_year(p) for p in sorted(DATA.glob("spy_eod_*.parquet"))],
                   ignore_index=True)
print(f"{len(quotes):,} clean OTM quotes, 2010 to 2023")

That keeps a little over a quarter of the raw rows, nearly six million clean quotes spanning the fourteen years. The surface is massively oversampled, so later we work with a subsample rather than the full pool.

The features

The network sees geometry, not prices. From log-moneyness x = ln(K / S) and time to expiry τ it builds nine features. Five polynomial terms give it the basic bowl. Two smile terms, x / √τ and 1 / √τ, are the standardized-moneyness coordinate that every parametric vol model leans on. And two wing terms, x² / τ and x / τ, blow up at short maturities, exactly where the smile steepens into a sharp corner.

def make_features(df):
    x, tau = df["log_m"].to_numpy(), df["tau"].to_numpy()
    rt, it = np.sqrt(np.maximum(tau, 1e-6)), 1.0 / np.maximum(tau, 1e-6)
    return np.stack([
        x, tau, x**2, tau**2, x*tau,     # polynomial base
        x / rt, 1.0 / rt,                # smile: standardized moneyness, 1/sqrt(T)
        x**2 * it, x * it,               # wing: 1/T blow-up for the short-dated corner
    ], axis=1)

Why geometry, and not the state of the market. Because moneyness and maturity mean the same thing in every regime. A ten percent out-of-the-money put is a ten percent out-of-the-money put whether it is 2011 or 2022. That stationarity is what lets a model trained on old data say something about a year it has never seen. It matters later, when the split gets strict.

The network, from scratch

No PyTorch. A plain numpy class holds the weights, runs a forward pass, and updates with Adam, first and second moment estimates of the gradient, bias-corrected for the warmup. The hidden layers pass through an activation function. The output layer stays linear, because implied vol is a real number, not a probability. The one line that is the whole experiment is the choice of that activation: ReLU, tanh, or softplus.

class MLP:
    """Feed-forward net with Adam, pure numpy. activation = relu | tanh | softplus."""
    def __init__(self, sizes, seed=42, activation="tanh"):
        rng = np.random.default_rng(seed)
        self.activation, self.W, self.b, self.t = activation, [], [], 0
        scale = 2.0 if activation in ("relu", "softplus") else 1.0   # He vs Xavier init
        for i in range(len(sizes) - 1):
            self.W.append(rng.standard_normal((sizes[i], sizes[i+1])) * math.sqrt(scale / sizes[i]))
            self.b.append(np.zeros((1, sizes[i+1])))
        self.mW = [np.zeros_like(w) for w in self.W]; self.vW = [np.zeros_like(w) for w in self.W]
        self.mb = [np.zeros_like(b) for b in self.b]; self.vb = [np.zeros_like(b) for b in self.b]

    def _act(self, z):
        if self.activation == "relu":     return np.maximum(0, z)
        if self.activation == "softplus": return np.logaddexp(0.0, z)      # smooth log(1 + e^z)
        return np.tanh(z)

    def _dact(self, z):
        if self.activation == "relu":     return (z > 0).astype(z.dtype)
        if self.activation == "softplus": return 1 / (1 + np.exp(-np.clip(z, -30, 30)))
        return 1 - np.tanh(z)**2

    def forward(self, X):
        a, acts, zs = X, [X], []
        for i, (W, b) in enumerate(zip(self.W, self.b)):
            z = a @ W + b; zs.append(z)
            a = self._act(z) if i < len(self.W) - 1 else z    # hidden activated, output linear
            acts.append(a)
        return a, acts, zs

    def predict(self, X):
        return self.forward(X)[0]

    def backward(self, y, acts, zs, lr, b1=0.9, b2=0.999, eps=1e-8):
        self.t += 1
        delta = 2 * (acts[-1] - y) / y.shape[0]               # MSE gradient
        for i in reversed(range(len(self.W))):
            dW = acts[i].T @ delta
            db = delta.sum(0, keepdims=True)
            if i > 0:
                delta = (delta @ self.W[i].T) * self._dact(zs[i-1])
            self.mW[i] = b1*self.mW[i] + (1-b1)*dW;  self.vW[i] = b2*self.vW[i] + (1-b2)*dW*dW
            self.mb[i] = b1*self.mb[i] + (1-b1)*db;  self.vb[i] = b2*self.vb[i] + (1-b2)*db*db
            mW = self.mW[i]/(1-b1**self.t); vW = self.vW[i]/(1-b2**self.t)
            mb = self.mb[i]/(1-b1**self.t); vb = self.vb[i]/(1-b2**self.t)
            self.W[i] -= lr * mW / (np.sqrt(vW) + eps)
            self.b[i] -= lr * mb / (np.sqrt(vb) + eps)

ReLU is the default in nearly every tutorial, a hard hinge, the larger of zero and the input. tanh is a smooth S-curve that flattens toward plus and minus one. softplus is a rounded ReLU. Swap between them and nothing else about the training changes: same data, same optimiser, same learning rate. That is what makes the comparison fair.

Train, validate, test

This is where a lot of surface fitting quietly cheats. Fit and score on the same quotes and you are measuring memory, not learning. So the quotes get split three ways. The network trains on one slice, gets scored on a held-out validation slice it never fits, and a test slice is touched once at the very end.

There are two honest ways to draw that split, and they answer different questions. A random split, pooling all fourteen years and dealing the quotes out at random, measures how well the model interpolates a smooth surface. That is the setting the live dashboard uses, eighty thousand quotes to train, twenty thousand held out. A chronological split, training on the early years and testing on a later one the model has never seen, measures something harder: whether the fit survives a new regime. I ran both, and the gap between them is the real lesson of the project.

One rule holds either way. Standardize the features with the training set's own mean and standard deviation, then apply those same numbers to validation and test. Fit the scaler on everything and you have leaked the future into the past.

# pool every clean quote, then a random held-out split (the dashboard setting)
idx = np.random.default_rng(0).permutation(len(quotes))
train = quotes.iloc[idx[:80_000]]
val   = quotes.iloc[idx[80_000:100_000]]     # held out, never trained on

Xtr, Xva = make_features(train), make_features(val)
mean, std = Xtr.mean(0), Xtr.std(0) + 1e-9   # standardize on TRAIN only, no leakage
Xtr, Xva = (Xtr - mean) / std, (Xva - mean) / std

iv_mean = train["iv"].mean()                 # center the target around its mean
ytr = (train["iv"].to_numpy() - iv_mean).reshape(-1, 1)
yva = (val["iv"].to_numpy()   - iv_mean).reshape(-1, 1)

# the champion architecture: three hidden layers of 16 units, tanh
net = MLP([Xtr.shape[1], 16, 16, 16, 1], activation="tanh")
for epoch in range(100):
    perm = np.random.default_rng(epoch).permutation(len(Xtr))
    for s in range(0, len(Xtr), 8192):
        b = perm[s:s + 8192]
        _, acts, zs = net.forward(Xtr[b])
        net.backward(ytr[b], acts, zs, lr=0.01)

rmse = np.sqrt(np.mean((net.predict(Xva) - yva)**2)) * 100   # vol points, held out
print(f"held-out RMSE: {rmse:.2f} vol points")

Run exactly that, geometry alone, and the network lands around five to six vol points of held-out error. Good enough to trace the shape of the surface, not good enough to price it. Getting from there to the numbers you see in the dashboard takes one more idea.

The activation function decides

Geometry fixes the shape of the surface. What it cannot fix is the level, which drifts every day with the mood of the market. So the full script feeds the network a second block of features, the state of the market on the quote date: realized vol, the VIX, the rate curve, credit spreads, thirteen series in all. Under the random split those features carry the level directly, and the held-out error collapses from around six vol points to just over one. The dashboard below is trained this way.

Now hold the features fixed and sweep architecture against activation. tanh wins every single cell.

Read down the columns and two things stand out. tanh beats ReLU at every size, and depth beats width, the 929-parameter deep net edges out the wide one with twenty times as many weights. Capacity, the knob everyone reaches for first, barely moves the number once the activation is right.

The reason is not mysterious. A vol surface is smooth, and tanh is smooth, so it fits the gentle curvature with a handful of units. ReLU tiles the same surface out of little flat planes, so it needs far more of them to round off one corner, and it still leaves faint creases along the wings.

Smoothness pays off most exactly where it is hardest to see, at the edge of the data. This is where the chronological split earns its keep. Train on the near-zero rates of the 2010s and test on the five percent rates of 2023, and those market-state features drift into territory the model never trained on. ReLU keeps extrapolating in straight lines, off toward infinity, and its error explodes past ten vol points. tanh saturates, flattens out, and degrades gently. That 1.31 percent is interpolation skill on a random split. The activation function is what decides whether the model merely interpolates or actually holds up when the world moves.

Where it still misses

No model nails the whole surface. Break the error down by region and it looks like a bathtub, close to one vol point across the body, climbing steeply into the short-dated deep out-of-the-money wings. Seven-day puts far below spot, the crash skew, are the hardest part, and the short-dated far out-of-the-money call tip is the single worst cell for every one of the networks.

That corner is a feature problem, not a size problem. A polynomial cannot express the way the wing steepens like one over the square root of maturity, which is why the wing terms exist, and why the full script also leans the training weight toward the tail. Together they lift the surface into both wings and roughly halve the worst gap. The very tip still resists. It is genuinely thin data, a few hundred quotes on a good day, priced by the market with real fear that a smooth surface wants to average away.

The other two tabs

The dashboard carries two more views built on the same idea, a model opinion of fair implied vol.

Correlation map

Eight daily features drawn from the surface and the tape: ATM vol, the put minus call skew, the term slope, the curvature, the overall vol level, realized vol, the SPY return, and quote breadth. The map animates their ninety-day rolling correlations, stepping forward one month at a time from 2010, so you watch the structure tighten in every stress and loosen in the calm. Next to it sits a plain vol-timing idea, go long vol when the market prices it below the model and short when it prices it above, shown as an illustration rather than a track record.

Option screen

A fair-value network trained only on quotes up to 2016, then pointed at every month-end from 2017 to 2023. For each day it ranks the chain by the gap between market implied vol and model implied vol. The cheap ones, where the market sits below the model, flag as buys. The rich ones flag as sells. It turns the surface into something concrete, a list of specific options the model thinks are mispriced.

What holds up

Three things survive the whole exercise. A small numpy network with no pricing theory built in learns the SPY surface to about one vol point, so the structure really is sitting in the data. The activation function, not the parameter count, decides how well, because a smooth surface wants a smooth function and ReLU fights that with flat planes. And the honest test is the one that hurts: a model that interpolates a random split at one percent can still fall apart out of time, which is the whole reason to draw the split with care and report both numbers.

Drive it below. Switch the activation to ReLU and watch the wings go blocky, add depth and watch the corner fill in, then move to the screen and see the same model hand you a list of cheap and rich options. If you want the parametric cousin of this, the one that says the smile must come from five stochastic-vol parameters rather than a stack of weights, see the Heston surface notebook. The one-day, single-snapshot version of this idea is here.