本文研究求解一类时间分数阶B-S期权定价模型的MQ拟插值新方法.

Black-Scholes期权定价理论是期权方法的基石, 这一突破性的成果对计算金融衍生品内在价值具有非常重要的指导作用^[2]. B-S期权定价模型基于标的资产符合几何布朗运动的假设, 也就是在效率市场假说成立的条件下才能够实现, 但真实的市场上标的物价格往往不能完全反映历史信息, 即效率市场假说并不成立. 为了解决这个问题, 学者将经典的B-S方程推广到分数阶, 分数阶模型具有记忆性, 能够度量期权价值对价格变化的历史信息的依赖程度^[6].

在分数阶B-S期权定价理论方面, 许多学者做了大量工作, 例如文献[3, 5] 提出了空间分数阶模型, 文献[4, 15] 发展了时间分数阶模型, 文献[11–12] 则研究时空分数阶模型.然而, 各种分数阶B-S模型的解析解都很复杂, 一般很难求出, 因此需要发展数值解. 许多数值求解分数阶B-S方程的算法已经被发展起来了, 例如文献[4, 16] 给出了一些关于分数阶B-S模型的数值算法, 但是这些算法都是基于有限差分等网格类方法, 对网格具有很强的依赖性, 难以处理期权价格因在行权价附近变化剧烈、需要维持较多的配置点以达到准确模拟的要求.因此, 科研界亟需探索构造无网格的求解分数阶B-S模型的新途径.

MQ拟插值是一种具有保形性的高精度无网格逼近方法, 可以求解复杂的不规则区域、移动边界和高维问题. 而且MQ函数是无穷连续的, 便于求解函数的导数. 此外, 该方法不需要求解大尺度矩阵, 可以避免过大条件数问题.因而, MQ拟插值方法作为无网格算法的一类, 在科学和工程中得到了广泛的应用^{[1, 10, 13-14, 17]}. Beatson^[1]甚至将MQ拟插值用在电影《指环王3》的道具设计上, 以增强动画效果.

本文将借助三阶MQ拟插值, 构造求解时间分数阶B-S模型的新算法, 算法不仅是无网格的, 而且计算稳定、具有较高精度.

给定序列

$ {\bf x}=\cdots <x_{i-1}<x_{i}<x_{i+1}<\cdots , \quad \, h\triangleq \max\limits_{i}(x_i-x_{i-1}) . $

其中当$ i\longrightarrow\pm\infty $时, $ x_{\pm i}\longrightarrow\pm\infty $. 则函数$ f(x) $的三阶MQ拟插值具有形式

(1.1) $ \begin{equation} ({\cal L}f)=\sum f(x_i)\Psi_i(x). \end{equation} $

其中

$ \Psi_{i}(x)=\Psi(x-x_i)=\frac{\psi_{i+1}-\psi_i(x)} {2(x_{i+2}-x_{i-1})}-\frac{\psi_i(x)-\psi_{i-1}(x)}{2(x_{i+1}-x_{i-2})}, $

这里$ \psi_i(x) $是三阶MQ函数的线性组合, 具有如下形式

$ \psi_i(x)=\psi(x-x_i)=\frac{(\phi_{i+1}-\phi_i)/(x_{i+1}-x_i)-(\phi_i-\phi_{i-1})/(x_i-x_{i-1})}{x_{i+1}-x_{i-1}}, $

这里$ \phi(x)=\sqrt{(x^2+c^2)^3} $ ($ c $是形状参数), $ \phi_i=\phi(x-x_i) $. 如文献[1]中所述, 上述拟插值具有保形性和高阶逼近的特性. 矩阵$ {{\bf \Psi}}=\{\Psi_s(x_i)\} $, $ {{\bf \Psi}_1}=\{\frac{\partial \Psi_s(x_i)}{\partial x}\} $和$ {{\bf \Psi}_2}=\{\frac{\partial^2\Psi_s(x_i)}{\partial x^2}\} $是有限带宽矩阵, 因而在数值计算中是高效的.

