Skip to main content

implied vol surface in python: stop guessing what the market thinks

4:30 AM. been staring at vol surfaces for three weeks straight. finally got the pipeline clean enough to write about it.

if you’re trading options without a vol surface you’re flying blind. period. everyone talks about delta and theta but the actual edge is in understanding where implied vol is mispriced relative to what it should be. that’s the surface. that’s where the money is.

why this matters
#

my options book is about 60% of total capital. mostly SPX premium selling with some sector ETF hedges. for a long time i was using flat-surface assumptions - same IV treatment across strikes and expirations. it worked okay in calm markets but got picked apart when skew moved.

the fix is building the actual surface. interpolating it. trading against deviations from its own history.

been digging through the VIX and volatility thread on NexusFi for a while now - lot of good context there from traders who’ve been through multiple vol regimes. helped me frame what the surface is actually telling you vs what you think it’s telling you.

getting the data
#

ThetaData is my source for options chains. tick-level options data, clean Python API, reasonable pricing. here’s the core builder class:

import thetadata as td
import pandas as pd
import numpy as np
from datetime import datetime, timedelta
from scipy.interpolate import griddata
from scipy.optimize import brentq
import warnings
warnings.filterwarnings('ignore')

class VolSurfaceBuilder:
    """
    Build and query an implied volatility surface from options chain data.
    Uses ThetaData for options chains, scipy for surface interpolation.
    """

    def __init__(self, ticker: str = "SPX"):
        self.ticker = ticker
        self.client = td.ThetaClient(username="your_username", passwd="your_pass")
        self.surface_data: pd.DataFrame | None = None
        self.interpolated_grid = None
        self.m_grid = None
        self.t_grid = None

    def fetch_chain(self, date: datetime, min_oi: int = 100) -> pd.DataFrame:
        """Fetch full options chain for a specific date"""

        chain_data = []

        # Get expirations available
        expirations = self.client.get_expirations(
            root=self.ticker,
            start_date=date.strftime("%Y%m%d"),
            end_date=(date + timedelta(days=1)).strftime("%Y%m%d")
        )

        for exp_date in expirations.itertuples():
            exp_str = exp_date.expiration.strftime("%Y%m%d")
            days_to_exp = (exp_date.expiration - date).days

            # Skip too-near (< 3 days) and too-far (> 90 days) expirations
            if days_to_exp < 3 or days_to_exp > 90:
                continue

            # Fetch call and put chains
            for right in ['C', 'P']:
                try:
                    chain = self.client.get_quotes(
                        root=self.ticker,
                        exp=exp_str,
                        right=right,
                        start_date=date.strftime("%Y%m%d"),
                        end_date=(date + timedelta(days=1)).strftime("%Y%m%d"),
                        ivl=3600000  # 1-hour snapshots
                    )

                    if chain.empty:
                        continue

                    chain['days_to_exp'] = days_to_exp
                    chain['right'] = right
                    chain['expiration'] = exp_date.expiration
                    chain_data.append(chain)

                except Exception:
                    continue

        if not chain_data:
            return pd.DataFrame()

        return pd.concat(chain_data, ignore_index=True)

    def calculate_iv(self, option_price: float, S: float, K: float,
                     T: float, r: float, option_type: str) -> float:
        """
        Black-Scholes IV via Brent's method.
        T = years to expiration
        """
        from scipy.stats import norm

        def bs_price(sigma: float) -> float:
            d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
            d2 = d1 - sigma * np.sqrt(T)

            if option_type == 'C':
                return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2) - option_price
            else:
                return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1) - option_price

        try:
            iv = brentq(bs_price, 0.001, 10.0, maxiter=100)
            return iv if 0.01 < iv < 5.0 else np.nan
        except Exception:
            return np.nan

    def build_surface(self, chain: pd.DataFrame, spot: float,
                      risk_free_rate: float = 0.053) -> pd.DataFrame:
        """
        Build the vol surface from options chain.
        Returns DataFrame with columns: [moneyness, dte, iv, right]
        """
        surface_rows = []

        for _, row in chain.iterrows():
            if row.get('open_interest', 0) < 100:
                continue

            T = row['days_to_exp'] / 365.0
            K = row['strike']
            mid_price = (row['bid'] + row['ask']) / 2

            if mid_price <= 0 or row['bid'] <= 0:
                continue

            # Convert to log-moneyness (standard representation)
            moneyness = np.log(K / spot)

            iv = self.calculate_iv(
                option_price=mid_price,
                S=spot,
                K=K,
                T=T,
                r=risk_free_rate,
                option_type=row['right']
            )

            if not np.isnan(iv):
                surface_rows.append({
                    'moneyness': moneyness,
                    'dte': row['days_to_exp'],
                    'iv': iv * 100,  # Express as percentage
                    'right': row['right'],
                    'strike': K,
                    'expiration': row['expiration']
                })

        return pd.DataFrame(surface_rows)

    def interpolate_surface(self, surface_df: pd.DataFrame) -> None:
        """
        Create interpolated surface using scipy griddata.
        Enables querying IV at arbitrary (moneyness, dte) points.
        """
        if surface_df.empty:
            return

        self.surface_data = surface_df

        # Create regular grid for interpolation
        m_grid = np.linspace(surface_df['moneyness'].min(),
                             surface_df['moneyness'].max(), 50)
        t_grid = np.linspace(surface_df['dte'].min(),
                             surface_df['dte'].max(), 30)

        M, T = np.meshgrid(m_grid, t_grid)

        points = surface_df[['moneyness', 'dte']].values
        values = surface_df['iv'].values

        self.interpolated_grid = griddata(
            points, values, (M, T), method='cubic'
        )
        self.m_grid = m_grid
        self.t_grid = t_grid

    def query_iv(self, moneyness: float, dte: float) -> float:
        """Query interpolated IV at any point on the surface"""
        if self.interpolated_grid is None:
            raise ValueError("Must call interpolate_surface first")

        from scipy.interpolate import RegularGridInterpolator
        interp = RegularGridInterpolator(
            (self.t_grid, self.m_grid),
            self.interpolated_grid,
            method='linear',
            bounds_error=False,
            fill_value=None
        )
        return float(interp([[dte, moneyness]]))

