第 21c 章 数值方法之三：有限差分法与方法比较

树方法从"价格怎样随机变动"出发，蒙特卡洛从"路径的期望"出发，有限差分法（finite difference methods）则从衍生品满足的偏微分方程出发：把微分方程中的导数替换成网格上的差商，得到一组差分方程，然后从到期日开始逐层向后求解。

有限差分法是工程和物理中求解 PDE 的标准工具，移植到金融后有几个特别之处：边界条件来自期权条款（到期收益、深度实值/虚值时的极限行为），美式期权在每一层都要与内在价值比较，而方程的系数依赖于标的价格本身。本章还会揭示一个很漂亮的联系：显式有限差分法就是三叉树，隐式有限差分法相当于每个节点有很多分支的"多叉树"。理解这一点，就把 21a 章的树方法和本章的 PDE 方法统一起来了。

本章最后综合比较树、蒙特卡洛、有限差分三类方法，作为原书第 21 章的总结。

学习目标

1. 把 BSM 偏微分方程离散到 \((t,S)\) 网格上，写出隐式和显式差分格式的系数，设置美式看跌期权的边界条件。
2. 理解隐式法每一步要解一个三对角线性方程组，会用追赶法（或 solve_banded）以 \(O(M)\) 的代价求解。
3. 解释显式法与三叉树的等价性，并说明显式法为什么会出现负"概率"和不收敛，以及如何用 \(\ln S\) 网格修正。
4. 了解跳格法和 Crank–Nicolson 方法，知道 Crank–Nicolson 在时间方向上是二阶精度。
5. 能根据衍生品的特征（美式/欧式、路径依赖、标的个数）选择合适的数值方法。

读前导读

这一章在解决什么问题。 期权价值 \(f\) 同时依赖时间 \(t\) 和股价 \(S\)。BSM 偏微分方程（PDE）说的是：在任何一个 \((t,S)\) 点上，\(f\) 对时间的变化率（theta）、对股价的一阶变化率（delta）和二阶变化率（gamma）之间必须满足一个固定关系，否则就有套利。你在 CFA 里见过"delta 对冲后的组合应赚无风险利率"，这个方程就是把那句话写成了公式。有限差分法的思路很朴素：既然方程是关于 delta、gamma、theta 的，那就在一张时间×股价的表格上，用相邻格子的差值去近似这些导数，把一个微分方程变成很多个代数方程，再从到期日一层一层往回解。

本章的第二个重点是把三种方法串起来。显式有限差分其实就是 21a 的三叉树，只是换了一种说法；隐式法则相当于"每个节点通向下一层所有节点"的树。最后一节给出三种方法的选型表，是整个第 21 章的总结，实务中挑定价方法时可以直接用。

需要先想起来的数学。

• 偏导数。 \(\frac{\partial f}{\partial S}\) 是"时间不变、只动股价时 \(f\) 的变化率"，即 delta；\(\frac{\partial^2f}{\partial S^2}\) 是 gamma；\(\frac{\partial f}{\partial t}\) 是"股价不变、只动时间时的变化率"，即 theta。下标写法 \(f_S\)、\(f_{SS}\)、\(f_t\) 是同一个意思。见 第 00 册第 05 章 多元微积分与优化。
• 差商与泰勒展开。 \(f(S+h)\approx f(S)+f'(S)h+\frac12f''(S)h^2\)。用它能看出前向差分 \(\frac{f(S+h)-f(S)}{h}\) 的误差约为 \(\frac12f''h\)，而中心差分 \(\frac{f(S+h)-f(S-h)}{2h}\) 的误差约为 \(\frac16f'''h^2\)，小得多。例：\(f=e^S\) 在 \(S=0\)，\(h=0.1\)：前向差分 1.0517，中心差分 1.0017，真值 1。见 第 00 册第 02 章 导数与泰勒展开。
• 线性方程组与三对角矩阵。 三对角矩阵只有主对角线和紧挨着的上下两条对角线非零。这种方程组可以像"多米诺骨牌"一样逐个消元，计算量与未知数个数成正比，而一般方程组要 \(O(M^3)\)。见 第 00 册第 06 章 线性代数速成。
• 链式法则。 换元 \(Z=\ln S\) 时，\(\frac{\partial f}{\partial S}=\frac{1}{S}\frac{\partial f}{\partial Z}\)，二阶导数还会多出一项。这是把 BSM 方程变成常系数方程的关键。见 第 00 册第 02 章。

怎么读这一章。 核心必读是"微分方程与网格""隐式有限差分法"（重点看差分近似、边界条件和求解思路）、"与三叉树的关系"和"负概率与不收敛"。系数 \(a_j,b_j,c_j\) 的代数整理可以跟着下面的讲解框走一遍，之后只记结构。\(\ln S\) 网格那一大段系数第一次可以只看结论"系数与 \(j\) 无关、等同三叉树"。跳格法只需知道名字；Crank–Nicolson 要知道它是"隐式和显式的平均、精度更高、实务最常用"。最后的三方法比较表务必读。

21.8 有限差分法

微分方程与网格

以支付股息率 \(q\) 的股票上的美式看跌为例，期权价值 \(f(t,S)\) 满足（原书式 17.6）

