分数阶BS方程中自由边界的数值计算

Numerical Calculation of Free Boundaries in the Fractional BS Equations

摘要: 本文研究一类分数阶美式期权定价模型的数值计算问题。由于分数阶Brown运动能够更好地反映风险资产价格的自相关性和长期记忆性，本文将研究由分数阶Brown运动推导得到的期权价格模型，即分数阶Black-Scholes方程。由于美式期权可以提前执行，该模型的边界条件可归结为分数阶抛物型方程的自由边值问题。自由边界问题的困难之处在于需要同时求出期权的最佳执行边界及分数阶BS方程的数值解，因此本文在美式期权的持有区域使用θ-差分格式预估数值解，再利用边界条件及Taylor公式校正得到收敛的方程数值解和自由边界。

Abstract: This paper studies the numerical computation of a class of fractional-order American option pricing models. Since the fractional-order Brownian motion can better reflect the autocorrelation and long-term memory of risky asset prices, this paper will study the option price model derived from the fractional-order Brownian motion, i.e., the fractional-order Black-Scholes equation. And since American options can be executed in advance, the boundary conditions of this model can be reduced to the free boundary problem of the fractional order parabolic equation. The difficulty of the free boundary problem is that the optimal execution boundary of the option and the numerical solution of the fractional-order BS equation need to be found simultaneously, so in this paper, the numerical solution is predicted in the holding region of the American option using the θ-difference formats, and then the converged numerical solution of the equation and the free boundary are obtained by using the boundary conditions and the correction of Taylor’s formula.

1. 引言

Black-Scholes (BS)期权定价模型为期权定价、风险管理以及其他衍生金融产品的定价提供了一个重要的基础框架。1900年，Louis Bachelier的《投机理论》 [1] 率先将金融理论和数学思想融合到一起，提出利用布朗运动刻画股票价格变化的随机模型。Sprekle [2] 指出了股票价格遵循对数正态分布的假设，极大地推进Louis Bachelier [1] 的模型，避免了股票价格为负数的情形。1963年，Boness [3] 提出投资者具有风险偏好，并在Sprekle模型的基础之上添加了风险因素。1964年，Paul Samnelson [4] 按照期权相比原生标的资产具有更高风险等性质证实了更具实用性的欧式看涨期权模型。Black和Scholes (1973)在有效市场假说和股票价格遵循几何布朗运动等一系列的假设下，推导出了著名的BS期权定价模型 [5] 。但是由于资本市场的本质特征和状态都是随机波动的，与传统BS期权定价模型的假设并不完全吻合，研究人员凭借对股票市场的观察和探索发现，该模型定价与实际市场价格有着较大的差别。Bojdecki [6] 和Tudor [7] 先后研究了修正后的次分数Brown运动，提出了基于次分数Brown运动的BS方程能够更好地刻画原生标的资产价格的变动趋势。伴随微分方程的分形结构在金融领域内被发现，大部分研究人员开始关注金融领域中的分数阶布朗运动，这种模型能更好地反映金融市场中存在的长记忆效应和非正态分布特征。Wyss [8] (2000)最早得出时间分数BS方程，2006年，Cartea等 [9] 给出分数阶跳跃扩散期权定价模型和分数阶障碍期权定价模型，2008年，Jumarie [10] 基于分数阶Taylor公式以及整数阶BS方程的推导过程，得出改进的时间分数阶BS方程。

与欧式期权不同，美式期权的持有人可以在期权到期之前任何时候行使期权，因此期权持有人在整个期权的有效期内必须时刻注意是继续持有还是立即实施。根据期权的内在价值是大于还是等于实施期权的收益，美式期权定价分成两个区域，继续持有区域和终止持有区域，它们之间的交界面称为最佳实施边界。在美式期权中，自由边界是一个重要的概念，因为它可以影响期权的价格和交易策略。

本文针对分数阶美式期权定价方程自由边值问题，构造出一种自由边界的预估–校正方法。这种方法的基本思想是：对每一行，综合利用分数阶方程及自由边界满足的条件，对期权价格在自由边界附近做泰勒展开，从而得到自由边界点的近似估计，然后调整网格进行下一行计算。与现有的数值方法相比，本文方法有两个优点：一是求解区域限定在自由边界内，从而减少了计算量；二是在求出期权值的同时也求出了期权的最佳执行边界。

2. 时间分数阶期权定价模型

期权定价理论最早由Black和Scholes提出，在假设原生资产价格服从几何Brown运动，基于随机微分方程理论推导出了整数阶期权定价模型 [5] ：

分数阶布朗运动可以被用于更准确地描述市场中的非局部和非马尔可夫性质，例如长记忆、长程相关性和非对称性等，其模型更应用于现实世界。由此得到期权价格所满足的欧式期权定价模型，即时间分数阶欧式BS方程 [11] [12] ：

