KoBol分数阶期权定价模型的数值方法

摘要

自Black⁃Scholes期权定价模型提出以来， 大量的期权定价模型被陆续提出并加以研究，成为国内外金融工程和金融数学的研究热点。由于列维过程能够很好地描述资产运动的动力学特征， 近年来基于列维过程的期权定价模型吸引了广泛关注， 如FMLS（finite moment log stable）、CGMY和KoBol模型。这些模型最终归结为数值求解一类分数阶偏微分方程。为此提出了求解这类分数阶偏微分方程的数值离散格式， 理论分析给出了数值格式稳定的充分条件。数值实验验证数值格式和算法的可行性和有效性。基于上证50与沪深300的股指期权实际交易数据， 利用KoBol分数阶模型进行定价并反演计算波动率曲线， 进一步验证了KoBol模型在真实市场中的有效性。

关键词

股指期权; 期权定价; 列维过程; 分数阶方程; 有限差分

经典的Black⁃Scholes（BS）模型^［1］假设标的资产的收益变化服从几何布朗运动， 也就是资产的价格服从对数正态分布， 设$S_t$为标的资产在$t$时刻的价格， $r$为无风险利率， $\sigma$为波动率， $V\left(x,t\right)$为$t$时刻时 $x=\mathrm{l}\mathrm{n}(S_t)$下的期权价格， 那么$V\left(x,t\right)$满足以下方程：

$$\frac{\partial V\left(x,t\right)}{\partial t}+\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}\frac{{\partial }^{\mathrm{2}}V\left(x,t\right)}{\partial x^{\mathrm{2}}}+$$

$$\left(r-\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}\right)\frac{\partial V\left(x,t\right)}{\partial x}=rV\left(x,t\right)$$

Black⁃Scholes期权定价模型是金融工程中一个重大突破， 为风险中性条件下的期权定价理论提供了有力的数学支撑。但是， 该模型成立的前提条件非常严格， 例如市场交易无摩擦， 市场不存在套利机会， 资产收益服从正态分布， 在这种情况下由该模型计算出的结果往往与市场真实的情况不相吻合。实际市场中的隐含分布与正态分布相比存在尖峰肥尾的现象， 因此在期权市场中会存在波动率微笑^［2］。

为了弱化其模型假设对定价带来的影响， 很多学者都在此基础上做了深入的研究。Merton^［3］提出了跳跃过程服从对数正态分布的跳⁃扩散模型。Kou^［4］假设资产跳跃服从对数双指数分布， 从而构造了相应的跳⁃扩散模型。 还有随机波动率模型^［5- 6］， 考虑交易费用的期权定价模型^［7］等。对于跳⁃扩散模型^［8-10］以及美式期权^［11］的求解也有许多不同的数值方法。

随着研究的深入， Mandelbrot^［12］发现股票的收益分布有长尾的特点， 在此基础上推导出用$\alpha$稳定的列维过程代替标准的几何布朗运动， 用纯无限跳跃模型来描述金融资产价格的变化， 其中$\alpha <\mathrm{2}$。在此基础上，许多学者利用列维过程来代替Black⁃Scholes模型中的布朗运动， 从而提出了相应基于列维过程的期权定价模型。Carr等^［13］提出FMLS(finite moment log stable)模型， 可以表示标的资产对数收益的倾斜密度特征。 Koponen^［14］、Boyarchenko等^［15］将修正$\alpha$稳定的列维过程应用在模拟标的资产的动力学特征上， 称为KoBoL模型。Carr等^［16］提出允许标的资产价格出现有限或者无限跳跃幅度的过程， 称之为CGMY模型。近年来， 这些由列维过程推导出的金融衍生品定价模型由于可以准确描述标的运动而得到广泛且深入的研究^［17］。

利用分数阶模型进行期权定价时， 需要求解一个分数阶对流扩散方程。Marom等^［18］比较了上述3种分数阶期权定价模型定价欧式期权的数值结果， 但并未给出相应的稳定性与收敛性条件。Wang等^［19］给出了具有一阶精度的差分离散格式， 并在障碍期权上进行定价。Meng等^［20］利用CGNR算法对欧式看涨期权进行定价， 并与Black⁃Scholes模型进行了比较。Zhang等^［21］利用BiCGSTAB算法求解了单边分数阶欧式看涨定价。

本文主要研究一类双边分数阶期权定价方程的数值解法。首先对于这类双边分数阶方程的一般形式， 利用带位移的Grünwald格式， 给出每个时间层上的离散格式， 并分析了迭代格式的数值稳定性。然后结合KoBol模型下迭代矩阵的特殊结构， 构造了预处理Krylov子空间方法进行求解。最后， 在数值实验中对欧式看涨期权进行定价， 并使用国内上证50ETF与沪深300ETF场内期权数据进行实证研究， 验证算法有效性的同时也体现了KoBol模型在真实市场中的有效性。

1 基于分数阶列维过程期权定价模型

对于列维过程$X_t$， $t\ge \mathrm{0}$， 其为增量独立固定的随机过程， 且路径依概率连续。不失一般性， 假设$X_{\mathrm{0}}\equiv \mathrm{0}$。$X_t$的对数特征函数有以下Lévy⁃Khintchine表达式：

$$\mathrm{l}\mathrm{n}\left[{\mathrm{e}}^{\mathrm{i}\xi X_t}\right]\equiv t\Psi \left(\xi \right)=m\mathrm{i}t\xi -\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}t{\xi }^{\mathrm{2}}+$$

$$t{\int }_{\mathbf{R}\backslash \{\mathrm{0}\}}\left({\mathrm{e}}^{\mathrm{i}\xi x-\mathrm{1}-\mathrm{i}\xi h\left(x\right)}\right)W(\mathrm{d}x)$$

式中：$m\in \mathbf{R}$；$\sigma \ge \mathrm{0}$； $\mathrm{i}=\sqrt[]{-\mathrm{1}}$； $h\left(x\right)$为截断函数。列维测度$W$满足以下公式：

$${\int }_{\mathbf{R}}^{}\mathrm{m}\mathrm{i}\mathrm{n}\left\{\mathrm{1},x^{\mathrm{2}}\right\}W(\mathrm{d}x)<\mathrm{\infty }$$

$\Psi (\xi )$为列维过程的特征指数， 列维测度$W(\mathrm{d}x)$还可以写成$W\left(\mathrm{d}x\right)=w(x)\mathrm{d}x$， $w(x)$又称作列维密度函数。

特别地， KoBol模型的列维密度函数为如下形式：

$$w\left(x\right)=\left\{\begin{array}{l}Dq{\left|x\right|}^{-\mathrm{1}-\alpha }{\mathrm{e}}^{-\lambda |x|},x<\mathrm{0}\\ Dp{x}^{-\mathrm{1}-\alpha }{\mathrm{e}}^{-\lambda x},x>\mathrm{0}\end{array}\right.$$

式中：$D>\mathrm{0}$， $\lambda >\mathrm{0}$， $p$、$q\in \left[-\mathrm{1,1}\right]$， $p+q=\mathrm{1}$， $\mathrm{0}<\alpha \le \mathrm{2}$。

除了KoBol模型外， 分数阶期权定价模型还有FMLS与CGMY模型， 这3类模型都可以表示为如下的分数阶微分方程：

$$\begin{array}{l}\frac{\partial V\left(x,t\right)}{\partial t}+a\frac{\partial V\left(x,t\right)}{\partial x}+b\left(x\right)\frac{{\partial }^{\alpha }\left(\tilde{f}\left(x\right)V\left(x,t\right)\right)}{{\partial }_{+}x}+\\ c\left(x\right)\frac{{\partial }^{\alpha }\left(\tilde{h}\left(x\right)V\left(x,t\right)\right)}{{\partial }_{-}x}=dV\left(x,t\right)\end{array}$$

（1）

式中：$x\in \left(-\mathrm{\infty },+\mathrm{\infty }\right)$， $t\in \left(\mathrm{0},T\right)$， $\mathrm{0}<\alpha <\mathrm{2}$， $a$和$d$是非负常数， 函数$b(x)$和$c(x)$充分光滑， 函数$\tilde{f}$$(x)$和$\tilde{h}$$(x)$是连续的， 且均为非负函数。${\partial }^{\alpha }/{\partial }_{+}{x}^{\alpha }$与${\partial }^{\alpha }/{\partial }_{-}{x}^{\alpha }$分别为Riemann⁃Liouville左分数阶微分算子与右分数阶微分算子^［22］， 有如下形式：

$$\frac{{\partial }^{\alpha }f\left(x\right)}{{\partial }_{+}x}=\frac{\mathrm{1}}{\Gamma \left(n-\alpha \right)}\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}$$$${\int }_{-\mathrm{\infty }}^x\frac{f(\xi )}{{(x-\xi )}^{\alpha +\mathrm{1}-n}}\mathrm{d}\xi$$

