量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】

来源:https://uqer.io/community/share/5534ad3ff9f06c8f33904689

这是量化分析师的偏微分方程系列的第二篇,在这一篇中我们将解决上一篇显式格式留下的稳定性问题。本篇将引入隐式差分算法,读者可以学到:

  1. 隐式差分格式描述
  2. 三对角矩阵求解
  3. 如何使用scipy加速算法实现

在完成两天的基础学习之后,在下一天中,我们将把已经学到的知识运用到金融定价领域最重要的方程之一:Black - Shcoles - Merton 偏微分方差

  1. from matplotlib import pylab
  2. import seaborn as sns
  3. import numpy as np
  4. np.set_printoptions(precision = 4)
  5. font.set_size(20)
  6. def initialCondition(x):
  7. return 4.0*(1.0 - x) * x

1. 隐式差分格式

像上一天一样,我们从差分格式的数学表述开始。隐式格式与显式格式的区别,在于我们时间方向选择的基准点。显式格式使用k,而隐式格式选择k+1:

量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】 - 图1

剩下的推到过程我完全一样,我们看到无论隐式格式还是显式格式,它们的截断误差是一样的:

量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】 - 图2

用离散值Uj,k替换uj,k,我们得到差分方程:

量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】 - 图3

最后,到这里我们得到一个迭代方程组:

量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】 - 图4

其中量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】 - 图5

  1. N = 500 # x方向网格数
  2. M = 500 # t方向网格数
  3. T = 1.0
  4. X = 1.0
  5. xArray = np.linspace(0,X,N+1)
  6. yArray = map(initialCondition, xArray)
  7. starValues = yArray
  8. U = np.zeros((N+1,M+1))
  9. U[:,0] = starValues
  1. dx = X / N
  2. dt = T / M
  3. kappa = 1.0
  4. rho = kappa * dt / dx / dx

1.1 矩阵求解(TridiagonalSystem)

虽然看上去形式只是变了一点,但是求解的问题有很大的变化。在每个时间点上,我们需要求解如下的一个线性方程组:

量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】 - 图6

这里A为:

量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】 - 图7

幸运的是,这个是个三对角矩阵,可以很简单的利用Gauss消去法求解。我们这里不会详细讨论算法的描述,细节都可以在下面的python类TridiagonalSystem中了解到:

  1. class TridiagonalSystem:
  2. def __init__(self, udiag, cdiag, ldiag):
  3. '''
  4. 三对角矩阵:
  5. udiag -- 上对角线
  6. cdiag -- 对角线
  7. ldiag -- 下对角线
  8. '''
  9. assert len(udiag) == len(cdiag)
  10. assert len(cdiag) == len(ldiag)
  11. self.udiag = udiag
  12. self.cdiag = cdiag
  13. self.ldiag = ldiag
  14. self.length = len(self.cdiag)
  15. def solve(self, rhs):
  16. '''
  17. 求解以下方程组
  18. A \ dot x = rhs
  19. '''
  20. assert len(rhs) == len(self.cdiag)
  21. udiag = self.udiag.copy()
  22. cdiag = self.cdiag.copy()
  23. ldiag = self.ldiag.copy()
  24. b = rhs.copy()
  25. # 消去下对角元
  26. for i in range(1, self.length):
  27. cdiag[i] -= udiag[i-1] * ldiag[i] / cdiag[i-1]
  28. b[i] -= b[i-1] * ldiag[i] / cdiag[i-1]
  29. # 从最后一个方程开始求解
  30. x = np.zeros(self.length)
  31. x[self.length-1] = b[self.length - 1] / cdiag[self.length - 1]
  32. for i in range(self.length - 2, -1, -1):
  33. x[i] = (b[i] - udiag[i]*x[i+1]) / cdiag[i]
  34. return x
  35. def multiply(self, x):
  36. '''
  37. 矩阵乘法:
  38. rhs = A \dot x
  39. '''
  40. assert len(x) == len(self.cdiag)
  41. rhs = np.zeros(self.length)
  42. rhs[0] = x[0] * self.cdiag[0] + x[1] * self.udiag[0]
  43. for i in range(1, self.length - 1):
  44. rhs[i] = x[i-1] * self.ldiag[i] + x[i] * self.cdiag[i] + x[i+1] * self.udiag[i]
  45. rhs[self.length - 1] = x[self.length - 2] * self.ldiag[self.length - 1] + x[self.length - 1] * self.cdiag[self.length - 1]
  46. return rhs

1.2 隐式格式求解

  1. for k in range(0, M):
  2. udiag = - np.ones(N-1) * rho
  3. ldiag = - np.ones(N-1) * rho
  4. cdiag = np.ones(N-1) * (1.0 + 2. * rho)
  5. mat = TridiagonalSystem(udiag, cdiag, ldiag)
  6. rhs = U[1:N,k]
  7. x = mat.solve(rhs)
  8. U[1:N, k+1] = x
  9. U[0][k+1] = 0.
  10. U[N][k+1] = 0.
  1. from lib.utilities import plotLines
  2. plotLines([U[:,0], U[:, int(0.10/ dt)], U[:, int(0.20/ dt)], U[:, int(0.50/ dt)]], xArray, title = u'一维热传导方程', xlabel = '$x$',
  3. ylabel = r'$U(\dot, \tau)$', legend = [r'$\tau = 0.$', r'$\tau = 0.10$', r'$\tau = 0.20$', r'$\tau = 0.50$'])