$V_t^{\left(\alpha \right)}\left(S,t\right)+\left(\left(r-D\right)S{V}_S-rV\right)\frac{t^{1-\alpha }}{\Gamma \left(2-\alpha \right)}+\frac{\Gamma \left(1+\alpha \right)}{2}{\sigma }^2{S}^2{V}_{SS}=0.$

其中相关符号定义如下(后续章节继续沿用)：

S：风险资产的价格(以简单为计，本文都认为是股票价格)；

D：期权在有效期内标的资产的红利(或回报率)；

t：时间变量，且 $t>0$ ；

r：投资的无风险利率；

σ：标的资产的价格波动率；

V：欧式看涨期权的价格；

$V_t^{\left(\alpha \right)}\left(S,t\right)$ ： $V\left(S,t\right)$ 对t的Riemann-Liouville分数阶导数，且 $0<\alpha \le 1$ 。

定理(时间分数阶Black-Scholes方程)假设股票价格S的运动过程为白噪声下的分数阶布朗运动，则：

$\text{d}S=rS\text{d}t+\sigma Sw\left(t\right){\left(\text{d}t\right)}^{\alpha /2}.$ (1)

式中，r为无风险利率(令股票价格的漂移率 $\mu =r$ 以便于计算)，σ为期权价格波动率， $w\left(t\right)$ 是与分数阶布朗运动 $\left\{B_H\left(t\right)=B\left(t,\alpha /2\right),0<\alpha \le 1\right\}$ 相关的Guassian白噪声。设 $V\left(S,t\right)$ 是关于股票价格S和时间t的欧式看涨期权价格函数，且 $V\left(S,t\right)$ 关于t有α阶导数，关于S可导。则有：

$V_t^{\left(\alpha \right)}\left(S,t\right)+\left(\left(r-D\right)S{V}_S-rV\right)\frac{t^{1-\alpha }}{\Gamma \left(2-\alpha \right)}+\frac{\Gamma \left(1+\alpha \right)}{2}{\sigma }^2{S}^2{V}_{SS}=0.$

证明

$V\left(S,t\right)$ 作泰勒展开，有

$\text{d}V=\frac{1}{\Gamma \left(1+\alpha \right)}\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}{\left(\text{d}t\right)}^{\alpha }+\frac{\partial V}{\partial S}\text{d}S+\frac{1}{2}\frac{{\partial }^{2}V}{{\partial }^{2}S}{\left(\text{d}S\right)}^2+o\left({\left(\text{d}t\right)}^{\alpha }\text{d}S\right).$ (2)

式中，

$\begin{array}{l}{\left(\text{d}S\right)}^2={\left(rS\text{d}t+\sigma Sw\left(t\right){\left(\text{d}t\right)}^{\alpha /2}\right)}^2\\ =r^2{S}^2{\left(\text{d}t\right)}^2+2r\sigma S^{2}w\left(t\right){\left(\text{d}t\right)}^{1+\alpha /2}+{\sigma }^2{S}^2{w}^2\left(t\right){\left(\text{d}t\right)}^{\alpha }.\end{array}$ (3)

当 $\text{d}t\to 0$ 时，忽略等式右边阶数高于 $\text{d}t$ 的项，即有

${\left(\text{d}S\right)}^2={\sigma }^2{S}^2{w}^2\left(t\right){\left(\text{d}t\right)}^{\alpha }.$ (4)

将(1)和(4)代入(2)，则

$\text{d}V=\frac{1}{\Gamma \left(1+\alpha \right)}\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}{\left(\text{d}t\right)}^{\alpha }+\frac{\partial V}{\partial S}\left(rS\text{d}t+\sigma Sw\left(t\right){\left(\text{d}t\right)}^{\alpha /2}\right)+\frac{1}{2}{\sigma }^2{S}^2{w}^2\left(t\right)\frac{{\partial }^{2}V}{{\partial }^{2}S}{\left(\text{d}t\right)}^{\alpha }.$ (5)

对上式取均值，得到

$\text{d}V=\left(\frac{1}{\Gamma \left(1+\alpha \right)}\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}+\frac{1}{2}{\sigma }^2{S}^2\frac{{\partial }^{2}V}{{\partial }^{2}S}\right){\left(\text{d}t\right)}^{\alpha }+rS\frac{\partial V}{\partial S}\text{d}t.$ (6)

利用∆对冲建立相应的期权定价连续模型及计算公式。

形成投资组合

$\Pi =V-\Delta S$

选取∆，使得 $\Pi$ 在 $\left(t,t+\text{d}t\right)$ 内是无风险的。