$$\frac{{\partial }^{\alpha }f\left(x\right)}{{\partial }_{-}x}=\frac{{(-\mathrm{1})}^n}{\Gamma \left(n-\alpha \right)}\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}$$$${\int }_x^{+\mathrm{\infty }}\frac{f(\xi )}{{(\xi -x)}^{\alpha +\mathrm{1}-n}}\mathrm{d}\xi$$

式中：$\alpha \in \left(n-\mathrm{1},n\right)$， $n$为整数， $\Gamma (\cdot )$表示Gamma函数， 在本文中只考虑$\alpha \in \left(\mathrm{1,2}\right)$的情况。

KoBol模型中的相关参数可以作如下表示：

$$\begin{array}{l}a=r-v-{\lambda }^{\alpha }\left(q-p\right),b\left(x\right)=\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\alpha }q{\mathrm{e}}^{-\lambda x},\\ \tilde{f}\left(x\right)={\mathrm{e}}^{\lambda x},d=r+\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\alpha }{\lambda }^{\alpha },c\left(x\right)=\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\alpha }p{\mathrm{e}}^{\lambda x},\\ \tilde{h}(x)={\mathrm{e}}^{-\lambda x},\\ v=\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\alpha }[p{(\lambda -\mathrm{1})}^{\alpha }+q{(\lambda +\mathrm{1})}^{\alpha }-\\ {\lambda }^{\alpha }-\alpha {\lambda }^{\alpha -\mathrm{1}}(q-p)]\end{array}$$

相比于BS模型中的对数正态分布， 在分数阶期权模型中， 可以通过调整参数的取值， 使隐含分布更接近市场的实际分布。例如在KoBol模型中， 可以通过参数$\alpha$和$\lambda$调整隐含分布的峰度， 参数$p$调整隐含分布的偏度， 这样在一定程度上可以消除波动率微笑对期权定价带来的影响。

期权定价问题是一个终值问题， 自变量$x$定义在无界区域$\left(-\mathrm{\infty },\mathrm{\infty }\right)$上， 为了能够使用数值方法求解该问题， 需要用合理的方法截断为$x\in [L,R]$， 参考文献［23］。在欧式期权中， 终值条件与边值条件为

