5. 样条插值基础：三次样条插值的数学原理与Python实现

各位好，我是老张。今天咱们聊聊三次样条插值。

说实话，在波动率曲面构建这个领域，样条插值是个绕不开的基础工具。我最早接触它是在做期权做市商系统的时候，当时需要把离散的波动率报价平滑成一条连续的曲线。试过线性插值，结果曲面看着像锯齿；试过多项式插值，又出现了龙格现象——两端剧烈震荡。后来一位前辈跟我说：「试试三次样条吧。」一试，效果出奇的好。

嗯，今天我就把这块内容掰开揉碎了讲给你听。

5.1 为什么需要样条插值？

先问个问题：给你一组离散点，比如不同行权价的隐含波动率，你怎么得到任意行权价上的波动率？

最简单的办法是线性插值。但线性插值有个毛病——它不光滑。在金融里，不光滑意味着导数不连续，而导数的变化直接关系到希腊字母的计算。你想想看，如果vega或者gamma突然跳变，交易员会疯掉的。

高次多项式插值呢？理论上可以做到光滑，但实际效果很差。我记得有一次用10次多项式去拟合一组波动率数据，结果在两端直接飞到了天上——这就是著名的龙格现象。

所以，我们需要一种折中方案：分段低次多项式，且保证连接处光滑。这就是样条插值的核心思想。

核心要点：三次样条插值用分段的三次多项式来拟合数据，在节点处保证函数值、一阶导数、二阶导数都连续。说白了，就是既光滑又稳定。

5.2 三次样条的数学原理

数学上，三次样条的定义是这样的：

假设我们有 n+1 个数据点 (x₀, y₀), (x₁, y₁), ..., (xₙ, yₙ)，且 x₀ < x₁ < ... < xₙ。

我们要构造 n 个三次多项式 S₀(x), S₁(x), ..., Sₙ₋₁(x)，每个多项式定义在区间 [xᵢ, xᵢ₊₁] 上。

每个三次多项式可以写成：

Sᵢ(x) = aᵢ + bᵢ(x - xᵢ) + cᵢ(x - xᵢ)² + dᵢ(x - xᵢ)³

这里有 4n 个未知数（每个多项式有 a, b, c, d 四个系数）。我们需要 4n 个方程才能解出来。

这些方程来自哪里？

  1. 插值条件：每个多项式必须通过左右两个端点。这给出 2n 个方程。
  2. 一阶导数连续：在内部节点处，左右多项式的一阶导数相等。这给出 n-1 个方程。
  3. 二阶导数连续：在内部节点处，左右多项式的二阶导数相等。这给出 n-1 个方程。
  4. 边界条件：还需要 2 个额外条件。最常用的是自然边界条件——两端二阶导数为 0。这给出 2 个方程。

加起来：2n + (n-1) + (n-1) + 2 = 4n。完美，方程数等于未知数个数。

个人经验：我在实际项目中几乎只用自然边界条件。为什么？因为它能保证曲线在两端不会「翘起来」。有一次我试了固定一阶导数的边界条件，结果曲线在两端出现了不自然的弯曲，交易员直接说「这曲线看着不对劲」。

5.3 求解三对角方程组

解这个方程组，最终会归结到一个三对角线性系统。具体推导过程我就不展开了，直接给结论：

设 Mᵢ = S''(xᵢ) 为二阶导数在节点处的值。对于自然边界条件，M₀ = Mₙ = 0。

对于内部节点，有：

μᵢ Mᵢ₋₁ + 2 Mᵢ + λᵢ Mᵢ₊₁ = dᵢ

其中：

hᵢ = xᵢ₊₁ - xᵢ
λᵢ = hᵢ / (hᵢ₋₁ + hᵢ)
μᵢ = 1 - λᵢ
dᵢ = 6 / (hᵢ₋₁ + hᵢ) * [(yᵢ₊₁ - yᵢ)/hᵢ - (yᵢ - yᵢ₋₁)/hᵢ₋₁]

这个方程组可以用追赶法高效求解，时间复杂度 O(n)。

避坑指南：我曾经在求解时直接用 numpy.linalg.solve 去解这个三对角矩阵，结果数据量一大（比如1000个点），速度慢得让人抓狂。后来改用 scipy.linalg.solve_banded，速度提升了两个数量级。记住：三对角矩阵一定要用专门算法。

5.4 Python 实现

好了，理论说完了，咱们直接上代码。我习惯从零实现一遍，这样能真正理解原理。当然，实际工作中直接用 scipy 的 CubicSpline 就行。

import numpy as np
import matplotlib.pyplot as plt

