分数阶Black-Scholes模型下的波动率校准问题

Volatility Calibration Problems in A Fractional Black-Scholes Models

摘要

对分数阶Black-Scholes模型中的波动率函数进行反演， 提出了一种准确、稳健的数值算法. 首先， 对于正问题， 考虑到收益函数的奇异性会影响 L1格式的收敛速度， 提出了一种基于改进 L1格式的有限差分方法. 此数值方法能有效恢复 L1格式的收敛性， 且在计算中只需求解稀疏的三对角线性系统. 其次， 对于反问题， 考虑了随时间变化的波动率函数， 波动率反演问题可以表述为求解损失函数的最小值. 构造了一个连续、分段线性的波动率函数， 并采用预测-校正方法来抑制可能的振荡. 数值计算和实证分析的结果表明了所提方法的准确性和可靠性.

Abstract

In this paper， an accurate and robust numerical algorithm is proposed to invert the volatility function in the fractional Black-Scholes model. First， for the direct problem， considering that the singularity of the payoff function affects the convergence speed of the L1 method， a finite difference method based on the improved L1 method is proposed. This numerical method can effectively recover the convergence of the L1 method， and only sparse tridiagonal linear systems need to be solved during the computation. Moreover， for the inverse problem， considering the time-dependent volatility function， the volatility inversion problem can be formulated as minimizing the loss function. A continuous and piecewise linear volatility function is constructed and a predictor-corrector approach to mitigate potential oscillations is employed. The results of numerical simulations and empirical analyses demonstrate the accuracy and reliability of the proposed method.

Graphical abstract

关键词

分数阶BS模型 / 反问题 / 改进 L1格式 / 欧式期权 / 波动率校准

Key words

fractional Black-Scholes model / inverse problem / improved L1 method / European options / volatility calibration

期权定价理论的重大突破始于Black等^［1］于1973年提出的Black-Scholes（BS）模型， 在过去几十年中取得了前所未有的成功. 该模型在标的资产服从几何布朗运动以及波动率是常数的假设下， 给出了一个简洁的欧式期权定价公式. 然而， 近年来的金融实证表明， BS 模型的一些假设是比较理想的， 与市场的数据存在一些偏差. 因此， 许多学者提出了 BS 模型的改进模型，包括跳-扩散模型^［2-3］、随机波动率模型^［4-6］和分数阶模型^［7-12］等. 近年来， 分数阶微积分方程引起了学者们的浓厚兴趣， 分数阶导数的非局部特征与金融市场中的分形特征相契合，使得分数阶微积分方程能更灵活地刻画资产变化的自相关性和基于时间的长期记忆性等， 为金融市场的价格波动提供更准确的描述. 2000年， Wyss^［7］首次将分形的思想运用到金融领域， 提出了时间分数阶 BS 期权定价模型， 并给出了欧式期权价格的闭式解. Cartea等^［8］结合FMLS 过程、CGMY 过程和 KoBoL 过程， 推导了3种空间分数阶 BS 定价模型. Jumarie^［9-10］利用分数阶 Taylor 级数提出了时空分数阶 BS 模型. Liang等^［11］建立了双分数 BS 模型， 并获得了欧式期权的显式定价公式. Zhang等^［12］讨论了时空分数阶期权定价模型.

期权定价模型的校准问题是数理金融学的一个经典问题. Crépey^［13］在广义 BS 模型下研究了局部波动率校准的 Tikhonov 正则化方法. Geng等^［14］利用二阶 Tikhonov 正则化校准局部波动率曲面. Zhang等^［15］结合基于罚函数的数值方法和双三次样条正则化方法， 通过美式期权的价格对波动率函数进行校准. Jin等^［16］在广义 BS 模型下，基于预测-校正方法提出了一种简单而稳健的数值算法. Jiang等^［17］在时间分数阶BS 模型下利用线性化方法反演了隐含波动率.

本文考虑时间分数阶 BS 模型下的波动率反演问题， 这是一个富有挑战性的反问题. 对期权价格高效的数值计算是波动率反演的一个重要部分， 分数阶导数的非局部性增加了数值计算的复杂度. 基于上述考虑， 文中采用基于改进 L1 格式的有限差分方法计算期权价格. 此外， 还构造了一个连续、分段线性的波动率函数， 并采用预测-校正算法抑制可能产生的振荡， 提出了一种准确、稳健的数值算法.

1 模型介绍

考虑时间分数阶BS模型

\left\{\begin{array}{l}\begin{array}{l}\frac{{\partial }^{\alpha }C(S,t)}{\partial t^{\alpha }}+\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}(t)S^{\mathrm{2}}\frac{{\partial }^{\mathrm{2}}C(S,t)}{\partial S^{\mathrm{2}}}+\\ r(t)S\frac{\partial C(S,t)}{\partial S}-r(t)C(S,t)=\mathrm{0},\end{array}\\ C(\mathrm{0},t)=p(t),C(\mathrm{\infty },t)=q(t),\\ C(S,T)=g(S),\end{array}\right.（1）

其中：C(S,t)表示执行价格为K， 到期时间为T的欧式期权价格；\sigma (t)和r(t)分别表示波动率函数和无风险利率；\frac{{\partial }^{\alpha }C(S,t)}{\partial t^{\alpha }}表示右修正Riemann-Liouville 分数阶导数， 定义为