量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】 - 图8

  1. from lib.utilities import plotSurface
  2. tArray = np.linspace(0, 0.2, int(0.2 / dt) + 1)
  3. tGrids, xGrids = np.meshgrid(tArray, xArray)
  4. plotSurface(xGrids, tGrids, U[:,:int(0.2 / dt) + 1], title = u"热传导方程 $u_\\tau = u_{xx}$,隐式格式($\\rho = 50$)", xlabel = "$x$", ylabel = r"$\tau$", zlabel = r"$U$")

量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】 - 图9

2. 继续组装

像我们在显示格式那一节介绍的同样做法,我们把之前的代码整合起来,归集与一个完整的类ImplicitEulerScheme中:

  1. from lib.utilities import HeatEquation

上面的代码(使用library功能,关于该功能的具体介绍请见帮助 — Library是干什么的)导入我们在上一期中已经定义过的类HeatEquation,避免代码重复。

  1. class ImplicitEulerScheme:
  2. def __init__(self, M, N, equation):
  3. self.eq = equation
  4. self.dt = self.eq.T / M
  5. self.dx = self.eq.X / N
  6. self.U = np.zeros((N+1, M+1))
  7. self.xArray = np.linspace(0,self.eq.X,N+1)
  8. self.U[:,0] = map(self.eq.ic, self.xArray)
  9. self.rho = self.eq.kappa * self.dt / self.dx / self.dx
  10. self.M = M
  11. self.N = N
  12. def roll_back(self):
  13. for k in range(0, self.M):
  14. udiag = - np.ones(self.N-1) * self.rho
  15. ldiag = - np.ones(self.N-1) * self.rho
  16. cdiag = np.ones(self.N-1) * (1.0 + 2. * self.rho)
  17. mat = TridiagonalSystem(udiag, cdiag, ldiag)
  18. rhs = self.U[1:self.N,k]
  19. x = mat.solve(rhs)
  20. self.U[1:self.N, k+1] = x
  21. self.U[0][k+1] = self.eq.bcl(self.xArray[0])
  22. self.U[self.N][k+1] = self.eq.bcr(self.xArray[-1])
  23. def mesh_grids(self):
  24. tArray = np.linspace(0, self.eq.T, M+1)
  25. tGrids, xGrids = np.meshgrid(tArray, self.xArray)
  26. return tGrids, xGrids

然后我们可以使用下面的三行简单调用完成功能:

  1. ht = HeatEquation(1.,X, T)
  2. scheme = ImplicitEulerScheme(M,N, ht)
  3. scheme.roll_back()
  4. scheme.U
  5. array([[ 0.0000e+00, 0.0000e+00, 0.0000e+00, ..., 0.0000e+00,
  6. 0.0000e+00, 0.0000e+00],
  7. [ 7.9840e-03, 7.2843e-03, 6.9266e-03, ..., 3.8398e-07,
  8. 3.7655e-07, 3.6926e-07],
  9. [ 1.5936e-02, 1.4567e-02, 1.3852e-02, ..., 7.6795e-07,
  10. 7.5308e-07, 7.3851e-07],
  11. ...,
  12. [ 1.5936e-02, 1.4567e-02, 1.3852e-02, ..., 7.6795e-07,
  13. 7.5308e-07, 7.3851e-07],
  14. [ 7.9840e-03, 7.2843e-03, 6.9266e-03, ..., 3.8398e-07,
  15. 3.7655e-07, 3.6926e-07],
  16. [ 0.0000e+00, 0.0000e+00, 0.0000e+00, ..., 0.0000e+00,
  17. 0.0000e+00, 0.0000e+00]])

3. 使用 scipy加速

软件工程行业里有句老话,叫做:“不要重复发明轮子!”。实际上,之前的代码里面,我们就造了自己的轮子:TridiagonalSystem。三对角矩阵作为最最常见的稀疏矩阵,关于它的线性方程组求解算法实际上早已为业界熟知,也已经有很多库内置了工业级别强度实现。这里我们取scipy作为例子,来展示使用外源库实现的好处:

  • 更加稳健的算法: 知名库算法由于使用者广泛,有更大的概率发现一些极端情形下的bug。库作者可以根据用户反馈,及时调整算法;
  • 更高的性能: 由于库的使用更为广泛,库作者有更大的动力去使用各种技术去提高算法的性能:例如使用更高效的语言实现,例如C。scipy中的情形就是一例。
  • 持续的维护: 库的受众范围广,社区的力量会推动库作者持续维护。 下面的代码展示,如何使用scipy中的solve_banded算法求解三对角矩阵:
  1. import scipy as sp
  2. from scipy.linalg import solve_banded
  3. A = np.zeros((3, 5))
  4. A[0, :] = np.ones(5) * 1. # 上对角线
  5. A[1, :] = np.ones(5) * 3. # 对角线
  6. A[2, :] = np.ones(5) * (-1.) # 下对角线
  7. b = [1.,2.,3.,4.,5.]
  8. x = solve_banded ((1,1), A,b)
  9. print 'x = A^-1b = ',x
  10. x = A^-1b = [ 0.1833 0.45 0.8333 0.95 1.9833]

