第九章-期权定价的有限差分方法

1 第九章第九章 期权定价的有限差分方法期权定价的有限差分方法 在本章中 我们将给出几个简单的例子来说明基于偏微分方程 PDE 框架的期权定价方法 具体的方法的是利用第五章中讲述的 有限差分方法来解决 Black scholes 偏微分方程 在 9 1 节中 我们 会回顾衍生品定价的数值解法以及指出如何利用适当的边界条件来 模拟一个特定的期权 在 9 2 节中我们将会应用简单的显式 差分 方法来求解一个简单的欧式期权 正如你已熟知的那样 这种方法 只能解出一些可以从金融角度来解释的不稳定的数值解 在 9 3 节 中我们将可以看到使用完全的隐式方法可以解决这种不稳定问题 在 9 4 节中我们将介绍 Crank Nicolson 方法在障碍期权定价中的应 用 它可以看做是一种显式与完全隐式方法的混合 最后 在 9 5 节中 我们会看到迭代松弛方法可以用于解决使用全隐式方法来解 决美式期权定价时由于存在提前执行的可能性而导致的自由边界问 题 9 19 1 使用有限差分法解使用有限差分法解 BSBS 方程方程 在 2 6 2 节中 我们给出了一个标的资产在时间 的价格为t 的期权 该期权的价格是一个函数 且满足偏微分 tS tSf tSf 方程 9 1 通过不同的边界条件可以让这个方程刻画不同的期权的特征 在某些地方可能因为假设的改变或者对路径依赖的改变而导致方程 式的具体形式改变 但是此处仅仅作为一个起点 帮助读者了解如 何应用基于有限差分方法来解决期权定价的问题 2 正如我们在第五章中遇到的情况那样 要用有限差分方法来解 偏微分方程 在此处我们必须建立资产价格和时间的离散网格 设 T 是期权的到期日 而 Smax 是一个足够大的资产价格 在我们所考 虑的时间范围内 的数值不能超过 Smax 设定 Smax 是因为偏微 tS 分方程的区域关于资产价格是无边界的 但是为了达到计算的目的 必须要求它是有界的 Smax 相当于 网格通过点 S t 取得 其 中 S t 满足 0 SS S 2maxSSM 2 0 tt t 2TtN 本章中使用网格符号为 我们回顾一下 9 1 方 程式的几种不同解法 向前差分 向后差分 中心 或对称 差分 对于第二个差分式子 有 至于究竟采用哪种方法进行离散化 我们将在后面的实际操作 3 过程中对显式和隐式的方法作出详细的阐述说明 另一个值得我们注意的问题是如何设置边界条件 对于执行价 格为的欧式看涨期权 其终值条件为K 对于执行价格为 K 的欧式看跌期权 其终值条件为 但是当我们涉及到资产价格的边界问题时 这个问题就变得复 杂了 因为在数学上要求必须在一个有界区域里来解这个方程 而 资产价格的区域是无界的 我们可以通过一些例子来说明这个问题 例例 9 19 1 我们首先考虑一支简单的欧式看跌期权 当资产的价 格非常大时 看跌期权几乎是一文不值 因为我们几乎可以肯定 tS 的是其执行价低于现价 Smax 的值必须足够大才能使得假设的边界条件是合理的 当资 产的价格时 鉴于资产的价格符合几何布朗运动的动力学模0 tS 型 我们可以认为资产的价格会保持为 0 因此到期收益为 K 在时 间 贴现值为t 其边界条件可以写为 例例 9 29 2 我们可以从 9 1 推导出处理普通欧式看涨期权的方法 不妨设资产的价格 那么在任何时刻 我们都可以知道该期0 tSt 4 权价值为 0 而对于一个足够大的资产价格 我们可以肯定的是该期权一 tS 定为价内期权 即具有大于零的价值的期权 那么期满行权我们可 以得到收益为 将终值贴现到 时刻 并且考虑标的资产KTS Kt 的价格是一个函数 那么我们根据无套利定价原理可以得到其边 tS 界条件如下 具体条件可以写为 对于资产价格 较大时 需要一个替代的边界条件 但这时候S 必须要求期权的 为 1 在这种情况下 我们所提出的边界条件一 般针对未知函数的倒数 而非函数本身 这便是数学物理方程中所 谓的 Neumann 边界条件 我们不会采取这种做法 因为这将使得数 值解法复杂化了 当我们面对障碍式期权时 事情可能就变得简单了 对敲出期 权来说 例如一个向下触及失效看跌期权 在标的资产价格触碰到 障碍时期权作废 价值为 0 向上触及失效看涨期权的情况也是类 似的 它们的优点都在于我们计算必须涉及到的边界即域是自然有 界的 显然存在的 美式期权的定价就显得更加复杂难以处理了 因为这涉及到了提前执行边界的问题 我们必须考虑什么水平的资 产价格和在什么时间是期权的最佳实施点 因此 在我们解决问题 的过程中必然产生一个自由边界 这就需要为不同的奇异期权设定 不同的边界条件 如何计算出正确的边界条件以及使用数值方法逼 5 近正确的结果 这依赖于对具体期权的选择 利用显式法求解欧式期权利用显式法求解欧式期权 由于初次尝试解方程 9 1 我们不妨考虑一个简单的欧式看 跌期权 我们分别用中心差分对 S 进行求导 用向后差分对时间进 行求导来逼近求解 这并不是唯一的方法 但是无论采取什么样的 方法都必须保证边界条件是符合特定期权模型的 例 9 1 中的边界 条件差分后的结果可以用以下的方程式来表示 9 2 需要特别指出说明的是 由于存在终点条件的假设 因此方程 必须适当地用回推法求解 不妨设 9 2 式中有 j N 由于存在给 定的终点条件 我们有一个未知量 可以认为它是关于三个已 1i N f 知量的函数 对其他的每个时间层 i 作同样的考虑 整理方程 我们可以得到 9 3 此时有 7 function price EuPutExpl S0 K r T sigma Smax dS dt 建立网格 并且在必要时调整增量 M round Smax dS dS Smax M N round T dt dt T N matval zeros M 1 N 1 vetS linspace 0 Smax M 1 veti 0 M vetj 0 N 建立边界条件 matval N 1 max K vetS 0 matval 1 K exp r dt N vetj matval M 1 0 建立三对角矩阵 a 0 5 dt sigma 2 veti r veti b 1 dt sigma 2 veti 2 r c 0 5 dt sigma 2 veti r veti 求解方程 for j N 1 1 for i 2 M matval i j a i matval i 1 j 1 b i matval i j 1 8 c i matval i 1 j 1 end end 返回价格 有可能在网格外线性插值生成 price interp1 vetS matval 1 S0 图图 9 19 1 简易欧式期权显式差分方法定价的简易欧式期权显式差分方法定价的 MATLABMATLAB 代码代码 定价不是直接利用 MATLAB 来完成的 图 9 1 中是定价所使用的 MATLAB 的程序代码 另外计算还需要给出 Smax 的值以及上文提到 的两个离散步骤 唯一一点需要注意的是 在数学符号中幂指数可 以从 0 开始 但是在 MATLAB 的矩阵计算中 幂指数是从 1 开始的 另外 如果资产的初始价格不在网格上 那么我们只能够寻找两个 相邻点进行插值 一般我们采取比较粗略的线性插值 但是如果有 需要计算分析期权价格的敏感性 在实际案例中经常需要这样做 则需要进行更加复杂的插值计算 c p blsprice 50 50 0 1 5 12 0 4 P P 4 076 EuPutExpl 50 50 0 1 5 12 0 4 100 2 5 1200 ans 4 0669 c p blsprice 50 50 0 1 5 12 0 3 P P 2 8446 EuPutExpl 50 50 0 1 5 12 0 3 100 2 5 1200 9 ans 2 8288 我们可以看到 数值方法能够提供比较精确的结果 为了提高 精度 我们也可以采用更加精细的网格进行计算 EuPutExpl 50 50 0 1 5 12 0 3 100 1 5 5 1200 ans 2 8597 EuPutExpl1 50 50 0 1 5 12 0 3 100 1 5 1200 ans 2 8271e 022 在此处我们看到在第五章中分析过的不稳定数值解的例子 有 一种解决方案就是使用隐式方法进行求解 另外一种方法则是进行 稳定性分析 然后在从离散化过程中导出边界 此处我们不会再使 用第二种方案 这与第五章中对热传导方程进行的简易变换处理方 法十分相似 接着 本章在下一节中将从金融的角度来阐释不稳定 性 这暗示我们需要通过改变变量来改变方程式 9 29 2 显式法不稳定性的金融学解释显式法不稳定性的金融学解释 i j f 1i j f 1 1ij f a 1i j f 1 1ij f i j f b 图图9 29 2 显式法 显式法 a a 与隐式法 与隐式法 b b 的求解 的求解BSBS方程的直观图示方程的直观图示 在显式法中 我们通过 和 ttSSf ttSf 三个值得到 这看起来与我们在 7 4 这部分 见 ttSSf tSf 图 9 2 a 介绍过的三叉树方法有点类似 这样方便我们从显式 1 1ij f 1 1ij f 10 法推导出新的方法 如下 参考文献 1 中第 18 章 假设在点 i j 和点 i j 1 对 S 求一阶导数和二阶导数是相等的 另外一种方法可以取得同样的效果 用代替 9 2 式右侧的 1 ji f 这种方法引入的误差是有界的 并且这种方法的网格细化至趋 ji f 于零 改变过后的替换后的差分方程变成了 也可写成 当 i 1 2 M 1 且 j 1 2 N 1 时 同时有 这个方法仍然是显式的而且仍然受到不稳定性的影响 系数 ia 以及 本身也对此作出了很好的解释 回顾之前我们使用过的二叉 i b ic 树和三叉树方法 我们是在风险中性的假设下将后继节点对期权的 预期价值进行贴现得到期权的价格 事实上 上述系数 包括 可以看做是长度为的时间间隔的贴现系数 1 1tr t 此外我们还有 这也说明了上述系数等于概率乘以一个折现系数 他们是否是 11 风险中性的概率呢 我们应当首先测试资产价格在时间间隔内增t 长的值 这正是我们所期望的风险中性世界 对于其方差的增加 有 此外 当极小时有 t 这也是服从风险中性假设的几何布朗运动的 因此我们可以认 为排除一些小问题后 显式法是可以和三叉树方法取得同样的效果 的 即 概率 和可能是负数 细心的读者可能会发现一个循 d 0 环 在第五章中我们学习了一个线性组合的系数的稳定性条件 在 热传导方程中 我们必须确定这个组合是一个凸组合 即系数是正 的并且相加等于 1 就像一个概率函数 有一种解决这个问题的方法 就是 参考文献 1 中通过改变 变量来实现的 通过 Z lnS 来代换重写 B S 方程 我们得到了实现 稳定性的简单条件 然而 变量的改变对某些奇异期权来说却也未 必是好事 在下一节中我们将使用隐式法来避免数值不稳定的问题 9 39 3 使用全隐式法对欧式期权定价使用全隐式法对欧式期权定价 为了克服显式法带来的不稳定性问题 我们采用全隐式法 首 对时间 t 采用向前差分求其偏导数 我们得到一个网格方程 也可写成 当 i 1 2 M 1 且 j 1 2 N 1 时 12 此时有 此时我们有需要使用 3 个未知量去计算出一个未知量 见图 9 2b 首先要注意的是 在每一个时间层 我们有 M 1 个方程式和 N 1 个未知量 边界条件决定两个未知量 终点条件也给出了最后 时间层的值 正如显式法的案例那样 我们必须用回推法 解决当 j N 1 1 情况下的所有线性方程组 在时间层 j 得到的方程组 如下 值得我们注意的是 上述的矩阵是个三对角矩阵 并且相对于 时间层 i 来说它是一个不变的常量 因此我们可以使用 LU 因式分解 来解决问题 1 所有工作的 MATLAB 代码在图 9 3 中给出 1 由于该矩阵具有稀疏结构 编写特定的代码来解决这类线性系统更加容易 此处我们只使用现成的 MATLAB 功能 13 function price EuPutImpl S0 K r T sigma Smax dS dt 建立网格并在必要时调整增量 M round Smax dS dS Smax M N round T dt dt T N matval zeros M 1 N 1 vetS linspace 0 Smax M 1 veti 0 M vetj 0 N 建立边界条件 matval N 1 max K vetS 0 matval 1 K exp r dt N vetj matval M 1 0 建立三对角矩阵系数 a 0 5 r dt veti sigma 2 dt veti 2 b 1 sigma 2 dt veti 2 r dt c 0 5 r dt veti sigma 2 dt veti 2 coeff diag a 3 M 1 diag b 2 M diag c 2 M 1 1 L U lu coeff 求解线性方程组 14 aux zeros M 1 1 for j N 1 1 aux 1 a 2 matval 1 j other term from BC is zero matval 2 M j U L matval 2 M j 1 aux end 返回价格 如果初始价格在网格外则可能进行性线性插值 price interp1 vetS matval 1 S0 图图9 39 3 全隐式方法求解简易欧式期权的全隐式方法求解简易欧式期权的MATLABMATLAB代码代码 c p blsprice 50 50 0 1 5 12 0 4 P P 4 0760 EuPutlmpl 50 50 0 1 5 12 0 4 100 0 5 5 2400 ans 4 0718 这个计算结果是非常精确的 这是因为使用了排除不稳定因素 存在的更加精细的网格的结果 另外一种方法就是通过 Crank Nicolson 方法来提高精度 这将在后面的章节中做进一步的介绍 9 49 4 使用使用 Crank NicolsonCrank Nicolson 方法对障碍式期权定价方法对障碍式期权定价 Crank Nicolson 方法已经在 5 3 3 节中提到过 它通过结合使用 显式法和隐式法来提高数值解的精确度 在 B S 模型的求解中应用 此方法可以得到 15 上述方程可以改写为 此时有 我们首先考虑在 2 7 1 中介绍过的向下触及失效期权 假设该 障碍期权是处于连续监测下的 在这种情况下我们需要考虑 S Smax 边界条件是 b S 综合考虑上述边界条件 我们可以将方程 9 5 改写为如下矩 阵形式 同时 有 16 相应的 MATLAB 代码在图 9 4 中列出 这个结果可以同 2 7 1 中通过定价分析公式得到的结果进行比 较 DOPut 50 50 0 1 5 12 0 4 40 ans 0 5424 DOPutCK 50 50 0 1 5 12 0 4 40 100 0 5 1 1200 ans 0 5414 障碍期权还有很多种 可以在 参考文献 9 中找到更多的关 于障碍式期权定价的偏微分方程应用 function price DOPutCK S0 K r T sigma Sb Smax dS dt 建立网格并且在必要的时候调整增量 M round Smax Sb dS dS Smax Sb M N round T dt dt T N 17 matval zeros M 1 N 1 vetS linspace Sb Smax M 1 veti vetS dS vetj 0 N 建立边界条件 matval N 1 max K vetS 0 matval 1 0 matval M 1 0 建立系数矩阵 alpha 0 25 dt sigma 2 veti 2 r veti beta dt 0 5 sigma 2 veti 2 r gamma 0 25 dt sigma 2 veti 2 r veti M1 diag alpha 3 M 1 diag 1 beta 2 M diag gamma 2 M 1 1 L U lu M1 M2 diag alpha 3 M 1 diag 1 beta 2 M diag gamma 2 M 1 1 求解线性方程组 for j N 1 1 matval 2 M j U L M2 matval 2 M j 1 end 返回价格 这个结果可能因为初始资产价格在网格外而有线性 18 插值产生 price interp1 vetS matval 1 S0 图图9 49 4 Crank NicolsonCrank Nicolson方法对向下触及失效期权定价的方法对向下触及失效期权定价的MATLABMATLAB代码 代码 9 59 5 美式期权定价美式期权定价 有限差分方法在对简单的欧式期权进行定价时的确是十分方便 的 但是这在美式期权定价中却显得并不是很实用 我们将这个方 法应用于美式期权定价的时候 精确的公式就不大适用了 这是因 为美式期权定价时必须考虑存在提前行权的可能性 根据无套利原 理 在 S t 空间中任一点的期权的价格不能低于其内在价值 即立 刻行权的所得收益 对于一个一般的美式期权 这意味着 从一个实际操作来看 在显式法中要充分考虑这种情况是很困 难的 但是我们可以做一些修改以便能够简单的应用 9 2 节中提到 的方法 在计算得到后 我们可以检测提前行权的概率 并且可 i j f 以像我们处理二叉树模型那样设 由于存在着不稳定性问题 我们可能更加倾向于使用隐式法 但是这样处理后又导致了另一个复杂的问题 即根据上述关系要求 必须是已知的 而在隐式法中却并不能已知要解决这个问题 i j f i j f 我们可以采取一个迭代的方法来求解这个线性方程组 而不是基于 LU 分解直接求解 在 3 2 5 节中我们介绍了 Gauss Seidel 松弛迭代 算法 为了方便我们先回顾一下这个方法 给出一个线性方程组 如下 19 我们将使用下面的迭代方法 初始点为 0 X 其中 k 是迭代计数器 而 是松弛变量 该迭代将会一直进行 到满足收敛条件 如 其中 是容限参数 现在我们将尝试使用 Crank Nicolson 方法对美式期权进行定价 我们要解决的是和 9 6 式大致相同的一个线性方程组 但是这里 的边界条件有一点不同 价值为 0 的期权不用再考虑是否有障碍 我们要用回推法求解的线性方程组如下 其中右边有 式子中额外增加的一项是为了将看跌期权独有的边界条件带入 模型方程 在使用松弛方法的时候应当考虑矩阵的三对角性质 1 M 并且这能够适应美式期权提前执行的特征 设为内在价值 i g i 1 2 M 1 同时 对于每个时间层 j 我们进行如下SiS 迭代 20 当计算从一个时间层进行到下一个时间层的时候 将上一个时 间层的计算结果作为下一个时间层迭代计算的初始向量可能是更加 合理的 完成上述步骤需要的代码在图 9 5 中将列出 使用该代码 有点麻烦 原因在于 MATLAB 的向量幂运算是从 1 开始的 而现实却 要求从 0 开始 这样我们便不能建立一个包含所有的向量矩阵 i j f 同时系数矩阵也无法储存 上述迭代最好是直接通过 和 1 M 计算得出 该代码的实际效果可能要比 MATLAB 金融工具箱中附带的使用二 叉树方法对美式期权定价的 binprice 函数 见 7 1 节 要好 tic pr opt binprice 50 50 0 1 5 12 1 1200 0 4 0 toc Elapsed time is 0 seconds opt 1 l ans 4 2830 tic AmPutCK 50 50 0 1 5 12 0 4 100 1 1 600 1 5 0 001 toc ans 21 4 2815 Elapsed time is 0 seconds tic AmPutCK 50 50 0 1 5 12 0 4 100 1 1 600 1 8 0 001 toc ans 4 2794 Elapsed time is 0 seconds tic AmPutCK 50 50 0 1 5 12 0 4 100 1 1 600 1 2 0 001 toc ans 4 2800 Elapsed time is 0 seconds tic AmPutCK 50 50 0 1 5 12 0 4 100 1 1 1200 1 2 0 001 toc ans 4 2828 Elapsed time is 0 seconds tic AmPutCK 50 50 0 1 5 12 0 4 100 1 1 100 1 2 0 001 toc ans 4 2778 Elapsed time is 0 seconds function price AmPutCK S0 K r T sigma Smax dS dt omega tol M round Smax dS dS Smax M 建立网格 23 N round T dt dt T N oldval zeros M 1 1 Gauss Seidel更新向量 newval zeros M 1 1 vetS linspace 0 Smax M 1 veti 0 M vetj 0 N 建立边界条件 payoff max K vetS 2 M 0 pastval payoff values for the last layer boundval K exp r dt N vetj 边界值 建立系数矩阵和式子右边矩阵 alpha 0 25 dt sigma 2 veti 2 r veti beta dt 0 5 sigma 2 veti 2 r gamma 0 25 dt sigma 2 veti 2 r veti M2 diag alpha 3 M 1 diag 1 beta 2 M diag gamma 2 M 1 1 使用SOR方法求解线性方程组 aux zeros M 1 1 for j N 1 1 aux 1 alpha 2 boundval 1 j boundval 1 j 1 建立右端矩阵并进行初始化 rhs M2 pastval aux oldval pastval 24 error realmax while tol error newval 1 max payoff 1 oldval 1 omega 1 beta 2 rhs 1 1 beta 2 oldval 1 gamma 2 oldval 2 for k 2 M 2 newval k max payoff k oldval k omega 1 beta k 1 rhs k alpha k 1 newval k 1 1 beta k 1 oldval k gamma k 1 oldval k 1 end newval M 1 max payoff M 1 oldval M 1 omega 1 beta M rhs M 1 alpha M newval M 2 1 beta M oldval M 1 error norm newval oldval oldval newval end pastval newval end 25 newval boundval 1 newval 0 加入缺少的值 返回价格 这个价格可能因为初始资产价格在网格外而由线性 插值生成 price interp1 vetS newval S0 图图9 59 5 Crank NicolsonCrank Nicolson方法求解美式看跌期权方法求解美式看跌期权MATLABMATLAB代码代码 通过以上例子我们可以看到 松弛变量 在迭代方法中对结果 的收敛性有十分显著的影响 就运算速度而言 上述有限差分方法 CK 似乎要比二叉树方法快 但是我们必须要注意一些问题 当 我们将其与隐式法比较时 两者是不相上下的 而且 CPU 对 MATLAB 代码的解码方式也可能影响到计算效率 不管怎样 一个足够大而 2 且精确的网格要比二叉树的点能够提供更加有效的敏感度估计 这 在求解 B S 方程中会用到 此外 在处理一些奇异期权的定价问题 的时候 上述提到的有限差分方法可能是个更好的选择 扩展阅读扩展阅读 许多金融工程应用偏微分方程的例子可以