Numerical Methods and Algorithms in Finance › Interpolation Techniques

Michael BrenndoerferNovember 1, 202512 min read
Quantitative FinanceAll chapters
Part 1 of 15Mathematical Foundations
Chapter 8Numerical Methods and Algorithms in Finance › Interpolation Techniques
Reading section 3 of 10 · Interpolation Techniques · 12 minSection 3 of 10 · 12 min

Book contents

Section 3 of 10 · 12 min read

Numerical Methods and Algorithms in Finance › Interpolation TechniquesOne section · Sections in this chapter3 / 10 ↓
Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

The linear and cubic-spline interpolation methods covered here construct continuous functions from discrete data points. In finance, you observe prices or rates at specific points, but need values at other locations. A yield curve might have observed rates at 3-month, 6-month, 1-year, 2-year, 5-year, and 10-year maturities, but you need to price a bond maturing in 3.5 years. The market provides data at certain locations, and interpolation fills in the values between them.

The challenge goes beyond convenience. Financial instruments exist at arbitrary maturities and strikes, not just at the discrete points where we have direct market observations. A corporate treasurer hedging a loan with a specific maturity, a portfolio manager pricing an off-the-run bond, or a derivatives trader valuing an exotic option at an unusual strike: all require rates or volatilities at points not directly quoted. Interpolation bridges this gap, but the choice of interpolation method affects pricing, hedging, and risk management in subtle but important ways.

Linear Interpolation

Linear interpolation connects adjacent data points with straight lines. Given points (x0,y0)(x_0, y_0) and (x1,y1)(x_1, y_1), the interpolated value at xx is:

y=y0+y1−y0x1−x0(x−x0)y = y_0 + \frac{y_1 - y_0}{x_1 - x_0}(x - x_0)

where:

  • (x0,y0)(x_0, y_0) and (x1,y1)(x_1, y_1): the two known data points bracketing xx
  • xx: the point at which we want to interpolate
  • y1−y0x1−x0\frac{y_1 - y_0}{x_1 - x_0}: the slope of the line connecting the two points
  • (x−x0)(x - x_0): horizontal distance from the left point

The formula starts at y0y_0 and adds the slope times the horizontal distance traveled. This geometric interpretation makes linear interpolation intuitive: we're simply walking along the straight line connecting two known points. The rate of change (slope) is constant within each segment. At a knot, it jumps only when the adjacent segment slopes differ.

This can be written as a weighted average:

y=(1−t)y0+ty1y = (1 - t)y_0 + ty_1

where:

  • t=x−x0x1−x0t = \frac{x - x_0}{x_1 - x_0}: interpolation parameter ranging from 0 to 1
  • When t=0t = 0: x=x0x = x_0 and y=y0y = y_0 (at the left endpoint)
  • When t=1t = 1: x=x1x = x_1 and y=y1y = y_1 (at the right endpoint)
  • When t=0.5t = 0.5: xx is the midpoint and yy is the average of y0y_0 and y1y_1

This weighted average form shows how linear interpolation works. The result is a convex combination of the two endpoint values, with weights determined by relative position. The parameter tt measures progress from the left point to the right point, normalized to the unit interval. This formulation generalizes naturally to higher dimensions and forms the basis for barycentric coordinates in computer graphics.

Linear interpolation is simple and computationally efficient, but its first derivative is discontinuous at a knot only when the adjacent segment slopes differ. If zero rates are interpolated linearly and then differentiated to obtain instantaneous forward rates, any such derivative discontinuities appear as jumps. Whether that representation is acceptable depends on the curve construction and its intended use.

