Implied Volatility Surface in C++

Setup

The Problem

A vanilla option's Black-Scholes price is a strictly increasing function of volatility. Given a market price, there exists a unique $\sigma_{\text{impl}}$ such that $\text{BS}(S, K, T, r, \sigma_{\text{impl}}) = \text{market price}$. This is the implied volatility.

If the Black-Scholes model were correct, $\sigma_{\text{impl}}$ would be the same for every strike $K$ and maturity $T$. In practice, it is not. The surface $\sigma^*(T, K)$ of market-implied volatilities exhibits:

• A smile or skew along the strike axis (more pronounced for equity index options)
• A term structure along the maturity axis (implied vol changes with expiry)

For a Monte Carlo pricer to match all market prices simultaneously, it needs to read from this surface at arbitrary $(t, S_t)$ during path simulation. This is the Dupire local volatility model, and the ImpliedVolatilitySurface class is its data structure.

What this module builds:

1. A ThomasSolver for the tridiagonal system arising from natural cubic spline interpolation.
2. An ImpliedVolatilitySurface that interpolates across strikes (cubic spline, per maturity) and across maturities (linear in total variance).
3. The numerical differentiation formulae that feed the Dupire formula.

Conventions:

• Volatilities are annualised, in decimal (0.20 = 20%).
• Maturities are in years from today.
• The risk-free rate $r$ is continuously compounded and used for forward-moneyness adjustments.
• No dividends.

Theory

1. Cubic Spline Interpolation

Given $N$ strike nodes $K_1 < K_2 < \cdots < K_N$ with corresponding implied vols $\sigma_1^*, \ldots, \sigma_N^*$ (at a fixed maturity $T_i$), we want a function $s(K)$ that:

• Passes exactly through each node: $s(K_j) = \sigma_j^*$
• Is a cubic polynomial on each interval $[K_j, K_{j+1}]$
• Has continuous first and second derivatives everywhere

Let $\Delta_j = K_{j+1} - K_j$. On interval $[K_j, K_{j+1}]$, write:

$s(K) = \alpha_j (K - K_j)^3 + \beta_j (K - K_j)^2 + \gamma_j (K - K_j) + \sigma_j^*$

The continuity of the second derivative at each interior knot $K_j$ yields the system of equations:

$\Delta_{j-1} \beta_{j-1} + 2(\Delta_{j-1} + \Delta_j) \beta_j + \Delta_j \beta_{j+1} = 3\!\left(\frac{\sigma_{j+2}^* - \sigma_{j+1}^*}{\Delta_{j+1}} - \frac{\sigma_{j+1}^* - \sigma_j^*}{\Delta_j}\right), \quad j = 1, \ldots, N-2$

with natural spline boundary conditions $\beta_1 = \beta_N = 0$ (zero second derivative at the endpoints).

This is a tridiagonal linear system for $\beta_2, \ldots, \beta_{N-1}$. Once the $\beta_j$ are known:

$\alpha_j = \frac{\beta_{j+1} - \beta_j}{3\Delta_j}, \qquad \gamma_j = \frac{\sigma_{j+1}^* - \sigma_j^*}{\Delta_j} - \alpha_j \Delta_j^2 - \beta_j \Delta_j$

2. The Thomas Algorithm

The tridiagonal system $A\mathbf{x} = \mathbf{b}$ has the form:

$\begin{pmatrix} b_1 & c_1 & & \\ a_2 & b_2 & c_2 & \\ & \ddots & \ddots & \ddots \\ & & a_{N-1} & b_{N-1} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_{N-1} \end{pmatrix} = \begin{pmatrix} R_1 \\ R_2 \\ \vdots \\ R_{N-1} \end{pmatrix}$

The Thomas algorithm solves this in $O(N)$ operations (versus $O(N^3)$ for Gaussian elimination on a general matrix):