本文的结构如下. 在第二节, 我们构造了具有空间$ \ell $阶和时间$ 2-\alpha $阶精度的离散格式. 在第三节中, 研究了该格式的稳定性和收敛性. 第四节给出数值结果.

本文考虑如下B-S模型^{[6, 12, 16]}

(2.1) $ \begin{equation} \left\{\begin{array} {l} { }\frac{\partial^\alpha V(S, t)}{\partial t^\alpha}+\frac{1}{2}\sigma^2S^2\frac{\partial^2V(S, t)}{\partial s^2}+rS\frac{\partial V(S, t)}{\partial S}-rV(S, t)=0, (S, t)\in(0, +\infty)\times(0, T), \\ V(0, t)=p(t), V(\infty, t)=q(t), \\ V(S, T)=v(s), \\ \end{array} \right. \end{equation} $

其中$ V(S, t) $, $ T $, $ r $和$ \sigma $分别是欧式期权的价格, 到期时间, 无风险利率和波动率. 这里$ 0<\alpha\leq1 $, 且有

(2.2) $ \begin{equation} \frac{\partial^\alpha V(S, t)}{\partial t^\alpha}= \left\{\begin{array}{ll} { }\frac{1}{\Gamma(1-\alpha)}\frac{\rm d}{{\rm d}t}\int_t^T\frac{V(S, \xi)-V(S, T)}{(\xi-t)^{\alpha}}{\rm d}\xi, \quad 0<\alpha<1, \\ { }\frac{\partial V(S, t)}{\partial t}, \quad \quad\quad \quad\quad \quad\quad \quad\quad\quad\alpha=1. \end{array}\right. \end{equation} $

令$ \tau=T-t $, $ x=\ln S $, $ U(x, \tau)=V(e^x, T-t) $, 并将原无界区域截断为有限区间$ (L_1, L_2) $, 则(2.1)式转换为

(2.3) $ \begin{equation} \left\{\begin{array} {l} { }{}^C_{0}D_{\tau}^{\alpha}U(x, \tau)=a\frac{\partial^2U(x, \tau)}{\partial x^2}+b\frac{\partial U(x, \tau)}{\partial x}-rU(x, \tau)+f(x, \tau), (x, \tau)\in(L_1, L_2)\times(0, T), \\ U(L_1, \tau)=p(\tau), U(L_2, \tau)=q(\tau), \\ U(x, 0)=u(x), \end{array} \right. \end{equation} $

这里$ ^C_{0}D_{\tau}^{\alpha} $是Caputo导数, $ a=\frac{1}{2}\sigma^2>0 $, $ b=r-a $, $ f(x, \tau) $是一个附加项(方便计算起见).

将区间$ (0, T) $分为$ N $个子区间, 其中$ \Delta t=T/N $, $ \tau_j=j\Delta t $$ (j=0, 1, 2, \cdots , N) $. 根据文献[9]中的经典$ L1 $离散式, 可以得到

(2.4) $ \begin{equation} ^C_{0}D_{\tau}^{\alpha}U(x_i, \tau_{k+1})=\frac{\Delta t^{-\alpha}}{\Gamma(2-\alpha)}\sum\limits^k_{j=0}\Big(U(x_i, \tau_{k+1-j})-U(x_i, \tau_{k-j}) \Big)b_j+{\cal O}(\Delta t^{2-\alpha}). \end{equation} $

其中$ b_j=(j+1)^{1-\alpha}-j^{1-\alpha} $. 取数据点

$ L_1=x_0<\cdots <x_{i-1}<x_i<x_{i+1}<\cdots <x_M=L_2. $

根据三阶MQ拟插值(0.1)式, 有

$ U(x_i, \tau_{k+1})\approx({\cal L}U)(x_i, \tau_{k+1})=\sum\limits_sU(x_s, \tau_{k+1})\Psi_s(x_i), i=0, 1, \cdots , M. $

进一步, 根据文献[13], 函数的一阶和二阶导数具有逼近形式

$ \begin{eqnarray*} U_x(x_i, \tau_{k+1})&\approx&\sum\limits_sU(x_i, \tau_{k+1})\Psi_s'(x_i)+{\cal O}(h^\ell), U_{xx}(x_i, \tau_{k+1})\\ &\approx&\sum\limits_sU(x_i, \tau_{k+1})\Psi_s''(x_i)+{\cal O}(h^\ell). \end{eqnarray*} $

将上述$ L1 $逼近式和三阶MQ拟插值的逼近式代入(2.3)式, 可得

(2.5) $ \begin{eqnarray} &&\frac{\Delta t^{-\alpha}}{\Gamma(2-\alpha)}\sum\limits^k_{j=0}(U^{k+1-j}_i-U^{k-j}_i)b_j {}\\ & =& a\sum\limits_sU^{k+1}_s\Psi''_s(x_i)+b\sum\limits_sU^{k+1}_s\Psi'_s(x_i)-rU^{k+1}_{i}+f^{k+1}_i+{\cal O}(\Delta t^{2-\alpha}+h^\ell), \end{eqnarray} $

其中$ U^k_i=U(x_i, \tau_k) $, $ f^k_i=f(x_i, \tau_k) $. 在(2.5)式的两边乘以$ d=\Gamma(2-\alpha)\Delta t^\alpha $, 有

$ \left\{\begin{array} {l} { } -ad\sum\limits_sU^1_s\Psi''_s(x_i)-bd\sum\limits_sU^1_s\Psi'_s(x_i)+(1+rd)U^1_i=U^0_i+df^1_i+{\cal O}(\Delta t^2+\Delta t^\alpha h^\ell), k=0, \\ { } -ad\sum\limits_sU^{k+1}_s\Psi''_s(x_i)-bd\sum\limits_sU^{k+1}_s\Psi'_s(x_i)+(1+rd)U^{k+1}_i\\ \ \ \ \ \ \ \ = \sum\limits^{k-1}_{j=0}(b_j-b_j+1)U^{k-j}_i+b_kU^0_i+df^{k+1}_i+{\cal O}(\Delta t^2+\Delta t^\alpha h^\ell), k\geq1. \end{array} \right. $

舍掉截断误差, 取$ \widetilde{U}^k_i $作为$ U^k_i $的逼近值, 得

(2.6) $ \begin{equation} \left\{\begin{array} {l} { } -ad\sum\limits_s\widetilde{U}^1_s\Psi''_s(x_i)-bd\sum\limits_s\widetilde{U}^1_s\Psi'_s(x_i)+(1+rd)\widetilde{U}^1_i=\widetilde{U}^0_i+df^1_i, k=0, \\ { } -ad\sum\limits_s\widetilde{U}^{k+1}_s\Psi''_s(x_i)-bd\sum\limits_s\widetilde{U}^{k+1}_s\Psi'_s(x_i)+(1+rd)\widetilde{U}^{k+1}_i\\ { } = \sum\limits^{k-1}_{j=0}(b_j-b_{j+1})\widetilde{U}^{k-j}_i+b_k\widetilde{U}^0_i+df^{k+1}_i, k\geq1. \end{array} \right. \end{equation} $

设$ {\bf A}=-ad{{\bf \Psi}}_2-bd{{\bf \Psi}}_1+(1+rd){\bf E} $, $ \widetilde{U}^k=[\widetilde{U}^k_0, \cdots , \widetilde{U}^k_M]^T $, $ u=[u_0, \cdots , u_M]^T $, $ f^k=[f_0, \cdots , f_M]^T $, $ G^1=[(1+rd)p^1-df^1_0-u_0, 0, \cdots , 0, (1+rd)q^1-df^1_M-u_M]^T $以及

$ \begin{eqnarray*} G^k&=& \bigg[(1+rd)p^k-df^k_0-\sum\limits^{k-2}_{j=0}(b_j-b_{j+1})p^{k-1-j}-b_{k-1}u_0, 0, \cdots , 0, \\ && (1+rd)q^k-df^k_M-\sum\limits^{k-2}_{j=0}(b_j-b_{j+1})q^{k-1-j}-b_{k-1}u_M\bigg]^T, \end{eqnarray*} $

这里$ {\bf E} $是单位矩阵,

$ {{\bf \Psi}}_1\triangleq\{\Psi'_{is}\}=\{\Psi'_s(x_i):2\leq i\leq M-1, 1\leq s\leq M, 0\quad elsewhere\}, $

以及

$ {{\bf \Psi}}_2\triangleq\{\Psi''_{is}\}=\{\Psi''_s(x_i):2\leq i\leq M-1, 1\leq s\leq M, 0\quad elsewhere\}. $

则, (2.6)式具有矩阵形式

(2.7) $ \begin{equation} \left\{\begin{array} {l} {\bf A}\widetilde{U}^1=\widetilde{U}^0+df^1+G^1, k=0, \\ { } {\bf A}\widetilde{U}^{k+1}=\sum\limits^{k-1}_{j=0}(b_j-b_{j+1})\widetilde{U}^{k-j}+b_k\widetilde{U}^0+df^{k+1}+G^{k+1}, k\geq1. \end{array} \right. \end{equation} $

设$ \widetilde{U}^k_i $为(2.6)式的逼近值, 令$ \varepsilon^k_i=\widetilde{U}^k_i-\widehat{U}^k_i $$ (i=0, 1, 2, \cdots , M, k=0, 1, 2, \cdots , N) $, 则

(3.1) $ \begin{equation} \left\{\begin{array} {l} { } -ad\sum\limits_s\varepsilon^1_s\Psi''_s(x_i)-bd\sum\limits_s\varepsilon^1_s\Psi'_s(x_i)+(1+rd)\varepsilon^1_i=\varepsilon^0_i+df^1_i, k=0, \\ { } -ad\sum\limits_s\varepsilon^{k+1}_s\Psi''_s(x_i)-bd\sum\limits_s\varepsilon^{k+1}_s\Psi'_s(x_i)+(1+rd)\varepsilon^{k+1}_i\\ \ \ \ \ \ \ \ \ =\sum\limits^{k-1}_{j=0}(b_j-b_{j+1})\varepsilon^{k-j}_i+b_k\varepsilon^0_i+df^{k+1}_i, k\geq1, \end{array} \right. \end{equation} $

这里$ \varepsilon^k_0=\varepsilon^k_M=0 $. 定义函数

$ \varepsilon^k(x)= \left\{\begin{array}{ll} { }\varepsilon^k_i, x_i-\frac{h}{2}<x\leq x_i+\frac{h}{2}, i=1, 2, \cdots , M-1, \\ { } 0, L_1\leq x\leq L_1+\frac{h}{2}\quad {\rm or} \quad L_2-\frac{h}{2}<x\leq L_2. \end{array}\right. $

设$ I^2=-1 $, $ L=L_1-L_2 $. 则

$ \varepsilon^k(x)=\sum\limits^\infty_{j=-\infty}\beta^k_je^{\frac{2j\pi x}{L}I}, $

这里$ \beta^k_j=\frac{1}{L}\int^L_0\varepsilon^k(x)e^{-\frac{2j\pi x}{L}I}{\rm d}x $是傅里叶展开式系数.

设$ \|\varepsilon^k\|=\sqrt{\sum\limits^{M-1}_{j=1}h_j|\varepsilon^k_j|^2} $, 则有$ \|\varepsilon^k\|^2=\int^L_0|\varepsilon^k(x)|^2{\rm d}x $. 根据Parseval等式, 有

$ \sum\limits^\infty_{j=-\infty}|\beta^k_j|^2=\frac{1}{L}\int^L_0{|\varepsilon^k(x)|^2{\rm d}x}, $

进一步可得

(3.2) $ \begin{equation} \|\varepsilon^k\|^2=L\sum\limits^\infty_{j=-\infty}|\beta^k_j|. \end{equation} $

假设$ x_i\approx L_1+{\rm i}h $, 基于以上分析, 设(3.1)式的解为

(3.3) $ \begin{equation} \varepsilon^k_i=\beta^ke^{\lambda(L_1+{\rm i}h)I}, \quad \lambda=\frac{2\pi j}{L}. \end{equation} $

将(3.3)式带入(3.1)式, 可得

(Ⅰ) 当$ k=0 $时, 有

$ \begin{eqnarray*} \left (-ad\sum\limits_se^{\lambda(L_1+sh)I}\Psi''_s(x_i)-bd\sum\limits_se^{\lambda(L_1+sh)I}\Psi' _s(x_i)+(1+rd)e^{\lambda (L_1+{\rm i}h)I} \right )\beta^1=\beta ^0e^{\lambda (L_1+{\rm i}h)I}, \end{eqnarray*} $

等价地

(3.4) $ \begin{equation} \left (-ad\sum\limits_se^{\lambda(s-i)hI}\Psi''_s(x_i)-bd\sum\limits_se^{\lambda(s-i)hI}\Psi' _s(x_i)+1+rd \right )\beta^1=\beta ^0. \end{equation} $

注意到函数$ g(x)=e^{I\lambda(x-x_i)} $的MQ拟插值的逼近形式为

$ \begin{eqnarray*} \sum\limits_se^{\lambda(x_s-x_i)I}\Psi_s(x)=\sum\limits_se^{\lambda(L_1+sh)I}\cdot e^{-\lambda (L_1+{\rm i}h)I}\Psi_s(x)=\sum\limits_se^{\lambda(s-i)hI}\Psi_s(x). \end{eqnarray*} $

因此, $ g^{\prime}(x)=\lambda Ie^{I\lambda (x-x_i)} $和$ g^{\prime\prime}(x)=-\lambda ^2 e^{I\lambda (x-x_i)} $的逼近形式分别为$ \sum\limits_se^{I\lambda (s-i)h} \Psi' _s(x) $和$ \sum\limits_se^{I\lambda (s-i)h} \Psi''_s(x) $, 因而(3.4)式等价于

(3.5) $ \begin{equation} (1+rd+ad\lambda^2-bd\lambda I)\beta^1=\beta^0. \end{equation} $

(Ⅱ) 当$ k\geq 1 $时, 同理可得

(3.6) $ \begin{equation} (1+rd+ad\lambda^2-bd\lambda I)\beta^{k+1}=\sum\limits^{k+1}_{j=0}(b_j-b_{j+1})\beta^{k-j}+b_k\beta^0. \end{equation} $

引理3.1  根据文献[9, 16], $ b_j=(j+1)^{1-\alpha}-j^{1-\alpha} $满足

(1) $ 1=b_0>b_1>b_2>\cdots >b_k>0 $, 当$ k\rightarrow \infty $时, $ b_k\rightarrow 0 $;

(2) $ \sum\limits^{k=1}_{j=0}(b_j-b_{j+1})+b_k=1 $.

命题3.1  对任何$ k>0 $, $ |\beta^k|\leq |\beta^0| $.

证  由$ a, r, d>0 $, 则有$ |1+rd+ad\lambda^2-bd\lambda I|\geq 1 $, 从(3.5)式我们很容易推出$ |\beta^1|\leq |\beta^0| $. 对于$ j=2, 3\cdots k $, 假设$ |\beta^j|\leq |\beta^0| $. 根据(3.6)式和引理3.1, 得

$ \begin{eqnarray*} |(1+rd+ad\lambda^2-bd\lambda I)\beta^{k+1}| &=&|(1+rd+ad\lambda^2-bd\lambda I)||\beta^{k+1}|\\ & =&\bigg|\sum\limits^{k-1}_{j=0}(b_j-b_{j+1})\beta^{k-j}+b_k\beta^0\bigg| \\ &\leq &\bigg(\sum\limits^{k-1}_{j=0}(b_j-b_{j+1})+b_k\bigg)|\beta^0|=|\beta^0|. \end{eqnarray*} $

这就意味着$ |\beta^{k+1}|\leq |\beta^0| $. 证毕.

定理3.1  MQ拟插值格式(2.6)是无条件稳定的.

证  应用(3.2)式和命题3.1, 得

$ \|\varepsilon^k\|\leq \|\varepsilon^0\|, \quad k=1, 2, \cdots , N. $

即三阶MQ拟插值格式是无条件稳定的. 证毕.

由文献[13, 定理2]简单修改可得出下面定理.

定理3.2  三阶MQ拟插值格式是$ L_2 $收敛的, 且存在一个常数$ C>0 $使得

$ \|U^k-\widetilde{U}_k\|\leq C(\Delta t^{2-\alpha}+h^\ell), \quad k=1, 2, \cdots , N. $

注3.1  与经典有限差分法的误差项$ {\cal O}(\Delta t^{2-\alpha}+h^2) $相比(参见文献[16]), 本文所构造格式的误差项为$ {\cal O}(\Delta t^{2-\alpha}+h^\ell) $, 且$ \ell $大于2 (参见文献[14]).

对于给定数据点$ \{x_i\}^M_{i=0} $, 添加四个辅助点, 满足

$ x_{-2}<x_{-1}<x_0<x_1<\cdots <x_M<x_{M+1}<x_{M+2}. $

选择径向基函数为

$ \phi_i(x)= \left\{\begin{array}{ll} (x-x_i)^3, & -2\leq i\leq 1, \\ \sqrt{((x-x_i)^2+c^2)^3}, {\quad} & 2\leq i \leq M-2, \\ (x_i-x)^3, & M-1 \leq i \leq M+2. \end{array}\right. $

以下是Matlab 2010的运算结果, 实验是在一个装有8 GB RAM和2.70 GHz速度的Inter Core i5笔记本电脑上完成的.

例4.1  考虑模型

(4.1) $ \begin{equation} \left\{\begin{array} {l} { }{}^C_{0}D_{\tau}^{\alpha}U(x, \tau)=a\frac{\partial^2U(x, \tau)}{\partial x^2}+b\frac{\partial U(x, \tau)}{\partial x}-rU(x, \tau)+f(x, \tau), (x, \tau)\in(-1, 1)\times(0, T), \\ U(-1, \tau)=0, U(1, \tau)=0, \\ U(x, 0)=x^2(1-x), \end{array} \right. \end{equation} $

这里

$ \begin{eqnarray*} f(x, \tau)&=&\Big(\frac{2\tau^{2-\alpha}}{\Gamma(3-\alpha)}+\frac{2\tau^{1-\alpha}}{\Gamma(2-\alpha)}\Big)x^2(1-x)\\ &&-(\tau+1)^2\Big(a(2-6x)+b(2x-3x^2)-rx^2(1-x)\Big), \end{eqnarray*} $

满足精确解表达式$ U(x, \tau)=(\tau+1)^2x^2(1-x^2) $.

设$ r=0.05, \sigma=0.25, a=\sigma^2/2, b=r-a $, $ T=1 $. 考虑非均匀数据点$ \{x_i:x_i=\cos\frac{{\rm i}\pi}{M}\}^M_{i=0} $上的求解问题, 根据文献[13], 选取$ c=0.2h^{1/3} $, $ h=1/(M-1) $. 定义均方误差$ (L_2) $和最大误差$ (L_\infty) $分别为

$ L_2=\sqrt{\frac{1}{M+1}\sum\limits^M_{i=0}\Big(U(x_i, T)-U_{\rm appr}(x_i, T)\Big)^2}, $

$ L_\infty=\max\limits_{0\leq i\leq M}|U(x_i, T)-U_{\rm appr}(x_i, T)|. $

它们相应的收敛速度分别列于表 1和表 2. 这里$ {\rm rate=}\frac{\log(\frac{E_i}{E_{i+1}})}{\log(\frac{h_i}{h_{i+1}})}, $$ E_i $, $ E_{i+1} $分别是相应于$ h_i $, $ h_{i+1} $的均方误差(或者最大误差). 如表格所示, 通过三阶MQ拟插值法获得的数值解与精确解高度一致. 在运用拟插值求解微分方程时, 不需要求解大尺度矩阵方程, 因而算法对较多的形状参数都比较稳定, 见表 3. 此外, 由于MQ函数是无限连续且可微的, 因此可以直接用于计算导数. 图 1给出了一阶和二阶导数的图形, 如图所示, 导数精确解与数值解高度一致.

例4.2  考虑如下系统

(4.2) $ \begin{equation} \left\{\begin{array} {l} { }{} ^C_{0}D_{\tau}^{\alpha}U(x, \tau)=a\frac{\partial^2U(x, \tau)}{\partial x^2}+b\frac{\partial U(x, \tau)}{\partial x}-rU(x, \tau)+f(x, \tau), (x, \tau)\in(0, 1)\times(0, T), \\ U(0, \tau)=(1+\tau)^2, U(1, \tau)=3(1+\tau)^2, \\ U(x, 0)=1+x^2+x^3, \end{array} \right. \end{equation} $

其中

$ \begin{eqnarray*} f(x, \tau)&=&\Big(\frac{2\tau^{2-\alpha}}{\Gamma(3-\alpha)}+\frac{2\tau^{1-\alpha}}{\Gamma(2-\alpha)}\Big) (1+x^2+x^3)\\ &&-(1+\tau)^2\Big(2a(1+3x)+b(2x+3x^2)-r(1+x^2+x^3)\Big), \end{eqnarray*} $

满足精确解表达式

(4.3) $ \begin{equation} U(x, \tau)=(1+\tau)^2(1+x^2+x^3). \end{equation} $

选取$ \alpha=0.7, r=0.5, a=1, b=r-a $, 做非均匀数据点$ \{x_i:x_i=1/2\cos\frac{{\rm i}\pi}{M}+1/2\}^M_{i=0} $$ (M=21) $上的模拟. 图 2给出了$ T=0.5 $, $ T=1.5 $, $ T=2.5 $时的三个截图. 如图所示, 在非均匀数据点的情况下, 依然可以提供很好的数值模拟. 当计算到$ T=2.5 $时, 在经过近$ 2500 $个时间步时$ (\Delta t=0.001) $后, CPU时间小于$ 2 $秒, 而运用有限差分计算相同的时间步需要47秒, 这也说明这种方法是高效的.

取$ T=1.5 $, $ M=21 $, $ N=1000 $, 对形状参数$ c $做敏感性分析(见表 4), 本例中采用留一交叉验证法(leave-one-out cross validation, 参见文献[7–8]选取最优的$ c=0.63 $. 由表 4可看出, 最小误差集中在"最优的$ c $" 附近, 随着$ c $的取值偏离, 误差也变得很大.

本文运用三阶MQ拟插值方法求解时间分数阶B-S期权模型. 借助傅立叶分析, 验证了算法的稳定性和收敛性. 无论是均匀点还是非均匀点, 数值结果都能有较好的模拟效果. 将这种方法应用于时间分数阶的障碍期权、回望期权和其他奇异期权的求解是我们以后要研究的课题.