实践中的 Heston 模型之有限差分

实践中的 Heston 模型之有限差分

Heston PDE

资产价格遵循以下随机过程：

\[\begin{align*} dS_\tau &= (\mu - q)S d\tau + \sqrt{v} S dB^1_\tau\\ dv_\tau &= \kappa(\theta - v)d\tau + \sigma \sqrt{v}dB^2_\tau\\ [dB^1_t, dB^2_t] &= \rho d\tau \end{align*}\]

通过构造无风险投资组合可以得到 Heston PDE：

\[\frac{\partial U}{\partial \tau} + S(r-q)\frac{\partial U}{\partial S} + \frac{1}{2}S^2v\frac{\partial^2 U}{\partial S^2} + \kappa(\theta - v)\frac{\partial U}{\partial v} + \frac{1}{2}\sigma^2v\frac{\partial^2 U}{\partial v^2} + \rho\sigma Sv\frac{\partial^2 U}{\partial S \partial v} - rU = 0\]

令 $t = T-\tau$，则原 PDE 变成：

\[\frac{\partial U}{\partial t} = S(r-q)\frac{\partial U}{\partial S} + \frac{1}{2}S^2v\frac{\partial^2 U}{\partial S^2} + [\kappa(\theta - v) + \sigma\Lambda]\frac{\partial U}{\partial v} + \frac{1}{2}\sigma^2v\frac{\partial^2 U}{\partial v^2} + \rho\sigma Sv\frac{\partial^2 U}{\partial S \partial v} - rU\]

Heston PDE 实质上是一个对流扩散反应方程（Convection-Diffusion-Reaction Equation），等号右边的微分算子 $L$ 可以分成三部分：

• 扩散项，$\frac{1}{2}S^2v\frac{\partial^2 U}{\partial S^2} + \frac{1}{2}\sigma^2v\frac{\partial^2 U}{\partial v^2} + \rho\sigma Sv\frac{\partial^2 U}{\partial S \partial v}$
• 对流项，$S(r-q)\frac{\partial U}{\partial S} + \kappa(\theta - v)\frac{\partial U}{\partial v}$
• 线性的反应项，$-rU$

Heston 模型下的衍生品定价问题将归结为求解上述 PDE 的初边值问题（Initial Boundary Value Problem，IBVP）。求解这个 IBVP 的一般步骤是先确定求解域的大小和网格划分，在求解域内用差分近似导数，进而将 PDE 转成“半离散”形式，也就是一个大型 ODE，最后求解这个 ODE 的初值问题（IVP）。因此，需要依次解决四个方面的问题：

• 将求解范围限定为一个充分大的矩形，并设定价格和方差维度上边界条件；
• 确定价格和方差维度上的网格划分；
• 用合适的差分格式近似微分算子 $L$ 中的导数；
• 选择时间维度上的演化方式求解半离散形式。

下面将分别讨论这几个问题。

初值与边界条件

IVP 的初值条件通常简洁明了，就是期权在到期日的内在价值函数。

与 BS PDE 相比，Heston PDE 多了一个变量——随机方差，并且方差是不可交易的对象，这让边界条件的设定变得复杂。价格维度 $S$ 上的边界条件通常根据期权的支付特点来设定 Dirichlet 条件或 Neumann 条件，例如：

• 欧式看涨期权：$S=0$ 时 $U=0$；$S = S_{max}$ 时 $\frac{\partial U}{\partial S}=1$；
• 欧式看跌期权：$S=0$ 时 $U=K$；$S = S_{max}$ 时 $U=0$；
• Up-Out 看涨期权：$S=0$ 时 $U=0$；$S = B$ 时 $U=rebate$

$v=v_{max}$ 处的边界条件

$v_{max}$ 处的边界条件设定见仁见智，既可能是启发式的又可能要考虑 PDE 在极限处的渐进性，常见的边界条件有：

• 针对欧式看涨期权，$v = v_{max}$ 时设定 $U=S$，用到了波动率趋于无穷时 BS-PDE 的解收敛于 $S$ 的结论；
• 一般的，$v = v_{max}$ 时设定 $\frac{\partial U}{\partial v}=0$，即认为期权价格关于方差的边际变化随着方差增长而趋近于 0。

另一种处理方式是在 $v = v_{max}$ 处不设定 Dirichlet 条件或 Neumann 条件，而是要求 Heston PDE 退化成 BS PDE（文献¹），即将 $v_{max}$ 代入到 Heston PDE 中并去掉关于 $v$ 的导数项，相当于将 $\frac{\partial U}{\partial v}=0$ “内嵌”到 Heston PDE 中得到新的 PDE。

$v=0$ 处的边界条件

$v=0$ 处边界条件的设定更加复杂，需要依赖于 Fichera 理论的指导（文献²、文献³）。

记求解域 $\Omega$ 上有微分算子 $L$：

\[LU = \sum_{i,j=1}^{n} a_{i,j}\frac{\partial^2 U}{\partial x_i \partial x_j} + \sum_{i=1}^{n} b_{i}\frac{\partial U}{\partial x_i} + cU\]

$\Sigma$ 是 $\Omega$ 的边界，并且矩阵 $\vert a_{i,j} \vert$ 在 $\Sigma \cup \Omega$ 上半正定的，

\[\sum_{i,j=1}^{n} a_{i,j} x_i x_j \ge 0,x \in \Sigma \cup \Omega\]

令

\[\Sigma_3 =\left \{\sum_{i,j=1}^n a_{i,j} x_i x_j > 0,x \in \Sigma \right \}\]

在 $\Sigma - \Sigma_3$ 上定义 Fichera 函数：

\[f(x) = \sum_{i=1}^n\left(b_i - \sum_{k=1}^n \frac{\partial a_{i,k}}{\partial x_k} \right)\vec{v}_i\]

其中 $\vec{v}$ 是 $\Sigma$ 上向内的单位法向量。

