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");
}