实践中的 Heston 模型之半解析解

目录

• 实践中的 Heston 模型之半解析解
  • 引言
  • 简单复习复变函数
  • 计算机中的复数计算
    • 造成的问题
  • 不连续性的解决办法
    • 直接积分
    • Rotation Count
    • 追踪前值
    • 修正公式
  • 展开
  • 参考文献
  • 附录

实践中的 Heston 模型之半解析解

引言

对于随机波动率驱动的资产价格过程，

\[dS = (r-q)Sdt + \sqrt{v}Sdz_1(t)\\ dv=\kappa[\theta - v]dt+\sigma\sqrt{v}dz_2(t)\\ d[z_1(t),z_2(t)] = \rho dt \]

Heston 在文献【1】中给出了对应欧式期权的“半解析解”，借助特征函数方法将期权价格表达为一个数值积分问题，以看涨期权（执行价格 \(K\)，期限 \(\tau\)）为例：

\[C(S_0,v_0,\tau) = S_0e^{-q\tau}P_1 - Ke^{-r\tau}P_2 \]

其中

\[P_{j}=\frac{1}{2} + \frac{1}{\pi} \int_{0}^{\infty} \operatorname{Re}\left[\frac{e^{-i \phi \ln K} f_{j}(x, v_0, \tau ; \phi)}{i \phi}\right] d \phi \]

或者

\[P_{j}=\frac{1}{2} + \frac{1}{\pi} \int_{0}^{\infty} \operatorname{Im}\left[\frac{e^{-i \phi \ln K} f_{j}(x, v_0, \tau ; \phi)}{\phi}\right] d \phi \]

而

\[\begin{aligned} &f_j(x,v_0,\tau;\phi) = \exp\left(C(\tau,\phi)+D(\tau,\phi)v_0+i\phi x\right)\\ &C(\tau ; \phi) = i(r-q)\phi\tau + \frac{a}{\sigma^{2}}\left\{\left(b_j-i\rho\sigma\phi + d\right) \tau-2 \ln \left[\frac{1-g e^{d \tau}}{1-g}\right]\right\} \\ &D(\tau ; \phi)=\frac{b_{j}-i\rho\sigma\phi + d}{\sigma^{2}}\left[\frac{1-e^{d \tau}}{1-g e^{d \tau}}\right]\\ &u_1 = 1/2 ; u_2 =-1/2\\ &a=\kappa\theta\\ &b_1 = \kappa-\rho\sigma;b_2 = \kappa\\ &d=\sqrt{ (i\rho\sigma\phi-b_j)^2 - \sigma^2(2iu_j\phi-\phi^2)}\\ &g=\frac{b_j-i\rho\sigma\phi+d}{b_j-i\rho\sigma\phi-d}\\ &x=\ln S_0\\ &j = 1,2 \end{aligned} \]

然而 Heston 所给出的公式仅是理论上正确的，在工程实践中却存在着缺陷和隐患（文献【2】）。问题就出在 \(f_j\) 身上，更具体地讲，是因为 \(f_j\) 涉及了复数的开方和对数计算，通常的标准库无法完美地表达复数的平方根和对数。

简单复习复变函数

欧拉公式将复数、三角函数、极坐标和指数函数四者联系了起来，

\[z = r(\cos\theta + i\sin\theta) = re^{i\theta},r \ge 0 \]

其中 \(r\) 是复数 \(z\) 的绝对值（模长），\(\theta\) 则称为 \(z\) 的辐角（argument 或 phase），若强制要求 \(-\pi < \theta \le \pi\)，则称 \(\theta\) 为 \(z\) 的主辐角（principal argument）。

对于代数计算和某些函数来说，实数复数之间在形式上没有区别，然而对于某些函数计算来说，实数复数之间存在明显的差异。

根据欧拉公式，复数和三角函数之间存在着联系，这种联系使得某些复变量函数存在着“多值性”。不幸的是，开方和对数函数就是典型的“多值函数”。根据定义，\(\sqrt{z}\) 有两个值，而 \(\ln{z}\) 则有无限个值，

\[\sqrt{z} = \sqrt{r} (\cos \frac{\theta + 2k\pi}{2} + i\sin \frac{\theta + 2k\pi}{2}), k=0,1\\ \ln{z} = \ln{r} + i(\theta + 2k\pi),k\in Z \]

有关复变函数更详尽的介绍请浏览在线教科书 Complex Analysis。

计算机中的复数计算

计算机中“函数”的概念更接近实变函数，即函数只能返回“一个”确定的值。通常的标准库无法处理复变函数中的这种多值性，一旦遇到多值函数，返回值会被强制调整到一个特定范围内。

以 C++ 的 STL 和 Python 的 cmath 为例，复数的对数结果被要求虚部落在 \(-\pi\) 和 \(\pi\) 之间（也就是变量的主辐角上），复数的开方结果被要求实部是非负数。这种强制措施会导致计算结果和人的直觉预期大相径庭。

以 Python 为例，

from cmath import log, sqrt, exp, pi

c0 = complex(-1.0, 0.01)
c1 = complex(-1.0, -0.01)

print(sqrt(c0))
print(sqrt(c1))

结果：

