引言

这篇文章记录了如何使用特征函数及 FFT 为期权定价。整篇文章的起点是特征函数 $\phi_X(u)$。它包含了随机变量 $X_T$ 的全部分布信息，因此即使概率密度或期权价格写不出来，也可以从特征函数把需要的量反演回来。模型换成 Black–Scholes、Merton 或 Heston 时，主要改变的是 $\phi_X$；后面的定价步骤基本不变。下文统一记

$$ \phi_X(u)=\mathbb E[e^{iuX}],\qquad X_T=\log S_T. $$

$r$、$q$ 和 $T$ 分别表示无风险利率、连续股息率和到期时间。

从特征函数到价格有两条相连的路。第一条先用 Gil–Pelaez 公式求 $X_T$ 的尾概率，再把看涨期权拆成两个尾概率。第二条直接处理价格：先用 Carr–Madan 的阻尼让价格关于对数行权价可积，再做 Fourier 反演。只算少数行权价时，可以逐点做数值积分；需要同一期限下的一整条价格曲线时，则把同一组离散求和写成 FFT。直接积分和 FFT 因而不是两套定价原理，只是同一反演公式的两种计算方式。

第一条路可以压缩成

$$ \phi_X(u) \xrightarrow{\text{Gil--Pelaez}} \left\{ \mathbb Q(X_T>\log K),\, \mathbb Q^S(X_T>\log K) \right\}, $$

再代入

$$ C(K,T) =S_0e^{-qT}\mathbb Q^S(X_T>\log K) -Ke^{-rT}\mathbb Q(X_T>\log K). $$

第二条路先由特征函数构造 Carr–Madan 的被积函数

$$ \psi(v)= \frac{e^{-rT}\phi_X\!\left(v-i(\alpha+1)\right)} {\alpha^2+\alpha-v^2+i(2\alpha+1)v}, $$

然后直接反演价格：

$$ C(k)=\frac{e^{-\alpha k}}{\pi} \int_0^\infty \operatorname{Re}\!\left(e^{-ivk}\psi(v)\right)\,\mathrm dv, \qquad k=\log K. $$

把最后这个积分离散成求和后，对一个 $k$ 逐项相加就是直接积分；把 $k$ 排成等距网格并同时计算整组求和，就是 FFT。

从特征函数反推概率

设 $F(x)=\mathbb P(X\le x)$。如果 $x$ 是 $F$ 的连续点，Gil–Pelaez 反演公式为 1

$$ F(x)=\frac12-\frac1\pi\int_0^\infty \frac{\operatorname{Im}\!\left(e^{-iux}\phi_X(u)\right)}{u}\,\mathrm du. $$

相应地，

$$ \mathbb P(X>x)=\frac12+\frac1\pi\int_0^\infty \frac{\operatorname{Im}\!\left(e^{-iux}\phi_X(u)\right)}{u}\,\mathrm du. $$

原文忽略了一个细节：如果 $x$ 恰好是跳点，积分给出的不是 $F(x)$，而是

$$ \frac{F(x^-)+F(x)}2. $$

在 Black–Scholes、Merton 和 Heston 这类连续终值分布中不会碰到这个区别，但反演公式本身并不要求分布连续，所以在此还是写清楚。

Gil–Pelaez 反演公式的推导

取 $0\lt\varepsilon\lt R\lt\infty$，先考虑截断后的积分

$$ I_{\varepsilon,R}(x) =\frac1\pi\int_\varepsilon^R \frac{e^{iux}\phi_X(-u)-e^{-iux}\phi_X(u)}{iu}\,\mathrm du. $$

将 $\phi_X(u)=\int e^{iuz}\,\mathrm dF(z)$ 代进去，可得

$$ \begin{aligned} I_{\varepsilon,R}(x) &=\frac2\pi\int_{-\infty}^{\infty} \left[\int_\varepsilon^R\frac{\sin(u(x-z))}{u}\,\mathrm du\right]\mathrm dF(z). \end{aligned} $$

这里 $u$ 只在有限区间 $[\varepsilon,R]$ 上变化，被积函数绝对可积，所以交换积分次序没有问题。真正需要多考虑一步的是后面的 $\varepsilon\downarrow0$ 和 $R\uparrow\infty$。

Dirichlet 积分告诉我们

$$ \lim_{\varepsilon\downarrow0,\,R\uparrow\infty} \frac2\pi\int_\varepsilon^R\frac{\sin(uy)}u\,\mathrm du =\operatorname{sign}(y). $$

截断正弦积分是一致有界的，因此可以在概率测度下取极限。于是