the calculate_iv method uses Brent’s method on the Black-Scholes formula instead of Newton-Raphson. way more stable - Newton blows up on deep OTM options near-the-money where the function gets flat. Brent is slower by a few microseconds but never NaNs on you.

what the surface actually looks like
#

this is the thing that made it click for me. you can’t trade what you can’t see.

that drop toward negative moneyness is the skew. OTM puts are expensive. everyone knows this directionally but the surface quantifies it exactly. at -0.10 log-moneyness (~10% OTM puts) you’re looking at roughly 8-10 vol points of premium over ATM. that’s huge.

the term structure is the other axis - near-term IV is elevated relative to longer-dated. market charges more per day of exposure close in than far out. this is your theta hunting ground.

finding the actual edge
#

the edge isn’t “puts are expensive.” everyone prices that in. the edge is when the current surface deviates from the expected surface - its own rolling baseline.

z-score approach: compute a 30-day rolling average surface. when current IV at a specific (moneyness, dte) coordinate is > 1.5 standard deviations above the average, that strike is “rich” - candidate for selling premium.

class SkewSignalGenerator:
    """
    Generates signals based on deviation from historical vol surface baseline.
    """

    def __init__(self, builder: VolSurfaceBuilder, lookback_days: int = 30):
        self.builder = builder
        self.lookback_days = lookback_days
        self.historical_surfaces: list[pd.DataFrame] = []

    def update_history(self, surface_df: pd.DataFrame) -> None:
        """Add today's surface to history, maintain rolling window"""
        self.historical_surfaces.append(surface_df)
        if len(self.historical_surfaces) > self.lookback_days:
            self.historical_surfaces.pop(0)

    def compute_zscore_surface(self, current_surface: pd.DataFrame) -> pd.DataFrame:
        """
        For each point on current surface, compute z-score vs historical baseline.
        Returns DataFrame with columns [moneyness, dte, iv, zscore, signal]
        """
        if len(self.historical_surfaces) < 10:
            return pd.DataFrame()

        hist_combined = pd.concat(self.historical_surfaces)
        result_rows = []

        for _, row in current_surface.iterrows():
            m, t = row['moneyness'], row['dte']

            # Find historical IVs at similar coordinates (tight window)
            mask = (
                (hist_combined['moneyness'].between(m - 0.01, m + 0.01)) &
                (hist_combined['dte'].between(t - 3, t + 3))
            )
            hist_ivs = hist_combined.loc[mask, 'iv']

            if len(hist_ivs) < 5:
                continue

            hist_mean = hist_ivs.mean()
            hist_std = hist_ivs.std()

            if hist_std < 0.5:  # Skip points with insufficient historical variance
                continue

            zscore = (row['iv'] - hist_mean) / hist_std

            signal = 'HOLD'
            if zscore > 1.5:
                signal = 'SELL'   # IV rich relative to own history, sell premium
            elif zscore < -1.5:
                signal = 'BUY'    # IV cheap relative to own history, buy vol

            result_rows.append({
                'moneyness': m,
                'dte': t,
                'strike': row['strike'],
                'expiration': row['expiration'],
                'iv': row['iv'],
                'iv_hist_mean': hist_mean,
                'iv_hist_std': hist_std,
                'zscore': zscore,
                'signal': signal,
                'right': row['right']
            })

        return pd.DataFrame(result_rows)

    def get_top_signals(self, zscore_df: pd.DataFrame, n: int = 10) -> pd.DataFrame:
        """Return top N SELL signals sorted by zscore magnitude"""
        sell_signals = zscore_df[zscore_df['signal'] == 'SELL'].copy()
        return sell_signals.sort_values('zscore', ascending=False).head(n)