设在时刻t形成投资组合，并在时段 $\left(t,t+\text{d}t\right)$ 内不改变份额∆。那么由于 $\Pi$ 是无风险的，因此在时刻 $t+\text{d}t$ ，投资组合的回报是

$\frac{{\Pi }_{t+\text{d}t}-{\Pi }_t}{{\Pi }_t}=r\text{d}t$

由于考虑到支付股息，因此 ${\Pi }_{t+\text{d}t}=V_{t+\text{d}t}-{\Delta }_t{S}_{t}D\text{d}t-\Delta S_{t+\text{d}t}$ ，代入上式，有

$\text{d}{V}_t-{\Delta }_t\text{d}{S}_t=r_t{\Pi }_t\text{d}t+{\Delta }_t{S}_{t}D\text{d}t$

利用伊藤公式，同时取 ${\Delta }_t=\frac{\partial V}{\partial S}$ ，得到

$\text{d}V=\left(rV-\left(r-D\right)S\frac{\partial V}{\partial S}\right)\text{d}t+\frac{\partial V}{\partial S}\text{d}S$ (7)

联立(6)和(7)整理可得考虑红利情况下，时间分数阶Black-Scholes微分方程：

$V_t^{\left(\alpha \right)}+\frac{\Gamma \left(1+\alpha \right)}{2}{\sigma }^2{S}^2{V}_{SS}+\left(\left(r-D\right)S{V}_S-rV\right)\frac{t^{1-\alpha }}{\Gamma \left(2-\alpha \right)}=0,0<S<S^{\ast }\left(t\right),0<t<T.$

美式期权定价问题的数学模型一般可归结为抛物型方程自由边值问题。自由边界问题 [13] 是指微分方程求解区域的某些边界是未知的(甚至是随时间变化的)，需要和微分方程的解共同确定。Holmes [14] 和Yang [15] 提出有限元方法，利用Front-fixing变换，将自由边界消除，得到一个矩形求解域，由此化简计算。但一旦网格剖分极其细密，也会致使运算量增加，进而出现计算数据溢出计算机内存的情况，降低计算效率。Janson和Tyskc [16] 提出神经网络法，这种方法是将回归方法和神经网络相结合，在某些情况下会出现较大误差或不符合实际情况的结果。Muthuraman [17] 提出移动边界法，只要求解一个通常的线性初边值问题，用新的初始边界代替自由边界，重复上述步骤即可求出自由边界并解决期权问题。张铁和祝丹梅 [18] 提出变网格差分方法，首先建立一个自由边界所满足的方程，再利用变网格技术求解此问题。Khaliq，Voss和Kazm [19] [20] 提出惩罚函数法，但这种方法的缺点是无法求出自由边界。此外还有蒙特卡罗模拟方法 [21] ，紧致差分算法 [22] [23] 等许多算法 [24] - [29] 。

至今对时间分数阶BS方程数值解的探讨不多，SongLina等 [30] (2013)结合欧式看跌期权提出时间分数阶BS方程的隐式差分算法，2017年，张瑜等 [31] 提出时间分数阶BS方程的奇数层显性差分与偶数层隐性差分交替的并行差分的方法。

假设市场是风险中性的，在期权有效期内有红利支付。此时，美式看涨期权的分数阶BS方程为：

$V_t^{\left(\alpha \right)}+\frac{\Gamma \left(1+\alpha \right)}{2}{\sigma }^2{S}^2{V}_{SS}+\left(\left(r-D\right)S{V}_S-rV\right)\frac{t^{1-\alpha }}{\Gamma \left(2-\alpha \right)}=0,0<S<S^{\ast }\left(t\right),0<t<T.$ (8)

其边界条件为：

$V\left(S,T\right)=\max \left(S-E,0\right),0<S\le S^{\ast }\left(T\right),$ (9)

$V\left(0,t\right)=0,0\le t<T,$ (10)

$V\left(S^{\ast }\left(t\right),t\right)=S^{\ast }\left(t\right)-E,V_S\left(S^{\ast }\left(t\right),t\right)=1,0\le t<T,$ (11)

$V\left(S,t\right)=S-E,\text{}{S}^{\ast }\left(t\right)\le S<+\infty .$ (12)

其中T为期权到期日，E为期权的执行价格， $S^{\ast }\left(t\right)$ 是期权提前执行的自由边界，在期权执行日有 $S^{*}\left(T\right)=\max \left(E,\frac{r}{D}E\right)$ 。

分数阶美式期权边界条件的物理意义体现在投资者在期权价格、时间和市场价格等方面做出行权或继续持有的决策，反映了投资者对期权价值和市场情况的判断，对于期权定价和投资决策具有重要的参考价值。为了便于构造数值方法，首先对问题进行变换，目的是将变系数方程化为常系数方程，将反向时间问题化为正向时间问题，引进变量变换