$$ \begin{aligned} \lim I_{\varepsilon,R}(x) &=\mathbb E[\operatorname{sign}(x-X)]\\ &=\mathbb P(X\lt x)-\mathbb P(X\gt x)\\ &=F(x^-)+F(x)-1. \end{aligned} $$

对实值随机变量，$\phi_X(-u)=\overline{\phi_X(u)}$。令 $z=e^{-iux}\phi_X(u)$，那么

$$ \frac{\overline z-z}{iu}=-\frac{2\operatorname{Im}z}{u}. $$

代回 $I_{\varepsilon,R}$ 就得到前面的反演公式。

我在旧版本里用 $1/\varepsilon$ 控制被积函数。这个估计可以说明固定 $\varepsilon,R$ 时能够换序，却不能直接处理 $\varepsilon\downarrow0$ 的极限；这里应当用正弦积分的一致有界性。

从计算角度看，特征函数好用的原因很简单。若 $X=Y+Z$ 且 $Y,Z$ 独立，则

$$ \phi_X(u)=\phi_Y(u)\phi_Z(u). $$

通俗解释：实空间里的卷积在频率空间变成了乘法。很多跳扩散或仿射模型的密度很难写，特征函数却还能保持比较规整的指数形式。

从概率分布到期权价格

风险中性测度 $\mathbb Q$ 下，看涨期权价格是

$$ C(K,T)=e^{-rT}\mathbb E^{\mathbb Q}[(S_T-K)^+]. $$

如果已经知道 $X_T=\log S_T$ 的密度 $f_X$，当然可以直接算

$$ C(K,T)=e^{-rT}\int_{\log K}^{\infty}(e^x-K)f_X(x)\,\mathrm dx. $$

Black–Scholes 的 $X_T$ 是正态分布，积分甚至有闭式解，此时没必要特意改用特征函数。麻烦出现在更复杂的模型里：例如加入跳跃或随机波动率后，密度和期权价格未必容易写出，$X_T$ 的特征函数却可能从独立增量或仿射结构中算出来。特征函数确定了整个终值分布，因而可以把上面的期望转成 Fourier 积分。这里换的是计算方式，并没有换一套定价原理。若同一期限需要算很多行权价，后面的 FFT 还能重复利用同一组频率值。

因此，这类方法通常被称为半解析（semi-analytical）或半闭式（semi-closed form）定价。Black–Scholes 价格是严格的闭式解；Heston 等模型则往往把价格写成一个或几个一维积分。本文的直接积分和 FFT 是在数值计算这种半闭式表示：它既不抽样模拟大量路径，也不在时间—状态网格上求解 PDE，但仍然要处理积分截断、数值求积和插值误差。所以，特征函数法避免了 Monte Carlo 和 PDE，却仍是一种数值计算方法，不能笼统地称为闭式解。

展开支付函数：

$$ C=e^{-rT}\mathbb E^{\mathbb Q}[S_T\mathbf1_{\{S_T>K\}}] -Ke^{-rT}\mathbb Q(S_T>K). $$

第二项已经是尾概率。第一项还多了一个 $S_T$，可以通过换 numéraire 处理。定义股票测度 $\mathbb Q^S$：

$$ \frac{\mathrm d\mathbb Q^S}{\mathrm d\mathbb Q} =\frac{e^{-rT}S_T}{S_0e^{-qT}}, $$

于是

$$ C=S_0e^{-qT}\mathbb Q^S(X_T>\log K) -Ke^{-rT}\mathbb Q(X_T>\log K). $$

如果 $\phi(u)=\mathbb E^{\mathbb Q}[e^{iuX_T}]$，股票测度下的特征函数为

$$ \phi_S(u) =\mathbb E^{\mathbb Q^S}[e^{iuX_T}] =\frac{\phi(u-i)}{\phi(-i)}. $$

分别对 $\phi$ 和 $\phi_S$ 使用 Gil–Pelaez 公式，就能算出上式中的两个概率。到这里其实已经有了一套可用的特征函数定价方法。

需要留意的是，$\phi_S$ 用到了复数点 $u-i$。这不是形式上的代换，它要求相应的指数矩存在。后面 Carr–Madan 公式中的复数位移也有同样的问题。

Carr–Madan 的阻尼

如果只算几个行权价，上一节的方法已经够用。但 FFT 希望把不同的行权价排在一个规则网格上，因此更适合直接变换价格关于对数行权价的函数。

令 $k=\log K$，把看涨价格记作 $C(k)$。当 $k\to-\infty$ 时，$C(k)$ 不会衰减到零，不能直接对它做 Fourier 变换。Carr 和 Madan 的做法是先乘上一个阻尼项 2 ：