我们使用上面的算法替代我们之前的TridiagonalSystem

  1. import scipy as sp
  2. from scipy.linalg import solve_banded
  3. for k in range(0, M):
  4. udiag = - np.ones(N-1) * rho
  5. ldiag = - np.ones(N-1) * rho
  6. cdiag = np.ones(N-1) * (1.0 + 2. * rho)
  7. mat = np.zeros((3,N-1))
  8. mat[0,:] = udiag
  9. mat[1,:] = cdiag
  10. mat[2,:] = ldiag
  11. rhs = U[1:N,k]
  12. x = solve_banded ((1,1), mat,rhs)
  13. U[1:N, k+1] = x
  14. U[0][k+1] = 0.
  15. U[N][k+1] = 0.
  1. plotLines([U[:,0], U[:, int(0.10/ dt)], U[:, int(0.20/ dt)], U[:, int(0.50/ dt)]], xArray, title = u'一维热传导方程,使用scipy', xlabel = '$x$',
  2. ylabel = r'$U(\dot, \tau)$', legend = [r'$\tau = 0.$', r'$\tau = 0.10$', r'$\tau = 0.20$', r'$\tau = 0.50$'])

量化分析师的Python日记【第11天 Q Quant兵器谱之偏微分方程2】 - 图10

同样的我们定义一个新类ImplicitEulerSchemeWithScipy使用scipy的算法:

  1. class ImplicitEulerSchemeWithScipy:
  2. def __init__(self, M, N, equation):
  3. self.eq = equation
  4. self.dt = self.eq.T / M
  5. self.dx = self.eq.X / N
  6. self.U = np.zeros((N+1, M+1))
  7. self.xArray = np.linspace(0,self.eq.X,N+1)
  8. self.U[:,0] = map(self.eq.ic, self.xArray)
  9. self.rho = self.eq.kappa * self.dt / self.dx / self.dx
  10. self.M = M
  11. self.N = N
  12. def roll_back(self):
  13. for k in range(0, self.M):
  14. udiag = - np.ones(self.N-1) * self.rho
  15. ldiag = - np.ones(self.N-1) * self.rho
  16. cdiag = np.ones(self.N-1) * (1.0 + 2. * self.rho)
  17. mat = np.zeros((3,self.N-1))
  18. mat[0,:] = udiag
  19. mat[1,:] = cdiag
  20. mat[2,:] = ldiag
  21. rhs = self.U[1:self.N,k]
  22. x = solve_banded((1,1), mat, rhs)
  23. self.U[1:self.N, k+1] = x
  24. self.U[0][k+1] = self.eq.bcl(self.xArray[0])
  25. self.U[self.N][k+1] = self.eq.bcr(self.xArray[-1])
  26. def mesh_grids(self):
  27. tArray = np.linspace(0, self.eq.T, M+1)
  28. tGrids, xGrids = np.meshgrid(tArray, self.xArray)
  29. return tGrids, xGrids

下面的代码,比较了两种做法的性能。可以看到仅仅简单的替代三对角矩阵算法,我们就获得了接近8倍的性能提升

  1. import time
  2. startTime = time.time()
  3. loop_round = 10
  4. # 不使用scipy
  5. for k in range(loop_round):
  6. ht = HeatEquation(1.,X, T)
  7. scheme = ImplicitEulerScheme(M,N, ht)
  8. scheme.roll_back()
  9. endTime = time.time()
  10. print '{0:<40}{1:.4f}'.format('执行时间(s) -- 不使用scipy.linalg: ', endTime - startTime)
  11. # 使用scipy
  12. startTime = time.time()
  13. for k in range(loop_round):
  14. ht = HeatEquation(1.,X, T)
  15. scheme = ImplicitEulerSchemeWithScipy(M,N, ht)
  16. scheme.roll_back()
  17. endTime = time.time()
  18. print '{0:<40}{1:.4f}'.format('执行时间(s) -- 使用scipy.linalg: ', endTime - startTime)
  19. 执行时间(s) -- 不使用scipy.linalg: 12.1589
  20. 执行时间(s) -- 使用scipy.linalg: 1.6224

4. 尾声

到这里为止,我们已经结束了偏微分方差差分格式的基础学习。这是一个很大的学科,这两天也只能做到“管中窥豹”。但是有了以上的基础知识,读者已经有了足够的积累,可以处理一些金融工程中会实际遇到的方程。在下一天中,我们将把这两天学习到的知识运用到金融工程史上最重要的方程:Black - Scholes - Merton 偏微分方程。