$\begin{array}{l}\tau =\frac{1}{2}{\sigma }^2\left(T-t\right),\\ x=\ln \left(\frac{S}{E}\right)+k\tau ,x_f\left(\tau \right)=\ln \left(\frac{S^{*}\left(t\right)}{E}\right)+k\tau ,k=2\left(r-D\right)/{\sigma }^2-1,\\ V\left(S,t\right)=Eu\left(x,\tau \right)+S-E.\end{array}$

在此变换下，自由边界问题(8)~(12)问题等同于如下问题：

$u_{\tau }^{\alpha }+W\left(\tau \right)u_x-h\left(\tau \right)u_{xx}+f\left(\tau \right)u=g\left(x,\tau \right),0<\tau \le T_1,-\infty <x<x_f\left(\tau \right),$ (13)

$u\left(x,0\right)=\max \left(1-{\text{e}}^x,0\right),-\infty <x<x_f\left(0\right),$ (14)

$\underset{x\to -\infty }{\lim }u\left(x,\tau \right)=1-{\mathrm{e}}^{x-k\tau },\text{}0<\tau \le T_1,$ (15)

$u\left(x_f\left(\tau \right),\tau \right)=0,u_x\left(x_f\left(\tau \right),\tau \right)=0,0<\tau \le T_1,$ (16)

$u\left(x,\tau \right)=0,x_f\left(\tau \right)\le x<+\infty .$ (17)

其中，

$\begin{array}{l}h\left(\tau \right)=\frac{\Gamma \left(1+\alpha \right){\tau }^{1-\alpha }}{{\left(T-\frac{2\tau }{\sigma }\right)}^{1-\alpha }},f\left(\tau \right)=\frac{2r{\tau }^{1-\alpha }}{{\sigma }^2\Gamma \left(2-\alpha \right)},W\left(\tau \right)=h\left(\tau \right)-\frac{2\left(r-D\right){\tau }^{1-\alpha }}{E{\sigma }^2\Gamma \left(2-\alpha \right)},\\ g\left(x,\tau \right)=\frac{\left[2r+\left(2r-4D-{\sigma }^2\right){\text{e}}^{x-k\tau }\right]{\tau }^{1-\alpha }}{{\sigma }^2\Gamma \left(2-\alpha \right)},T_1=\frac{{\sigma }^2}{2}T.\end{array}$

此时 $x_f\left(\tau \right)$ 为新的自由边界，并且 $x_f\left(0\right)=\max \left(\ln \left(r/D\right),0\right)$ 。另外，为方便计算将公式(13)、(14)中的左边界条件限定为 $x=-L$ 。

3. 分数阶期权定价方程的数值算法

3.1. 方程解的差分近似

由于B-S方程在时间上是分数阶求导而在空间上是整数阶求导，在此需要用不同的离散化方法。

3.1.1. 时间与空间上的差分近似

方程(13~17)数值求解的困难在于必须同时确定持有区域的函数值以及自由边界的位置，而具有特定性质和边界条件的分数阶抛物型偏微分方程，可以通过适当的离散化和迭代求解，由此得出分数阶抛物型偏微分方程的数值解。同时，由于BS方程关于空间仅有一维，有限元和差分法等效。差分方法是通过差商代替微商对方程以及定解问题离散化。通过建立与偏微分方程相应的差分方程有多种方式，从求解的方式可以分为两大类：一类是显式差分格式，求解的过程是显式，通过直接运算求出它的值；另一种是隐式差分格式，求解过程必须通过一个代数方程组才能得到它的值 [32] 。

众所周知，如果 $y=f\left(x\right)$ 充分光滑，那么 $f^{\prime }\left(x\right)$ 和 $f^{″}\left(x\right)$ 可以用以下差商形式代替 [33] ：

$f^{\prime }\left(x\right)~\frac{\Delta f}{\Delta x}=\frac{f\left(x+\Delta x\right)-f\left(x-\Delta x\right)}{2\Delta x},f^{″}\left(x\right)~\frac{{\Delta }^{2}f}{\Delta x^2}=\frac{f\left(x+\Delta x\right)-2f\left(x\right)+f\left(x-\Delta x\right)}{\Delta x^2}.$

而对于时间分数阶期权定价模型，其中分数阶导数使用Riemann-Liouville定义 [33] ，时间分数阶导数的离散格式为 [34] [35] ：