$$ c_\alpha(k)=e^{\alpha k}C(k),\qquad \alpha>0. $$

在 $\mathbb E[S_T^{\alpha+1}]\lt\infty$ 的条件下，阻尼后的价格有合适的可积性。按照本文的符号约定，其 Fourier 变换为

$$ \psi(v)= \frac{e^{-rT}\phi\!\left(v-i(\alpha+1)\right)} {\alpha^2+\alpha-v^2+i(2\alpha+1)v}. $$

再做反演：

$$ C(k)=\frac{e^{-\alpha k}}\pi \int_0^\infty\operatorname{Re}\!\left(e^{-ivk}\psi(v)\right)\,\mathrm dv. $$

$\alpha$ 并不是越大越好。首先，$\phi(v-i(\alpha+1))$ 必须存在；其次，$\alpha$ 太小时阻尼不够，太大时又会带来很悬殊的数值尺度。Black–Scholes 和正态跳跃的 Merton 模型有所有实数阶矩，因此这里比较宽松。换成厚尾模型或某些随机波动率模型时，需要先看特征函数的解析条带。

求解：直接积分与 FFT

先做直接积分

把积分截断在 $[0,U]$，令 $v_j=j\eta$、$\eta=U/N$。用梯形法离散：

$$ C(k)\approx\frac{e^{-\alpha k}}\pi \sum_{j=0}^{N}w_j \operatorname{Re}\!\left(e^{-iv_jk}\psi(v_j)\right)\eta, $$

其中 $w_0=w_N=1/2$，其余权重为 $1$。对应的 Python 实现如下：

def carr_madan_direct(strikes, cf, r, T, alpha=1.5, U=160.0, N=4096):
    v = np.linspace(0.0, U, N + 1)
    shifted = v - 1j * (alpha + 1.0)
    psi = np.exp(-r * T) * cf(shifted) / (
        alpha**2 + alpha - v**2 + 1j * (2 * alpha + 1) * v
    )
    k = np.log(np.atleast_1d(strikes))
    integrand = np.real(np.exp(-1j * np.outer(k, v)) * psi)
    return np.exp(-alpha * k) * np.trapezoid(integrand, v, axis=1) / np.pi

再把求和写成 DFT

若每个 $k$ 都单独求和，许多计算会重复。FFT 的做法是同时规定频率网格和对数行权价网格：

$$ \begin{aligned} v_j&=j\eta,\\ k_m&=-b+m\lambda,\\ j,m&=0,\ldots,N-1. \end{aligned} $$

取

$$ \lambda\eta=\frac{2\pi}{N}, $$

则

$$ e^{-iv_jk_m}=e^{ibv_j}e^{-i2\pi jm/N}. $$

右边第二项就是 DFT 的核。把 $\psi(v_j)$、梯形权重、$\eta$ 和 $e^{ibv_j}$ 预先乘在一起，调用一次 FFT，就能得到 $N$ 个行权价的价格。复杂度从逐点求和的 $O(N^2)$ 降到 $O(N\log N)$。

通常取

$$ b=\frac{N\lambda}{2}=\frac\pi\eta, $$

对数行权价网格于是覆盖 $[-b,b)$。又因为 $U=N\eta$，

$$ \lambda=\frac{2\pi}{U}. $$

这个关系对调参很重要。固定 $U$ 而只增加 $N$，会减小频率步长并扩大 $k$ 的覆盖范围，却不会让行权价网格变密。想减小 $\lambda$，还得增加 $U$。市场行权价一般也不会正好落在 FFT 网格上，所以最后还有一层插值误差。所以只算一个期权时，FFT 未必比直接积分划算。FFT 的优势在于同一期限下批量计算一串行权价。

用 Black–Scholes 和 Merton 检查结果

在 Black–Scholes 模型下，

$$ X_T=\log S_0+\left(r-q-\frac12\sigma^2\right)T+\sigma\sqrt T Z, \qquad Z\sim N(0,1), $$

所以

$$ \phi_{BS}(u)=\exp\!\left( iu\left[\log S_0+\left(r-q-\frac12\sigma^2\right)T\right] -\frac12\sigma^2Tu^2 \right). $$

它有闭式期权价格，适合检查 Fourier 公式中的正负号、贴现和股息。

再看 Merton 跳扩散。设对数跳幅 $Y\sim N(\mu_J,\delta_J^2)$，跳跃强度为 $\lambda_J$。记

$$ \kappa=\mathbb E[e^Y-1]=e^{\mu_J+\delta_J^2/2}-1. $$