根据 $f$ 的值 $\Sigma - \Sigma_3$ 被分为三部分：

• $\Sigma_0$：$f(x) = 0$
• $\Sigma_1$：$f(x) > 0$
• $\Sigma_2$：$f(x) < 0$

Fichera 理论指出 $\Sigma_0$ 和 $\Sigma_1$ 上不要求边界条件，PDE 在 $\Sigma_0$ 和 $\Sigma_1$ 上遵循退化的形式；在 $\Sigma_2$ 上要求设置边界条件。

回到 Heston 问题，Fichera 函数在 $v = 0$ 边界上的值等于 $\frac{1}{2} \kappa\theta - \sigma^2$，恰好是 Feller 条件（数学真美妙！），也就是说是否需要设置边界条件取决于随机方差过程是否满足 Feller 条件，

• 当满足 Feller 条件时，无需设置边界条件，PDE 在边界上变成退化形式，即将 $v=0$ 代入到原 PDE 中形成新的 PDE（文献⁴）；
• 当不满足 Feller 条件时，需要设置边界条件，通常设置为期权的支付函数（文献⁵）。

Fichera 理论对涉及到 CIR 利率模型的 PDE 同样适用（文献²）。

价格和方差维度上的网格划分

求解域的上下限设定是启发式的，例如对于普通的看涨和看跌期权来说，价格和方差的下限通常自然地设定为 0（方差过程不满足 Feller 条件时，方差有正的概率达到零），上限通常设定为充分大的数。Heston 模型下价格近似于对数正态分布，因此可以用近似分布的某个分位数作为上限，以保证价格能达到该上限的可能性充分小。方差的分布是非中心的卡方分布，同理可以用某个分位数设定方差的上限。

为减少计算量的同时保有足够的精度，价格和方差维度上的网格通常是非均匀的，价格维度上在当前价格和执行价格附近稠密，在边界处稀疏；方差维度上在 0 附近稠密，在上限处稀疏。非均匀网格的构造通常采用文献⁶的方法论。

价格和方差维度上的差分

二元函数 $f(x,y)$ 导数的差分近似格式存在多种形式，记 $\Delta x_i = x_i - x_{i-1}$，一阶导数 $\frac{\partial f}{\partial x}$ 可以用下面几种差分格式近似：

一阶精度格式，分别是向后和向前差分：

\[\begin{aligned} f^{\prime}\left(x_i\right) & \approx \alpha_{i,-1} f\left(x_{i-1}\right)+\alpha_{i, 0} f\left(x_i\right) \\ f^{\prime}\left(x_i\right) & \approx \beta_{i, 0} f\left(x_i\right)+\beta_{i, 1} f\left(x_{i+1}\right) \\ \end{aligned}\] \[\begin{gathered} \alpha_{i,-1}=\frac{-1}{\Delta x_i}, \alpha_{i, 0}=\frac{1}{\Delta x_i} \\ \beta_{i, 0}=\frac{-1}{\Delta x_{i+1}}, \beta_{i, 1}=\frac{1}{\Delta x_{i+1}} \end{gathered}\]

二阶精度格式，分别是向后、中心和向前差分

\[\begin{aligned} f^{\prime}\left(x_i\right) & \approx \alpha_{i,-2} f\left(x_{i-2}\right)+\alpha_{i,-1} f\left(x_{i-1}\right)+\alpha_{i, 0} f\left(x_i\right) \\ f^{\prime}\left(x_i\right) & \approx \beta_{i,-1} f\left(x_{i-1}\right)+\beta_{i, 0} f\left(x_i\right)+\beta_{i, 1} f\left(x_{i+1}\right) \\ f^{\prime}\left(x_i\right) & \approx \gamma_{i, 0} f\left(x_i\right)+\gamma_{i, 1} f\left(x_{i+1}\right)+\gamma_{i, 2} f\left(x_{i+2}\right) \end{aligned}\] \[\begin{gathered} \alpha_{i,-2}=\frac{\Delta x_i}{\Delta x_{i-1}\left(\Delta x_{i-1}+\Delta x_i\right)}, \alpha_{i,-1}=\frac{-\Delta x_{i-1}-\Delta x_i}{\Delta x_{i-1} \Delta x_i}, \alpha_{i, 0}=\frac{\Delta x_{i-1}+2 \Delta x_i}{\Delta x_i\left(\Delta x_{i-1}+\Delta x_i\right)} \\ \beta_{i,-1}=\frac{-\Delta x_{i+1}}{\Delta x_i\left(\Delta x_i+\Delta x_{i+1}\right)}, \quad \beta_{i, 0}=\frac{\Delta x_{i+1}-\Delta x_i}{\Delta x_i \Delta x_{i+1}}, \quad \beta_{i, 1}=\frac{\Delta x_i}{\Delta x_{i+1}\left(\Delta x_i+\Delta x_{i+1}\right)} \\ \gamma_{i, 0}=\frac{-2 \Delta x_{i+1}-\Delta x_{i+2}}{\Delta x_{i+1}\left(\Delta x_{i+1}+\Delta x_{i+2}\right)}, \quad \gamma_{i, 1}=\frac{\Delta x_{i+1}+\Delta x_{i+2}}{\Delta x_{i+1} \Delta x_{i+2}}, \quad \gamma_{i, 2}=\frac{-\Delta x_{i+1}}{\Delta x_{i+2}\left(\Delta x_{i+1}+\Delta x_{i+2}\right)} \end{gathered}\]

二阶导数 $\frac{\partial^2 f}{\partial x^2}$ 就简单多了，通常只采用中心差分格式：