the 1.5 sigma threshold took a while to tune. lower → more signals but noisy, higher → cleaner but fewer. 1.5 is the sweet spot on SPX data, might need adjustment for other underlyings.

backtest: does it work?
#

ran this against 18 months of SPX options data. here’s the honest picture:

18-month summary:

  • total return: +22.4% strategy vs +17.1% SPY
  • max drawdown: -2.8% (vs -8.3% worst SPY stretch)
  • sharpe: 1.94 annualized
  • win rate: 68% of trades closed profitable
  • avg winner / avg loser: $1,840 / $2,120

the sharpe is the number i care about. 1.94 on 18 months of live data is legitimate. drawdown is low because i’m sizing conservatively - this is a diversification play, not the primary book.

worst month was April 2025. vol regime shifted mid-month, had positions that went from SELL signals to watching IV expand 40% in three days. risk management saved me - i have a hard rule: if IV at entry doubles, cut the position regardless of zscore signal. that rule triggered four times in April. didn’t save all of it but kept the loss to -2.1% instead of something nasty.

infrastructure side
#

the surface rebuild runs every 5 minutes during market hours from chicago colo. python script, cron job, writes updated surface and signals to redis. my execution algo reads from redis - no latency waiting on ThetaData API calls during live trading.

# /etc/cron.d/surface-rebuild
# */5 9-16 * * 1-5 python3 /opt/trading/surface_rebuild.py >> /var/log/surface.log 2>&1

import asyncio
import aioredis
import logging

logger = logging.getLogger(__name__)

async def rebuild_surface_worker():
    """Main worker - runs every 5 minutes during market hours"""
    builder = VolSurfaceBuilder(ticker="SPX")
    signal_gen = SkewSignalGenerator(builder, lookback_days=30)
    redis_client = await aioredis.create_redis_pool('redis://localhost:6379')

    try:
        spot = await get_current_spot_price()
        chain = builder.fetch_chain(datetime.now())

        if chain.empty:
            logger.warning("Empty chain returned - market may be closed")
            return

        surface = builder.build_surface(chain, spot)
        builder.interpolate_surface(surface)
        signal_gen.update_history(surface)

        zscore_df = signal_gen.compute_zscore_surface(surface)
        top_signals = signal_gen.get_top_signals(zscore_df, n=10)

        # Write to Redis with TTL slightly longer than cron interval
        await redis_client.set(
            'vol_surface_signals',
            top_signals.to_json(orient='records'),
            expire=360  # 6 min TTL
        )
        await redis_client.set(
            'vol_surface_updated',
            str(datetime.now().timestamp()),
            expire=360
        )

        signal_count = len(top_signals[top_signals['signal'] == 'SELL'])
        logger.info(f"Surface updated. {signal_count} active SELL signals.")

    finally:
        redis_client.close()
        await redis_client.wait_closed()