Forward sweep — eliminate the lower diagonal: $c_j' = \frac{c_j}{b_j - a_j c_{j-1}'}, \qquad R_j' = \frac{R_j - a_j R_{j-1}'}{b_j - a_j c_{j-1}'}$

Back substitution: $x_{N-1} = R_{N-1}', \qquad x_j = R_j' - c_j' x_{j+1}$

The algorithm is numerically stable for diagonally dominant matrices (which the cubic spline system is, since the off-diagonal elements equal $\Delta_j$ and the diagonal equals $2(\Delta_{j-1}+\Delta_j) > \Delta_{j-1} + \Delta_j$).

3. Interpolation Along the Maturity Axis

For a query $(T, K)$ with $T_i \leq T < T_{i+1}$, we cannot simply interpolate implied vols — linear interpolation in vol introduces calendar spread arbitrage. The no-arbitrage condition requires total variance $w(T,K) = \sigma^*(T,K)^2 \cdot T$ to be non-decreasing in $T$.

Forward-moneyness adjustment: At maturity $T_i$, the strike equivalent to $K$ at time $T$ is: $K^{(i)} = K \cdot e^{r(T_i - T)}$ This aligns the forward moneyness $K/F(T)$ across maturities, where $F(T) = S_0 e^{rT}$ is the forward price.

Linear interpolation in total variance: $w(T, K) = w(T_i, K^{(i)}) + \frac{w(T_{i+1}, K^{(i+1)}) - w(T_i, K^{(i)})}{T_{i+1} - T_i} \cdot (T - T_i)$

$\sigma^*(T, K) = \sqrt{w(T, K) / T}$

4. Dupire Local Volatility

Given the implied volatility surface $\sigma^*(T,K)$, the Dupire (1994) formula extracts the local volatility $\sigma_D(t,S)$ such that the model: $dS_t = r S_t\,dt + \sigma_D(t, S_t) S_t\,dW_t$ prices all vanilla options consistently:

$\sigma_D^2(T, K) = \frac{\frac{\partial w}{\partial T}}{ 1 - \frac{y}{w}\frac{\partial w}{\partial y} + \frac{1}{4}\left(-\frac{1}{4} - \frac{1}{w} + \frac{y^2}{w^2}\right)\!\left(\frac{\partial w}{\partial y}\right)^{\!2} + \frac{1}{2}\frac{\partial^2 w}{\partial y^2} }$

where $y = \ln(K/F(T))$ is log-moneyness and $w(T,K) = \sigma^*(T,K)^2 T$.

In the implementation, the partial derivatives are computed by finite differences on the surface using central stencils:

$\frac{\partial \sigma^*}{\partial T}(T,K) \approx \frac{\sigma^*(T+\varepsilon_T, K) - \sigma^*(T-\varepsilon_T, K)}{2\varepsilon_T}$

$\frac{\partial \sigma^*}{\partial K}(T,K) \approx \frac{\sigma^*(T, K+\varepsilon_K) - \sigma^*(T, K-\varepsilon_K)}{2\varepsilon_K}$

$\frac{\partial^2 \sigma^*}{\partial K^2}(T,K) \approx \frac{\sigma^*(T, K+\varepsilon_K) - 2\sigma^*(T,K) + \sigma^*(T, K-\varepsilon_K)}{\varepsilon_K^2}$

The DupireLocalVolatilityModel wraps an ImpliedVolatilitySurface and computes these derivatives on the fly during path simulation.

Implementation

ThomasSolver.h

#pragma once
#include <vector>

using Vector = std::vector<double>;

class ThomasSolver {
public:
    // Solves the tridiagonal system A x = rhs where:
    //   lower_diag = [a_2, ..., a_N]   (size N-1)
    //   central_diag = [b_1, ..., b_N] (size N)
    //   upper_diag = [c_1, ..., c_{N-1}] (size N-1)
    //   rhs = [R_1, ..., R_N]           (size N)
    ThomasSolver(const Vector& lower_diag,
                 const Vector& central_diag,
                 const Vector& upper_diag,
                 const Vector& rhs);

    // Returns solution vector [x_1, ..., x_N]
    Vector solve() const;

private:
    Vector _lower_diagonal;
    Vector _central_diagonal;
    Vector _upper_diagonal;
    Vector _right_hand_side;
};

ThomasSolver.cpp

#include "ThomasSolver.h"
#include <stdexcept>

ThomasSolver::ThomasSolver(const Vector& lower,
                           const Vector& central,
                           const Vector& upper,
                           const Vector& rhs)
    : _lower_diagonal(lower),
      _central_diagonal(central),
      _upper_diagonal(upper),
      _right_hand_side(rhs)
{
    const size_t N = central.size();
    if (lower.size() != N - 1 || upper.size() != N - 1 || rhs.size() != N)
        throw std::invalid_argument("ThomasSolver: inconsistent vector sizes");
}

Vector ThomasSolver::solve() const {
    const size_t N = _central_diagonal.size();

    // Working copies — forward sweep modifies them
    Vector c_prime(N - 1);
    Vector R_prime(N);

    // Forward sweep: eliminate lower diagonal
    c_prime[0] = _upper_diagonal[0] / _central_diagonal[0];
    R_prime[0] = _right_hand_side[0] / _central_diagonal[0];

    for (size_t j = 1; j < N; ++j) {
        double denom = _central_diagonal[j]
                       - (j > 0 ? _lower_diagonal[j - 1] * c_prime[j - 1] : 0.0);
        if (j < N - 1)
            c_prime[j] = _upper_diagonal[j] / denom;
        R_prime[j] = (_right_hand_side[j]
                      - _lower_diagonal[j - 1] * R_prime[j - 1]) / denom;
    }

    // Back substitution
    Vector x(N);
    x[N - 1] = R_prime[N - 1];
    for (int j = static_cast<int>(N) - 2; j >= 0; --j)
        x[j] = R_prime[j] - c_prime[j] * x[j + 1];

    return x;
}

ImpliedVolatilitySurface.h

#pragma once
#include <vector>

using Vector = std::vector<double>;
using Matrix = std::vector<Vector>;

class ImpliedVolatilitySurface {
public:
    // maturities: strictly increasing, positive (years from today)
    // strikes:    non-decreasing, non-negative
    // market_implied_vols: matrix[i][j] = sigma*(T_i, K_j), all positive
    // risk_free_rate: continuously compounded, annualised
    ImpliedVolatilitySurface(const Vector& maturities,
                             const Vector& strikes,
                             const Matrix& market_implied_vols,
                             double risk_free_rate);

    // Interpolated/extrapolated sigma*(maturity, strike)
    double implied_volatility(double maturity, double strike) const;

    double risk_free_rate() const { return _risk_free_rate; }

private:
    bool ordered_vectors()    const;
    bool market_dimensions()  const;
    bool vols_positivity()    const;

    // Per-maturity cubic spline at a fixed maturity index
    double smile_implied_volatility(size_t maturity_index, double strike) const;

    // Populates _alpha, _beta, _gamma coefficient matrices
    void evaluate_cubic_spline_coefficients();

    // Cubic spline coefficients [maturity_index][strike_interval_index]
    Matrix _alpha_coefficients;
    Matrix _beta_coefficients;
    Matrix _gamma_coefficients;
    Vector _delta_strikes;          // delta_strikes[j] = K_{j+1} - K_j

    Vector _maturities;
    Vector _strikes;
    Matrix _market_implied_volatilities;
    double _risk_free_rate;
};

ImpliedVolatilitySurface.cpp (key methods)

#include "ImpliedVolatilitySurface.h"
#include "ThomasSolver.h"
#include <algorithm>
#include <cmath>
#include <stdexcept>

ImpliedVolatilitySurface::ImpliedVolatilitySurface(
    const Vector& maturities,
    const Vector& strikes,
    const Matrix& market_implied_vols,
    double risk_free_rate)
    : _maturities(maturities),
      _strikes(strikes),
      _market_implied_volatilities(market_implied_vols),
      _risk_free_rate(risk_free_rate)
{
    if (!ordered_vectors())
        throw std::invalid_argument("Maturities and strikes must be strictly ordered and positive");
    if (!market_dimensions())
        throw std::invalid_argument("market_implied_vols dimensions must be (M x N)");
    if (!vols_positivity())
        throw std::invalid_argument("All implied vols must be positive");

    for (size_t j = 0; j < _strikes.size() - 1; ++j)
        _delta_strikes.push_back(_strikes[j + 1] - _strikes[j]);

    evaluate_cubic_spline_coefficients();
}

double ImpliedVolatilitySurface::implied_volatility(double maturity, double strike) const {
    const size_t M = _maturities.size();

    // Flat extrapolation before first maturity (forward-moneyness adjusted)
    if (maturity <= _maturities[0]) {
        double K0 = strike * std::exp(_risk_free_rate * (_maturities[0] - maturity));
        return smile_implied_volatility(0, K0);
    }

    // Flat extrapolation after last maturity
    if (maturity >= _maturities[M - 1]) {
        double KM = strike * std::exp(_risk_free_rate * (_maturities[M - 1] - maturity));
        return smile_implied_volatility(M - 1, KM);
    }

    // Find the maturity bracket [T_i, T_{i+1}]
    size_t i = 0;
    while (maturity >= _maturities[i]) ++i;
    --i;  // now T_i <= maturity < T_{i+1}

    // Forward-moneyness adjusted strikes at T_i and T_{i+1}
    double K_left  = strike * std::exp(_risk_free_rate * (_maturities[i]     - maturity));
    double K_right = strike * std::exp(_risk_free_rate * (_maturities[i + 1] - maturity));

    // Total variances at bracket endpoints
    double sigma_l = smile_implied_volatility(i,     K_left);
    double sigma_r = smile_implied_volatility(i + 1, K_right);
    double w_left  = sigma_l * sigma_l * _maturities[i];
    double w_right = sigma_r * sigma_r * _maturities[i + 1];

    // Linear interpolation in total variance
    double alpha = (maturity - _maturities[i]) / (_maturities[i + 1] - _maturities[i]);
    double w     = w_left + alpha * (w_right - w_left);

    return std::sqrt(w / maturity);
}

double ImpliedVolatilitySurface::smile_implied_volatility(
    size_t mat_idx, double strike) const
{
    const size_t N = _strikes.size();

    // Left extrapolation: linear with slope = spline derivative at K_1
    if (strike <= _strikes[0]) {
        double deriv = _gamma_coefficients[mat_idx][0];
        return _market_implied_volatilities[mat_idx][0]
               + deriv * (strike - _strikes[0]);
    }

    // Right extrapolation: linear with slope = spline derivative at K_N
    if (strike >= _strikes[N - 1]) {
        size_t last = N - 2;
        double dK   = _delta_strikes[last];
        double deriv = 3.0 * _alpha_coefficients[mat_idx][last] * dK * dK
                     + 2.0 * _beta_coefficients[mat_idx][last]  * dK
                     +       _gamma_coefficients[mat_idx][last];
        return _market_implied_volatilities[mat_idx][N - 1]
               + deriv * (strike - _strikes[N - 1]);
    }

    // Cubic spline interpolation: find the interval K_j <= strike < K_{j+1}
    size_t j = 0;
    while (strike >= _strikes[j]) ++j;
    --j;

    double dx    = strike - _strikes[j];
    double alpha = _alpha_coefficients[mat_idx][j];
    double beta  = _beta_coefficients[mat_idx][j];
    double gamma = _gamma_coefficients[mat_idx][j];

    return alpha * dx * dx * dx
         + beta  * dx * dx
         + gamma * dx
         + _market_implied_volatilities[mat_idx][j];
}

void ImpliedVolatilitySurface::evaluate_cubic_spline_coefficients() {
    const size_t M = _maturities.size();
    const size_t N = _strikes.size();

    for (size_t i = 0; i < M; ++i) {
        const Vector& Y = _market_implied_volatilities[i];

        // Build the tridiagonal system for interior beta values
        // System size: (N-2) x (N-2), for j = 1, ..., N-2
        Vector lower, central, upper, rhs;
        for (size_t j = 0; j < N - 3; ++j) {
            lower.push_back(_delta_strikes[j + 1]);
            upper.push_back(_delta_strikes[j + 1]);
        }
        for (size_t j = 0; j < N - 2; ++j) {
            central.push_back(2.0 * (_delta_strikes[j] + _delta_strikes[j + 1]));
            rhs.push_back(3.0 * ((Y[j + 2] - Y[j + 1]) / _delta_strikes[j + 1]
                               - (Y[j + 1] - Y[j])     / _delta_strikes[j]));
        }

        ThomasSolver solver(lower, central, upper, rhs);
        Vector interior_beta = solver.solve();

        // Full beta vector: beta_1 = 0 (natural BC), interior, beta_N = 0
        Vector beta(N, 0.0);
        for (size_t j = 1; j < N - 1; ++j)
            beta[j] = interior_beta[j - 1];
        _beta_coefficients.push_back(beta);

        // Alpha and gamma from beta
        Vector alpha(N - 1), gamma(N - 1);
        for (size_t j = 0; j < N - 1; ++j) {
            alpha[j] = (beta[j + 1] - beta[j]) / (3.0 * _delta_strikes[j]);
            gamma[j] = (Y[j + 1] - Y[j]) / _delta_strikes[j]
                     - alpha[j] * _delta_strikes[j] * _delta_strikes[j]
                     - beta[j]  * _delta_strikes[j];
        }
        _alpha_coefficients.push_back(alpha);
        _gamma_coefficients.push_back(gamma);
    }
}

Validation

Exact Recovery at Nodes

The interpolated surface must reproduce input vols exactly:

sigma*(T_i, K_j) == market_implied_vols[i][j]   for all i, j

Arbitrage-Free Check: Calendar Spread

For any $K$ and $T_1 < T_2$:

$w(T_2, K) \geq w(T_1, K), \quad \text{i.e.,} \quad \sigma^*(T_2,K)^2 T_2 \geq \sigma^*(T_1,K)^2 T_1$

This is enforced automatically by linear interpolation in total variance.

Benchmark: Flat Surface

If the input surface is constant $\sigma^*(T_i, K_j) = \sigma_0$ for all $i, j$, then:

• implied_volatility(T, K) == sigma_0 for any $(T, K)$
• Local vol $\sigma_D(t, S) = \sigma_0$ everywhere (pure Black-Scholes)
• No numerical issues in Thomas solver (diagonal dominant, uniform spacing)

Numerical Example

Using the market surface below (equity index ATM smile, 3 maturities × 5 strikes):

At the query point $(T=0.375, K=95)$, the total-variance interpolation gives a vol between the 0.25 and 0.50 maturity smiles at the forward-adjusted strike.

Limitations

Arbitrage at smile boundary. Natural spline conditions ($\beta_1 = \beta_N = 0$) can produce negative second derivatives in the interior, which implies a negative butterfly spread — a static arbitrage. Arbitrage-free splines (e.g., Fengler 2009) add monotonicity and convexity constraints but require quadratic programming.

Forward-moneyness assumes constant $r$. The forward adjustment $K^{(i)} = K e^{r(T_i - T)}$ uses a flat risk-free rate. A term structure of rates requires replacing $r$ with the relevant forward rate for each interval.

Thomas solver is not parallelised. The per-maturity systems are independent and trivially parallelisable with std::for_each and std::execution::par. For large surfaces (100+ maturities × 500+ strikes), this matters.

The Dupire formula amplifies noise. Taking second derivatives of an interpolated surface amplifies any fitting error. For noisy market data, regularisation (e.g., Tikhonov) or SVI parametrisation of the input surface is required before computing local vols.

Interview Angle

Junior (L1): Why interpolate in total variance rather than in implied volatility directly?

Implied volatility is not a linear quantity in time. If you have vol = 20% at 6 months and vol = 22% at 12 months, a naïve linear interpolation gives 21% at 9 months — but this implies total variance $0.21^2 \times 0.75 = 0.033$, which is less than the 6-month total variance $0.20^2 \times 0.5 = 0.020$. Wait — that's fine. But in general, the no-arbitrage condition is $\partial w / \partial T \geq 0$, which is automatically satisfied when interpolating $w$ linearly but not when interpolating $\sigma$ linearly.

Senior (L2): Describe how the Thomas algorithm works and why it is $O(N)$.

A general dense $N \times N$ system requires $O(N^3)$ via LU decomposition. The tridiagonal structure means that Gaussian elimination only touches three diagonals — the forward sweep processes each row in constant time, and the back substitution does likewise. Total work: $O(N)$ multiplications and additions. The spline system is always diagonally dominant ($b_j = 2(\Delta_{j-1}+\Delta_j) > \Delta_{j-1} + \Delta_j$), so no pivoting is needed and the algorithm is numerically stable.

Researcher (L3): The Dupire formula uses implied vol derivatives computed by bump-reval. What are the accuracy requirements on the step sizes $\varepsilon_T$ and $\varepsilon_K$, and how do you balance truncation error against floating-point cancellation?

Central difference stencils give $O(h^2)$ truncation error. Floating-point cancellation begins to dominate when $h \lesssim \sqrt{\epsilon_{\text{mach}}} \approx 10^{-8}$. The optimal step for a central first derivative is $h_{\text{opt}} \approx \epsilon_{\text{mach}}^{1/3} \approx 6\times10^{-6}$. For a second derivative: $h_{\text{opt}} \approx \epsilon_{\text{mach}}^{1/4} \approx 10^{-4}$. In practice, these are set per-surface based on the smile curvature; typical values are $\varepsilon_K = 0.01S_0$ and $\varepsilon_T = 0.001$ years. A more robust approach uses complex-step differentiation ($\text{Im}[f(K+ih)] / h$), which eliminates cancellation entirely but requires complex arithmetic throughout the surface evaluation.