白话解释：把方程读成希腊字母：\(\Theta+(r-q)S\Delta+\frac12\sigma^2S^2\Gamma=rf\)。左边是持有期权一小段时间、股价按风险中性方式变动时期权价值的期望变化率：时间流逝带来 \(\Theta\)，股价漂移带来 \((r-q)S\Delta\)，股价波动通过凸性带来 \(\frac12\sigma^2S^2\Gamma\)（这一项和债券凸性带来的额外收益是同一个道理）。右边说这个期望变化率必须等于"期权价值 × 无风险利率"。所以方程就是"风险中性世界里期权的期望收益率等于 \(r\)"的微分版本。 注意这个方程对欧式和美式都一样，二者的区别只在边界条件和"每层与内在价值比较"这一步。

构造网格（原书图 21.15）：

• 时间方向：期限 \(T\) 分成 \(N\) 段，\(\Delta t=T/N\)，共 \(N+1\) 个时间点 \(0,\Delta t,2\Delta t,\dots,T\)；
• 价格方向：选一个足够高的 \(S_{max}\)，使看跌期权在这个价格上几乎没有价值；\(\Delta S=S_{max}/M\)，共 \(M+1\) 个价格点 \(0,\Delta S,2\Delta S,\dots,S_{max}\)。\(S_{max}\) 的选取应使当前股价恰好是网格点之一。

网格共 \((M+1)(N+1)\) 个点，点 \((i,j)\) 对应时间 \(i\Delta t\)、股价 \(j\Delta S\)，期权值记为 \(f_{i,j}\)。

隐式有限差分法

导数的差分近似。在内部点 \((i,j)\) 处，\(\partial f/\partial S\) 可以用

• 前向差分：\(\dfrac{f_{i,j+1}-f_{i,j}}{\Delta S}\) (21.22)
• 后向差分：\(\dfrac{f_{i,j}-f_{i,j-1}}{\Delta S}\) (21.23)

取二者的平均，得到更精确的中心差分：

时间导数用前向差分，联系 \(i\Delta t\) 与 \((i+1)\Delta t\)：

二阶导数：\((i,j+1)\) 处的后向差分 \(\frac{f_{i,j+1}-f_{i,j}}{\Delta S}\) 减去 \((i,j)\) 处的后向差分 \(\frac{f_{i,j}-f_{i,j-1}}{\Delta S}\)，再除以 \(\Delta S\)：

中心差分 (21.24) 和 (21.26) 的截断误差都是 \(O(\Delta S^2)\)，而前向/后向差分只有 \(O(\Delta S)\)——这就是取平均的好处（可以用泰勒展开验证）。