In[16]:
Code
def linear_interpolate(x, x_data, y_data):
    """Perform piecewise linear interpolation."""
    x = np.asarray(x, dtype=float)
    x_data = np.asarray(x_data, dtype=float)
    y_data = np.asarray(y_data, dtype=float)

    if x_data.ndim != 1 or y_data.ndim != 1 or len(x_data) != len(y_data):
        raise ValueError("x_data and y_data must be one-dimensional and equal length")
    if not (
        np.all(np.isfinite(x))
        and np.all(np.isfinite(x_data))
        and np.all(np.isfinite(y_data))
    ):
        raise ValueError("Queries, knots, and values must be finite")
    if len(x_data) < 2 or np.any(np.diff(x_data) <= 0):
        raise ValueError("x_data must contain at least two strictly increasing values")

    # Handle single value case
    scalar_input = x.ndim == 0
    x = np.atleast_1d(x)
    if np.any((x < x_data[0]) | (x > x_data[-1])):
        raise ValueError("x lies outside the interpolation range")

    result = np.zeros_like(x, dtype=float)

    for i, xi in enumerate(x):
        # Find the interval containing xi
        idx = np.searchsorted(x_data, xi) - 1
        idx = np.clip(idx, 0, len(x_data) - 2)

        # Linear interpolation
        x0, x1 = x_data[idx], x_data[idx + 1]
        y0, y1 = y_data[idx], y_data[idx + 1]
        t = (xi - x0) / (x1 - x0)
        result[i] = (1 - t) * y0 + t * y1

    return result[0] if scalar_input else result

Cubic Spline Interpolation

Cubic spline interpolation addresses the smoothness problem by fitting piecewise cubic polynomials that are continuous in both first and second derivatives. The bending-strip analogy specifically describes the natural cubic spline, which minimizes an integrated squared-curvature functional subject to interpolation. The example below instead uses SciPy's default not-a-knot boundary condition; clamped and other boundary conditions define different cubic splines.

Given nn data points, cubic splines fit n−1n-1 cubic polynomials, each of the form:

Si(x)=ai+bi(x−xi)+ci(x−xi)2+di(x−xi)3S_i(x) = a_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3

where:

  • Si(x)S_i(x): cubic polynomial on the interval [xi,xi+1][x_i, x_{i+1}]
  • aia_i: constant term, equals yiy_i (the data value at xix_i)
  • bib_i: linear coefficient, controls the slope at xix_i
  • cic_i: quadratic coefficient, related to curvature
  • did_i: cubic coefficient, allows the curvature to vary across the interval

A cubic polynomial has four coefficients, and we have n−1n-1 such polynomials, giving 4(n−1)4(n-1) unknowns in total. To determine all these coefficients uniquely, we need 4(n−1)4(n-1) equations. These come from the interpolation and smoothness requirements that make splines so useful.

The coefficients are determined by requiring:

  • Interpolation: each polynomial passes through its endpoint data
  • Continuity: adjacent polynomials match at data points
  • Smoothness: first and second derivatives match at interior data points
  • Boundary conditions: natural and clamped conditions are common alternatives; SciPy's CubicSpline uses not-a-knot conditions by default, as in the example below

The interpolation conditions provide 2(n−1)2(n-1) equations (each polynomial must match the data at both its left and right endpoints). Derivative matching at the n−2n-2 interior points adds 2(n−2)2(n-2) more equations (one each for first and second derivatives). This gives 2(n−1)+2(n−2)=4n−62(n-1) + 2(n-2) = 4n - 6 equations. The two boundary conditions supply the remaining two equations needed to determine all coefficients uniquely.

This creates a curve that passes exactly through all data points, with continuous first and second derivatives. That numerical smoothness can help when downstream calculations differentiate the curve. The CubicSpline call below requests an interpolating cubic spline; it does not request a finance-specific forward-rate or no-arbitrage constraint. Validate those properties separately when they matter.

In[17]:
Code
from scipy.interpolate import CubicSpline

# Synthetic illustrative yield curve (maturity in years, continuously compounded zero rates in %)
maturities = np.array([0.25, 0.5, 1, 2, 3, 5, 7, 10, 20, 30])
zero_rates = np.array([4.5, 4.6, 4.7, 4.5, 4.3, 4.2, 4.3, 4.4, 4.6, 4.7])

# Create interpolators
linear_interp = lambda x: linear_interpolate(x, maturities, zero_rates)
cubic_spline = CubicSpline(maturities, zero_rates)

# Generate fine grid for plotting
x_fine = np.linspace(0.25, 30, 200)
y_linear = linear_interp(x_fine)
y_spline = cubic_spline(x_fine)

Let's visualize the difference between linear and cubic spline interpolation:

Out[18]:
Visualization
Line chart comparing linear and cubic spline interpolation of yield curve data points.
Comparison of linear interpolation and cubic spline interpolation for yield curve construction. The cubic spline produces a smooth curve while linear interpolation creates visible kinks at data points.