\[f^{\prime\prime}\left(x_i\right) \approx \delta_{i,-1} f\left(x_{i-1}\right)+\delta_{i, 0} f\left(x_i\right)+\delta_{i, 1} f\left(x_{i+1}\right)\] \[\delta_{i,-1}=\frac{2}{\Delta x_i\left(\Delta x_i+\Delta x_{i+1}\right)}, \quad \delta_{i, 0}=\frac{-2}{\Delta x_i \Delta x_{i+1}}, \quad \delta_{i, 1}=\frac{2}{\Delta x_{i+1}\left(\Delta x_i+\Delta x_{i+1}\right)}\]

混合导数 $\frac{\partial^2 f}{\partial x \partial y}$ 通常用九点差分格式：

\[\frac{\partial^2f}{\partial x\partial y}(x_i,y_j) \approx \sum_{k,l=-1}^{1}\beta_{i,k}\hat{\beta}_{j,l}f(x_{i+k}, y_{j+l})\] \[\beta_{i,-1}=\frac{-\Delta x_{i+1}}{\Delta x_i\left(\Delta x_i+\Delta x_{i+1}\right)}, \quad \beta_{i, 0}=\frac{\Delta x_{i+1}-\Delta x_i}{\Delta x_i \Delta x_{i+1}}, \quad \beta_{i, 1}=\frac{\Delta x_i}{\Delta x_{i+1}\left(\Delta x_i+\Delta x_{i+1}\right)}\]

除了常见的九点格式，混合导数的近似还存在两个特殊的七点差分格式（文献⁷、文献⁸）。点 $(x_i,y_j)$ 处的混合导数可以由 $(x_{i-1},y_j)$、$(x_i,y_j)$、$(x_{i+1},y_j)$、$(x_i,y_{j-1})$、$(x_i,y_{j+1})$、$(x_{i+1},y_{j+1})$、$(x_{i-1},y_{j-1})$ 七个点的线性组合近似，这是因为

\[\begin{align*} \frac{\partial^2 u}{\partial x \partial y} &\approx \frac{1}{h_r h_u}\left(u\left(x+h_r, y+h_u\right)-u-h_r \frac{\partial u}{\partial x}-h_u \frac{\partial u}{\partial y}-\frac{1}{2} h_r^2 \frac{\partial^2 u}{\partial x^2}-\frac{1}{2} h_u^2 \frac{\partial^2 u}{\partial y^2}\right)\\ \frac{\partial^2 u}{\partial x \partial y} &\approx \frac{1}{h_l h_d}\left(u\left(x-h_l, y-h_d\right)-u+h_l \frac{\partial u}{\partial x}+h_d \frac{\partial u}{\partial y}-\frac{1}{2} h_l^2 \frac{\partial^2 u}{\partial x^2}-\frac{1}{2} h_d^2 \frac{\partial^2 u}{\partial y^2}\right) \end{align*}\] \[\begin{gathered} h_l = \Delta x_i, h_r = \Delta x_{i+1}\\ h_u = \Delta y_j, h_d = \Delta y_{j+1} \end{gathered}\]

两者加权（通常是等权）就可以得到混合导数的差分近似。同理，点 $(x_i,y_j)$ 处的混合导数也可以由 $(x_{i-1},y_j)$、$(x_i,y_j)$、$(x_{i+1},y_j)$、$(x_i,y_{j-1})$、$(x_i,y_{j+1})$、$(x_{i-1},y_{j+1})$、$(x_{i+1},y_{j-1})$ 七个点的线性组合近似。后面会提到基于七点格式可以创造出一种高效的数值算法。现实中选用哪种七点格式依赖于 $x$ 和 $y$ 维度的物理意义和关系。

PDE 中的扩散项通常采用中心差分格式，而对流项的差分格式选择比较复杂（文献⁹）：

• 当对流项系数相对于扩散项系数很小时可以采用中心差分格式；
• 当对流项系数相对于扩散项系数较大时（即对流项占优）中心差分可能会造成解出现振荡（文献¹⁰），此时对流项需要采用“迎风格式”，也就是说当对流项系数大于零时，采用向前差分；当对流项系数小于零时，采用向后差分。

时间维度上的演化方式

在 $[0, S_{max}] \times [0, v_{max}]$ 划分成的网格上，一旦确定好差分格式和边界条件，原始 PDE 的初边值问题（IBVP）就可以用一个大型的 ODE 的初值问题（IVP）近似，这就是 PDE 的半离散形式，写作

\[\frac{dU}{dt} = AU + R(t)\]

其中，$U_t$ 是状态价格的近似向量（按照某种索引规则将原先的状态矩阵拉平为向量），矩阵 $A$ 是 PDE 右侧微分算子 $L$ 的差分近似，$R(t)$ 是边界条件对应的修正项。

基础求解方法

求解 ODE 初值问题的经典方法分别是隐式 Euler 法（IE）和 Crank-Nicolson 法（CN），两种方法都是无条件稳定的，不要求时间差分间隔满足特定条件就能保证数值解的收敛。

隐式 Euler 法：

\[\frac{1}{\Delta t}(U_{t + \Delta t} - U_t) = A U_{t + \Delta t} + R_{t + \Delta t}\]

Crank-Nicolson 法：

\[\frac{1}{\Delta t}(U_{t + \Delta t} - U_t) = A(\frac{1}{2} U_{t + \Delta t} + \frac{1}{2}U_{t}) + \frac{1}{2}(R_t + R_{t + \Delta t})\]

以 IE 为例，要求解出 $U_{t + \Delta t}$ 不可避免地需要求解矩阵 $I - \Delta t A$ 确定的线性方程组。与一维的 BS PDE 的情况不同，此时的矩阵 $A$ 并不是一个高度结构化的矩阵。全部采用中心差分的话，BS PDE 对应的矩阵 $A$ 将是一个三对角矩阵，因此 $I - \Delta t A$ 也是三对角的，可以用追赶法（也称 Thomas 算法）这种高效的方法求解对应的线性方程组。同样采用中心差分，假定价格维度划分为 100 份，方差维度划分为 60 份，按照先价格后方差的字典序拉平状态价格矩阵，那么 Heston PDE 对应的矩阵 $A$ 将是一个 6000 x 6000 的矩阵，尽管 $A$ 将是高度稀疏的，但 $A$ 的带宽将达到 2 x 100 + 1（如果对流项用二阶的迎风格式，带宽还会翻倍），这种情况下通常用迭代法求解矩阵 $I - \Delta t A$ 确定的线性方程组，显然计算成本过高。