$\frac{\partial u\left(x_i,{\tau }_{m+1}\right)}{\partial {\tau }^{\alpha }}=\frac{{\left(d\tau \right)}^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{j=1}{\overset{m}{\sum }}\left[u\left(x_i,{\tau }_{m+2-j}\right)-u\left(x_i,{\tau }_{m+1-j}\right)\right]}\left[j^{1-\alpha }-{\left(j-1\right)}^{1-\alpha }\right]+O\left(k\right),$

忽略误差，可以得到

$\frac{\partial u\left(x_i,{\tau }_{m+1}\right)}{\partial {\tau }^{\alpha }}=\frac{{\left(d\tau \right)}^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left(u\left(x_i,{\tau }_{m+1}\right)-{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}{d}_{j}u\left(x_i,{\tau }_{m+1-j}\right)-l_{m}u}\left(x_i,{\tau }_1\right)\right).$

其中， $l_j=j^{1-\alpha }-{\left(j-1\right)}^{1-\alpha },d_j=2j^{1-\alpha }-{\left(j+1\right)}^{1-\alpha }-{\left(j-1\right)}^{1-\alpha }.$

将求解区域划分为均匀网格，将时间区间 $\left[0,T_1\right]$ 划分成M个间隔相等的小区间，记

${\tau }_m=\left(m-1\right)d\tau ,m=1,2,\cdots ,M,M+1.$

将区间 $\left[-L,x_f\left(0\right)\right]$ 划分成N个间隔相等的区间段，

$x_i=-L+\left(i-1\right)dx,i=1,2,\cdots ,N,N+1.$

不妨记 $u_i^m=u\left(x_i,{\tau }_m\right)$ ，边界条件转化为：

$u_i^0=\max \left(1-{\text{e}}^{-L+\left(i-1\right)dx},0\right),u_0^m=1-\exp \left(-L-k\left(m-1\right)d\tau \right),u_N^m=0.$

建立问题(6)~(10)的古典显式格式(见图1)：

$\frac{\partial u\left(x_i,{\tau }_{m+1}\right)}{\partial {\tau }^{\alpha }}=-W\left({\tau }_m\right)\frac{u_{i+1}^m-u_{i-1}^m}{2dx}+h\left({\tau }_m\right)\frac{u_{i+1}^m-2u_i^m+u_{i-1}^m}{{\left(dx\right)}^2}-f\left({\tau }_m\right)u_i^m+g_i^m.$

古典隐式格式(见图2)：

$\frac{\partial u\left(x_i,{\tau }_{m+1}\right)}{\partial {\tau }^{\alpha }}=-W\left({\tau }_{m+1}\right)\frac{u_{i+1}^{m+1}-u_{i-1}^{m+1}}{2dx}+h\left({\tau }_{m+1}\right)\frac{u_{i+1}^{m+1}-2u_i^{m+1}+u_{i-1}^{m+1}}{{\left(dx\right)}^2}-f\left({\tau }_{m+1}\right)u_i^{m+1}+g_i^{m+1}.$

Figure 1. Fractional order explicit difference scheme

图1. 分数阶显式差分格式

Figure 2. Fractional order implicit difference scheme

图2. 分数阶隐式差分格式

3.1.2. θ-差分格式

由于显示差分的条件稳定性和隐式差分的绝对稳定性，在此取加权差分格式 [36] [37] ：

$\begin{array}{l}\frac{{\left(d\tau \right)}^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\left[u\left(x_i,{\tau }_{m+2-j}\right)-u\left(x_i,{\tau }_{m+1-j}\right)\right]}\left[j^{1-\alpha }-{\left(j-1\right)}^{1-\alpha }\right]\\ =\left(1-\theta \right)\left[-W\left({\tau }_m\right)\frac{u_{i+1}^m-u_{i-1}^m}{2dx}+h\left({\tau }_m\right)\frac{u_{i+1}^m-2u_i^m+u_{i-1}^m}{{\left(dx\right)}^2}-f\left({\tau }_m\right)u_i^m+g_i^m\right]\\ +\theta \left[-W\left({\tau }_{m+1}\right)\frac{u_{i+1}^{m+1}-u_{i-1}^{m+1}}{2dx}+h\left({\tau }_{m+1}\right)\frac{u_{i+1}^{m+1}-2u_i^{m+1}+u_{i-1}^{m+1}}{{\left(dx\right)}^2}-f\left({\tau }_{m+1}\right)u_i^{m+1}+g_i^{m+1}\right],\end{array}$

其中 $0\le \theta \le 1$ ，当 $\theta =1/2$ 时，上式即为分数阶BS方程的C-N格式。整理得：