The cubic spline creates a smooth curve that follows the general shape of the data without the angular appearance of linear interpolation. Notice how the spline captures the inverted portion of the curve (rates dipping between 2-7 years) more naturally. The transition through the minimum is gradual, without the artificial corners that linear interpolation would produce at the 3-year and 5-year data points.

Out[19]:
Visualization
Line chart comparing zero rates from linear and cubic spline interpolation, including the synthetic input points.
Zero rates from linear and cubic spline interpolation appear similar.
Line chart comparing implied forward rates from linear versus cubic spline interpolation.
In this construction, linear interpolation produces derivative jumps while the cubic spline remains smooth.

Interpolation in Volatility Surfaces

Volatility surfaces present a harder interpolation problem because they are two-dimensional: volatility depends on both strike price and time to expiration. Options traders need to interpolate these surfaces to price options at strikes and expirations not directly quoted in the market. A yield curve requires fitting a function of maturity alone, but a volatility surface requires fitting a function of two variables simultaneously.

The two-dimensional nature adds complexity. The generic interpolator below fits the supplied grid; the code does not add positivity, calendar-arbitrage, or strike-arbitrage constraints. Treat those properties as separate validation checks, together with checks for unstable oscillations.

In[20]:
Code
from scipy.interpolate import RectBivariateSpline

# Synthetic illustrative volatility surface
# Strikes as moneyness (K/S)
strikes = np.array([0.80, 0.90, 0.95, 1.00, 1.05, 1.10, 1.20])
# Expirations in years
expirations = np.array([0.083, 0.25, 0.5, 1.0, 2.0])

# Implied volatilities (%) - constructed smile-like pattern with a higher high-strike wing
vol_surface = np.array(
    [
        [32, 28, 25, 24, 25, 28, 35],  # 1 month
        [30, 26, 24, 22, 23, 26, 32],  # 3 months
        [28, 24, 22, 21, 22, 24, 30],  # 6 months
        [26, 23, 21, 20, 21, 23, 28],  # 1 year
        [24, 22, 20, 19, 20, 22, 26],  # 2 years
    ]
)

# Create 2D interpolator
vol_interp = RectBivariateSpline(expirations, strikes, vol_surface)
In[21]:
Code
# Interpolate at specific points
test_expiry = 0.75  # 9 months
test_strike = 1.025  # 2.5% OTM call

# vol_interp was created in the previous code block
interpolated_vol = vol_interp(test_expiry, test_strike)[0, 0]
Out[22]:
Console
Interpolated volatility at T=0.75y, K/S=1.025:
  Implied volatility: 20.84%

The interpolated volatility of approximately 21% for this 9-month, slightly out-of-the-money option falls between the synthetic 6-month and 1-year input values, as expected. The bivariate spline smoothly interpolates both the strike and expiration dimensions simultaneously, producing a value consistent with the surrounding inputs in the synthetic volatility surface.

Let's visualize the full interpolated surface:

Out[23]:
Visualization
3D surface plot of implied volatility against strike moneyness and time to expiration.
Interpolated synthetic implied-volatility surface. The constructed inputs impose an asymmetric smile-like pattern: the high-strike wing is 2-3 volatility points above the low-strike wing, and the cross-strike range narrows at longer expirations.

The constructed surface illustrates an asymmetric smile-like pattern: the input volatilities rise away from at-the-money, the high-strike wing is 2-3 volatility points above the corresponding low-strike wing, and the cross-strike range narrows from 11 to 7 volatility points as expiration increases. The inputs are constructed in the chapter code rather than loaded from market observations, so the figure demonstrates interpolation behavior rather than an empirical law or a claim about market beliefs.

Choosing Interpolation Methods

The choice of interpolation method depends on the application.

Linear interpolation is appropriate when:

  • A simple interpolation rule is useful; benchmark a vectorized implementation for large lookup workloads rather than this teaching loop
  • The data is already dense enough that smoothness does not matter much
  • You want to avoid oscillation artifacts that can occur with higher-order methods

Cubic splines can be useful when:

  • C2C^2 smoothness is needed
  • Boundary behavior, overshoot, extrapolation, and economic constraints have been validated for the application
  • A smooth presentation is useful after those checks

Shape-preserving interpolation can help when the input data's monotonicity should be retained. For example, SciPy's PCHIP preserves that shape without overshooting the supplied values. This property alone does not enforce positive forwards or option-price no-arbitrage constraints; those require separate checks.