回想 BS PDE 的情况，当求解一维对流扩散反应方程时需要求解一个带状矩阵（三对角）决定的线性方程。如果求解二维对流扩散反应方程的过程能够转换成为求解两个一维对流扩散反应方程，就能将求解一个复杂方程组的问题变为求解两个简单方程组的问题，从而为减少计算成本创造可能性。

事实上，上述想法是可行的，并且已经独立地发展出了两套理论体系，分别是交替方向隐式方法（ADI）和算子分裂方法（OS）。特别是 ADI，已经广泛地应用到了 Heson 模型及其扩展的数值求解实践中，OS 则在物理学领域应用广泛。

ADI 和 OS 两套理论几乎同时诞生于冷战时代，分属苏联和美国大阵营。

高级求解方法

ADI 概述

由于矩阵 $A$ 是微分算子的差分近似，因此可以将矩阵 $A$ 按照维度分解为三部分，分别是

• $A_S$，对应 $S(r-q)\frac{\partial U}{\partial S} + \frac{1}{2}S^2v\frac{\partial^2 U}{\partial S^2} - \frac{r}{2}U$；
• $A_v$，对应 $\kappa(\theta - v) \frac{\partial U}{\partial v} + \frac{1}{2}\sigma^2v\frac{\partial^2 U}{\partial v^2} - \frac{r}{2}U$；
• $A_{Sv}$，对应 $\rho\sigma Sv\frac{\partial^2 U}{\partial S \partial v}$

同理，$R(t)$ 也做相似的分解：$R_S(t)$、$R_v(t)$ 和 $R_{Sv}(t)$

将传统的隐式 Euler 法拆解成两步，用隐式格式交替求解 $A_S$ 和 $A_v$ 决定的子问题，

\[\begin{align*} \frac{1}{\Delta t}(\bar{U}_t - U_t) &= A_S \bar{U}_t + (A_v + A_{Sv})U_t + R_S(t)\\ \frac{1}{\Delta t}(U_{t + \Delta t} - \bar{U}_t) &= A_v U_{t + \Delta t} + (A_v + A_{Sv})\bar{U}_t + R_v(t) \end{align*}\]

在第一步中采用先价格后方差的字典序拉平状态价格矩阵，保证 $A_S$ 是带状矩阵。在第二步中需要变换索引，采用先方差后价格的字典序，这样能保证 $A_v$ 也是带状矩阵。最终，求解一个复杂方程组变为求解两个结构化方程组。

上述两步便是 ADI 的雏形，下面介绍若干 ADI 的变种和一般理论。

若干 ADI 格式

对于 k 维的 PDE，矩阵 $A$ 和向量 $R$ 存在分解：

• $A = A_0 + A_1 + A_2 + \cdots + A_k$，$A_i$ 对应第 i 个空间维度，$A_0$ 对应混合倒数；
• $R = R_0 + R_1 + R_2 + \cdots + R_k$，$R_i$ 对应第 i 个空间维度，$R_0$ 对应混合倒数。

令 $F(t, U) = AU_t + R_t$、 $F_i(t, U) = A_i U(t) + R_i(t)$，将上述两步法推广并改写成递归形式就得到历史上第一个 ADI 格式——Douglas 格式（Do）的一般形式：