def natural_cubic_spline(x, y):
    """
    自然三次样条插值
    参数:
        x: 节点横坐标 (n+1 个点)
        y: 节点纵坐标
    返回:
        系数矩阵 [a, b, c, d]
    """
    n = len(x) - 1  # 区间个数
    h = np.diff(x)  # 区间长度
    
    # 构建三对角矩阵的系数
    A = np.zeros((n-1, n-1))
    rhs = np.zeros(n-1)
    
    for i in range(1, n):
        if i == 1:
            A[i-1, i-1] = 2 * (h[0] + h[1])
            A[i-1, i] = h[1]
            rhs[i-1] = 6 * ((y[2] - y[1])/h[1] - (y[1] - y[0])/h[0])
        elif i == n-1:
            A[i-1, i-2] = h[n-2]
            A[i-1, i-1] = 2 * (h[n-2] + h[n-1])
            rhs[i-1] = 6 * ((y[n] - y[n-1])/h[n-1] - (y[n-1] - y[n-2])/h[n-2])
        else:
            A[i-1, i-2] = h[i-1]
            A[i-1, i-1] = 2 * (h[i-1] + h[i])
            A[i-1, i] = h[i]
            rhs[i-1] = 6 * ((y[i+1] - y[i])/h[i] - (y[i] - y[i-1])/h[i-1])
    
    # 求解 M 值（二阶导数）
    M = np.zeros(n+1)
    M[1:-1] = np.linalg.solve(A, rhs)
    
    # 计算每个区间的系数
    a = y[:-1]
    b = np.zeros(n)
    c = M[:-1] / 2
    d = np.zeros(n)
    
    for i in range(n):
        b[i] = (y[i+1] - y[i]) / h[i] - h[i] * (M[i+1] + 2 * M[i]) / 6
        d[i] = (M[i+1] - M[i]) / (6 * h[i])
    
    return a, b, c, d

def evaluate_spline(x_eval, x_nodes, coeffs):
    """
    计算样条在给定点上的值
    """
    a, b, c, d = coeffs
    n = len(x_nodes) - 1
    
    result = np.zeros_like(x_eval)
    for i, xv in enumerate(x_eval):
        # 找到所在的区间
        idx = np.searchsorted(x_nodes, xv) - 1
        idx = max(0, min(idx, n-1))
        
        dx = xv - x_nodes[idx]
        result[i] = a[idx] + b[idx] * dx + c[idx] * dx**2 + d[idx] * dx**3
    
    return result

# 示例：用波动率数据测试
x_nodes = np.array([0.0, 0.25, 0.50, 0.75, 1.0])
y_nodes = np.array([0.20, 0.22, 0.25, 0.23, 0.21])  # 模拟的波动率

coeffs = natural_cubic_spline(x_nodes, y_nodes)

# 生成密集点用于绘图
x_dense = np.linspace(0, 1, 100)
y_dense = evaluate_spline(x_dense, x_nodes, coeffs)

# 绘图
plt.figure(figsize=(10, 6))
plt.plot(x_nodes, y_nodes, 'ro', label='原始数据点')
plt.plot(x_dense, y_dense, 'b-', label='三次样条插值')
plt.xlabel('行权价 (标准化)')
plt.ylabel('隐含波动率')
plt.title('三次样条插值示例')
plt.legend()
plt.grid(True)
plt.show()

实用建议：实际工作中，直接用 scipy.interpolate.CubicSpline 就好。但自己实现一遍有个好处——当遇到边界条件需要定制时，你知道底层在干什么。我去年做一个奇异期权的定价引擎时，就需要在边界处固定一阶导数为某个特定值，这时候自己实现的代码就派上用场了。

5.5 用 scipy 快速实现

当然，日常开发中我们直接用 scipy：

from scipy.interpolate import CubicSpline

# 同样的数据
x = np.array([0.0, 0.25, 0.50, 0.75, 1.0])
y = np.array([0.20, 0.22, 0.25, 0.23, 0.21])

# 创建三次样条对象（自然边界条件）
cs = CubicSpline(x, y, bc_type='natural')

# 计算任意点的值
x_new = np.linspace(0, 1, 100)
y_new = cs(x_new)

# 还可以计算一阶导数和二阶导数
y_deriv = cs(x_new, 1)  # 一阶导
y_deriv2 = cs(x_new, 2) # 二阶导

你看，scipy 一行代码就搞定了。但理解背后的原理，能让你在遇到问题时知道怎么调参。

5.6 知识结构总览

下面这张图总结了三次样条插值的核心逻辑：

三次样条插值知识结构 输入：离散数据点 (x₀, y₀), (x₁, y₁), ..., (xₙ, yₙ) 核心原理：分段三次多项式 每个区间 [xᵢ, xᵢ₊₁] 上定义一个三次多项式 节点处：函数值、一阶导、二阶导均连续 求解过程：构建三对角方程组 追赶法求解 → 得到各区间多项式系数 a, b, c, d 输出：光滑连续的插值曲线 数学原理 数值计算

5.7 实际应用中的注意事项

最后，分享几个我在实战中踩过的坑：

  • 节点分布不均匀：如果 x 轴上的点分布不均匀，样条可能会出现震荡。我建议在波动率曲面上，尽量让行权价的对数间距均匀一些。
  • 外推问题：三次样条不适合外推。一旦超出数据范围，曲线会按最后一段多项式的趋势延伸，结果可能完全不合理。所以，永远不要用样条做外推。
  • 数据噪声：如果原始数据有噪声，三次样条会完美地拟合噪声。这时候可以考虑用平滑样条（smoothing spline），它允许曲线不完全通过数据点，但更平滑。

一句话总结：三次样条插值是波动率曲面构建的基石。它用分段三次多项式在保证光滑的同时避免了高次多项式的震荡。理解它的原理，你就能在遇到问题时灵活调整——比如换边界条件、加单调性约束等。

好了，这一章就到这里。下一章我们会深入讨论如何用样条插值构建完整的波动率曲面，包括二维插值和曲面平滑的技巧。到时候见。

公众号:蓝海资料掘金营，微信deep3321

← 上一章 📖 返回目录 下一章 →