$\begin{array}{l}\left[-\theta \left(B_{m+1}-A_{m+1}\right)\right]u_{i+1}^{m+1}+\left[\theta \left(2B_{m+1}+C_{m+1}\right)+1\right]u_i^{m+1}-\theta \left(A_{m+1}+B_{m+1}\right)u_{i-1}^{m+1}\\ =\left(1-\theta \right)\left(B_m-A_m\right)u_{i+1}^m-\left(1-\theta \right)\left(2B_m+C_m\right)u_i^m+\left(1-\theta \right)\left(A_m+B_m\right)u_{i-1}^m\\ +\theta R{g}_i^{m+1}+{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}{d}_j{u}_i^{m+1-j}+l_m{u}_i^1}+\left(1-\theta \right)R{g}_i^m,\end{array}$

即

$\begin{array}{l}{c}_{m+1}{u}_{i+1}^{m+1}+a_{m+1}{u}_{i-1}^{m+1}+b_{m+1}{u}_i^{m+1}\\ ={c^{\prime }}_m{u}_{i+1}^m+{a^{\prime }}_m{u}_{i-1}^m+{b^{\prime }}_m{u}_i^m+{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}{d}_j{u}_i^{m+1-j}+u_i^1{l}_m}+\left(1-\theta \right)R{g}_i^m+\theta R{g}_i^{m+1}.\end{array}$

其矩阵形式的具体表达式

$\begin{array}{l}\left(\begin{array}{ccccc}{b}_{m+1}& c_{m+1}& & & \\ a_{m+1}& b_{m+1}& \ddots & & \\ & \ddots & \ddots & c_{m+1}& \\ & & a_{m+1}& b_{m+1}& c_{m+1}\\ & & & a_{m+1}& b_{m+1}\end{array}\right)\left(\begin{array}{c}{U}_2^{m+1}\\ U_3^{m+1}\\ ⋮\\ U_{N-1}^{m+1}\\ U_N^{m+1}\end{array}\right)=\left(\begin{array}{ccccc}{b^{\prime }}_m& {c^{\prime }}_m& & & \\ {a^{\prime }}_m& {b^{\prime }}_m& \ddots & & \\ & \ddots & \ddots & {c^{\prime }}_m& \\ & & {a^{\prime }}_m& {b^{\prime }}_m& {c^{\prime }}_m\\ & & & {a^{\prime }}_m& {b^{\prime }}_m\end{array}\right)\left(\begin{array}{c}{U}_2^m\\ U_3^m\\ ⋮\\ U_{N-1}^m\\ U_N^m\end{array}\right)\\ +\left(\begin{array}{c}{a^{\prime }}_m{U}_1^m-a_{m+1}{U}_1^{m+1}\\ 0\\ ⋮\\ 0\\ {c^{\prime }}_m{U}_{N+1}^m-c_{m+1}{U}_{N+!}^{m+1}\end{array}\right)+l_m\left(\begin{array}{c}{U}_2^1\\ U_3^1\\ ⋮\\ U_{N-1}^1\\ U_N^1\end{array}\right)+{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}{d}_j{U}^{m+1-j}}+R\left(1-\theta \right)\left(\begin{array}{c}{g}_2^m\\ g_3^m\\ ⋮\\ g_{N-1}^m\\ g_N^m\end{array}\right)+R\theta \left(\begin{array}{c}{g}_2^{m+1}\\ g_3^{m+1}\\ ⋮\\ g_{N-1}^{m+1}\\ g_N^{m+1}\end{array}\right),\end{array}$ ${\displaystyle \underset{j=1}{\overset{m-1}{\sum }}{d}_j\left(\begin{array}{c}{U}_2^{m+1-j}\\ U_3^{m+1-j}\\ ⋮\\ U_{N-1}^{m+1-j}\\ U_N^{m+1-j}\end{array}\right)}=d_1\left(\begin{array}{c}{U}_2^m\\ U_3^m\\ ⋮\\ U_{N-1}^m\\ U_N^m\end{array}\right)+d_2\left(\begin{array}{c}{U}_2^{m-1}\\ U_3^{m-1}\\ ⋮\\ U_{N-1}^{m-1}\\ U_N^{m-1}\end{array}\right)+d_3\left(\begin{array}{c}{U}_2^{m-2}\\ U_3^{m-2}\\ ⋮\\ U_{N-1}^{m-2}\\ U_N^{m-2}\end{array}\right)+\cdots +d_{m-1}\left(\begin{array}{c}{U}_2^2\\ U_3^2\\ ⋮\\ U_{N-1}^2\\ U_N^2\end{array}\right).$

其中，

