Black-Scholes基础模型公式：从理论到工程落地的实现与优化

1次阅读

共计 2628 个字符，预计需要花费 7 分钟才能阅读完成。

Black-Scholes 模型被誉为量化金融的 ”Hello World”，但真正把它从数学公式变成可靠的生产代码时，往往会遇到一些教科书不会告诉你的坑。最近在开发期权定价引擎时，我重新梳理了整个实现流程，总结出这份兼顾数学严谨性和工程实用性的指南。

先看几个实际案例：

某私募的定价系统在价外期权上频繁返回 NaN，最终发现是 d2 项出现极大负值导致 exp 函数溢出
做市商系统在开盘时出现响应延迟，溯源发现是 for 循环实现的希腊字母计算拖慢整体速度
回测结果与理论值偏差 3 个标准差，根源在于使用了精度不足的标准正态分布近似

这些问题的本质，是数学上的连续性与计算机的离散性之间的鸿沟。让我们从最基础的公式开始：

$$ C = S_0N(d_1) – Ke^{-rT}N(d_2) $$

其中：
$$ d_1 = \frac{\ln(S_0/K) + (r + \sigma^2/2)T}{\sigma\sqrt{T}} $$
$$ d_2 = d_1 – \sigma\sqrt{T} $$

import math
from scipy.stats import norm
import numpy as np

# 方法 1：scipy 官方实现（推荐）def N(x):
    return norm.cdf(x)

# 方法 2：Abramowitz 近似（无外部依赖）def N_approx(x):
    a = 0.2316419
    t = 1 / (1 + a * abs(x))
    b = 0.319381530
    c = -0.356563782
    d = 1.781477937
    e = -1.821255978
    f = 1.330274429
    y = (b*t + c*t**2 + d*t**3 + e*t**4 + f*t**5) * math.exp(-x*x/2)/math.sqrt(2*math.pi)
    return 1 - y if x < 0 else y

# 方法 3：直接积分（仅用于教学）def N_integral(x, steps=10000):
    dx = x/steps
    return sum(math.exp(-0.5*(i*dx)**2)*dx/math.sqrt(2*math.pi) for i in range(steps))

def black_scholes(S, K, T, r, sigma, option_type='call'):
    """
    S: 标的资产现价
    K: 行权价
    T: 到期时间 (年)
    r: 无风险利率
    sigma: 波动率
    option_type: 'call' 或 'put'
    """
    # 参数边界检查
    assert S > 0 and K > 0 and T >= 0 and sigma >= 0

    d1 = (math.log(S/K) + (r + 0.5*sigma**2)*T) / (sigma * math.sqrt(T))
    d2 = d1 - sigma * math.sqrt(T)

    if option_type == 'call':
        price = S * N(d1) - K * math.exp(-r*T) * N(d2)
    else:
        price = K * math.exp(-r*T) * N(-d2) - S * N(-d1)

    # 处理极端情况
    if math.isnan(price):
        return 0.0 if option_type == 'call' and S < K*0.01 else S

    return price

# 低效写法
prices = [black_scholes(s, K, T, r, sigma) for s in spot_prices]

# 高效向量化
S_array = np.array(spot_prices)
def black_scholes_vectorized(S, K, T, r, sigma):
    d1 = (np.log(S/K) + (r + 0.5*sigma**2)*T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S * norm.cdf(d1) - K * np.exp(-r*T) * norm.cdf(d2)

from numba import njit

@njit
def N_numba(x):
    a = 0.2316419
    t = 1 / (1 + a * abs(x))
    b = 0.319381530
    c = -0.356563782
    d = 1.781477937
    e = -1.821255978
    f = 1.330274429
    y = (b*t + c*t**2 + d*t**3 + e*t**4 + f*t**5) * math.exp(-x*x/2)/math.sqrt(2*math.pi)
    return 1 - y if x < 0 else y

@njit
def bs_numba(S, K, T, r, sigma):
    d1 = (math.log(S/K) + (r + 0.5*sigma**2)*T) / (sigma * math.sqrt(T))
    d2 = d1 - sigma * math.sqrt(T)
    return S * N_numba(d1) - K * math.exp(-r*T) * N_numba(d2)

浮点精度陷阱
比较价格时使用相对误差而非绝对误差
避免直接比较 a == b，应使用 math.isclose(a, b, rel_tol=1e-9)
蒙特卡洛替代场景
当遇到美式期权或路径依赖期权时，BS 模型不再适用

使用控制变量法减少模拟次数：

def monte_carlo(S0, K, T, r, sigma, n_sims=100000):
    # 用 BS 解作为控制变量
    bs_price = black_scholes(S0, K, T, r, sigma)
    z = np.random.normal(0, 1, n_sims)
    ST = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*z)
    payoff = np.maximum(ST - K, 0)
    mc_price = np.exp(-r*T) * payoff.mean()
    # 计算协方差调整
    cov = np.cov(payoff, ST)[0,1]
    var = np.var(ST)
    adj = cov/var * (np.mean(ST) - S0*np.exp(r*T))
    return mc_price - adj, bs_price