(0.004999937502734214+1.0000124996093955j)
(0.004999937502734214-1.0000124996093955j)

c0 和 c1 是两个接近的数，但开方的结果相去甚远，虚部的符号完全相反，破坏了连续性的预期。

c = complex(1.0, 4.0)

print(log(exp(c)))

结果：

(1-2.2831853071795867j)

c 和 log(exp(c)) 则变成了完全不相等的两个数，原本预期 log(exp(c)) 的虚部依然等于 4，但 4 大于 \(\pi\)，最终被强制减去 \(2\pi\)才能落在主辐角的范围内。

造成的问题

标准库对复数的强制处理将可能导致 Heston 公式中的积分项出现“间断”，进而得到错误的定价结果，文献【2】中给出了不连续的案例，并分析了错误定价的误差大小，通常来说期限越长误差越大。

不连续性的解决办法

在 Heston 的公式问世一段时间之后，业内人士才发现模型实践过程中隐藏的陷阱，并先后提出了几种办法解决积分项的不连续性。

直接积分

文献【1】中的（A7）显示，\(C\) 和 \(D\) 其实是一对特定 ODE 的解，若直接求 ODE 的数值解就可以避免遇到复数的多值函数。但代价是计算量过大，文献【3】指出，直接积分法本身耗时就长，而且随着期权期限的增长而增长。

Rotation Count

前文提到，复数的开方和对数函数是多值函数，而标准库会强制调整计算结果的辐角或虚部，正是这种强制调整带来了不连续性。所以，若已知调整的幅度，就可以避免调用标准库函数，自行计算正确的辐角或虚部。

幸运的是 Lord 和 Kahl 在文献【4】中指出 \(f_j\) 关于 \(d\) 是个偶函数，因此无须对开方得到结果进行调整，而 Jackel 和 Kahl 在文献【5】中提出 Rotation Count 算法来计算复数对数的正确结果。其大体思路是分别推测出标准库对 \(g-1\) 和 \(ge^\tau-1\) 辐角的调整幅度，再不借助标准库而直接计算出 \(\ln \frac{ge^\tau-1}{g-1}\)。文献【4】对 Rotation Count 做了更详细的分析。

追踪前值

另一个调整复数对数虚部的算法由 Kangro 等人在文献【6】中提出，通过持续追踪先前复数对数虚部的值来调整当前复数对数的虚部。

注：QuantLib 中 AnalyticHestonEngine 的枚举型 BranchCorrection 就是指追踪法。

修正公式

文献【7】、【8】和【9】提出了等价形式的定价公式，仅需要对 Heston 的结果稍加变化即可规避复对数的不连续性，非常巧妙。文献【2】详细论证了修正公式为何不会出现不连续性。

沿用 Heston 的符号，修正公式如下：

\[\begin{aligned} &f_j^{adj}(x,v_0,\tau;\phi) = \exp\left(C^{adj}(\tau,\phi)+D^{adj}(\tau,\phi)v+i\phi x\right)\\ &C^{adj}(\tau ; \phi) = i(r-q)\phi\tau + \frac{a}{\sigma^{2}}\left\{\left(b_j-i\rho\sigma\phi - d\right) \tau-2 \ln \left[\frac{1-g^{-1}e^{-d \tau}}{1-g^{-1}}\right]\right\} \\ &D^{adj}(\tau ; \phi)=\frac{b_{j}-i\rho\sigma\phi - d}{\sigma^{2}}\left[\frac{1-e^{-d \tau}}{1-g^{-1} e^{-d \tau}}\right]\\ \end{aligned} \]

也就是用 \(-d\) 替换 \(d\)，同时用 \(g^{-1}\) 替换 \(g\)。

注：QuantLib 中 AnalyticHestonEngine 的枚举型 Gatheral 就是指修正公式。

展开

不单单是半解析解，当使用傅立叶变换求解 Heston 模型时 \(f_2\) 也会出现在数值积分当中（文献【12】）。

此外，使用特征函数求随机波动率模型（半）解析解这一方法论影响了一大批相关的研究，例如时间依赖参数的 Heston 模型（文献【11】）、Variance Gamma 和 Schöbel-Zhu 模型（文献【4】）以及用控制变量求解 Heston 模型（文献【10】）等等。而快速求解这些模型通常都绕不开处理复变函数中的多值函数（幂函数或对数函数）。

正确处理多值函数，并确保积分项的连续性是实践随机波动率模型的关键。

参考文献

【1】A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options
【2】The Little Heston Trap
【3】Implementation Problems and Solutions in Stochastic Volatility Models of the Heston Type
【4】Complex Logarithms in Heston‐Like Models
【5】Not-so-complex logarithms in the Heston model
【6】Pricing European-style options under jump diffusion processes with stochastic volatility: Applications of Fourier transform
【7】A Note on the Discontinuity Problem in Heston's Stochastic Volatility Model
【8】A Perfect Calibration! Now What?
【9】The Volatility Surface: A Practitioner's Guide
【10】Interest Rate Modeling（Andersen 和 Piterbarg）
【11】Models with Time-Dependent Parameters Using Transform Methods: Application to Heston's Model
【12】A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions

附录

C++ STL 中模板类 complex 的说明。