$\begin{array}{l}R=\Gamma \left(2-\alpha \right){\left(d\tau \right)}^{\alpha },A_m=\frac{RW\left({\tau }_m\right)}{2dx},B_m=\frac{Rh\left({\tau }_m\right)}{{\left(dx\right)}^2},C_m=Rf\left({\tau }_m\right),g_i^m=g\left(x_i,{\tau }_m\right).\\ c_{m+1}=-\theta \left(B_{m+1}-A_{m+1}\right),a_{m+1}=-\theta \left(A_{m+1}+B_{m+1}\right),b_{m+1}=\theta \left(C_{m+1}+2B_{m+1}\right)+1.\\ {c^{\prime }}_m=\left(1-\theta \right)\left(B_m-A_m\right),{a^{\prime }}_m=\left(1-\theta \right)\left(A_m+B_m\right),{b^{\prime }}_m=-\left(1-\theta \right)\left(C_m+2B_m\right).\\ l_j=j^{1-\alpha }-{\left(j-1\right)}^{1-\alpha },d_j=2j^{1-\alpha }-{\left(j+1\right)}^{1-\alpha }-{\left(j-1\right)}^{1-\alpha }.\end{array}$

3.2. 自由边界的差分近似

自由边界的确定可以看作是分数阶期权定价模型在边界条件下的最优决策问题。在分数阶抛物型方程中，自由边界的确定涉及到期权价格在边界上的数值解。通过求解分数阶抛物型方程，可以得到期权价格在不同时间和价格水平下的数值解，从而确定自由边界。关于自由边界的数值求解产生了许多数值方法，包括焓方法(E. L. Albasiny, 1956)，前沿固定法，有限差分法(J. Crank, 1957)和变网格法等。本文主要使用前沿追踪法，建立了变网格，在每一个时间步长，离边界远处用正常的有限差分，在自由边界附近将差分公式进行修正，以适应不均匀的空变化的自由边界。上一节给出了期权定价方程在持有区域内部的离散算法，本节用“预估–校正”的方法同时确定自由边界及边界附近的函数值。

3.2.1. 预估

由于函数值和自由边界的位置必须同时确定，因此利用u在τ时刻各点的值来确定自由边界的位置 $x_f\left(\tau \right)$ 。

设第m行值已知，即自由边界 $x_f\left({\tau }_m\right)$ 的值已知，现确定m + 1行函数的值及自由边界的位置。首先假设 $x_f\left({\tau }_{m+1}\right)=x_f\left({\tau }_m\right)$ 由上节的θ-差分格式或显–隐差分格式可以求出执行区域上的函数值，在此建立方程估计 $x_f\left(\tau \right)$ 的值。已知 $x_{N-1}\le x_f\left(\tau \right)$ 是邻近 $x_f\left(\tau \right)$ 的一个网格节点，设 $x_f\left(\tau \right)=x_{N-1}+qdx$ ，其中dx为网格步长。我们要确定 $x_f\left(\tau \right)$ 的值，只需要求出q。由泰勒公式得，