；

g(S)表示期权的收益函数. 对于欧式看涨期权，其收益函数为g(S)=\mathrm{m}\mathrm{a}\mathrm{x}\{S-K,\mathrm{0}\}.对于欧式看跌期权，其收益函数为g(S)=\mathrm{m}\mathrm{a}\mathrm{x}\{K-S,\mathrm{0}\}.边界条件p(t)和q(t)取决于期权种类. 令x=\mathrm{l}\mathrm{n}S,\tau =T-t,U(x,\tau )=C({\mathrm{e}}^x,T-\tau ),则模型（1）变为

\left\{\begin{array}{l}\begin{array}{l}{}_{\mathrm{0}}{D}_{\tau }^{\alpha }U\left(x,\tau \right)=\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}\left(\tau \right)\frac{{\partial }^{\mathrm{2}}U\left(x,\tau \right)}{\partial x^{\mathrm{2}}}+\\ \left(r\left(\tau \right)\left.-\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}\left(\tau \right)\right)\right.\frac{\partial U\left(x,\tau \right)}{\partial x}-r\left(\tau \right)U\left(x,\tau \right),\end{array}\\ U\left(-\mathrm{\infty },\tau \right)=p\left(\tau \right),U\left(\mathrm{\infty },\tau \right)=q\left(\tau \right),\\ U\left(x,\mathrm{0}\right)=g\left(x\right),\end{array}\right.

其中：{}_{\mathrm{0}}{D}_{\tau }^{\alpha }U\left(x,\tau \right)表示左修正Riemann-Liouville分数阶导数， 其定义为

{}_{\mathrm{0}}{D}_{\tau }^{\alpha }U\left(x,\tau \right)=\frac{\mathrm{1}}{\mathrm{\Gamma }\left(\mathrm{1}-\alpha \right)}\frac{\partial }{\partial \tau }{\int }_{\mathrm{0}}^{\tau }\frac{U\left(x,\eta \right)-U\left(x,\mathrm{0}\right)}{{\left(\tau -\eta \right)}^{\alpha }}\mathrm{d}\eta.

为了对上述模型进行数值求解，需要将无界区域截断为有界区间\left[B_l,B_r\right]以便计算， 于是

\left\{\begin{array}{l}\begin{array}{l}{}_{\mathrm{0}}{D}_{\tau }^{\alpha }U\left(x,\tau \right)=\\ a\left(\tau \right)\frac{{\partial }^{\mathrm{2}}U\left(x,\tau \right)}{\partial x^{\mathrm{2}}}+b\left(\tau \right)\frac{\partial U\left(x,\tau \right)}{\partial x}-\\ r\left(\tau \right)U\left(x,\tau \right)\triangleq \mathrm{ℒ}U\left(x,\tau \right),\end{array}\\ U\left(B_l,\tau \right)=p\left(\tau \right),U\left(B_r,\tau \right)=q\left(\tau \right),\\ U\left(x,\mathrm{0}\right)=g\left(x\right),\end{array}\right.（2）

其中：a(\tau )=\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}(\tau ),b(\tau )=r(\tau )-\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}(\tau ).

2 数值解

本节给出欧式期权定价的数值方法. 令{\tau }_k=k\mathrm{\Delta }\tau ,k=\mathrm{0,1},\mathrm{2},\cdots ,N；x_i=B_l+ih,i=\mathrm{0,1},\mathrm{2},\cdots ,M,其中\mathrm{\Delta }\tau =T/N表示时间步长，h=(B_r-B_l)/M表示空间步长. 当时， 分数阶导数{}_{\mathrm{0}}{D}_{\tau }^{\alpha }U(x,\tau )与Caputo分数阶导数是一致的， 于是{}_{\mathrm{0}}{D}_{\tau }^{\alpha }U(x,\tau )在点(x_i,{\tau }_k)处可离散为

\begin{array}{c}{}_{\mathrm{0}}{D}_{\tau }^{\alpha }U(x_i,{\tau }_k)=\rho {\displaystyle \sum _{j=\mathrm{0}}^k}{\omega }_{k-j}U(x_i,{\tau }_j)+O(\mathrm{\Delta }{\tau }^{\mathrm{2}-\alpha }),\end{array}

其中\rho =\mathrm{\Delta }{\tau }^{-\alpha }/\mathrm{\Gamma }(\mathrm{2}-\alpha )且