$$\begin{array}{c}\left\{\begin{array}{c}V\left(x,T\right)=s\left(x\right),L\le x<R\\ V\left(L,t\right)=l\left(t\right),\mathrm{0}\le t<T,r>\mathrm{0}\\ V\left(R,t\right)=r\left(t\right),\mathrm{0}\le t<T,r>\mathrm{0}\end{array}\right.\end{array}$$

（2）

对于看涨期权有

$$\left\{\begin{array}{l}s(x)=\mathrm{m}\mathrm{a}\mathrm{x}\{{\mathrm{e}}^x-K,\mathrm{0}\}\\ l(t)=\mathrm{0}\\ r(t)={\mathrm{e}}^R-K{\mathrm{e}}^{-r(T-t)}\end{array}\right.$$

对于看跌期权有

$$\left\{\begin{array}{l}s(x)=\mathrm{m}\mathrm{a}\mathrm{x}\{K-{\mathrm{e}}^x,\mathrm{0}\}\\ l\left(t\right)=K{\mathrm{e}}^{-r\left(T-t\right)}-{\mathrm{e}}^L\\ r(t)=\mathrm{0}\end{array}\right.$$

首先对截断区域$[L,R]\times [\mathrm{0},T]$进行网格划分。 将空间层$N$等分， 步长$h=(R-L)/N$， 对应的节点为$x_n=L+nh$， $n=\mathrm{0,1},\dots ,N$； 将时间层$M$等分， 步长$\tau =T/M$， 对应的节点为$t_m=m\tau$， $m=\mathrm{0,1},\dots ,M$。函数$V(x,t)$在对应节点上取值简记为$V(x_n,t_m)=V_n^m$， 其余记号类似。

Meerschaert与Tadjeran证明了使用Grünwald格式离散分数阶扩散方程得到的迭代格式不稳定， 并提出了带位移的Grünwald⁃Letnikov格式^［24］如下：

$$\begin{array}{cc}\frac{{\partial }^{\alpha }(\tilde{f}V)}{{\partial }_{+}{x}^{\alpha }}(x_n,t_m)& =\frac{\mathrm{1}}{h^{\alpha }}\sum _{k=\mathrm{0}}^{(x_n-x_L)/h+\mathrm{1}}{g}_k^{(\alpha )}\tilde{f}(x_n-(k-\mathrm{1})h)V(x-(k-\mathrm{1})h,t)+O(h)=\\ & \frac{\mathrm{1}}{h^{\alpha }}\sum _{k=\mathrm{0}}^{n+\mathrm{1}}{g}_k^{(\alpha )}\tilde{f}{}_{n-k+\mathrm{1}}{V}_{n-k+\mathrm{1}}^m+O(h)\end{array}$$

$$\begin{array}{cc}\frac{{\partial }^{\alpha }(\tilde{h}V)}{{\partial }_{-}{x}^{\alpha }}(x_n,t_m)& =\frac{\mathrm{1}}{h^{\alpha }}\sum _{k=\mathrm{0}}^{(x_R-x_n)/h+\mathrm{1}}{g}_k^{\left(\alpha \right)}\tilde{h}(x_n+(k-\mathrm{1})h)V(x+(k-\mathrm{1})h,t)+O(h)=\\ & \frac{\mathrm{1}}{h^{\alpha }}\sum _{k=\mathrm{0}}^{N-n+\mathrm{1}}{g}_k^{(\alpha )}{\tilde{h}}_{n+k-\mathrm{1}}{V}_{n+k-\mathrm{1}}^m+O(h)\end{array}$$

其中，$g_k^{(\alpha )}$=（-1）^k $\left(\begin{array}{c}\alpha \\ k\end{array}\right)=$（-1）^kΓ（α+1）/（Γ（k+1）·Γ（α-k+1）） 称为Grünwald⁃Letnikov系数。 $g_k^{\alpha }={(-\mathrm{1})}^k\left(\begin{array}{c}\alpha \\ k\end{array}\right)$实际上是生成函数$g(z)={\left(\mathrm{1}-z\right)}^{\alpha }$的泰勒展开多项式系数， 该系数还可以更简单地由如下递推公式直接求出：

$$\begin{array}{l}{g}_{\mathrm{0}}^{\left(\alpha \right)}=\mathrm{1},\\ g_k^{\left(\alpha \right)}=\left(\mathrm{1}-\frac{\alpha +\mathrm{1}}{k}\right)g_{k-\mathrm{1}}^{\left(\alpha \right)},\end{array}$$$$k=\mathrm{1,2},\dots$$

当$\mathrm{1}<\alpha <\mathrm{2}$时， Grünwald⁃Letnikov系数$g_j^{\left(\alpha \right)}={\left(-\mathrm{1}\right)}^j\left(\begin{array}{c}\alpha \\ j\end{array}\right)$有如下性质：

$$\left\{\begin{array}{l}{g}_{\mathrm{0}}^{(\alpha )}=\mathrm{1},g_{\mathrm{1}}^{(\alpha )}=-\alpha <\mathrm{0},\mathrm{0}\le \cdots \le g_{\mathrm{3}}^{(\alpha )}\le g_{\mathrm{2}}^{(\alpha )}\le \mathrm{1}\\ \sum _{k=\mathrm{0}}^{+\mathrm{\infty }}{g}_k^{(\alpha )}=\mathrm{0},\sum _{k=\mathrm{0}}^l{g}_k^{(\alpha )}<\mathrm{0},\mathrm{0}<l<\mathrm{\infty }\end{array}\right.$$

对方程$(\mathrm{1})$使用中心差分格式离散对流项， 得到如下的半离散格式：

$$\frac{\partial V}{\partial t}+a\frac{V_{n+\mathrm{1}}^m-V_{n-\mathrm{1}}^m}{\mathrm{2}h}+$$

$$b_n\frac{\mathrm{1}}{h^{\alpha }}\sum _{k=\mathrm{0}}^{n+\mathrm{1}}{g}_k^{\left(\alpha \right)}{\tilde{f}}_{n-k+\mathrm{1}}{V}_{n-k+\mathrm{1}}^m+$$

$$c_n\frac{\mathrm{1}}{h^{\alpha }}\sum _{k=\mathrm{0}}^{N-n+\mathrm{1}}{g}_k^{(\alpha )}{\tilde{h}}_{n+k-\mathrm{1}}{V}_{n+k-\mathrm{1}}^m-d{V}_n^m=\mathrm{0}$$

式中：$n=\mathrm{1,2},\cdots ,N-\mathrm{1}$。记${\mathit{V}}^{\mathit{m}}={\left(V_{\mathrm{1}}^m,V_{\mathrm{2}}^m,\cdots ,V_{N-\mathrm{2}}^m,V_{N-\mathrm{1}}^m\right)}^{\mathrm{T}}$。当$n$取遍$\mathrm{1,2},\cdots ,N-\mathrm{1}$时， 可以得到半离散方程组， 写成如下的矩阵形式：

$$\frac{\partial V}{\partial t}+\mathit{S}{\mathit{V}}^{\mathit{m}}+{\mathit{q}}^{\mathit{m}}=\mathrm{0}$$

（3）

其中

$$\mathit{S}=\frac{a}{\mathrm{2}h}\mathit{J}+\frac{\mathrm{1}}{h^{\alpha }}\mathit{B}\cdot \mathit{G}\cdot \tilde{\mathit{F}}+\frac{\mathrm{1}}{h^{\alpha }}\mathit{C}\cdot {\mathit{G}}^{\mathit{*}}\cdot \tilde{H}-d\mathit{I}$$

（4）

$${\mathit{q}}^{\mathit{m}}={\mathit{q}}_{\mathrm{1}}^{\mathit{m}}+{\mathit{q}}_{\mathrm{2}}^{\mathit{m}}+{\mathit{q}}_{\mathrm{3}}^{\mathit{m}}+{\mathit{q}}_{\mathrm{4}}^{\mathit{m}}$$

$$\mathit{J}=$$$${\left(\begin{array}{cccccc}\mathrm{0}& \mathrm{1}& & & & \\ -\mathrm{1}& \mathrm{0}& \mathrm{1}& & \mathrm{0}& \\ & -\mathrm{1}& \mathrm{0}& \mathrm{1}& & \\ & & -\mathrm{1}& \ddots & \ddots & \\ & \mathrm{0}& & \ddots & \mathrm{0}& \mathrm{1}\\ & & & & -\mathrm{1}& \mathrm{0}\end{array}\right)}_{(N-\mathrm{1})\times (N-\mathrm{1})}$$

$$\mathit{G}=$$$${\left(\begin{array}{cccccc}{g}_{\mathrm{1}}^{(\alpha )}& g_{\mathrm{0}}^{(\alpha )}& \begin{array}{cc}\mathrm{0}& \cdots \end{array}& \mathrm{0}& \cdots & \mathrm{0}\\ g_{\mathrm{2}}^{(\alpha )}& g_{\mathrm{1}}^{(\alpha )}& \begin{array}{cc}{g}_{\mathrm{0}}^{(\alpha )}& \cdots \end{array}& \mathrm{0}& \cdots & \mathrm{0}\\ g_{\mathrm{3}}^{(\alpha )}& g_{\mathrm{2}}^{(\alpha )}& \begin{array}{cc}{g}_{\mathrm{1}}^{(\alpha )}& g_{\mathrm{0}}^{(\alpha )}\end{array}& \cdots & \cdots & \mathrm{0}\\ \begin{array}{c}⋮\\ ⋮\end{array}& \begin{array}{c}⋮\\ ⋮\end{array}& \begin{array}{c}\begin{array}{cc}\ddots & \ddots \end{array}\\ \begin{array}{cc}⋮& \ddots \end{array}\end{array}& \begin{array}{c}\ddots \\ \ddots \end{array}& \begin{array}{c}\cdots \\ \ddots \end{array}& \begin{array}{c}⋮\\ ⋮\end{array}\\ g_{N-\mathrm{2}}^{(\alpha )}& \cdots & \begin{array}{cc}\cdots & \cdots \end{array}& g_{\mathrm{2}}^{(\alpha )}& g_{\mathrm{1}}^{(\alpha )}& g_{\mathrm{0}}^{(\alpha )}\\ g_{N-\mathrm{1}}^{(\alpha )}& g_{N-\mathrm{2}}^{(\alpha )}& \begin{array}{cc}\cdots & \cdots \end{array}& g_{\mathrm{3}}^{(\alpha )}& g_{\mathrm{2}}^{(\alpha )}& g_{\mathrm{1}}^{(\alpha )}\end{array}\right)}_{(N-\mathrm{1})\times (N-\mathrm{1})}$$

$${\mathit{q}}_{\mathrm{1}}^{\mathit{m}}=\frac{a}{\mathrm{2}h}\left(\begin{array}{c}-V_{\mathrm{0}}^m\\ \mathrm{0}\\ ⋮\\ \mathrm{0}\\ V_N^m\end{array}\right)$$，$${\mathit{q}}_{\mathrm{2}}^{\mathit{m}}=\frac{{\tilde{f}}_{\mathrm{0}}{V}_{\mathrm{0}}^m}{h^{\alpha }}B\left(\begin{array}{c}{g}_{\mathrm{2}}^{(\alpha )}\\ g_{\mathrm{3}}^{(\alpha )}\\ ⋮\\ g_{N-\mathrm{1}}^{(\alpha )}\\ g_N^{(\alpha )}\end{array}\right)$$，

$${\mathit{q}}_{\mathrm{3}}^{\mathit{m}}=\frac{{\tilde{h}}_N{V}_N^m}{h^{\alpha }}C\left(\begin{array}{c}{g}_N^{(\alpha )}\\ g_{N-\mathrm{1}}^{(\alpha )}\\ ⋮\\ g_{\mathrm{3}}^{(\alpha )}\\ g_{\mathrm{2}}^{(\alpha )}\end{array}\right)$$， $${\mathit{q}}_{\mathrm{4}}^{\mathit{m}}=\frac{g_{\mathrm{0}}}{h^{\alpha }}\left(\begin{array}{c}{c}_{\mathrm{1}}{\tilde{h}}_{\mathrm{0}}{V}_{\mathrm{0}}^m\\ \mathrm{0}\\ ⋮\\ \mathrm{0}\\ b_{N-\mathrm{1}}{\tilde{f}}_N{V}_N^m\end{array}\right)$$

$$\mathit{B}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(b\left(x_{\mathrm{1}}\right),b\left(x_{\mathrm{2}}\right),\dots ,b\left(x_{N-\mathrm{1}}\right)\right)$$

$$\tilde{\mathit{F}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\tilde{f}(x_{\mathrm{1}}),\tilde{f}(x_{\mathrm{2}}),\dots ,\tilde{f}(x_{N-\mathrm{1}})\right)$$

$$\mathit{C}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(c\left(x_{\mathrm{1}}\right),c\left(x_{\mathrm{2}}\right),\dots ,c\left(x_{N-\mathrm{1}}\right)\right)$$

$$\tilde{H}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\tilde{h}\left(x_{\mathrm{1}}\right),\tilde{h}\left(x_{\mathrm{2}}\right),\dots ,\tilde{h}\left(x_{N-\mathrm{1}}\right)\right)$$

对于半离散格式$(\mathrm{3})$， 时间方向的采用加权隐式差分格式离散， 第$m$层与第$m+\mathrm{1}$层的加权平均如下：

$$\frac{{\mathit{V}}^{\mathit{m}+\mathrm{1}}-{\mathit{V}}^{\mathit{m}}}{\tau }+\theta \left(\mathit{S}{\mathit{V}}^{\mathit{m}+\mathrm{1}}+{\mathit{q}}^{\mathit{m}+\mathrm{1}}\right)+$$

$$\left(\mathrm{1}-\theta \right)\left(\mathit{S}{\mathit{V}}^{\mathit{m}}+{\mathit{q}}^{\mathit{m}}\right)=\mathrm{0}$$

整理得

$$\begin{array}{l}\left(\mathit{I}-\left(\mathrm{1}-\theta \right)\tau \mathit{S}\right){\mathit{V}}^{\mathit{m}}=\left(\mathit{I}+\theta \tau \mathit{S}\right){\mathit{V}}^{\mathit{m}+\mathrm{1}}+\\ \tau \left(\theta {\mathit{q}}^{\mathit{m}+\mathrm{1}}+(\mathrm{1}-\theta ){\mathit{q}}^{\mathit{m}}\right)\end{array}$$

（5）

其中$\theta \in [\mathrm{0,1}]$。 特别地， 当$\theta =\mathrm{1}$时， 为显式格式

$${\mathit{V}}^{\mathit{m}}=\left(\mathit{I}+\tau \mathit{S}\right){\mathit{V}}^{\mathit{m}+\mathrm{1}}+\tau {\mathit{q}}^{\mathit{m}+\mathrm{1}}$$

当$\theta =\mathrm{0}$时， 为隐式格式

$$\left(\mathit{I}-\tau \mathit{S}\right){\mathit{V}}^{\mathit{m}}={\mathit{V}}^{\mathit{m}+\mathrm{1}}+\tau {\mathit{q}}^{\mathit{m}}$$

当$\theta =\mathrm{1}/\mathrm{2}$时， 为Crank⁃Nicolson格式

$$\begin{array}{l}\left(\mathit{I}-\frac{\mathrm{1}}{\mathrm{2}}\tau \mathit{S}\right){\mathit{V}}^{\mathit{m}}=\left(\mathit{I}+\frac{\mathrm{1}}{\mathrm{2}}\tau \mathit{S}\right){\mathit{V}}^{\mathit{m}+\mathrm{1}}+\\ \frac{\mathrm{1}}{\mathrm{2}}\tau \left({\mathit{q}}^{\mathit{m}}+{\mathit{q}}^{\mathit{m}+\mathrm{1}}\right)\end{array}$$

（6）

Tadjeran^［25］证明了带位移的Grünwald格式在空间层上为一阶精度， 又因为Crank⁃Nicolson格式在时间层上为二阶精度， 那么差分格式$(\mathrm{6})$的截断误差为$O\left({\tau }^{\mathrm{2}}\right)+O(h)$， 所以该格式相容。在此基础上， 对式$(\mathrm{6})$的稳定性做以下分析。

引理1   （Gerschgorin圆盘定理） 设A=［a_ij］ ∈ ${\mathit{C}}^{n\times n}$， 令

$$G_i\left(\mathit{A}\right)=\left\{z\in \mathit{C}:\left|z-a_{ii}\right|\le \sum _{j\ne i}\left|a_{ij}\right|\right\},$$

$$i=\mathrm{1},\cdots ,n$$

则有

$$\lambda \left(\mathit{A}\right)\subset G_{\mathrm{1}}\left(\mathit{A}\right)\bigcup G_{\mathrm{2}}\left(\mathit{A}\right)\bigcup \cdots \bigcup G_n\left(\mathit{A}\right)$$

命题1   若$\mathrm{1}<\alpha \le \mathrm{2}$， 当下式成立时，

$$\frac{\left(b_M\left({\tilde{f}}_M-{\tilde{f}}_m\right)+c_M\left({\tilde{h}}_M-{\tilde{h}}_m\right)\right)\alpha }{h^{\alpha }}-d+\frac{a}{h}<\mathrm{0}$$

（7）

Crank⁃Nicolson离散格式$(\mathrm{6})$稳定。 其中$\tilde{f}$${}_M=\mathrm{m}\mathrm{a}\mathrm{x}\left(\tilde{f}\right)$， ${\tilde{h}}_M=\mathrm{m}\mathrm{a}\mathrm{x}\left(\tilde{h}\right)$， $\tilde{f}$${}_m=\mathrm{m}\mathrm{i}\mathrm{n}\left(\tilde{f}\right)$， ${\tilde{h}}_m=\mathrm{m}\mathrm{i}\mathrm{n}\left(\tilde{h}\right)$， $b_M=\mathrm{m}\mathrm{a}\mathrm{x}\left(b\right)$， $c_M=\mathrm{m}\mathrm{a}\mathrm{x}\left(c\right)$。

证明   记$\mathit{Z}={\left(\mathit{I}-\frac{\mathrm{1}}{\mathrm{2}}\tau \mathit{S}\right)}^{-\mathrm{1}}\left(\mathit{I}+\frac{\mathrm{1}}{\mathrm{2}}\tau \mathit{S}\right)$， 若要满足$\rho \left(\mathit{Z}\right)<\mathrm{1}$， 则矩阵$\mathit{S}$的特征值${\lambda }_S$需满足以下方程：

$$\left|\frac{\mathrm{1}+\frac{\mathrm{1}}{\mathrm{2}}\tau {\lambda }_S}{\mathrm{1}-\frac{\mathrm{1}}{\mathrm{2}}\tau {\lambda }_S}\right|<\mathrm{1}$$

解之得$\mathrm{R}\mathrm{e}\left({\lambda }_S\right)<\mathrm{0}$。记S_nn为矩阵S第n行的对角元，r_n为该行所对应的圆盘半径， 由式$(\mathrm{4})$ 可以得到

$$S_{nn}=\frac{b_n{g}_{\mathrm{1}}^{\left(\alpha \right)}{\tilde{f}}_n+c_n{g}_{\mathrm{1}}^{\left(\alpha \right)}{\tilde{h}}_n}{h^{\alpha }}-d=$$

$$-\frac{\left(b_n{\tilde{f}}_n+c_n{\tilde{h}}_n\right)\alpha }{h^{\alpha }}-d<\mathrm{0}$$

当$n=\mathrm{2},\cdots ,N-\mathrm{2}$时， 有

$$r_n=b_n\frac{\mathrm{1}}{h^{\alpha }}\sum _{k=\mathrm{0},k\ne \mathrm{1}}^n{g}_k^{\left(\alpha \right)}{\tilde{f}}_{n-k+\mathrm{1}}+$$

$$c_n\frac{\mathrm{1}}{h^{\alpha }}\sum _{k=\mathrm{0},k\ne \mathrm{1}}^{N-n}{g}_k^{(\alpha )}{\tilde{h}}_{n+k-\mathrm{1}}+\frac{a}{h}\le$$

$$\frac{b_n{\tilde{f}}_M+c_n{\tilde{h}}_M}{h^{\alpha }}\alpha +\frac{a}{h}$$

因为$\mathrm{R}\mathrm{e}\left({\lambda }_S\right)<\mathrm{0}$， 所以有$S_{nn}+r_n<\mathrm{0}$， 即

$$\frac{\left(b_n\left({\tilde{f}}_M-{\tilde{f}}_n\right)+c_n\left({\tilde{h}}_M-{\tilde{h}}_n\right)\right)\alpha }{h^{\alpha }}-d+\frac{a}{h}<\mathrm{0}$$

同理， 当$n=\mathrm{1}$或$n=N-\mathrm{1}$时， 有

$$r_n\le \frac{b_n{\tilde{f}}_M+c_n{\tilde{h}}_M}{h^{\alpha }}\alpha +\frac{a}{\mathrm{2}h}$$

同样有

$$\frac{\left(b_n\left({\tilde{f}}_M-{\tilde{f}}_n\right)+c_n\left({\tilde{h}}_M-{\tilde{h}}_n\right)\right)\alpha }{h^{\alpha }}-d+\frac{a}{\mathrm{2}h}<\mathrm{0}$$

合并之后有， 当下式成立时，

$$\frac{\left(b_M\left({\tilde{f}}_M-{\tilde{f}}_m\right)+c_M\left({\tilde{h}}_M-{\tilde{h}}_m\right)\right)\alpha }{h^{\alpha }}-d+\frac{a}{h}<\mathrm{0}$$

Crank⁃Nicolson离散格式$\left(\mathrm{6}\right)$稳定。

2 预处理技术

在求解大规模稀疏线性方程组时， 以CG、 GMRES、BiCGSTAB和CGNR为代表的Krylov子空间迭代法是目前广泛使用的方法^［26］， 并在金融领域有广泛的应用^［20-21］。这类子空间迭代法的收敛速度与迭代矩阵特征值的聚集程度有关， 使用预处理方法能有效改善原线性方程组系数矩阵的性质， 使得相应迭代方法的步数和求解时间大大减少， 提高计算稳定性和计算效率。

预处理方法是指对于线性方程组$\mathit{A}\mathit{x}=\mathit{b}$， 其中$\mathit{A}$为系数矩阵， $\mathit{b}$为右端向量， $\mathit{x}$为需要求解的向量， 寻找一个非奇异矩阵$\mathit{H}$， 然后应用Krylov子空间迭代法求解以下同解线性方程组：

$${\mathit{H}}^{-\mathrm{1}}\mathit{A}\mathit{x}={\mathit{H}}^{-\mathrm{1}}\mathit{b}$$

或

$$\mathit{A}{\mathit{H}}^{-\mathrm{1}}\mathit{y}=\mathit{b},\mathit{x}={\mathit{H}}^{-\mathrm{1}}\mathit{y}$$

相应得到原算法的左预处理格式与右预处理格式， 其中$\mathit{H}$称为预处理矩阵。考虑到右预处理方法不会改变GMRES算法中的残差， 在本文中使用右预处理格式。

在上节中离散得到的Crank⁃Nicolson格式如下：

$$\left(\mathit{I}-\frac{\mathrm{1}}{\mathrm{2}}\tau \mathit{S}\right){\mathit{V}}^{\mathit{m}}=\left(\mathit{I}+\frac{\mathrm{1}}{\mathrm{2}}\tau \mathit{S}\right){\mathit{V}}^{\mathit{m}+\mathrm{1}}+\frac{\mathrm{1}}{\mathrm{2}}\tau \left({\mathit{q}}^{\mathit{m}}+{\mathit{q}}^{\mathit{m}+\mathrm{1}}\right)$$

$$\mathit{S}=\frac{a}{\mathrm{2}h}\mathit{J}+\frac{\mathrm{1}}{h^{\alpha }}\mathit{B}\cdot \mathit{G}\cdot \tilde{\mathit{F}}+\frac{\mathrm{1}}{h^{\alpha }}\mathit{C}\cdot {\mathit{G}}^{\mathit{*}}\cdot \tilde{\mathit{H}}-d\mathit{I}$$

同时， 对于KoBol模型， 可以将方程$(\mathrm{1})$中的系数写成以下格式：

$$\left\{\begin{array}{c}b\left(x\right)=k_{\mathrm{1}}{\mathrm{e}}^{-\lambda x},\tilde{f}\left(x\right)={\mathrm{e}}^{\lambda x}\\ c\left(x\right)=k_{\mathrm{2}}{\mathrm{e}}^{\lambda x},\tilde{h}\left(x\right)={\mathrm{e}}^{-\lambda x}\end{array}\right.$$

那么对于矩阵$\mathit{B}$、 $\tilde{\mathit{F}}$， 就有如下形式：

$$\left\{\begin{array}{ll}\mathit{B}=k_{\mathrm{1}}\cdot \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathrm{e}}^{-\lambda x_{\mathrm{1}}},{\mathrm{e}}^{-\lambda x_{\mathrm{2}}},\dots ,{\mathrm{e}}^{-\lambda x_{N-\mathrm{1}}})& \\ \tilde{\mathit{F}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathrm{e}}^{\lambda x_{\mathrm{1}}},{\mathrm{e}}^{\lambda x_{\mathrm{2}}},\dots ,{\mathrm{e}}^{\lambda x_{N-\mathrm{1}}})& \end{array}\right.$$

所以可以得到如下形式：

$$\mathit{B}\cdot \mathit{G}\cdot \tilde{\mathit{F}}=k_{\mathrm{1}}\cdot \left(\begin{array}{cccccc}{g}_{\mathrm{1}}^{\left(\alpha \right)}& g_{\mathrm{0}}^{\left(\alpha \right)}{\mathrm{e}}^{\mathrm{\lambda }h}& \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}\\ g_{\mathrm{2}}^{\left(\alpha \right)}{\mathrm{e}}^{-\mathrm{\lambda }h}& g_{\mathrm{1}}^{\left(\alpha \right)}& g_{\mathrm{0}}^{\left(\alpha \right)}{\mathrm{e}}^{\mathrm{\lambda }h}& \mathrm{0}& \cdots & \mathrm{0}\\ g_{\mathrm{3}}^{\left(\alpha \right)}{\mathrm{e}}^{-\mathrm{2}\mathrm{\lambda }h}& g_{\mathrm{2}}^{\left(\alpha \right)}{\mathrm{e}}^{-\mathrm{\lambda }h}& g_{\mathrm{1}}^{\left(\alpha \right)}& \ddots & \ddots & ⋮\\ ⋮& ⋮& \ddots & \ddots & \ddots & \mathrm{0}\\ g_{N-\mathrm{2}}^{\left(\alpha \right)}{\mathrm{e}}^{-\mathrm{\lambda }\left(N-\mathrm{3}\right)h}& g_{N-\mathrm{3}}^{\left(\alpha \right)}{\mathrm{e}}^{-\mathrm{\lambda }\left(N-\mathrm{4}\right)h}& \cdots & \ddots & g_{\mathrm{1}}^{\left(\alpha \right)}& g_{\mathrm{0}}^{\left(\alpha \right)}{\mathrm{e}}^{\mathrm{\lambda }h}\\ g_{N-\mathrm{1}}^{\left(\alpha \right)}{\mathrm{e}}^{-\mathrm{\lambda }\left(N-\mathrm{2}\right)h}& g_{N-\mathrm{2}}^{\left(\alpha \right)}{\mathrm{e}}^{-\mathrm{\lambda }\left(N-\mathrm{3}\right)h}& g_{N-\mathrm{3}}^{\left(\alpha \right)}{\mathrm{e}}^{-\mathrm{\lambda }\left(N-\mathrm{4}\right)h}& \cdots & g_{\mathrm{2}}^{\left(\alpha \right)}{\mathrm{e}}^{-\mathrm{\lambda }h}& g_{\mathrm{1}}^{\left(\alpha \right)}\end{array}\right)$$

可以发现矩阵$\mathit{B}\cdot \mathit{G}\cdot \tilde{\mathit{F}}$为Toeplitz矩阵， 同理$\mathit{C}\cdot {\mathit{G}}^{\mathit{*}}\cdot \tilde{\mathit{H}}$也为Toeplitz矩阵， 从而矩阵S也同样为Toeplitz矩阵， 可以证明在Crank⁃Nicolson格式中的迭代矩阵也为Toeplitz矩阵。

可以采用循环预处理子来加速子空间方法^［27］， 如Strang循环预处理子^［28］与Chan循环预处理子^［29］。记${\mathit{T}}_{\mathit{n}}$为一个$n\times n$的Toeplitz矩阵， 那么Strang循环预处理矩阵$s({\mathit{T}}_{\mathit{n}})$是一个与${\mathit{T}}_{\mathit{n}}$阶数相同的Toeplitz矩阵， 其元素可由长度为$\mathrm{2}n-\mathrm{1}$的序列$s_k$所决定。 其中

$$s_k=\left\{\begin{array}{cc}{t}_k,& \mathrm{0}\le k<n/\mathrm{2}\\ \mathrm{0},& \mathrm{若}n\mathrm{为}\mathrm{偶}\mathrm{数}\mathrm{且}k=n/\mathrm{2}\\ t_{k-n},& n/\mathrm{2}<k<n\\ s_{k+n},& -n<k<\mathrm{0}\end{array}\right.$$

类似地，Chan循环预处理矩阵$c({\mathit{T}}_{\mathit{n}})$也是一个与${\mathit{T}}_{\mathit{n}}$阶数相同的Toeplitz矩阵， 其元素可由长度为$\mathrm{2}n-\mathrm{1}$的序列$s_k$所决定。 其中

$$c_k=\left\{\begin{array}{ll}\frac{(n-k)t_k+k{t}_{k-n}}{n},& \mathrm{0}\le k\le n-\mathrm{1}\\ c_{n+k},& \mathrm{1}-n\le k<\mathrm{0}\end{array}\right.$$

在后面的数值实验中将该预处理技术应用于GMRES、 BiCGSTAB和CGNR算法上并比较计算效果。

3 数值实验

首先对一个带精确解的双边分数阶扩散方程来验证该离散格式的精度与收敛阶。然后， 在KoBol模型下对欧式看涨期权进行定价。

例1   考虑如下终值问题：

$$\left\{\begin{array}{cc}\begin{array}{l}\frac{\partial u(x,t)}{\partial t}+c_{+}(x,t)\frac{{\partial }^{\alpha }u(x,t)}{{\partial }_{+}{x}^{\alpha }}+\\ c_{-}(x,t)\frac{{\partial }^{\alpha }u(x,t)}{{\partial }_{-x}{x}^{\alpha }}+s(x,t)=\mathrm{0}\end{array}& \\ \begin{array}{l}{c}_{+}\left(x,t\right)=\Gamma \left(\mathrm{2}-\alpha \right)x^{\alpha }\\ c_{-}(x,t)=\Gamma (\mathrm{2}-\alpha ){(\mathrm{2}-x)}^{\alpha }\end{array}& \\ \begin{array}{l}s(x,t)=-\mathrm{32}{\mathrm{e}}^{t-T}\left[\frac{x^{\mathrm{2}}+{\left(\mathrm{2}-x\right)}^{\mathrm{2}}}{\mathrm{2}-\alpha }-\right.\\ \frac{\mathrm{3}(x^{\mathrm{3}}+{(\mathrm{2}-x)}^{\mathrm{3}})}{(\mathrm{2}-\alpha )(\mathrm{3}-\alpha )}+\\ \left.\frac{\mathrm{3}(x^{\mathrm{4}}+{(\mathrm{2}-x)}^{\mathrm{4}})}{(\mathrm{2}-\alpha )(\mathrm{3}-\alpha )(\mathrm{4}-\alpha )}\right]-\\ \mathrm{4}{\mathrm{e}}^{t-T}{x}^{\mathrm{2}}{(\mathrm{2}-x)}^{\mathrm{2}}\end{array}& \\ u(x,T)=\mathrm{4}{x}^{\mathrm{2}}{(\mathrm{2}-x)}^{\mathrm{2}}& \\ u(\mathrm{0},t)=u(\mathrm{2},t)=\mathrm{0}& \\ \mathrm{0}\le x\le \mathrm{2},\mathrm{0}\le t\le T,\alpha =\mathrm{1.8},T=\mathrm{1}& \end{array}\right.$$

其中该方程的精确解为$U\left(x,t\right)=\mathrm{4}{\mathrm{e}}^{t-T}{x}^{\mathrm{2}}{\left(\mathrm{2}-x\right)}^{\mathrm{2}}$。

根据前文提到的离散格式对上述方程进行差分离散， 并将例1中的系数代入式$(\mathrm{7})$验证， 发现此时Crank⁃Nicolson格式无条件稳定。取$h=\frac{\mathrm{1}}{\mathrm{20}}$， $\tau =\frac{\mathrm{1}}{\mathrm{40}}$， $t=\mathrm{0}$时的数值解与真实解如图1， 并记精度误差为$E_{\mathrm{\infty }}=\underset{i,j}{\mathrm{m}\mathrm{a}\mathrm{x}}\{|u(i,j)-U(i,j)|\}$， 计算结果如表1。

图1 例1中数值解与真实解的比较(t=0)

Fig. 1 Comparison of exact solution with numerical solution for Example 1 (t=0)

表1 例1中Crank⁃Nicolson格式计算精度${\mathit{E}}_{\mathrm{\infty }}$与收敛阶

Tab. 1 Numerical accuracy ${\mathit{E}}_{\mathrm{\infty }}$ and convergence order of Crank⁃Nicolson scheme for Example 1

N	M	$E_{\mathrm{\infty }}$	$\mathrm{l}\mathrm{o}{\mathrm{g}}_2\frac{E_{\mathrm{\infty }}^i}{E_{\mathrm{\infty }}^{i+1}}$
$2^6$	$2^4$	0.013 7
$2^7$	$2^5$	0.005 7	1.264 9
$2^8$	$2^6$	0.002 7	1.055 7
$2^9$	$2^7$	0.001 3	1.009 6
$2^{10}$	$2^8$	0.000 7	1.004 7
$2^{11}$	$2^9$	0.000 3	1.002 2

例2   考虑如下欧式看涨期权KoBol模型：

$$\left\{\begin{array}{c}\begin{array}{c}\begin{array}{l}\frac{\partial V\left(x,t\right)}{\partial t}+\left(r-v-{\lambda }^{\alpha }\left(q-p\right)\right)\frac{\partial V\left(x,t\right)}{\partial x}+\\ \frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\alpha }\left[p{\mathrm{e}}^{\lambda x}\frac{{\partial }^{\alpha }\left({\mathrm{e}}^{-\lambda x}V(x,t)\right)}{{\partial }_{-}{x}^{\alpha }}+\right.\\ \left.q{\mathrm{e}}^{-\lambda x}\frac{{\partial }^{\alpha }\left({\mathrm{e}}^{\lambda x}V(x,t)\right)}{{\partial }_{+}{x}^{\alpha }}\right]=\left(r+\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\alpha }{\lambda }^{\alpha }\right)V(x,t)\end{array}\\ V(x,T)=\mathrm{m}\mathrm{a}\mathrm{x}\{{\mathrm{e}}^x-K,\mathrm{0}\},L<x<R\end{array}\\ V(L,t)=\mathrm{0,0}\le t<T\\ V(R,t)={\mathrm{e}}^R-K{\mathrm{e}}^{-r(T-t)},\mathrm{0}\le t<T\end{array}\right.$$

其中参数选取分别为：$r=\mathrm{0.05}$， $T=\mathrm{0.5}$， $p=\mathrm{0.5}$， $q=\mathrm{0.5}$， $\sigma =\mathrm{0.2}$， $\lambda =\mathrm{3}$， $\alpha =\mathrm{1.9}$， $K=\mathrm{80}$， $v=\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\alpha }[p{(\lambda -\mathrm{1})}^{\alpha }+q{(\lambda +\mathrm{1})}^{\alpha }-{\lambda }^{\alpha }-\alpha {\lambda }^{\alpha -\mathrm{1}}(q-p)]$。

例2的解曲面如图2a。将相同的期权参数利用Black⁃Scholes公式进行定价， 将得出的期权价格与Kobol模型下的做差比较如图2b， 其价格在一定程度上体现了列维分布与正态分布相比具有尖峰肥尾的特点。在期权交易中， 将期权分为实值期权、 平值期权与虚值期权， 那么， 尖峰意味着在平值附近其概率分布更靠近现价， 从而在相同波动率下其理论收益略低于Black⁃Scholes模型， 所以价格也略低。而肥尾意味着在深度虚值与深度实值部分的理论收益高于正态分布的估计， 其价格略高于Black⁃Scholes模型的价格， 这也与实验结果观测一致， 说明该模型更接近于实际分布。

图2 例2欧式看涨期权在KoBol模型下的解

Fig. 2 Solution for a European call option in KoBol model

将Strang和T.Chan循环预处理子用在GMRES、BiCGSTAB、 CGNR算法上并与未经预处理的算法作比较， 在计算例2的同时， 记录每个时间层上求解线性方程组的迭代步数IT与计算时间并取平均， 其中最大迭代步数$\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{I}\mathrm{T}=\mathrm{10}\mathrm{000}$， 停止准则为$||r_k||/||r_{\mathrm{0}}||\le \mathrm{1}\times {\mathrm{10}}^{-\mathrm{7}}$， 实验结果见表2。

表2 例2采用预处理子空间方法的计算结果对比

Tab. 2 Numerical results of preconditioned subspace methods for Example 2

算法	$N=2^{10}$	$M=2^5$	$N=2^{11}$	$M=2^6$	$N=2^{12}$	$M=2^7$
	迭代步数	时间/s	迭代步数	时间/s	迭代步数	时间/s
GMRES	50.09	0.018 3	60.59	0.116 2	71.28	0.491 4
GMRES⁃Strang	5.00	0.002 3	5.00	0.011 4	5.05	0.042 7
GMRES⁃T.Chan	6.00	0.002 6	6.00	0.013 3	6.00	0.049 3
BiCGSTAB	40.16	0.103 5	46.56	0.389 3	53.63	1.656 4
BiCGSTAB⁃Strang	4.00	0.010 9	4.00	0.041 6	4.01	0.154 4
BiCGSTAB⁃T. Chan	4.00	0.011 8	4.00	0.041 2	4.01	0.153 9
CGNR	581.63	0.272 5	940.44	2.459 9	1 482.13	15.660 9
CGNR⁃Strang	8.00	0.004 0	8.00	0.020 9	7.99	0.087 2
CGNR⁃T.Chan	9.00	0.004 3	9.00	0.023 5	9.00	0.096 6

实验结果表明， 无论是迭代步数还是计算时间， 使用预处理技术之后的计算效率都明显优于未经预处理的算法。同时， 可以发现相同条件下预处理GMRES算法是所有算法中计算时间最快的方法。预处理技术之所以可以降低子空间算法的迭代步数， 原因在于预处理之后系数矩阵的特征值较为聚集。将例2中$N={\mathrm{2}}^{\mathrm{10}}$时， 经两类预处理子预处理之后的系数矩阵特征值绘制如图3所示。可以发现， 经过预处理之后的特征值确实较为聚集。

图3 例2经预处理后的系数矩阵特征值分布

Fig. 3 Eigenvalue of preconditioned coefficient matrix

4 基于股指期权的实证研究

通过对中国股票市场进行研究发现， 在国内股票市场也存在尖峰肥尾的现象。图4为上证50指数与沪深300指数在2019年1月至12月每日对数收益率的归一直方图与对应的正态分布概率密度函数图， 可以发现在中国股票市场中， 这一特点十分明显。所以Black⁃Scholes模型在我国市场中进行定价时会有较大误差。接下来验证KoBol模型可以更好地描述国内期权市场。

图4 2019年上证50指数与沪深300指数的尖峰肥尾现象

Fig. 4 Skewed and fat tailed phenomenon of SSE 50 and CSI 300 Index in 2019

在国内的场内期权交易市场中， 交易最为活跃、成交量最大的是当月合约， 其买卖价差也更为接近， 所以考虑上交所2020年1月3日收盘时， 1月22日到期的上证50ETF与沪深300ETF场内期权收盘数据。在实际交易中， 市场上更多地采用Black⁃Scholes模型进行定价， 将比较KoBol模型及Black⁃Scholes模型的定价与市场价格之间的差距。对KoBol与Black⁃Scholes（BS）模型选取的参数如表3， 标的现价$S_{\mathrm{0}}$与时间$T$均为市场收盘数据， 无风险收益率$r$为1年期国债收益率。采用过去180 d的年化历史波动率作为定价模型中的波动率， 年化天数为252 d， 经计算得50ETF为$\sigma =\mathrm{16.74}\mathrm{\%}$， 300ETF为$\sigma =\mathrm{17.29}\mathrm{\%}$。

表3 KoBol模型与Black⁃Scholes模型参数选取

Tab. 3 Parameter selection for Kobol model and ( BS model )

参数	上证50ETF			沪深300ETF
	BS	KoBol		BS	KoBol
r/%	2.435 9	2.435 9		2.435 9	2.435 9
T	20/365	20/365		20/365	20/365
$S_0$	3.068	3.068		4.139	4.139
$\alpha$	—	1.86		—	1.86
$\lambda$	—	1.0		—	1.8
p	—	0.507		—	0.51

对于Black⁃Scholes模型直接利用欧式期权的解析解进行定价。用$V_{\mathrm{M}}$表示市场实际价格， 用均方误差$E_{\mathrm{1}}$、相对误差$E_{\mathrm{2}}$与最大误差$E_{\mathrm{\infty }}$来衡量误差大小， 分别有如下形式：

$$\begin{array}{l}{E}_{\mathrm{1}}=\frac{{\sum }_{i=\mathrm{1}}^n(V^i-V_{\mathrm{M}}^i{)}^{\mathrm{2}}}{n}\\ E_{\mathrm{2}}=\frac{{{\sum }_{i=\mathrm{1}}^n\left(\frac{V^i-V_{\mathrm{M}}^i}{V_{\mathrm{M}}^i}\right)}^{\mathrm{2}}}{n}\\ E_{\mathrm{\infty }}=\mathrm{m}\mathrm{a}\mathrm{x}\left(\left|V^i-V_{\mathrm{M}}^i\right|\right)\end{array}$$

两种模型的定价结果如表4与表5， 与市场价格作差取绝对值之后的结果见表6与表7。可以看出KoBol模型下的误差相比Black⁃Scholes模型都要更小， 说明KoBol模型下的价格更接近市场价格。

表4 50ETF场内期权两种模型定价

看涨期权/元			行权价K/元	看跌期权/元
$V_{\mathrm{M}}$	$V_{\mathrm{B}\mathrm{S}}$	$V_{\mathrm{K}\mathrm{o}\mathrm{B}\mathrm{o}l}$		$V_{\mathrm{M}}$	$V_{\mathrm{B}\mathrm{S}}$	$V_{\mathrm{K}\mathrm{o}\mathrm{B}\mathrm{o}\mathrm{l}}$
0.327 1	0.321 7	0.324 6	2.75	0.000 6	0.000 08	0.000 8
0.273 8	0.272 1	0.275 0	2.80	0.000 9	0.000 30	0.001 1
0.224 0	0.223 0	0.225 8	2.85	0.001 4	0.001 20	0.001 8
0.176 8	0.175 5	0.177 5	2.90	0.002 4	0.003 70	0.003 4
0.128 6	0.131 1	0.131 6	2.95	0.005 1	0.009 20	0.007 5
0.085 9	0.091 8	0.090 4	3.00	0.012 2	0.019 80	0.016 2
0.028 1	0.035 5	0.032 4	3.10	0.053 5	0.063 40	0.058 0
0.006 9	0.009 4	0.008 6	3.20	0.132 1	0.137 20	0.134 2
0.001 8	0.001 7	0.002 7	3.30	0.228 7	0.229 30	0.228 1
0.001 1	0.000 2	0.001 4	3.40	0.327 2	0.327 70	0.326 7

表5 300ETF场内期权两种模型定价

Tab. 5 Option prices of 300ETF in two models

看涨期权/元			行权价K/元	看跌期权/元
$V_{\mathrm{M}}$	$V_{\mathrm{B}\mathrm{S}}$	$V_{\mathrm{K}\mathrm{o}\mathrm{B}\mathrm{o}\mathrm{l}}$		$V_{\mathrm{M}}$	$V_{\mathrm{B}\mathrm{S}}$	$V_{\mathrm{K}\mathrm{o}\mathrm{B}\mathrm{o}\mathrm{l}}$
0.350 7	0.345 0	0.352 8	3.80	0.001 3	0.000 9	0.001 5
0.253 3	0.248 9	0.255 4	3.90	0.004 1	0.004 7	0.004 0
0.162 1	0.161 6	0.164 3	4.00	0.011 1	0.017 2	0.012 7
0.086 2	0.091 1	0.089 3	4.10	0.034 9	0.046 6	0.037 6
0.038 8	0.043 2	0.039 8	4.20	0.087 0	0.098 6	0.088 0
0.015 5	0.016 8	0.015 0	4.30	0.165 0	0.172 1	0.163 0
0.006 1	0.005 3	0.005 6	4.40	0.256 0	0.260 4	0.253 5
0.002 6	0.001 4	0.002 6	4.50	0.349 6	0.356 4	0.350 4

表6 50ETF场内期权两种模型与市场价格比较

Tab. 6 Comparisons of two models with market price of 50ETF options

看涨期权/元		行权价K/元	看跌期权/元
$|V_{\mathrm{M}}-V_{\mathrm{B}\mathrm{S}}|$	$|V_{\mathrm{M}}-V_{\mathrm{K}\mathrm{o}\mathrm{B}\mathrm{o}\mathrm{l}}|$		$|V_{\mathrm{M}}-V_{\mathrm{B}\mathrm{S}}|$	$|V_{\mathrm{M}}-V_{\mathrm{K}\mathrm{o}\mathrm{B}\mathrm{o}\mathrm{l}}|$
0.005 4	0.002 5	2.75	0.000 5	0.000 2
0.001 7	0.001 2	2.80	0.000 5	0.000 2
0.001 0	0.001 8	2.85	0.000 2	0.000 4
0.001 3	0.000 7	2.90	0.001 3	0.001 0
0.002 5	0.003 0	2.95	0.004 1	0.002 4
0.005 9	0.004 5	3.00	0.007 6	0.004 0
0.007 4	0.004 3	3.10	0.009 9	0.004 5
0.002 5	0.001 7	3.20	0.005 1	0.002 1
0.000 1	0.000 9	3.30	0.000 6	0.000 6
0.000 9	0.000 3	3.40	0.000 5	0.000 5
0.138 2	0.062 4	$E_1({10}^{-4})$	0.201 2	0.048 5
0.089 8	0.041 7	$E_2$	0.248 9	0.069 1
0.007 4	0.004 5	$E_{\mathrm{\infty }}$	0.009 9	0.004 5

表7 300ETF场内期权两种模型与市场价格比较

Tab. 7 Comparisons of two models with market price of 300ETF options

看涨期权/元		行权价K/元	看跌期权/元
$|V_{\mathrm{M}}-V_{\mathrm{B}\mathrm{S}}|$	$|V_{\mathrm{M}}-V_{\mathrm{K}\mathrm{o}\mathrm{B}\mathrm{o}\mathrm{l}}|$		$|V_{\mathrm{M}}-V_{\mathrm{B}\mathrm{S}}|$	$|V_{\mathrm{M}}-V_{\mathrm{K}\mathrm{o}\mathrm{B}\mathrm{o}\mathrm{l}}|$
0.005 7	0.002 10	3.8	0.000 4	0.000 2
0.004 4	0.002 10	3.9	0.000 6	0.000 1
0.000 5	0.002 20	4.0	0.006 1	0.001 6
0.004 9	0.003 10	4.1	0.011 7	0.002 7
0.004 4	0.001 00	4.2	0.011 6	0.000 9
0.001 3	0.000 50	4.3	0.007 1	0.002 0
0.000 7	0.000 50	4.4	0.004 4	0.002 5
0.001 2	0.000 02	4.5	0.006 8	0.000 8
0.123 2	0.030 70	$E_1({10}^{-4})$	0.530 2	0.026 9
0.033 7	0.001 20	$E_2$	0.068 6	0.007 4
0.005 7	0.003 10	$E_{\mathrm{\infty }}$	0.011 7	0.002 7

因场内交易为竞价交易， 考虑到实值期权在临近到期日时行权风险上升， 所以其买卖价差较大且成交量有限， 市场作用接近于期货， 此时其成交价格并不能精确体现期权实际价格， 而虚值期权相对来说买卖价差更小， 且不存在套利空间， 所以其价格能较好地反映期权的实际价格。因此考虑计算虚值期权的3种误差， 得到结果如表8至表10， 可以看出，与Black⁃Scholes模型相比， KoBol模型的定价更贴近市场实际价格。

表8 50ETF与300ETF虚值期权的均方误差E₁

Tab. 8 Error ${\mathit{E}}_{\mathrm{1}}$ of 50ETF and 300ETF OTM options

$E_1/{10}^{-4}$	50ETF			300ETF
	看涨	看跌		看涨	看跌
BS	0.156 0	0.128 6		0.057 5	0.436 8
KoBol	0.055 2	0.038 2		0.003 6	0.025 2

表9 50ETF与300ETF虚值期权的相对误差E₂

Tab. 9 Error ${\mathit{E}}_{\mathrm{2}}$ of 50ETF and 300ETF OTM options

$E_2$	50ETF			300ETF
	看涨	看跌		看涨	看跌
BS	0.223 2	0.408 8		0.066 4	0.132 0
KoBol	0.103 4	0.114 0		0.002 1	0.014 8

表10 50ETF与300ETF虚值期权的最大误差$E_{\mathrm{\infty }}$

Tab. 10 Error ${\mathit{E}}_{\mathrm{3}}$ of 50ETF and 300ETF OTM options

$E_{\mathrm{\infty }}$	50ETF			300ETF
	看涨	看跌		看涨	看跌
BS	0.007 4	0.007 6		0.004 4	0.011 7
KoBol	0.004 3	0.004 0		0.001 0	0.002 7

在前文提到了Black⁃Scholes模型不能完美地描述实际市场， 主要是因为Black⁃Scholes模型假设在市场风险中性下标的资产的预期收益分布服从正态分布， 即隐含收益分布与市场不符。这样的不足之处导致了有波动率微笑的存在， 所以在利用Black⁃Scholes模型进行场外期权定价时， 需要根据期权的期限与行权价来调整波动率， 从而会产生相应的误差， 如果能够将不同行权价期权的隐含波动率控制在更小的范围甚至接近一致， 那么对于期权的定价则具有重要意义。

接下来反演计算隐含波动率微笑曲线， KoBol模型和Black⁃Scholes模型的参数选取依然如表3， 考虑到波动率与期权价格的单调关系以及期权的内在价值， 利用二分法反推虚值期权的隐含波动率。 选取虚值期权， 即行权价高于标的价格的看涨期权与行权价低于标的价格的看跌期权数据进行计算后， KoBol模型和Black⁃Scholes模型下不同行权价对应的虚值期权隐含波动率曲线即为图5， 横轴为虚值期权行权价， 纵轴为期权的隐含波动率。表11给出了虚值期权在KoBol模型和Black⁃Scholes模型下隐含波动率的具体数据。

可以从图5看出，无论是50ETF还是300ETF，场内虚值期权的隐含波动率在KoBol模型下被控制在了一个更小的范围内，且波动率曲线更为平稳。由表11也可以看出，KoBol模型下的隐含波动率标准差更小。这都说明了KoBol模型在此参数下更好地描述了该时刻标的资产的隐含收益分布。

表11 2020年1月3日股指期权隐含波动率

Tab. 11 Implied volatility data of index options on 2020/1/3

上证50ETF				沪深300ETF
行权价 K/元	BS/%	KoBol/ %		行权价K/元	BS/%	KoBol/ %
2.75	21.01	15.11		3.80	18.21	16.33
2.80	19.07	15.51		3.90	16.75	17.38
2.85	17.11	15.46		4.00	14.76	16.39
2.90	15.27	14.84		4.10	14.11	16.46
2.95	13.99	14.70		4.20	16.10	17.05
3.00	13.34	14.67		4.30	16.77	17.56
3.10	14.08	15.12		4.40	17.86	17.69
3.20	15.15	15.26		4.50	19.25	17.18
3.30	16.98	14.42
3.40	20.90	15.06
标准差	2.68	0.33			1.61	0.51

图5 2020年1月3日股指期权隐含波动率图

Fig. 5 Implied volatility of index options on 2020/1/3

5 结语

本文利用带位移的Grünwald差分格式对一类基于列维过程的分数阶期权定价模型进行了离散，分析了数值格式的稳定性条件， 采用预处理Krylov子空间方法求解对应的线性代数方程组，数值结果验证了模型与算法的有效性。同时， 针对国内股指期权的实际交易数据，利用KoBol分数阶模型对股指期权进行定价并反演计算波动率微笑曲线， 通过实证分析说明该模型比Black⁃Scholes模型有更好的效果。