$u_{N-1}=u_{x_{{}_f}}-qdx{\left(\frac{\partial u}{\partial x}\right)}_{x_{{}_f}}+\frac{1}{2}{\left(qdx\right)}^2{\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{x_{{}_f}}+O\left({\left(qdx\right)}^3\right).$

利用边界条件(16)，且舍去高阶小项 $O\left({\left(dx\right)}^3\right)$ ，得到近似方程

$u_{N-1}-\frac{1}{2}{\left(qdx\right)}^2{\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{x_{{}_f}}=0.$ (18)

再次利用(16)式得到

$0={\left(\frac{\text{d}u}{\text{d}\tau }\right)}_{x_f}={\left(\frac{\partial u}{\partial \tau }\right)}_{x_f}+{\left(\frac{\partial u}{\partial x}\right)}_{x_f}\frac{\text{d}{x}_f}{\text{d}\tau }={\left(\frac{\partial u}{\partial \tau }\right)}_{x_f}.$

故由方程(13)得

$-h\left(\tau \right){\left(u_{xx}\right)}_{x_f}=g\left(x_f,\tau \right),$

代入(18)式有

$u_{N-1}+\frac{1}{2}{\left(qdx\right)}^2\frac{g\left(x_f,\tau \right)}{h\left(\tau \right)}=0.$

因此可得q所满足的非线性方程，

$u_{N-1}+\frac{q^2{\left(dx\right)}^2\left[2r+\left(2r-4D-{\sigma }^2\right){\text{e}}^{x_{N-1}+qdx-k\tau }\right]{\left(T-\frac{2\tau }{{\sigma }^2}\right)}^{1-\alpha }}{2{\sigma }^2\Gamma \left(1+\alpha \right)\Gamma \left(2-\alpha \right)}=0.$ (19)

在方程(19)中解q存在唯一的基础上，我们可利用牛顿法计算q：

$\begin{array}{l}{q}_{m+1}=q_m-\frac{H\left(q_m\right)}{H_1\left(q_m\right)},\\ \frac{H\left(q_m\right)}{H_1\left(q_m\right)}=\frac{u_{N_m-1}^{m+1}}{H_1\left(q_m\right)}-\frac{2q_{m}r+\left(2r-4D-{\sigma }^2\right)q_m{\text{e}}^{x_{N_m-1}+q_{m}dx-k\tau }}{4r+\left(4r-8D-2{\sigma }^2+2qdx-4D{q}_{m}dx-q_m{\sigma }^{2}dx\right){\text{e}}^{x_{N_m-1}+q_{m}dx-k\tau }}.\end{array}$

其中，

$\begin{array}{l}H\left(q\right)=u_{N-1}+\frac{q^2{\left(dx\right)}^2\left[2r+\left(2r-4D-{\sigma }^2\right){\text{e}}^{x_{N-1}+qh-k\tau }\right]{\left(T-\frac{2\tau }{{\sigma }^2}\right)}^{1-\alpha }}{2{\sigma }^2\Gamma \left(1+\alpha \right)\Gamma \left(2-\alpha \right)},\\ H_1\left(q\right)=H^{\prime }\left(q\right)=\frac{q{\left(dx\right)}^2\left[4r+\left(4r-8D-2{\sigma }^2+2qdx-4Dqdx-q{\sigma }^{2}dx\right){\text{e}}^{x_{N-1}+qh-k\tau }\right]}{2{\sigma }^2\Gamma \left(1+\alpha \right)\Gamma \left(2-\alpha \right){\left(T-\frac{2\tau }{{\sigma }^2}\right)}^{\alpha -1}}.\end{array}$

由于执行区域内 $u_{N-1}\ne 0$ ，由牛顿迭代可以由方程(19)求出q，即求出 $x_f\left({\tau }_{m+1}\right)$ 。

如果 $0\le q<1$ ，则取 $N_{m+1}=N_m$ ，m + l步计算完成；如果 $q\ge 1$ ，则需扩充网格，设

$N_{m+1}=N_m+\left[q\right],$

由上节中的差分格式补充计算节点 $u_i^{m+1},i=N_m,\cdots ,N_{m+1}-1$ ，从而完成m + 1步计算。在计算中用到 $U_i^m,i>N_m$ 的值时，可根据边界条件设

$U_i^m=0,i\ge N_m,m=0,1,\cdots ,M.$

下图(见图3)是从到 ${\tau }_m$ 到 ${\tau }_{m+1}$ 步网格扩充示意图，其中黑色小方框点表示扩充的节点。

Figure 3. Grid expansion diagram

图3. 网格扩充示意图

3.2.2. 校正

当计算第m + 1行 $x_f\left(\tau \right)$ 的值时，上节首先假设 $x_f^1\left({\tau }_{m+1}\right)=x_f\left({\tau }_m\right)$ ，再由方程得到边界的预估值 $x_f^2\left({\tau }_{m+1}\right)=x_{N_m-1}+qdx$ 。

对新的预估值 $x_f^2\left({\tau }_{m+1}\right)$ 重复预估过程，得到数列 $\left\{x_f^k\left({\tau }_{m+1}\right)\right\},k\in Z$ ，其极限即为自由边界的精确值。给定精度ε，当 $\left|x_f^k\left({\tau }_{m+1}\right)-x_f^{k-1}\left({\tau }_{m+1}\right)\right|\le \epsilon$ 时，停止计算，以 $x_f^k\left({\tau }_{m+1}\right)$ 表示精确值。重复以上过程即可得到执行区域上分数阶BS方程的数值解及自由边界。

4. 总结

由于资本市场的本质特征和状态都是随机波动的，与传统BS期权定价模型的假设并不完全吻合，而基于次分数Brown运动的BS方程能够更好地刻画原生标的资产价格的变动趋势，基于修正后的布朗运动能够更多地反映自相关性、长期记忆性和增量相关性等众多性质，最终得出时间分数阶BS方程。

本文研究了一种包含分数阶微分方程的美式期权定价方程。该模型更适用于实际金融市场是相对易于求解的期权定价模型。本文给出了该模型在持有区域上的θ-差分格式，其中离散化的精度和稳定性可以通过选择合适的网格大小和差分近似方法来提高。此外，根据该格式得到了分数阶BS微分方程数值解的算法。

由于自由边界的问题是一个非线性问题，同样在一般情况下也不容易得不到其显式表达式。本文根据自由边界所满足的条件，通过构建求解美式期权定价自由边值问题的预估–校正网格差分方法，同时求出期权的最佳执行边界及分数阶BS方程的数值解。

基金项目

国家自然科学基金面上项目支持：NSFC12071215。