\[\left\{\begin{array}{l} Y_0=U_{n-1}+\Delta t F\left(t_{n-1}, U_{n-1}\right), \\ Y_j=Y_{j-1}+\theta \Delta t\left(F_j\left(t_n, Y_j\right)-F_j\left(t_{n-1}, U_{n-1}\right)\right), \quad j=1,2, \ldots, k, \\ U_n=Y_k \end{array}\right.\]

$\theta = 1$ 时就是刚才的情况，若 $\theta = \frac{1}{2}$ 则表示用了 CN 法拆解。

观察 Douglas 格式的一般形式可以发现，第一步是用简单的显式 Euler 格式得到初步结果，随后再递归修正各个维度上的误差。Douglas 格式的推广包括 Craig-Sneyd 格式（CS）、修正 Craig-Sneyd 格式（MCS）和 Hundsdorfer-Verwer 格式（HV）也遵循相同的计算流程，只是递归修正的部分更加复杂。

Craig-Sneyd 格式的形式如下：

\[\left\{\begin{array}{l} Y_0=U_{n-1}+\Delta t F\left(t_{n-1}, U_{n-1}\right), \\ Y_j=Y_{j-1}+\theta \Delta t\left(F_j\left(t_n, Y_j\right)-F_j\left(t_{n-1}, U_{n-1}\right)\right), \quad j=1,2, \ldots, k, \\ \widehat{Y}_0=Y_0+\frac{1}{2} \Delta t\left(F_0\left(t_n, Y_k\right)-F_0\left(t_{n-1}, U_{n-1}\right)\right), \\ \widehat{Y}_j=\widehat{Y}_{j-1}+\theta \Delta t\left(F_j\left(t_n, \widehat{Y}_j\right)-F_j\left(t_{n-1}, U_{n-1}\right)\right), \quad j=1,2, \ldots, k, \\ U_n=\widehat{Y}_k \end{array}\right.\]

修正 Craig-Sneyd 格式的形式如下：

\[\left\{\begin{array}{l} Y_0=U_{n-1}+\Delta t F\left(t_{n-1}, U_{n-1}\right), \\ Y_j=Y_{j-1}+\theta \Delta t\left(F_j\left(t_n, Y_j\right)-F_j\left(t_{n-1}, U_{n-1}\right)\right), \quad j=1,2, \ldots, k, \\ \widehat{Y}_0=Y_0+\theta \Delta t\left(F_0\left(t_n, Y_k\right)-F_0\left(t_{n-1}, U_{n-1}\right)\right), \\ \widetilde{Y}_0=\widehat{Y}_0+(\frac{1}{2} - \theta) \Delta t\left(F\left(t_n, Y_k\right)-F\left(t_{n-1}, U_{n-1}\right)\right), \\ \widetilde{Y}_j=\widetilde{Y}_{j-1}+\theta \Delta t\left(F_j\left(t_n, \widetilde{Y}_j\right)-F_j\left(t_{n-1}, U_{n-1}\right)\right), \quad j=1,2, \ldots, k, \\ U_n=\widetilde{Y}_k \end{array}\right.\]

Hundsdorfer-Verwer 格式的形式如下：

\[\left\{\begin{array}{l} Y_0=U_{n-1}+\Delta t F\left(t_{n-1}, U_{n-1}\right), \\ Y_j=Y_{j-1}+\theta \Delta t\left(F_j\left(t_n, Y_j\right)-F_j\left(t_{n-1}, U_{n-1}\right)\right), \quad j=1,2, \ldots, k, \\ \widetilde{Y}_0=Y_0+\frac{1}{2} \Delta t\left(F\left(t_n, Y_k\right)-F\left(t_{n-1}, U_{n-1}\right)\right), \\ \widetilde{Y}_j=\widetilde{Y}_{j-1}+\theta \Delta t\left(F_j\left(t_n, \widetilde{Y}_j\right)-F_j\left(t_{n-1}, Y_k\right)\right), \quad j=1,2, \ldots, k, \\ U_n=\widetilde{Y}_k \end{array}\right.\]

ADI 的稳定性

假设 $k$ 维对流扩散方程中空间变量的微分算子 L 如下：

\[L U = \sum_{i=1}^k d_{i,i} \frac{\partial^2 U}{\partial x_i^2} + \sum_{i=1}^{k}\sum_{j=1}^k d_{i,j} \frac{\partial^2 U}{\partial x_i \partial x_j} + \sum_{i=1}^k c_i \frac{\partial U}{\partial x_i}, d_{i,j}=d_{j,i}\]

并且存在 $\gamma$ 使得 $\vert d_{i,j}\vert \le \gamma \sqrt{d_{i,i}, d_{j,j}}$。

在金融领域，$\gamma$ 就是相关性绝对值的最大值。

通常上述 ADI 格式的无条件稳定性由 $\theta$ 的大小控制，同时 $\theta$ 的选择也受到微分算子参数的影响（文献¹¹、文献¹²、文献¹³）。

Douglas 格式

对于扩散问题，Do 无条件稳定的充分条件是二维情况下，$\theta \ge 1/2$；三维情况下，$\theta \ge \max(\frac{1}{2}, \frac{2(2\gamma+1)}{9})$。Do 无条件稳定的必要条件是 $\theta$ 要满足

\[\theta \ge \max\left(\frac{1}{2}, \frac{1}{2} d_k [(k-1)\gamma+1]\right)\] \[d_k = \left(1-\frac{1}{k}\right)^{k-1}\]

对于对流扩散问题，二维情况下，$\theta \ge 1/2$ 时 Do 格式是无条件稳定的。

Craig-Sneyd 格式

对于扩散问题，CS 无条件稳定的充分条件是在二维、三维的情况下，$\theta \ge 1/2$。CS 无条件稳定的必要条件是 $\theta$ 要满足

\[\theta \ge \max \left(\frac{1}{2}, \frac{1}{2} c_k k \gamma\right)\] \[c_k = \left(1-\frac{1}{k}\right)^{k}\]

对于对流扩散问题，二维情况下，$\theta \ge 1/2$ 时 CS 格式是无条件稳定的。

修正 Craig-Sneyd 格式

对于扩散问题，MCS 无条件稳定的充分条件是二维情况下，$\theta \ge \max(\frac{1}{4}, \frac{\gamma+1}{6})$；三维情况下，$\theta \ge \max(\frac{1}{4}, \frac{2(2\gamma+1)}{13})$。MCS 无条件稳定的必要条件是 $\theta$ 要满足 \(\theta \ge \max \left(\frac{1}{4}, \frac{1}{2} b_k[(k-1)\gamma+1]\right)\\ b_k = \frac{1}{1+\left(1-\frac{1}{k}\right)^{k-1}}\)

对于对流扩散问题，二维情况下，$1 \ge \theta \ge 1/2$ 时 MCS 格式是无条件稳定的。

此外，MCS 格式的无条件稳定性存在更细致的结果（文献¹⁴，文献¹⁵ 做了充分的数值实验），令 \(\begin{aligned} & p_1(\theta)=\theta \sqrt{\frac{5-2 \theta}{1-2 \theta}}-\theta, \\ & q_1(\theta)=\frac{\theta}{2}\left(\frac{2 \theta+\sqrt{4 \theta^2-16 \theta+6}}{1-2 \theta}\right), \\ & p_2(\theta)=\frac{4 \theta^2}{\sqrt{4 \theta^4-24 \theta^3+24 \theta^2-8 \theta+1}-2 \theta^2+4 \theta-1} \\ & q_2(\theta)=\frac{-4 \theta^3+11 \theta^2-6 \theta+1+\sqrt{16 \theta^6-128 \theta^5+216 \theta^4-160 \theta^3+61 \theta^2-12 \theta+1}}{5 \theta^2-4 \theta+1} . \end{aligned}\)

• 当 $\theta = 1$ 时，若 $\gamma < \frac{2+\sqrt{10}}{6}$，则 MCS 无条件稳定。
• 当 $\frac{1}{4} \le \theta \le \frac{1}{3}$ 时，若 $\gamma < \min(p_1(\theta), q_1(\theta))$，则 MCS 无条件稳定
• 当 $\frac{1}{3} \le \theta \le \frac{1}{3}$ 时，若 $\gamma < \min(p_2(\theta), q_2(\theta))$，则 MCS 无条件稳定
• 若 $\gamma \in [0,\frac{1}{2}]$，当且仅当 $\theta \ge \frac{1}{4}$ 时 MCS 无条件稳定

Hundsdorfer-Verwer 格式

\[\left\{\begin{array}{l} Y_0=U_{n-1}+\Delta t F\left(t_{n-1}, U_{n-1}\right), \\ Y_j=Y_{j-1}+\theta \Delta t\left(F_j\left(t_n, Y_j\right)-F_j\left(t_{n-1}, U_{n-1}\right)\right), \quad j=1,2, \ldots, k, \\ \widetilde{Y}_0=Y_0+\frac{1}{2} \Delta t\left(F\left(t_n, Y_k\right)-F\left(t_{n-1}, U_{n-1}\right)\right), \\ \widetilde{Y}_j=\widetilde{Y}_{j-1}+\theta \Delta t\left(F_j\left(t_n, \widetilde{Y}_j\right)-F_j\left(t_{n-1}, Y_k\right)\right), \quad j=1,2, \ldots, k, \\ U_n=\widetilde{Y}_k \end{array}\right.\]

对于扩散问题，HV 无条件稳定的充分条件是二维情况下，$\theta \ge \max(\frac{1}{4}, \frac{\gamma+1}{4+2\sqrt{2}})$；三维情况下，$\theta \ge \max(\frac{1}{4}, \frac{2（2\gamma+1)}{4+2\sqrt{3}})$。HV 无条件稳定的必要条件是 $\theta$ 要满足

\[\theta \ge \max \left(\frac{1}{4}, \frac{1}{2} a_k[(k-1)\gamma+1]\right)\]

其中 $a_k \in (0, 1/2)$ ，并且是

\[2a\left(1+\frac{1-a}{k-1}\right)^{k-1} - 1 = 0\]

的解。

对于对流扩散问题，二维情况下，当 $\theta \ge \frac{1}{2} + \frac{1}{6}\sqrt{3}$ 时 HV 是无条件稳定的。

一般来说 $\theta$ 越小误差越小收敛速度越快，要尽量选择小的 $\theta$（文献¹⁶）。

OS 概述

OS 的其他叫法分别是分步法（Fractional Steps）和局部一维化法（LOD），由于是前苏联学者发展起来的的，也被称为苏联分裂法（真是很不吉利！）。其思想与 ADI 类似，即将多维 PDE 问题的求解转变为求解若干一维问题。与 ADI 相比，OS 的降维过程更为彻底，因为 OS 按照维度分解了原始的微分算子，直接从微分算子的层面降维。

目前 OS 在金融领域的应用不及 ADI 广泛，对于 Heston PDE 来说，最常用到的是 Yanenko 格式（文献¹⁷）。

Yanenko 格式

将 $L$ 分解为两部分，$L = L_S + L_v$，

\[\begin{align*} L_S U &= S(r-q)\frac{\partial U}{\partial S} + \frac{1}{2}S^2v\frac{\partial^2 U}{\partial S^2} + \frac{1}{2}\rho\sigma Sv\frac{\partial^2 U}{\partial S \partial v} - \frac{1}{2}rU\\ L_v U &= \frac{1}{2}\sigma^2v\frac{\partial^2 U}{\partial v^2} + \kappa(\theta - v) \frac{\partial U}{\partial v} + \frac{1}{2}\rho\sigma Sv\frac{\partial^2 U}{\partial S \partial v} - \frac{1}{2}rU \end{align*}\]

由于 $L_S + L_v \approx (I + L_S)(I + L_v) - I$，因此

\[\frac{\partial U}{\partial t} + U = (I + L_S)(I + L_v)U\]

原先二维的 PDE 问题现在被分解为两个 PDE 子问题的联合，

\[\begin{align*} \frac{\partial U}{\partial t} + U &= (I + L_S)U\\ \frac{\partial U}{\partial t} + U &= (I + L_v)U \end{align*}\]

依次求解两个子问题就得到原问题的解。沿用前面的记号，并用 $\theta$ 法求解两个子问题中非混合导数部分，混合导数的部分采用显示 Euler 格式，

\[\begin{align*} \frac{1}{\Delta t}(\bar{U}_t - U_t) &= A_S [\theta \bar{U}_t + (1-\theta)U_t] + \frac{1}{2} A_{Sv}U_t + R_S(t)\\ \frac{1}{\Delta t}(U_{t + \Delta t} - \bar{U}_t) &= A_v [\theta U_{t + \Delta t} + (1-\theta)\bar{U}_t] + \frac{1}{2}A_{Sv}\bar{U}_t + R_v(t) \end{align*}\]

和 Do 一样，在第一步中采用先价格后方差的字典序拉平状态价格矩阵，保证 $A_S$ 是带状矩阵。在第二步中采用先方差后价格的字典序，保证 $A_v$ 是带状矩阵。最终实现了类似 Do 的效果，该格式是一阶精度的，并且 $\frac{1}{2} \le \theta \le 1$ 的时是无条件稳定的。

OS 在物理领域的应用更加广泛，特别是不存在混合导数的情况，这在物理学中比较常见。

Lie-Trotter 格式和 Strang 格式

当不存在混合导数时，上述 Yanenko 格式也被称为 Lie-Trotter 格式（文献¹⁸、¹⁹、²⁰），该格式同样是一阶精度的。此外，针对不存在混合导数的情况还有一种 Strang 格式（文献¹⁸、¹⁹、²⁰），该格式是二阶精度的， \(\begin{align*} \frac{1}{2\Delta t}(\bar{U}_t - U_t) &= A_S [\theta \bar{U}_t + (1-\theta)U_t] + R_S(t)\\ \frac{1}{\Delta t}(\hat{U}_t - \bar{U}_t) &= A_v [\theta \hat{U}_t + (1-\theta)\bar{U}_t] + R_v(t)\\ \frac{1}{2\Delta t}(U_{t + \Delta t} - \hat{U}_t) &= A_S [\theta U_{t + \Delta t} + (1-\theta)\hat{U}_t] + R_S(t) \end{align*}\)

Strang 格式很像跳舞，$S$ 方向上跳半步，$v$ 方向上跳一步，最后在 $S$ 方向上再跳半步。Yanenko 和 Lie-Trotter 格式的舞步更简单一些，分别跳两步就好了。

分量形式分裂格式

OS 的思想和七点格式结合能产生更有趣的分量形式（Componentwise）分裂格式。现实中 $S$ 和 $v$ 通常有较强的负相关性（文献¹⁴、⁸），因此 $(S_i,v_j)$ 点上的混合导数用下列七个点的线性组合即可足够精确的表示：

• $(S_{i-1},v_j)$、
• $(S_i,v_j)$、$(S_{i+1},v_j)$
• $(S_i,v_{j-1})$、$(S_i,v_{j+1})$
• $(S_{i-1},v_{j+1})$、$(S_{i+1},v_{j-1})$

从概率角度直观理解七点格式：由于较强的负相关性，$(S_i,v_j)$ 点处的转移概率密度函数在上述七个点出的取值较 $(S_{i+1},v_{j+1})$ 和 $(S_{i-1},v_{j-1})$ 两个点要高得多，因为 $S$ 和 $v$ 同方向变动的可能性非常小。因此这七个点上的信息已经足够多，可以忽略 $(S_{i+1},v_{j+1})$ 和 $(S_{i-1},v_{j-1})$ 这两个点，采用七点格式而不是九点格式。

若混合导数采用七点格式，矩阵 $A$ 依然可以被分解为 3 部分，分别是

• $A_S’$，对应 $S(r-q)\frac{\partial U}{\partial S} + \frac{1}{2}S^2v\frac{\partial^2 U}{\partial S^2} - \frac{r}{3}U$，以及混合导数七点格式中引入的 $\frac{\partial U}{\partial S}$ 和 $\frac{\partial^2 U}{\partial S^2}$；
• $A_v’$，对应 $\kappa(\theta - v)\frac{\partial U}{\partial v} + \frac{1}{2}\sigma^2v\frac{\partial^2 U}{\partial v^2} - \frac{r}{3}U$，以及混合导数七点格式中引入的 $\frac{\partial U}{\partial v}$ 和 $\frac{\partial^2 U}{\partial v^2}$；
• $A_{Sv}’$，对应 $\rho\sigma Sv\frac{\partial^2 U}{\partial S \partial v} - \frac{r}{3}U$，但要去掉七点格式中引入的 $\frac{\partial U}{\partial S}$、$\frac{\partial^2 U}{\partial S^2}$、$\frac{\partial U}{\partial v}$ 和 $\frac{\partial^2 U}{\partial v^2}$

同理，$R(t)$ 也做相似的分解：$R_S’(t)$、$R_v’(t)$ 和 $R_{Sv}’(t)$。

有了上述分解，就可以直接用 Lie-Trotter 格式和 Strang 格式求解

\[\frac{dU}{dt} = AU + R(t)\]

Lie-Trotter 格式 \(\begin{align*} \frac{1}{\Delta t}(\bar{U}_t - U_t) &= A'_S [\theta \bar{U}_t + (1-\theta)U_t] + R'_S(t)\\ \frac{1}{\Delta t}(\hat{U}_t - \bar{U}_t) &= A'_{Sv} [\theta \hat{U}_t + (1-\theta)\bar{U}_t] + R'_{Sv}(t)\\ \frac{1}{\Delta t}(U_{t + \Delta t} - \hat{U}_t) &= A_v [\theta U_{t + \Delta t} + (1-\theta)\hat{U}_t] + R'_v(t) \end{align*}\)

Strang 格式 \(\begin{align*} \frac{1}{2\Delta t}(U^1_t - U_t) &= A'_S [\theta U^1_t + (1-\theta)U_t] + R'_S(t)\\ \frac{1}{2\Delta t}(U^2_t - U^1_t) &= A'_v [\theta U^2_t + (1-\theta)U^1_t] + R'_v(t)\\ \frac{1}{\Delta t}(U^3_t - U^2_t) &= A'_{Sv} [\theta U^3_t + (1-\theta)U^2_t] + R'_{Sv}(t)\\ \frac{1}{2\Delta t}(U^4_t - U^3_t) &= A'_v [\theta U^4_t + (1-\theta)U^3_t] + R'_v(t)\\ \frac{1}{2\Delta t}(U_{t+\Delta t} - U^4_t) &= A'_S [\theta U_{t+\Delta t} + (1-\theta)U^4_t] + R'_S(t) \end{align*}\)

和 Yanenko 格式一样，要保证 $A_S’$ 是带状矩阵，需要采用先价格后方差的字典序拉平状态价格矩阵。要保证 $A_v’$ 是带状矩阵，需要采用先方差后价格的字典序。要保证 $A_{Sv}’$ 是带状矩阵，需要按照对角线排序拉平状态价格矩阵，即按照下面的顺序：

1. $(S_0, v_0)$
2. $(S_1, v_0)$、$(S_0, v_1)$
3. $(S_2, v_0)$、$(S_1, v_1)$、$(S_0, v_2)$
4. $(S_3, v_0)$、$(S_2, v_1)$、$(S_1, v_2)$、$(S_0, v_3)$
5. …

其他话题

求解方程组

无论是 ADI 还是 OS，中间步骤都要求解带状矩阵决定的方程组，

• 如果迎风格式采用一阶精度的差分格式，矩阵是三对角（Tridiagonal）矩阵，带宽是 3；
• 如果迎风格式采用二阶精度的差分格式，矩阵是五对角（Pentadiagonal）矩阵，带宽是 5。

前面提到过，三对角方程组可以用 Thomas 算法，五对角方程组的解法可以参考文献⁹、²¹。

OS 与美式期权

OS 的含义很广泛，可以像上述 Yanenko、Lie-Thotter 和 Strang 格式一样按维度将算子分解，而在有些物理问题中则按照扩散项和对流项将算子分解。如果将美式期权定价问题看作是一种 Linear Complementary Problem（LCP），文献²²、²³应用 OS 的思想，将原问题的算子分解为两部分：PDE 演化部分和 Linear Complementary 部分分别求解。此外，PDE 演化部分和前述一样，也可以用 ADI 和 OS 格式求解（文献²⁴、²⁵）。

去除混合导数项

有文献指出混合导数项的存在会降低收敛速。此外，若没有混合导数还使得 Strang 格式能够使用。文献²⁶通过一个线性变换去除了混合导数，但代价是边界条件开始变得复杂。

不平滑的初始状态

应对不平滑的初始状态通常用 Luskin-Rannacher 阻尼技术（文献²⁷），即初始的两步采用隐式 Euler 格式。此外，也可以用 $\theta = 1$ 的 Do 格式实现 Luskin-Rannacher 阻尼（文献¹⁶）。

实证案例与模型扩展

Heston 模型最常见的扩展有两种，一是结合跳变成 Bates 模型，二是结合随机利率，这两方面的实践经验可以参考文献¹⁶、²⁸、²⁹、³⁰、³¹。

其他形式的随机波动率模型，例如和 Heston 非常相似的 Beta 波动率模型（文献³²）也可以采用 ADI 或 OS 求解（文献³³）。

比较 ADI 和 OS 的文献不多，目前只有文献³⁴、³⁵。

参考文献

1. Pricing Equity Derivatives under Stochastic Volatility A Partial Differential Equation Approach ↩︎

2. Fichera Theory and its Application in Finance ↩︎ ↩︎²

3. Unconditionally Stable and Second-Order Accurate Explicit Finite Difference Schemes Using Domain Transformation ↩︎

4. Artificial Boundary Method for the Solution of Pricing European Options Under the Heston Model ↩︎

5. American Option of Stochastic Volatility Model with Negative Fichera Function on Degenerate Boundary ↩︎

6. Pricing Financial Instruments: The Finite Difference Method ↩︎

7. Componentwise Splitting Methods For Pricing American Options Under Stochastic Volatility ↩︎

8. A Componentwise Splitting Method for Pricing American Options Under the Bates Model ↩︎ ↩︎²

9. A Parallel Cyclic Reduction Algorithm for Pentadiagonal Systems with Application to a Convection-Dominated Heston PDE ↩︎ ↩︎²

10. Finite Difference Methods in Derivatives Pricing under Stochastic Volatility Models ↩︎

11. Stability of ADI Schemes for Multidimensional Diffusion Equations with Mixed Derivative Terms ↩︎

12. Stability of ADI Schemes Applied to Convection–Diffusion Equations with Mixed Derivative Terms ↩︎

13. Unconditional Stability of Second-Order ADI Schemes Applied to Multi-Dimensional Diffusion Equations with Mixed Derivative Terms ↩︎

14. A New Stability Result for the Modified Craig–Sneyd Scheme Applied to Two-Dimensional Convection–Diffusion Equations with Mixed Derivatives ↩︎ ↩︎²

15. A Case Study on Pricing Foreign Exchange Options Using the Modified Craig–Sneyd ADI Scheme ↩︎

16. ADI Finite Difference Schemes for the Heston-Hull-White PDE ↩︎ ↩︎² ↩︎³

17. The Fractional Step Method versus the Radial Basis Functions for Option Pricing with Correlated Stochastic Processes ↩︎

18. Some Facts about Operator-Splitting and Alternating Direction Methods ↩︎ ↩︎²

19. Operator Splitting(Splitting Methods in Communication, Imaging, Science, and Engineering) ↩︎ ↩︎²

20. Analysis of Operator Splitting for Advection–Diffusion–Reaction Problems from Air Pollution Modelling ↩︎ ↩︎²

21. On Solving Pentadiagonal Linear Systems via Transformations ↩︎

22. Operator Splitting Methods for American Option Pricing ↩︎

23. Operator Splitting Methods for Pricing American Options Under Stochastic Volatility ↩︎

24. ADI Schemes with Ikonen-Toivanen Splitting for Pricing American Put Options in the Heston Model ↩︎

25. Efficient Numerical Methods for Pricing American Options Under Stochastic Volatility ↩︎

26. A Quick Operator Splitting Method for Option Pricing ↩︎

27. On the Smoothing Property of the Crank-Nicolson Scheme ↩︎

28. ADI FD Schemes for the Numerical Solution of the Three-Dimensional Heston-Cox-Ingersoll-Ross PDE ↩︎

29. ADI Finite Difference Schemes for Option Pricing in the Heston Model with Correlation ↩︎

30. ADI Schemes for Valuing European Options Under the Bates Model ↩︎

31. Efficient and Stable Numerical Solution of the Heston–Cox–Ingersoll–Ross Partial Differential Equation by Alternating Direction Implicit Finite Difference Schemes ↩︎

32. The Beta Stochastic Volatility Model ↩︎

33. Derivatives Pricing Under Beta Stochastic Volatility Model Using ADI Schemes ↩︎

34. A Comparison Study of ADI and Operator Splitting Methods on Option Pricing Models ↩︎

35. A Comparison Study of ADI and LOD Methods on Option Pricing Models ↩︎