推导拆解：用泰勒展开验证（\(h=\Delta S\)，\(f',f'',f'''\) 都在 \(S\) 处取值）。 展开两侧：\(f(S+h)=f+f'h+\frac12f''h^2+\frac16f'''h^3+\cdots\)，\(f(S-h)=f-f'h+\frac12f''h^2-\frac16f'''h^3+\cdots\)。 前向差分：\(\frac{f(S+h)-f(S)}{h}=f'+\frac12f''h+\cdots\)，误差首项正比于 \(h\)，即 \(O(\Delta S)\)。 中心差分：两式相减，偶数次项（\(f\)、\(f''h^2\)）抵消，得 \(\frac{f(S+h)-f(S-h)}{2h}=f'+\frac16f'''h^2+\cdots\)，误差 \(O(\Delta S^2)\)。 二阶差分：两式相加，奇数次项抵消，得 \(f(S+h)+f(S-h)-2f=f''h^2+\frac{1}{12}f''''h^4+\cdots\)，除以 \(h^2\) 后误差 \(O(\Delta S^2)\)。 实际含义：网格加密一倍，中心差分的误差约降为四分之一，前向差分只降为一半。固定收益里算有效久期时"收益率上下各移一次"而不是"只往上移"，用的正是中心差分的这个优点。

差分方程。把这些近似代入 (21.21)，并令 \(S=j\Delta S\)：

其中 \(j=1,\dots,M-1\)，\(i=0,\dots,N-1\)。注意 \(\Delta S\) 全部约掉了，系数只依赖 \(j\)。整理得

推导拆解：从差分方程到 (21.27) 的整理过程。 第一步，约去 \(\Delta S\)：第二项 \(j\Delta S\cdot\frac{1}{2\Delta S}=\frac j2\)，第三项 \(j^2\Delta S^2\cdot\frac{1}{\Delta S^2}=j^2\)。所以方程只剩 \(j\)，不再含 \(\Delta S\)。 第二步，两边同乘 \(\Delta t\)：\(f_{i+1,j}-f_{i,j}+\frac12(r-q)j\Delta t(f_{i,j+1}-f_{i,j-1})+\frac12\sigma^2j^2\Delta t(f_{i,j+1}+f_{i,j-1}-2f_{i,j})=r\Delta t\,f_{i,j}\)。 第三步，把已知的 \(f_{i+1,j}\) 留在一边，未知的 \(f_{i,\cdot}\) 按 \(j-1,j,j+1\) 归类到另一边：\(f_{i,j-1}\) 的系数是 \(\frac12(r-q)j\Delta t-\frac12\sigma^2j^2\Delta t\)（即 \(a_j\)）；\(f_{i,j}\) 的系数是 \(1+\sigma^2j^2\Delta t+r\Delta t\)（即 \(b_j\)）；\(f_{i,j+1}\) 的系数是 \(-\frac12(r-q)j\Delta t-\frac12\sigma^2j^2\Delta t\)（即 \(c_j\)）。

这个方程把时刻 \(i\Delta t\) 的三个未知值和时刻 \((i+1)\Delta t\) 的一个已知值联系起来。由于空间导数是在未知的 \(i\) 层上取的，所以叫"隐式"。

边界条件。网格的三条边：

(21.28) 是到期收益；(21.29) 说股价为 0 时美式看跌立即行权，价值为 \(K\)；(21.30) 说股价很高时看跌期权一文不值。

求解。先处理 \(T-\Delta t\) 这一层：

右边由 (21.28) 已知，又有 \(f_{N-1,0}=K\)（21.32）、\(f_{N-1,M}=0\)（21.33），于是得到 \(M-1\) 个方程、\(M-1\) 个未知数 \(f_{N-1,1},\dots,f_{N-1,M-1}\)。

这个方程组的系数矩阵是三对角的，无须求逆矩阵：由 \(j=1\) 的方程把 \(f_{N-1,2}\) 用 \(f_{N-1,1}\) 表示，代入 \(j=2\) 的方程把 \(f_{N-1,3}\) 用 \(f_{N-1,1}\) 表示……直到用 \(j=M-1\) 的方程解出 \(f_{N-1,1}\)，再回代得到其余各值。这就是追赶法（Thomas 算法），每层计算量为 \(O(M)\)。

白话解释：为什么要"解方程组"，以及为什么不贵。隐式法的每个方程都含三个未知数，单看一个方程解不出来，但 \(M-1\) 个方程首尾相连，加上两端边界已知，整体就能解。写成矩阵形式 \(A\mathbf f_i=\mathbf f_{i+1}\)，\(A\) 的每一行只有三个非零数，挨在对角线附近。追赶法先从上往下依次消去每行的"左邻"（追），再从下往上回代（赶），每个未知数只处理常数次，所以总量与 \(M\) 成正比。如果用一般的高斯消元或求逆矩阵，计算量是 \(O(M^3)\)，\(M=400\) 时相差上万倍。

解出后，把每个 \(f_{N-1,j}\) 与内在价值 \(K-j\Delta S\) 比较：若小于内在价值，则在该点提前行权是最优的，令 \(f_{N-1,j}=K-j\Delta S\)。\(T-2\Delta t\) 等各层依此类推，最终得到 \(f_{0,1},\dots,f_{0,M-1}\)，其中对应 \(S_0\) 的那个就是所求期权价格。总计算量为 \(O(MN)\)。

与控制变量结合。和树一样，可以用同一网格为一个有解析解的相似期权（如欧式看跌）估值，再用 (21.20) 修正。

例 21.10：隐式法为美式看跌定价

用隐式法为例 21.1 的美式看跌定价（\(S_0=K=50\)，\(r=10\%\)，\(\sigma=40\%\)，\(T=5/12\)）。取 \(M=20\)，\(N=10\)，\(\Delta S=5\)：股价从 0 到 100 每 5 美元一格，时间每半个月一格。原书表 21.4 给出了完整网格。

• 美式期权价格 4.07；
• 同一网格上的欧式期权价格 3.91；
• BSM 欧式真值 4.08；
• 控制变量估计 \(4.07+(4.08-3.91)=4.24\)。

表 21.4 中，股价不超过 40 的行在各时点都等于内在价值——这是提前行权区域。例如股价为 40 时，剩余期限在 3 个月以内的各列都是 10.00。

显式有限差分法

隐式法非常稳健：\(\Delta S\) 和 \(\Delta t\) 趋于零时，它总是收敛到微分方程的解（一般规则是让 \(\Delta S\) 与 \(\sqrt{\Delta t}\) 成比例地趋于零）。缺点是每一步都要解 \(M-1\) 个联立方程。

如果假设点 \((i,j)\) 处的 \(\partial f/\partial S\) 和 \(\partial^2f/\partial S^2\) 与点 \((i+1,j)\) 处相同：

差分方程就变成

（若对 \(\partial f/\partial t\) 改用后向差分而非前向差分，同样得到显式法。）

原书图 21.16 对比了两种方法的结构：隐式法把 \(i\Delta t\) 层的三个值与 \((i+1)\Delta t\) 层的一个值联系起来；显式法把 \(i\Delta t\) 层的一个值与 \((i+1)\Delta t\) 层的三个值联系起来，可以逐点直接计算，不需要解方程组。

例 21.11：显式法为美式看跌定价

用显式法、同样的 \(M=20\)、\(N=10\)、\(\Delta S=5\)，原书表 21.5 给出期权价格 4.26。但表的左上角（高股价、早时刻）出现了负数和不一致——一个看跌期权的价值不可能为负。原因见下文。

变量替换：对 \(\ln S\) 建网格

标的服从几何布朗运动时，以 \(Z=\ln S\) 为变量更高效。由链式法则，(21.21) 变为常系数方程

网格对 \(Z\) 等距（对 \(S\) 是等比的）。

推导拆解：链式法则如何把系数变成常数。 令 \(S=e^Z\)。一阶：\(\frac{\partial f}{\partial Z}=\frac{\partial f}{\partial S}\cdot\frac{dS}{dZ}=S\frac{\partial f}{\partial S}\)，所以 \(S\frac{\partial f}{\partial S}=\frac{\partial f}{\partial Z}\)。 二阶：对上式再对 \(Z\) 求导，\(\frac{\partial^2f}{\partial Z^2}=\frac{\partial}{\partial Z}\left(S\frac{\partial f}{\partial S}\right)=S\frac{\partial f}{\partial S}+S^2\frac{\partial^2f}{\partial S^2}\)（乘积法则，并再用一次 \(\frac{dS}{dZ}=S\)）。于是 \(S^2\frac{\partial^2f}{\partial S^2}=\frac{\partial^2f}{\partial Z^2}-\frac{\partial f}{\partial Z}\)。 代回 (21.21)：\((r-q)\frac{\partial f}{\partial Z}+\frac12\sigma^2\left(\frac{\partial^2f}{\partial Z^2}-\frac{\partial f}{\partial Z}\right)\)，合并一阶项得 \(\left(r-q-\frac{\sigma^2}{2}\right)\frac{\partial f}{\partial Z}\)。原来系数里的 \(S\) 和 \(S^2\) 都被吸收掉了。 这里又出现了 \(-\sigma^2/2\)，与 21b 章模拟 \(\ln S\) 时的修正项是同一个来源。

隐式法：

显式法：

这些系数与 \(j\) 无关，这是变量替换的主要好处。多数情况下 \(\Delta Z=\sigma\sqrt{3\Delta t}\) 是好的选择。

与三叉树的关系

显式有限差分法等价于三叉树（而隐式法等价于每个节点有 \(M+1\) 个分支的多叉树）。

把 (21.34) 中括号里的三项（去掉贴现因子 \(1/(1+r\Delta t)\)）解释为概率：

（原书图 21.17。）验证：

• 三者之和为 1；
• 股价的期望增量为 \((r-q)j\Delta t\cdot\Delta S=(r-q)S\Delta t\)，正是风险中性世界中的期望变化；
• \(\Delta t\) 很小时，增量的方差为 \(\sigma^2j^2\Delta t\cdot\Delta S^2=\sigma^2S^2\Delta t\)，与 \(S\) 的随机过程一致。

所以 (21.34) 的含义就是：\(f_{i,j}\) 等于下一步三个可能值的风险中性期望，按无风险利率贴现（贴现因子 \(1/(1+r\Delta t)\approx e^{-r\Delta t}\)）。这和三叉树的倒推完全一样。

推导拆解：表中三条验证逐条算一遍（记 \(p_d,p_m,p_u\) 为表中三项）。 和为 1：\(p_d+p_u=\sigma^2j^2\Delta t\)（漂移项一正一负抵消），再加 \(p_m=1-\sigma^2j^2\Delta t\)，正好为 1。 期望增量：\((+\Delta S)p_u+(-\Delta S)p_d=\Delta S(p_u-p_d)=\Delta S\cdot(r-q)j\Delta t\)，而 \(j\Delta S=S\)，所以等于 \((r-q)S\Delta t\)。 方差：\(E[(\text{增量})^2]=\Delta S^2(p_u+p_d)=\Delta S^2\sigma^2j^2\Delta t=\sigma^2S^2\Delta t\)；方差等于它减去期望增量的平方，后者是 \(\Delta t^2\) 量级，\(\Delta t\) 小时可以忽略。 结论：显式法不是一个"碰巧形似"树的公式，而是在网格上精确匹配了股价变化前两阶矩的三叉树，和 21a 章构造树的原则完全一致。

负概率与不收敛。三叉树要表现良好，三个"概率"都必须为正。在例 21.11 中，中间概率 \(1-\sigma^2j^2\Delta t=1-0.16\times j^2\times0.04167\) 在 \(j\ge13\)（即 \(S\ge65\)）时变为负数。这就解释了表 21.5 左上角的负价格和不一致。这是显式法的主要问题：概率可能为负，结果不一定收敛。数值分析的语言是：显式法是条件稳定的，需要 \(\sigma^2j^2\Delta t\le1\)，即 \(\Delta t\le1/(\sigma^2M^2)\)——价格网格越密，时间步必须越小（Hull & White 1990 给出了解决方法）。

白话解释：为什么高股价处先出问题。在 \(S\) 等距网格上，每格宽度 \(\Delta S\) 固定，但股价在 \(\Delta t\) 内的典型波动是 \(\sigma S\sqrt{\Delta t}\)，股价越高波动越大。当这个波动超过一格宽度太多（\(\sigma j\sqrt{\Delta t}>1\)），"只能走到相邻格"的三叉树就装不下这么大的方差，只好让中间概率变成负数来凑。\(\ln S\) 网格让格子宽度随股价按比例放大，每个位置的"波动 / 格宽"都一样，问题就消失了。 "条件稳定"与"无条件稳定"：前者指只有满足步长条件时误差才不会被逐层放大；后者指任何步长都不会放大误差（但精度仍取决于步长）。隐式法和 Crank–Nicolson 属于后者。

\(\ln S\) 网格解决问题。变量替换后，\(Z\) 下降 \(\Delta Z\)、不变、上升 \(\Delta Z\) 的概率就是 (21.37)–(21.39) 括号里的项，对应股价 \(S\to Se^{-\Delta Z}\)、\(S\)、\(Se^{\Delta Z}\)。这些概率与 \(j\) 无关。令 \(\Delta Z=\sigma\sqrt{3\Delta t}\)，中间概率为 \(1-\sigma^2\Delta t/(3\sigma^2\Delta t)=2/3\)，上下概率为 \(\frac16\pm\frac{\Delta t}{2\Delta Z}(r-q-\frac{\sigma^2}{2})\)，和 21a 章的三叉树完全相同，而且在 \(\Delta t\) 足够小时总是正的。

其他有限差分法

跳格法（hopscotch method）：在网格节点间交替使用显式和隐式计算（原书图 21.18）。每个时刻先在"显式节点"（E）上按常规计算；然后处理"隐式节点"（I），此时它们的相邻节点值已经算出，隐式方程只含一个未知数，无须解联立方程组。

Crank–Nicolson 方法：对 \(\frac{f_{i+1,j}-f_{i,j}}{\Delta t}\) 的估计取隐式法与显式法的平均——即空间导数取 \(i\) 层与 \(i+1\) 层的平均。这相当于在 \((i+\frac12)\Delta t\) 处做中心差分，时间方向的截断误差从 \(O(\Delta t)\) 降到 \(O(\Delta t^2)\)，同时保持无条件稳定。每步仍需解一个三对角方程组。Crank–Nicolson 是实务中最常用的 PDE 定价格式。它的一个已知缺点是：在收益函数有折点（如执行价处）时，最初几步可能产生小幅振荡，常见的补救是在开始时先走几步隐式法（Rannacher 平滑）。

白话解释：为什么"取平均"就能把精度提高一阶。隐式法在第 \(i\) 层取空间导数，显式法在第 \(i+1\) 层取，两者都偏向区间的一端，误差都是 \(O(\Delta t)\) 且方向相反。取平均相当于在区间中点 \((i+\frac12)\Delta t\) 估计，这和上面"中心差分比前向差分准"是同一个道理，只是这次用在时间方向上。可以类比为：估计一年的平均利率，用年中的利率比用年初或年末的利率更准。

有限差分法的应用

• 可以处理与树相同类型的问题，适用于美式与欧式期权；
• 难以处理收益依赖标的历史路径的情形；
• 以大幅增加计算时间为代价可以处理多个状态变量（网格变成多维）；
• 希腊字母：delta、gamma、theta 可以直接由网格上的 \(f_{i,j}\) 计算（与树类似）；vega 需要微调波动率并在同一网格上重算。

三种方法的综合比较（原书第 21 章小结）

原书第 21 章介绍了三种不能用解析公式时的估值方法：

• 树：每个 \(\Delta t\) 股价乘以 \(u\) 或 \(d\)，参数使风险中性世界中价格变化的均值和标准差正确；从到期日向后倒推；美式期权在每个节点取立即行权与继续持有的较大者。
• 蒙特卡洛：在风险中性世界中随机抽样路径，计算收益并按无风险利率贴现，样本均值即估值。
• 有限差分：把衍生品满足的 PDE 转为差分方程，从到期日向后求解。显式法在功能上等同于三叉树；隐式法更复杂，但无须特别措施就能保证收敛。

选择取决于衍生品的特征和所需的精度。一个经验法则是：能写出 PDE 且维度 ≤ 2 时，用有限差分或树；有路径依赖或维度 ≥ 3 时，用蒙特卡洛。

金融直觉：可以把三种方法和你熟悉的估值方法对应起来。树和有限差分像"从终值往回贴现的逐期 DCF"，每期都能插入管理层决策（行权、赎回、转换），所以擅长带选择权的产品；代价是要为每个可能状态都留一个格子，标的一多格子数就指数爆炸。蒙特卡洛像"情景分析"，生成大量未来情景、各自算收益再平均，情景里装什么变量都行，但每条情景只知道过去、不知道未来，所以难以做"现在要不要行权"这种需要未来信息的决策。表中"收敛速度"一行的 \(O(\cdot)\) 表示误差随网格或样本量缩小的速度，例如 \(O(1/N)\) 指步数翻倍、误差约减半。

量化实战

本章内容在量化交易中的用途

1. PDE 定价引擎。生产级的期权定价库通常用 Crank–Nicolson 在 \(\ln S\) 网格上求解，加上非均匀网格（在执行价和当前价附近加密）和 Rannacher 平滑。障碍期权（在障碍处设边界条件）、可转换债券（转换、赎回、回售构成的复杂边界）是有限差分的典型应用。
2. 三对角求解器。隐式法和 Crank–Nicolson 每步解一个三对角方程组，复杂度 \(O(M)\)。这类求解器在量化系统中还用于三次样条插值（收益率曲线、波动率曲面）、局部波动率模型校准等。
3. 稳定性检查。显式法的条件 \(\Delta t\le1/(\sigma^2M^2)\) 是一个可以写进代码的断言。任何显式时间推进的模拟（包括一些风险中性密度的演化）都应检查它。
4. 方法一致性。同一个美式期权，用树、显式 \(\ln S\) 网格、隐式网格应当收敛到同一个值。生产系统常用一种方法定价、另一种方法做独立验证（模型验证团队的标准做法）。

Python 示例：隐式、显式与 Crank–Nicolson

下面的代码复算例 21.10 和例 21.11，展示显式法的负价格和负概率，检查隐式法随网格加密的收敛，验证 \(\ln S\) 网格上的显式法（即三叉树）的稳定性，最后比较隐式法与 Crank–Nicolson 在欧式看跌上的精度。三对角方程组用 scipy.linalg.solve_banded 求解。

import numpy as np
from scipy.linalg import solve_banded
from scipy.stats import norm

def bs_put(S, K, T, r, sig):
    d1 = (np.log(S/K)+(r+0.5*sig**2)*T)/(sig*np.sqrt(T)); d2 = d1-sig*np.sqrt(T)
    return K*np.exp(-r*T)*norm.cdf(-d2) - S*norm.cdf(-d1)

def fd_put(S0, K, T, r, sig, Smax, M, N, q=0.0, method="implicit", american=True):
    """S 等距网格上的有限差分（原书 21.8 节）。返回 S0 处的期权价格和 t=0 一列。"""
    dS, dt = Smax/M, T/N
    j = np.arange(M+1); S = j*dS
    f = np.maximum(K - S, 0.0)                      # (21.28) 到期边界
    jj = j[1:M]
    if method == "implicit":                        # (21.27) 的系数
        a = 0.5*(r-q)*jj*dt - 0.5*sig**2*jj**2*dt
        b = 1 + sig**2*jj**2*dt + r*dt
        c = -0.5*(r-q)*jj*dt - 0.5*sig**2*jj**2*dt
        ab = np.zeros((3, M-1))                     # 三对角矩阵的带状存储
        ab[0, 1:] = c[:-1]; ab[1] = b; ab[2, :-1] = a[1:]
    else:                                           # (21.34) 的系数
        a = (-0.5*(r-q)*jj*dt + 0.5*sig**2*jj**2*dt)/(1+r*dt)
        b = (1 - sig**2*jj**2*dt)/(1+r*dt)
        c = (0.5*(r-q)*jj*dt + 0.5*sig**2*jj**2*dt)/(1+r*dt)
    for i in range(N-1, -1, -1):
        new = np.empty_like(f)
        new[0], new[M] = K, 0.0                      # (21.29)(21.30)
        if method == "implicit":
            rhs = f[1:M].copy()
            rhs[0] -= a[0]*K                        # 已知边界移到右边
            new[1:M] = solve_banded((1, 1), ab, rhs)
        else:
            new[1:M] = a*f[0:M-1] + b*f[1:M] + c*f[2:M+1]
        if american:
            new = np.maximum(new, K - S)            # 与内在价值比较
        f = new
    return np.interp(S0, S, f), S, f

args = dict(S0=50, K=50, T=5/12, r=0.10, sig=0.40, Smax=100, M=20, N=10)
fA, S, col = fd_put(**args)
fE, _, _ = fd_put(**args, american=False)
fBS = bs_put(50, 50, 5/12, 0.10, 0.40)
print("例21.10 隐式: 美式 %.4f  欧式 %.4f  BSM %.4f  控制变量 %.4f" % (fA, fE, fBS, fA + fBS - fE))
fX, S, colX = fd_put(**args, method="explicit")
print("例21.11 显式: 美式 %.2f" % fX)
print("  显式 t=0 列中 S=75..100:", np.round(colX[15:], 2))
print("  中间概率 1-sig^2 j^2 dt <0 的 j:", [j for j in range(1, 20) if 1 - 0.16*j*j*(5/12/10) < 0])

# 网格加密：隐式法收敛
print("\n隐式法网格加密（Smax=100）:")
for M, N in [(20, 10), (40, 40), (100, 200), (200, 800), (400, 3200)]:
    print(f"  M={M:4d} N={N:5d}  {fd_put(50,50,5/12,0.10,0.40,100,M,N)[0]:.4f}")

# ln S 网格上的显式法，dZ = sig*sqrt(3 dt)：等价于三叉树，概率恒正
def fd_logS_explicit(S0, K, T, r, sig, N, q=0.0):
    dt = T/N; dZ = sig*np.sqrt(3*dt); nu = r - q - 0.5*sig**2
    pd = -dt/(2*dZ)*nu + dt/(2*dZ**2)*sig**2      # (21.37) 括号内
    pm = 1 - dt/dZ**2*sig**2                       # (21.38) 括号内 = 2/3
    pu = dt/(2*dZ)*nu + dt/(2*dZ**2)*sig**2       # (21.39) 括号内
    disc = 1/(1 + r*dt)
    Z = np.log(S0) + dZ*np.arange(-N, N+1)
    f = np.maximum(K - np.exp(Z), 0)
    for i in range(N-1, -1, -1):
        Z = Z[1:-1]
        f = np.maximum(disc*(pd*f[:-2] + pm*f[1:-1] + pu*f[2:]), K - np.exp(Z))
    return f[0], (pd, pm, pu)
for N in (50, 200, 1000):
    v, probs = fd_logS_explicit(50, 50, 5/12, 0.10, 0.40, N)
    print(f"ln S 显式 N={N:4d}: {v:.4f}   概率 {np.round(probs, 4)}")

# θ 格式：θ=1 隐式，θ=0.5 Crank–Nicolson；欧式看跌与 BSM 比较精度
def fd_theta(S0, K, T, r, sig, Smax, M, N, theta=0.5, american=False):
    dS, dt = Smax/M, T/N
    S = np.arange(M+1)*dS; jj = np.arange(1, M)
    lo = 0.5*sig**2*jj**2 - 0.5*r*jj            # 空间算子 L f = lo*f[j-1] + di*f[j] + up*f[j+1]
    di = -sig**2*jj**2 - r
    up = 0.5*sig**2*jj**2 + 0.5*r*jj
    ab = np.zeros((3, M-1))
    ab[0, 1:] = -theta*dt*up[:-1]; ab[1] = 1 - theta*dt*di; ab[2, :-1] = -theta*dt*lo[1:]
    f = np.maximum(K - S, 0.0)
    for i in range(N-1, -1, -1):
        b_new = K if american else K*np.exp(-r*(T - i*dt))   # S=0 边界：欧式看跌为 K 的贴现值
        Lf = lo*f[:-2] + di*f[1:-1] + up*f[2:]
        rhs = f[1:-1] + (1-theta)*dt*Lf
        rhs[0] += theta*dt*lo[0]*b_new
        new = np.empty_like(f); new[0], new[M] = b_new, 0.0
        new[1:-1] = solve_banded((1, 1), ab, rhs)
        if american: new = np.maximum(new, K - S)
        f = new
    return np.interp(S0, S, f)

print("\n欧式看跌误差（M=200, Smax=150）: N   隐式误差    CN误差")
ex = bs_put(50, 50, 5/12, 0.10, 0.40)
for N in (10, 20, 40, 80):
    e1 = fd_theta(50, 50, 5/12, 0.10, 0.40, 150, 200, N, theta=1.0) - ex
    e2 = fd_theta(50, 50, 5/12, 0.10, 0.40, 150, 200, N, theta=0.5) - ex
    print(f"   N={N:3d}  {e1:+.5f}   {e2:+.5f}")

关键输出：

例21.10 隐式: 美式 4.0672  欧式 3.9112  BSM 4.0760  控制变量 4.2320
例21.11 显式: 美式 4.26
  显式 t=0 列中 S=75..100: [ 0.45 -0.13  0.28 -0.11  0.06  0.  ]
  中间概率 1-sig^2 j^2 dt <0 的 j: [13, 14, 15, 16, 17, 18, 19]

隐式法网格加密（Smax=100）:
  M=  20 N=   10  4.0672
  M=  40 N=   40  4.2274
  M= 100 N=  200  4.2731
  M= 200 N=  800  4.2813
  M= 400 N= 3200  4.2835
ln S 显式 N=  50: 4.2599   概率 [0.1653 0.6667 0.168 ]
ln S 显式 N= 200: 4.2783   概率 [0.166  0.6667 0.1673]
ln S 显式 N=1000: 4.2830   概率 [0.1664 0.6667 0.167 ]

欧式看跌误差（M=200, Smax=150）: N   隐式误差    CN误差
   N= 10  -0.06213   -0.02147
   N= 20  -0.03037   +0.00100
   N= 40  -0.01441   +0.00171
   N= 80  -0.00641   +0.00164

读输出：

• 例 21.10 的美式 4.07、欧式 3.91、显式法 4.26 与原书一致。控制变量用四位小数计算为 4.232，原书用两位小数的中间结果得 4.24。
• 显式法 \(t=0\) 那一列在 \(S=80\)、\(90\) 处出现 \(-0.13\)、\(-0.11\) 这样的负价格，正好落在中间概率为负的区域（\(j\ge13\)）。
• 隐式法在网格加密时单调收敛到约 4.284，与 21a 章二叉树的收敛值一致——三种方法殊途同归。
• \(\ln S\) 网格上的显式法中间概率恒为 2/3，上下概率都为正；\(N=50\) 时的 4.2599 与 21a 章三叉树 50 步的 4.2598 几乎相同，印证了二者的等价性（差别来自贴现因子 \(1/(1+r\Delta t)\) 与 \(e^{-r\Delta t}\)）。
• 在空间网格固定为 \(M=200\) 时，隐式法的误差随 \(N\) 加倍大约减半（一阶收敛）；Crank–Nicolson 在 \(N=20\) 时误差已降到 0.001 量级，之后被空间误差（约 0.0016）主导，不再随 \(N\) 下降。\(N=10\) 时 CN 误差较大，是收益在执行价处的折点引起的初始振荡。

本章小结

有限差分法把衍生品满足的偏微分方程离散到时间—价格网格上，用差商代替导数，从到期日开始逐层向后求解。隐式法在未知层上取空间导数，每层要解一个三对角方程组（追赶法，\(O(M)\)），无条件稳定，网格加密时总能收敛；显式法在已知层上取空间导数，可以逐点直接计算，但它等价于一棵三叉树，当中间"概率" \(1-\sigma^2j^2\Delta t\) 为负时会出现负价格和不收敛。把变量换成 \(Z=\ln S\) 后系数与 \(j\) 无关，取 \(\Delta Z=\sigma\sqrt{3\Delta t}\) 时显式法与 21a 章的三叉树完全相同且概率恒正。跳格法交替使用显式和隐式节点以避免解方程组，Crank–Nicolson 取二者平均，在时间方向上达到二阶精度，是实务中最常用的格式。综合原书第 21 章：树和有限差分向后计算，擅长美式期权，不擅长路径依赖和高维；蒙特卡洛向前计算，擅长路径依赖和高维，不擅长提前行权。

练习

基础

1. 用泰勒展开证明中心差分 (21.24) 的截断误差是 \(O(\Delta S^2)\)，而前向差分 (21.22) 是 \(O(\Delta S)\)。 提示：\(f(S\pm\Delta S)=f\pm f'\Delta S+\frac12f''\Delta S^2\pm\frac16f'''\Delta S^3+\cdots\)，相减后偶数阶项抵消。

2. 在例 21.11 的网格中（\(\sigma=0.4\)，\(M=20\)，\(T=5/12\)），要让所有内部点（\(j=1,\dots,19\)）的中间概率都非负，\(N\) 至少要取多少？ 提示：需要 \(\sigma^2\times19^2\times\Delta t\le1\)，即 \(\Delta t\le1/57.76\)，\(N\ge57.76\times5/12=24.07\)，取 \(N\ge25\)。

3. 用显式有限差分法为美式看跌期权定价（原书习题 21.18 的设定），说明你如何设置边界条件，以及在每一步如何处理提前行权。 提示：边界同 (21.28)–(21.30)；每算出一层，与内在价值取大。

4. 在显式有限差分法中，\(S=0\) 和 \(S=S_{max}\) 处的边界条件什么时候会影响结果？ 提示：显式法每步只向内传播一格，\(N\) 步后边界影响只能到达距离边界 \(N\) 格以内的点。若 \(S_0\) 距离边界超过 \(N\) 格，边界条件不影响 \(S_0\) 处的价格（原书习题 21.21）。

5. 说明隐式法为什么可以看作每个节点有 \(M+1\) 个分支的"多叉树"。 提示：解三对角方程组后，\(f_{i,j}\) 是 \(f_{i+1,\cdot}\) 所有 \(M+1\) 个值的线性组合（三对角矩阵的逆是满矩阵）。

进阶

6. 用隐式法为美式货币看涨期权定价时，(21.27)–(21.30) 需要做哪些修改？ 提示：\(q\) 换成 \(r_f\)；到期边界改为 \(\max(j\Delta S-K,0)\)；\(S=0\) 处价值为 0；\(S=S_{max}\) 处价值约为 \(S_{max}-K\)（美式看涨深度实值时立即行权）（原书习题 21.17）。

7. 为可转换债券设计一个有限差分方案：说明需要哪些边界条件（到期、转换、发行人赎回、持有人回售），每一层如何比较。 提示：每层取 \(\max(\min(\text{继续持有值},\ \text{赎回价}),\ \text{转换价值},\ \text{回售价})\) 的适当组合；到期为 \(\max(\text{面值},\ \text{转换价值})\)（原书习题 21.23）。

8. 修改"量化实战"中的 fd_theta 函数，加入 Rannacher 平滑：前两个时间步用隐式法（\(\theta=1\)），之后用 Crank–Nicolson。观察 \(N=10\) 时的误差是否改善。 提示：Rannacher 平滑通常能消除折点引起的振荡。

9. 用本章代码，分别用二叉树（21a 章 crr，需先运行第 21a 章量化实战代码或把其中的 crr 函数复制过来）、隐式有限差分（加密网格）、\(\ln S\) 显式法计算例 21.5（离散美元股息）的美式看跌价格。有限差分法如何处理离散美元股息？ 提示：在除息日处，令除息前的 \(f(S)=f_{\text{除息后}}(S-D)\)，需要在网格上插值；这是有限差分法处理离散事件的通用手法（跳跃条件）。

10. 构造一个两资产的问题（例如两只股票中较大者的欧式看涨期权），比较二维有限差分网格（每维 100 点、100 个时间步）与蒙特卡洛（\(10^5\) 条路径）的计算量。维度增加到 5 时呢？ 提示：二维网格每步 \(10^4\) 个点，五维 \(10^{10}\) 个点——"维度灾难"；蒙特卡洛的计算量只随维度线性增加。

原书推荐习题：21.18、21.28（显式有限差分实现，观察稳定性问题）；21.17（隐式法用于美式货币看涨时的修改）；21.21（显式法边界的影响范围）；21.23（可转换债券的边界条件与有限差分设计）。

原书对照

延伸阅读：Hull & White (1990) 显式有限差分法的改进；Wilmott (1998)《Derivatives: The Theory and Practice of Financial Engineering》；Clewlow & Strickland《Implementing Derivatives Models》；Press 等《Numerical Recipes》（三对角求解与 PDE 数值方法）。