{\omega }_j=\left\{\begin{array}{l}\mathrm{1},j=\mathrm{0},\\ {\left(j+\mathrm{1}\right)}^{\mathrm{1}-\alpha }+{\left(j-\mathrm{1}\right)}^{\mathrm{1}-\alpha }-\mathrm{2}{j}^{\mathrm{1}-\alpha },\\ j=\mathrm{1,2},\cdots ,k-\mathrm{1},\\ {\left(j-\mathrm{1}\right)}^{\mathrm{1}-\alpha }-j^{\mathrm{1}-\alpha },j=k.\end{array}\right.

由于收益函数的奇异性会影响数值格式的收敛速度， 考虑对L1格式^［18］的起始步骤进行修正， 这就是改进L1格式， 详见文献［19］. 令V(x,\tau )=U(x,\tau )-g(x)， 则模型（2）等价于

\left\{\begin{array}{l}{}_{\mathrm{0}}{D}_{\tau }^{\alpha }V\left(x,\tau \right)=\mathrm{ℒ}V\left(x,\tau \right)+\mathrm{ℒ}g\left(x\right),\\ \begin{array}{l}V\left(B_l,\tau \right)=p\left(\tau \right)-g\left(B_l\right),\\ V\left(B_r,\tau \right)=q\left(\tau \right)-g\left(B_r\right),\end{array}\\ V\left(x,\mathrm{0}\right)=\mathrm{0}.\end{array}\right.（3）

记V^k(x)=V(x,{\tau }_k)，k=\mathrm{0,1},\cdots ,N,则V^{\mathrm{0}}(x)=\mathrm{0}，于是模型（3）可近似为

\left\{\begin{array}{l}\rho V^{\mathrm{1}}(x)=\mathrm{ℒ}{V}^{\mathrm{1}}(x)+(\mathrm{1}+c_{\mathrm{0}})\mathrm{ℒ}g(x),k=\mathrm{1},\\ \rho (V^{\mathrm{2}}(x)+{\omega }_{\mathrm{1}}{V}^{\mathrm{1}}(x))=\\ \mathrm{ℒ}{V}^{\mathrm{2}}(x)+(\mathrm{1}+c_{\mathrm{1}})\mathrm{ℒ}g(x),k=\mathrm{2},\\ \rho {\displaystyle \sum _{j=\mathrm{1}}^k}{\omega }_{k-j}{V}^j(x)=\mathrm{ℒ}{V}^k(x)+\mathrm{ℒ}g(x),\\ k=\mathrm{3,4},\cdots ,N,\end{array}\right.

其中：c_{\mathrm{0}}=\mathrm{11}/\mathrm{12},c_{\mathrm{1}}=-\mathrm{5}/\mathrm{12}. 边界条件为V^k(B_l)=p({\tau }_k)-g(B_l),V^k(B_r)=q({\tau }_k)-g(B_r), 对于一阶和二阶空间导数， 其二阶有限差分近似为

\frac{\partial V(x_i,{\tau }_k)}{\partial x}=\frac{V(x_{i+\mathrm{1}},{\tau }_k)-V(x_{i-\mathrm{1}},{\tau }_k)}{\mathrm{2}h}+O(h^{\mathrm{2}}),

\begin{array}{l}\frac{{\partial }^{\mathrm{2}}V(x_i,{\tau }_k)}{\partial x^{\mathrm{2}}}=\\ \frac{V(x_{i+\mathrm{1}},{\tau }_k)-\mathrm{2}V(x_i,{\tau }_k)+V(x_{i-\mathrm{1}},{\tau }_k)}{h^{\mathrm{2}}}+O(h^{\mathrm{2}}).\end{array}

令V_i^k表示V(x_i,{\tau }_k)的近似值， 可得差分格式为

\begin{array}{l}{\tilde{\lambda }}^{(\mathrm{1})}{V}_{i-\mathrm{1}}^{\mathrm{1}}+{\tilde{\kappa }}^{(\mathrm{1})}{V}_i^{\mathrm{1}}+{\tilde{\theta }}^{(\mathrm{1})}{V}_{i+\mathrm{1}}^{\mathrm{1}}=\\ (\mathrm{1}+c_{\mathrm{0}})\left({\lambda }^{(\mathrm{1})}{g}_{i-\mathrm{1}}+{\kappa }^{(\mathrm{1})}{g}_i\right.\left.+{\theta }^{(\mathrm{1})}{g}_{i+\mathrm{1}}\right),\end{array}

\begin{array}{l}{\tilde{\lambda }}^{(\mathrm{2})}{V}_{i-\mathrm{1}}^{\mathrm{2}}+{\tilde{\kappa }}^{(\mathrm{2})}{V}_i^{\mathrm{2}}+{\tilde{\theta }}^{(\mathrm{2})}{V}_{i+\mathrm{1}}^{\mathrm{2}}=\\ (\mathrm{1}+c_{\mathrm{1}})\left({\lambda }^{(\mathrm{2})}{g}_{i-\mathrm{1}}+{\kappa }^{(\mathrm{2})}{g}_i\right.\left.+{\theta }^{(\mathrm{2})}{g}_{i+\mathrm{1}}\right)-\rho {\omega }_{\mathrm{1}}{V}_i^{\mathrm{1}},\end{array}

{\tilde{\lambda }}^{(k)}{V}_{i-\mathrm{1}}^k+{\tilde{\kappa }}^{(k)}{V}_i^k+{\tilde{\theta }}^{(k)}{V}_{i+\mathrm{1}}^k={\lambda }^{(k)}{g}_{i-\mathrm{1}}+{\kappa }^{(k)}{g}_i+{\theta }^{(k)}{g}_{i+\mathrm{1}}-\rho {\displaystyle \sum _{j=\mathrm{1}}^{k-\mathrm{1}}}{\omega }_{k-j}{V}_i^j,k\ge \mathrm{3},(\mathrm{4})

其中：

\begin{array}{c}{\tilde{\lambda }}^{(k)}=\frac{b({\tau }_k)}{\mathrm{2}h}-\frac{a({\tau }_k)}{h^{\mathrm{2}}},{\tilde{\kappa }}^{(k)}=\frac{\mathrm{2}a({\tau }_k)}{h^{\mathrm{2}}}+r({\tau }_k)+\rho ,\\ {\tilde{\theta }}^{(k)}=-\frac{a({\tau }_k)}{h^{\mathrm{2}}}-\frac{b({\tau }_k)}{\mathrm{2}h},{\lambda }^{(k)}=\frac{a({\tau }_k)}{h^{\mathrm{2}}}-\frac{b({\tau }_k)}{\mathrm{2}h},\\ {\kappa }^{(k)}=-\frac{\mathrm{2}a({\tau }_k)}{h^{\mathrm{2}}}-r({\tau }_k),{\theta }^{(k)}=\frac{a({\tau }_k)}{h^{\mathrm{2}}}+\frac{b({\tau }_k)}{\mathrm{2}h}.\end{array}

边界条件和初值条件可离散为

\begin{array}{c}{V}_{\mathrm{0}}^k=p({\tau }_k)-g(B_l),V_M^k=q({\tau }_k)-g(B_r),k\ge \mathrm{1},\\ V_i^{\mathrm{0}}=\mathrm{0},i=\mathrm{0,1},\mathrm{2},\cdots ,M.\end{array}

利用差分格式（4）计算V_i^k， 于是期权价格的数值解为U(x_i,{\tau }_k)=V_i^k+g_i. 令S_i={\mathrm{e}}^{x_i},t_k=k\mathrm{\Delta }\tau ,则C(S_i,t_k)=U(x_i,{{\tau }_{N-}}_k),i=\mathrm{0,1},\mathrm{2},\cdots ,M；k=\mathrm{0,1},

\mathrm{2},\cdots ,N.

3 反问题

本节给出波动率反演问题的表述， 并介绍预测-校正数值算法^［16］. 考虑随时间变化的波动率函数\sigma (t)， 并假设\sigma (t)是连续、分段线性的. 令U_{ij}表示到期日为T_i(i=\mathrm{1,2},\cdots ,M_t)， 执行价格为K_j(j=\mathrm{1,2},\cdots ,M_k)的期权市场价格， T_{\mathrm{1}}\le T_{\mathrm{2}}\le \cdots \leT_{M_t}，K_{\mathrm{1}}\le K_{\mathrm{2}}\le \cdots \le K_{M_k}. 令U_{ij}^{\mathrm{a}}和U_{ij}^{\mathrm{b}}分别表示期权的卖价和买价.

校准问题 给定一组期权的市场价格U_{ij}(i=\mathrm{1,2},\cdots ,M_t,j=\mathrm{1,2},\cdots ,M_k)， 寻找波动率函数\sigma (t)，使得下列均方误差最小：

\begin{array}{c}\begin{array}{l}{\Gamma }_i(\sigma )=\frac{\mathrm{1}}{M_k}{\displaystyle \sum _{j=\mathrm{1}}^{M_k}}{\left[U(S_{\mathrm{0}},\mathrm{0},K_j,T_i,\sigma (t))-U_{ij}\right]}^{\mathrm{2}}{\omega }_{ij},\\ i=\mathrm{1,2},\cdots ,M_t,\end{array}\end{array}

其中：U(S_{\mathrm{0}},\mathrm{0},K_j,T_i,\sigma (t))表示S=S_{\mathrm{0}}处欧式期权的理论价格，这里以期权价格的数值解近似理论价格，采用上节提出的改进L1方法计算； U_{ij}=(U_{ij}^{\mathrm{a}}+U_{ij}^{\mathrm{b}})/\mathrm{2} 是期权买价与卖价的平均值； {\omega }_{ij}表示权重. 采用最速下降法^［20］求解下面的优化问题.令T_{i+\mathrm{1}/\mathrm{2}}=(T_i+T_{i+\mathrm{1}})/\mathrm{2},i=\mathrm{0,1},\cdots ,M_t-\mathrm{1},T_{\mathrm{0}}=\mathrm{0}.下面给出预测-校正数值算法的详细过程.

步骤1 假设波动率\sigma (t)是[\mathrm{0},T_{\mathrm{1}}]上的常值函数， 即\sigma (t)=y_{\mathrm{1}},t\in [\mathrm{0},T_{\mathrm{1}}]. 下面只需确定参数y_{\mathrm{1}}的值. 寻求最优值y_{\mathrm{1}}以最小化下列目标函数：

\begin{array}{c}{\Gamma }_{\mathrm{1}}(\sigma )=\frac{\mathrm{1}}{M_k}{\displaystyle \sum _{j=\mathrm{1}}^{M_k}}{\left[U(S_{\mathrm{0}},\mathrm{0},K_j,T_{\mathrm{1}},\sigma (t))-U_{\mathrm{1}j}\right]}^{\mathrm{2}}{\omega }_{\mathrm{1}j},\end{array}

并记{\sigma }_{\mathrm{1}}(t)=\sigma (t).

步骤2 由步骤1已获得参数y_{\mathrm{1}}的值，进一步假设波动率\sigma (t)是[\mathrm{0},T_{\mathrm{2}}]上通过点(T_{\mathrm{1}/\mathrm{2}},y_{\mathrm{1}})的线性函数， 定义为

\begin{array}{c}\sigma (t)=\frac{y_{\mathrm{2}}-y_{\mathrm{1}}}{T_{\mathrm{2}}-T_{\mathrm{1}/\mathrm{2}}}(t-T_{\mathrm{2}})+y_{\mathrm{2}},t\in [\mathrm{0},T_{\mathrm{2}}].\end{array}

下面只需确定参数y_{\mathrm{2}}的值.寻求最优值y_{\mathrm{2}}以最小化下列目标函数：

\begin{array}{c}{\Gamma }_{\mathrm{2}}(\sigma )=\frac{\mathrm{1}}{M_k}{\displaystyle \sum _{j=\mathrm{1}}^{M_k}}{\left[U(S_{\mathrm{0}},\mathrm{0},K_j,T_{\mathrm{2}},\sigma (t))-U_{\mathrm{2}j}\right]}^{\mathrm{2}}{\omega }_{\mathrm{2}j},\end{array}

并记{\sigma }_{\mathrm{2}}(t)=\sigma (t).

步骤3 在上述步骤中已获得参数y_{\mathrm{1}},y_{\mathrm{2}}的值，进一步假定波动率\sigma (t)是[\mathrm{0},T_{\mathrm{3}}]上的分段线性函数， 定义为

\sigma (t)=\left\{\begin{array}{l}{\sigma }_{\mathrm{2}}(t),t\in [\mathrm{0},T_{\mathrm{3}/\mathrm{2}}],\\ {\sigma }_{\mathrm{3}}(t),t\in [T_{\mathrm{3}/\mathrm{2}},T_{\mathrm{3}}],\end{array}\right.

其中{\sigma }_{\mathrm{3}}(t)为通过点\left(T_{\mathrm{3}/\mathrm{2}},{\sigma }_{\mathrm{2}}(T_{\mathrm{3}/\mathrm{2}})\right)的线性函数， 即

\begin{array}{c}{\sigma }_{\mathrm{3}}(t)=\frac{y_{\mathrm{3}}-{\sigma }_{\mathrm{2}}(T_{\mathrm{3}/\mathrm{2}})}{T_{\mathrm{3}}-T_{\mathrm{3}/\mathrm{2}}}(t-T_{\mathrm{3}})+y_{\mathrm{3}}.\end{array}

为了获得波动率函数， 需要确定最优值y_{\mathrm{3}}以最小化下列目标函数：

\begin{array}{c}{\Gamma }_{\mathrm{3}}(\sigma )=\frac{\mathrm{1}}{M_k}{\displaystyle \sum _{j=\mathrm{1}}^{M_k}}{\left[U(S_{\mathrm{0}},\mathrm{0},K_j,T_{\mathrm{3}},\sigma (t))-U_{\mathrm{3}j}\right]}^{\mathrm{2}}{\omega }_{\mathrm{3}j}.\end{array}

步骤4 对于\mathrm{4}\le i\le M_t， 重复执行以下过程. 假设波动率\sigma (t)是[\mathrm{0},T_i] 上的分段线性函数：

\sigma (t)=\left\{\begin{array}{l}{\sigma }_{\mathrm{2}}(t),t\in [\mathrm{0},T_{\mathrm{3}/\mathrm{2}}],\\ {\sigma }_m(t),t\in [T_{m-\mathrm{3}/\mathrm{2}},T_{m-\mathrm{1}/\mathrm{2}}],\mathrm{3}\le m\le i-\mathrm{1},\\ {\sigma }_i(t),t\in [T_{i-\mathrm{3}/\mathrm{2}},T_i],\end{array}\right.

其中{\sigma }_i(t) 定义为

\begin{array}{c}{\sigma }_i(t)=\frac{y_i-{\sigma }_{i-\mathrm{1}}(T_{i-\mathrm{3}/\mathrm{2}})}{T_i-T_{i-\mathrm{3}/\mathrm{2}}}(t-T_i)+y_i.\end{array}

为了获得波动率函数， 需要确定最优值y_i以最小化下列目标函数：

\begin{array}{c}{\Gamma }_i(\sigma )=\frac{\mathrm{1}}{M_k}{\displaystyle \sum _{j=\mathrm{1}}^{M_k}}{\left[U(S_{\mathrm{0}},\mathrm{0},K_j,T_i,\sigma (t))-U_{ij}\right]}^{\mathrm{2}}{\omega }_{ij}.\end{array}

经过上述步骤，可以得到一个[\mathrm{0},T_{M_t}]上连续、分段线性的波动率函数\sigma (t).

4 数值实验

本节首先通过数值实验检验改进L1方法在欧式期权定价中的准确性和收敛性. 随后， 在最后两小节中，分别采用合成数据和真实市场数据对所提校准算法的效率进行进一步测试. 所有实验都是在Dell Inspiron 14 Intel Core i5-7200U 2.70 GHz上使用MATLAB R2017b软件进行的.

4.1 欧式期权价格的数值结果

本小节通过欧式看涨期权和欧式看跌期权这两个数值例子验证改进L1方法的有效性. 同时， 为了进行对比，也给出了L1方法的计算结果. 此外，收敛阶的定义为

\begin{array}{c}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}=\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{2}}\frac{{‖e_{(\mathrm{\Delta }\tau )}‖}_{\mathrm{\infty }}}{{‖e_{(\mathrm{\Delta }\tau /\mathrm{2})}‖}_{\mathrm{\infty }}},\end{array}

其中：e_{(\mathrm{\Delta }\tau )}=U_{(\mathrm{\Delta }\tau )}-U_{(\mathrm{\Delta }\tau /\mathrm{2})}， U_{(\mathrm{\Delta }\tau )}表示步长为\mathrm{\Delta }\tau时期权价格的数值解.

例1 考虑关于欧式看涨期权的时间分数阶BS方程^［21］：

\left\{\begin{array}{l}\frac{{\partial }^{\alpha }C(S,t)}{\partial t^{\alpha }}+\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}{S}^{\mathrm{2}}\frac{{\partial }^{\mathrm{2}}C(S,t)}{\partial S^{\mathrm{2}}}+\\ (r-D)S\frac{\partial C(S,t)}{\partial S}-rC(S,t)=\mathrm{0},\\ C(\mathrm{0.1},t)=\mathrm{0},\hspace{1em}\\ C(\mathrm{100},t)=\mathrm{100}-K\mathrm{e}\mathrm{x}\mathrm{p}(-r(\mathrm{1}-t)),\\ C(S,T)=\mathrm{m}\mathrm{a}\mathrm{x}\{S-K,\mathrm{0}\}.\end{array}\right.

参数值为r=\mathrm{0.05},\sigma =\mathrm{0.25}， T=\mathrm{1}， K=\mathrm{50}， D=\mathrm{0}.空间节点数设定为M=\mathrm{100}.N\mathrm{和}\alpha的取值见表1， 根据差分公式（4）计算得到V_i^k， 进而可得期权价格的数值解为

\begin{array}{c}C(S_i,t_k)=U(x_i,{{\tau }_{N-}}_k)=V_i^{N-k}+g_i,\\ i=\mathrm{0,1},\mathrm{2},\cdots ,M;k=\mathrm{0,1},\mathrm{2},\cdots ,N.\end{array}

表1和表2分别展示了不同参数\alpha下，改进L1方法与L1方法的数值误差和收敛阶的对比. 通过观察表1和表2可以发现，相较于原始L1方法，本文所提方法的误差值更小， 尤其是在参数\alpha值较小的情况下. 同时， 在收敛速度方面， 改进L1方法在时间上具有\mathrm{2}-\alpha阶精度， 而L1方法只有一阶精度. 这些对比结果证明了所提方法是精确且有效的.

例2 考虑关于欧式看跌期权的时间分数阶BS方程^［21］：

\left\{\begin{array}{l}\frac{{\partial }^{\alpha }C(S,t)}{\partial t^{\alpha }}+\frac{\mathrm{1}}{\mathrm{2}}{\sigma }^{\mathrm{2}}{S}^{\mathrm{2}}\frac{{\partial }^{\mathrm{2}}C(S,t)}{\partial S^{\mathrm{2}}}+\\ (r-D)S\frac{\partial C(S,t)}{\partial S}-rC(S,t)=\mathrm{0},\\ C(\mathrm{0.1},t)=K\mathrm{e}\mathrm{x}\mathrm{p}(-r(\mathrm{1}-t))-\mathrm{0.1},\\ C(\mathrm{100},t)=\mathrm{0},\\ C(S,T)=\mathrm{m}\mathrm{a}\mathrm{x}\{K-S,\mathrm{0}\}.\end{array}\right.

参数值为r=\mathrm{0.05},\sigma =\mathrm{0.25}，T=\mathrm{1}， K=\mathrm{50}， D=\mathrm{0}. 空间节点数M=\mathrm{100}. N\mathrm{和}\alpha的取值见表3，表3和表4 分别展示了不同参数\alpha下，L1方法与改进L1方法的数值误差和收敛阶的对比. 可以观察到， 本文所提方法误差值比L1方法误差值更小， 尤其对较小的\alpha值.此外， 改进L1方法的收敛阶在时间上达到了\mathrm{2}-\alpha. 结果表明，对于欧式期权，所提方法比原始L1方法表现更好.

4.2 波动率校准的数值模拟

本小节将通过两个具有已知准确解的数值实验验证算法的准确性和稳定性. 采用第3节提出的预测-校正数值算法对波动率进行重构，其中欧式期权的理论价格用改进L1方法计算的期权价格的数值解近似.

模型参数取值为r=\mathrm{0.015},\alpha =\mathrm{0.9},B_l=\mathrm{l}\mathrm{n}(\mathrm{0.1}K),B_r=\mathrm{l}\mathrm{n}(\mathrm{2}K). 初始标的资产价格S_{\mathrm{0}}=\mathrm{100}. 到期日为T_i=\mathrm{60}i/\mathrm{360}，i=\mathrm{1,2},\cdots ,\mathrm{6}. 执行价格为K_j=\mathrm{94}+\mathrm{2}j，j=\mathrm{1,2},\cdots ,\mathrm{5}. 离散参数M=\mathrm{100}，\mathrm{\Delta }\tau =\mathrm{1}/\mathrm{360}. 此外， 假设所有的期权价格具有相同的权重. 对于最优值y_i(i=\mathrm{1,2},\cdots ,M_t)，其初始猜测为0.2.

例3 考虑波动率函数^［22］：

\begin{array}{c}\sigma (t)=\mathrm{0.4}+\frac{\mathrm{0.1}}{\mathrm{1}+t}.\end{array}（5）

首先，对于给定的波动率函数（5），通过差分格式（4）计算欧式看涨期权价格， 将此作为期权的市场价格，不同执行价格和到期日的期权价格计算结果见表5. 其次，波动率的重构采用预测-校正算法，目标函数{\Gamma }_i(\sigma )中欧式期权的理论价格U(S_{\mathrm{0}},\mathrm{0},K_j,T_i,\sigma (t))采用差分格式（4）计算. 图1展示了重构波动率与准确波动率的对比. 可以看到，它们是非常一致的.

为了进一步验证算法的稳健性， 在期权价格上添加噪声， 加入噪声后的期权价格\begin{array}{c}{U}_{ij}^{\delta }=U_{ij}+\delta \xi ,\end{array} 其中\xi是区间(-\mathrm{1,1})上的随机数. 用下列误差衡量波动率反演的结果：

\begin{array}{c}\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{1}=\underset{\mathrm{1}\le i\le N}{\mathrm{m}\mathrm{a}\mathrm{x}}|{\sigma }_i-{\widehat{\sigma }}_i|,\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{2}=\sqrt[]{\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}|{\sigma }_i-{\widehat{\sigma }}_i{|}^{\mathrm{2}}},\end{array}

其中{\sigma }_i和{\widehat{\sigma }}_i分别表示真实波动率和重构波动率. 图2展示了添加噪音的校准结果. 观察图2可知， 除了一些细微的、可忽略的波动外， 重构波动率与真实波动率几乎是吻合的. 此外， 表6列出了不同\delta值下波动率反演的误差值. 这些结果表明了预测-校正数值算法的准确性和稳定性.

例4 考虑波动率函数^［23］：

\begin{array}{c}\sigma (t)=\mathrm{0.2}-\mathrm{0.1}\mathrm{l}\mathrm{n}(\mathrm{1.5}+\mathrm{3}t).\end{array}

波动率校准结果见图3.从图3可以看到重构获得的分段线性波动率与真实波动率高度一致. 此外， 误差值为\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{1}=\mathrm{0.001}\mathrm{670}\mathrm{46},\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{2}=\mathrm{0.000}\mathrm{784}\mathrm{21}. 所得结果证实了所提算法能为波动率提供准确的估计. 求解最优化问题的数值结果见表7. 从表7可知，只需几次迭代即可求得最优值，充分表明所提算法在计算角度方面同样高效.

例5 考虑波动率函数：

\begin{array}{c}\sigma (t)=\mathrm{0.2}+\frac{\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{10}t)}{\mathrm{100}}.\end{array}

由于此函数波动相对较大，选取到期日为T_i=\mathrm{20}i/\mathrm{360}，i=\mathrm{1,2},\cdots ,\mathrm{18}.图4对比了重构波动率与真实波动率. 此外，校准的误差值为\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{1}=\mathrm{0.001}\mathrm{265}\mathrm{236},\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{2}=\mathrm{0.000}\mathrm{254}\mathrm{875}. 可以看到所提算法能有效重构波动率函数.

4.3 实证分析

本小节将利用真实的市场数据展示所提方法的实用性. 基于2022年6月13日上海证券交易所50ETF欧式看涨期权价格，进行波动率函数估计. 将期权卖价与买价的平均值作为市场价格. 表8列出了具有不同执行价格和到期日的欧式看涨期权价格. 考虑了4种到期日：T_{\mathrm{1}}=\mathrm{9}/\mathrm{365},T_{\mathrm{2}}=\mathrm{44}/\mathrm{365},T_{\mathrm{3}}=\mathrm{107}/\mathrm{365},T_{\mathrm{4}}=\mathrm{198}/\mathrm{365}和8种执行价格：K_j=\mathrm{2.600}\mathrm{0}+\mathrm{0.05}j,j=\mathrm{1},\cdots ,\mathrm{8}.上证50ETF的现值为S_{\mathrm{0}}=\mathrm{2.835}\mathrm{0}， 无风险利率r=\mathrm{1.408}\mathrm{0}\mathrm{\%} （以上海银行间同业拆借利率为近似值）. 参数值取为B_l=\mathrm{l}\mathrm{n}(\mathrm{0.1}K),B_r=\mathrm{l}\mathrm{n}(\mathrm{2}K),M=\mathrm{100},\mathrm{\Delta }\tau =\mathrm{1}/\mathrm{365}. 校准结果可以通过理论期权价格与市场价格间的均方根误差 （RMSE）衡量.不同\alpha值下的均方根误差见表9.结果显示，所提方法的结果令人满意，并且当\alpha =\mathrm{0.9}时均方根误差取得最小值.\alpha =\mathrm{0.9}时由上证50ETF市场数据重构的波动率见图5. 图6展示了当\alpha =\mathrm{0.9}时SSE 50ETF市场价格与重构价格的对比.由图6可知对于4种不同期限日，重构的期权价格与市场价格是非常一致的. 除此之外， 求解最优化问题的迭代次数和计算时间（表10）表明了所提算法的效率.

5 结语

本文考虑分数阶BS模型下的波动率函数反演问题，提出了一种准确且稳健的数值算法.为了高效地求解该问题，需要一个既准确又快速的期权定价数值方法.为此，采用了基于改进L1格式的有限差分方法计算期权价格.该方法在计算过程中只需求解稀疏的三对角线性系统. 数值实验结果表明， 原始L1 方法对欧式期权仅具有一阶精度， 而改进方法的收敛阶达到了\mathrm{2}-\alpha.此外， 还给出了波动率校准问题的表述， 这是个不适定反问题. 为求解该问题， 构造了一个连续、分段线性的波动率函数， 并采用预测-校正算法抑制可能出现的振荡. 算法过程中的优化问题采用最速下降法进行求解. 数值模拟和市场数据结果表明，所提方法是准确且稳健的.未来将进一步对数值方法的截断误差、稳定性等进行研究.

参考文献

[1] BLACK F， SCHOLES M. The pricing of options and corporate liabilities［J］. Journal of Political Economy， 1973， 81（3）： 637-654.

[2] MERTON R C. Option pricing when underlying stock returns are discontinuous［J］. Journal of Financial Economics， 1976， 3（1/2）： 125-144.

[3] KOU S G. A jump-diffusion model for option pricing［J］. Management Science， 2002， 48（8）： 1086-1101.

[4] WHITE A D， HULL J C. The pricing of options on assets with stochastic volatilities［J］. The Journal of Finance， 1987， 42（2）：281-300.DOI：10.1111/j.1540-6261.1987.tb02568.x .

[5] STEIN E M， STEIN J C. Stock price distributions with stochastic volatility： an analytic approach［J］. The Review of Financial Studies， 1991， 4（4）： 727-752.

[6] HESTON S L. A closed-form solution for options with stochastic volatility with applications to bond and currency options［J］. The Review of Financial Studies， 1993， 6（2）： 327-343.

[7] WYSS W. The fractional Black-Scholes equation［J］. Fractional Calculus and Applied Analysis，2000，3（1）：51-61.

[8] CARTEA Á， DEL-CASTILLO-NEGRETE D. Fractional diffusion models of option prices in markets with jumps［J］. Physica A： Statistical Mechanics and Its Applications， 2007， 374（2）： 749-763.

[9] JUMARIE G. Derivation and solutions of some fractional Black-Scholes equations in coarse-grained space and time：Application to Merton’s optimal portfolio［J］. Computers & Mathematics with Applications， 2010， 59（3）： 1142-1164.

[10] JUMARIE G. Stock exchange fractional dynamics defined as fractional exponential growth driven by （usual） Gaussian white noise：Application to fractional Black-Scholes equations［J］. Insurance： Mathematics and Economics， 2008， 42（1）： 271-287.

[11] LIANG Jinrong， WANG Jun， ZHANG Wenjun， et al. The solution to a bi-fractional Black-Scholes-Merton differential equation［J］. International Journal of Pure and Applied Mathematics， 2010，58（1）：99-112.

[12] ZHANG Hongmei， LIU Fawang， CHEN Shanzhen， et al. Fast numerical simulation of a new time-space fractional option pricing model governing European call option［J］. Applied Mathematics and Computation， 2018， 339： 186-198.

[13] CRÉPEY S. Calibration of the local volatility in a generalized Black-Scholes model using Tikhonov regularization［J］. SIAM Journal on Mathematical Analysis， 2003， 34（5）： 1183-1206.

[14] GENG Jian， NAVON I M， CHEN Xiao. Non-parametric calibration of the local volatility surface for European options using a second-order Tikhonov regularization［J］. Quantitative Finance， 2014， 14（1）： 73-85.

[15] ZHANG Kai， TEO K L. A penalty-based method from reconstructing smooth local volatility surface from American options［J］. Journal of Industrial & Management Optimization， 2015， 11（2）： 631-644.

[16] JIN Yuzi， WANG Jian， KIM S， et al. Reconstruction of the time-dependent volatility function using the Black-Scholes model［J］. Discrete Dynamics in Nature and Society， 2018， 2018： 3093708.

[17] JIANG Xiaoying， XU Xiang. On implied volatility recovery of a time-fractional Black-Scholes equation for double barrier options［J］. Applicable Analysis， 2021， 100（15）： 3145-3160.

[18] ZHANG Hongmei， LIU Fawang， TURNER I， et al. Numerical solution of the time fractional Black-Scholes model governing European options［J］. Computers & Mathematics with Applications， 2016， 71（9）： 1772-1783.

[19] FORD N J， YAN Yubin. An approach to construct higher order time discretisation schemes for time fractional partial differential equations with nonsmooth data［J］. Fractional Calculus and Applied Analysis， 2017， 20（5）： 1076-1105.

[20] BURDEN R L， FAIRES J D. Numerical analysis， 9th Edition［M］. Boston： Brooks/Cole， 2011.

[21] ROUL P， GOURA V M K P. A compact finite difference scheme for fractional Black-Scholes option pricing model［J］. Applied Numerical Mathematics， 2021， 166： 40-60.

[22] HE Xinjiang， ZHU Songping. On full calibration of hybrid local volatility and regime-switching models［J］. Journal of Futures Markets， 2018， 38（5）： 586-606.

[23] ZHAO Jiajia， XU Zuoliang. Calibration of time-dependent volatility for European options under the fractional Vasicek model［J］. AIMS Mathematics， 2022， 7（6）： 11053-11069.

基金资助

国家自然科学基金资助项目(12071479)