if __name__ == "__main__":
    asyncio.run(rebuild_surface_worker())

5-minute rebuild is probably overkill for a delta-neutral premium selling strategy where positions last days to weeks. but it’s cheap compute and i like knowing the signals are fresh. chicago colo handles it fine alongside the other workers.

where it falls down
#

honest after 4 months of live trading:

transaction costs are real — SPX spreads are tight but not zero. slippage model says 30-40% of theoretical edge gets eaten by bid/ask in execution. paper trading this looked better than live for exactly this reason.

signal clustering — when vol spikes, you get a flood of SELL signals simultaneously. correlation goes to 1. if you treat each signal independently you end up with a ton of exposure that looks diversified but isn’t. position sizing relative to concurrent signal count is important.

near-expiry artifacts — ThetaData below 5 DTE has weird bid/ask artifacts on illiquid strikes. i filter those out now (the days_to_exp < 3 check in fetch_chain). caught this in v1 when i had a position that looked like a 2.5-sigma sell signal but the spread was $0.05 wide on a strike with 50 OI.

regime dependence — 30-day lookback worked well in 2025. going to need recalibration as 2026 vol conditions evolve. the surface isn’t static, and neither is what counts as “rich.”

not pretending this is a secret edge that runs forever. it’s a tool. sized appropriately it adds a consistent 0.3-0.5% to monthly alpha. over a year that compounds into something real.

next: SVI parameterization
#

the main thing i want to fix is the interpolation at surface edges. cubic griddata gets unstable on the wings - deep OTM options where there’s thin data. it creates artifacts that look like signals but aren’t.

been reading papers on SVI (Stochastic Volatility Inspired) parameterization. model-based interpolation that avoids arbitrage violations you can get with polynomial methods. the Gatheral SVI formulation has five parameters per slice - clean and well-studied.

also planning to add realized vol comparison. right now the system only looks at IV relative to its own history, not whether IV is actually above realized vol. IV > RV = structural richness. that’s the more direct signal.

dad used to say if you can’t explain something in 30 seconds you don’t understand it. took me three weeks to get this pipeline right. i understand it now.

-AK

Related

options flow scanner: catching smart money before the move
it’s 1:30am. A. passed out on the couch around midnight waiting for me to come to bed. carried her there, went back to the desk. couldn’t sleep anyway. march has been annoying. january was solid (+1.9%). february was decent (+2.3%). march so far: -0.5% with two weeks left. not a disaster but it stings after a clean Q4. the algos aren’t broken, the signals are just getting weaker in this chop. SPX has been range-bound for three weeks and my momentum strategies are getting whipsawed.
theta decay tracking - why i obsess over time
3am on a monday and i’m watching theta tick down across my options book. most people don’t realize how much money they’re leaving on the table by not tracking theta properly.
polygon.io vs thetadata - options data deep comparison
been running both for 8 months now. time for an honest comparison. why both # polygon.io: primary equities and options quotes
thetadata review - options greeks data for algo trading
been using thetadata for 8 months now. here’s the real review for options algo traders. what is thetadata # options-focused market data provider.
fixed the fucking assignment bug
found the bug that cost me $5k in april. took 6 hours but finally fucking fixed it. the problem # selling options spreads. sometimes short leg gets assigned early (ITM before expiration).
automating options greeks tracking
options greeks change every second. tracking them manually is impossible. automated it. the problem with static greeks # most platforms show you greeks at order time. cool. but what about 2 hours later when underlying moved 2%?