风险中性特征函数是

$$ \begin{aligned} \phi_M(u)=\exp\Bigg\{& iu\left[\log S_0+ \left(r-q-\lambda_J\kappa-\frac12\sigma^2\right)T\right] -\frac12\sigma^2Tu^2\\ &+\lambda_JT\left(e^{iu\mu_J-\delta_J^2u^2/2}-1\right) \Bigg\}. \end{aligned} $$

$-\lambda_J\kappa$ 是风险中性漂移补偿，漏掉它以后折现股票价格就不再是鞅。Merton 模型还能对跳跃次数做条件化，将价格写成 Poisson 加权的对数正态价格之和。我在代码里用这个和式作为参考值。

Heston 一类仿射模型也可以接到同一个定价函数上，只要替换 $\phi(u)$。不过，实际写 Heston 特征函数时还要处理复对数的分支连续性，这里没有展开。它与波动率曲面校准的关系可参见 市场报价如何变成模型标价？ 。

数值结果

下图的参数为 $S_0=100$、$r=2\%$、$q=1\%$、$T=1$、$\sigma=22\%$；Carr–Madan 积分使用 $\alpha=1.5$、$U=160$、$N=4096$。Merton 模型另取 $\lambda_J=0.45$、$\mu_J=-0.12$、$\delta_J=0.25$。

Black–Scholes 与 Merton 模型下直接 Fourier 积分、FFT 和参考价格的比较，以及特征函数的频率衰减

Black–Scholes 直接积分相对闭式解的最大绝对误差约为 $7.8\times10^{-14}$，已经接近浮点精度；FFT 加上线性插值后，最大误差约为 $2.7\times10^{-2}$。在 Merton 模型下，FFT 相对 Poisson 混合参考值的最大误差约为 $2.3\times10^{-2}$。

单看左边的价格曲线，这些差异几乎看不出来。右上角的误差图则显示得很清楚：这里 FFT 的主要误差不在积分本身，而在较疏的行权价网格和随后的线性插值。单纯把 $N$ 调大不一定有用，因为固定 $U$ 时 $\lambda=2\pi/U$ 不变。

Fourier 定价数值实验 参考值、直接积分与 FFT 的计算都在浏览器本地完成。

左图是价格，右图是绝对误差的 $\log_{10}$。FFT 先在等距对数行权价 网格上计算，再插值到图中的行权价。

容易踩的坑

我在实现时遇到的问题主要有下面几个：

  • Fourier 符号没有统一。 有的资料用 $e^{iuX}$，有的用 $e^{-iuX}$。只要定义一致都可以，但特征函数、反演相位、复数位移和 FFT 方向必须一起改。
  • 把 $\log S_T$ 和 $\log(S_T/S_0)$ 混在一起。 二者差一个 $\log S_0$，漏掉后整个相位都会偏。
  • 漏掉贴现、股息或跳跃漂移补偿。 这类错误仍可能画出很光滑的价格曲线，肉眼并不好发现。
  • 只调 $N$。 $N,U,\eta,\lambda$ 同时决定频率截断、积分步长、行权价间距和覆盖范围，不能把 $N$ 当成唯一的精度参数。更系统的误差分析可见 Roger Lee 的文章 3 。
  • 没有检查价格边界。 至少应检查
$$ \max(S_0e^{-qT}-Ke^{-rT},0)\le C(K,T)\le S_0e^{-qT}, $$

以及看涨价格是否随 $K$ 递减、是否关于 $K$ 凸。深度虚值区域偶尔出现 $-10^{-12}$ 量级的负数通常是浮点误差；若偏差更大，就应该回头检查公式和网格。

现在我对这套方法的用法大致是：先用直接积分和一个有参考解的模型核对公式；若只算少量 $(K,T)$，直接积分已经足够；若要在同一期限上生成一整条行权价曲线，再使用 FFT。对于路径依赖期权，只有终值特征函数并不够，还是要回到 Monte Carlo 或其他数值方法。有关随机模拟的部分见 随机模拟与随机微分方程数值解 。

  1. J. Gil-Pelaez, “Note on the inversion theorem,” Biometrika, 38(3/4), 1951, 481–482. DOI . ↩︎

  2. Peter Carr and Dilip Madan, “Option valuation using the fast Fourier transform,” Journal of Computational Finance, 2(4), 1999. DOI . ↩︎

  3. Roger W. Lee, “Option pricing by transform methods: extensions, unification, and error control,” Journal of Computational Finance, 7(3), 2004. DOI ; corrected manuscript . ↩︎