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

来源:https://uqer.io/community/share/555dc9e8f9f06c6c7404f96e

欢迎来到 Black - Scholes — Merton 的世界!本篇中我们将使用在第11天学习到的知识应用到这个金融学的具体方程之上!

  1. import numpy as np
  2. import math
  3. import seaborn as sns
  4. from matplotlib import pylab
  5. font.set_size(15)

1. 问题的提出

BSM模型可以设置为如下的偏微分方差初值问题:

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

做变量替换τ=T−t,并且设置上下边界条件:

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

2. 算法

按照之前介绍的隐式差分格式的方法,用离散差分格式代替连续微分:

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

其中

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

以上即为差分方程组。

这里有些细节需要处理,就是左右边界条件,我们这里使用Dirichlet边界条件,则:

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

3.实现

  1. import scipy as sp
  2. from scipy.linalg import solve_banded

描述BSM方程结构的类:BSMModel

  1. class BSMModel:
  2. def __init__(self, s0, r, sigma):
  3. self.s0 = s0
  4. self.x0 = math.log(s0)
  5. self.r = r
  6. self.sigma = sigma
  7. def log_expectation(self, T):
  8. return self.x0 + (self.r - 0.5 * self.sigma * self.sigma) * T
  9. def expectation(self, T):
  10. return math.exp(self.log_expectation(T))
  11. def x_max(self, T):
  12. return self.log_expectation(T) + 4.0 * self.sigma * math.sqrt(T)
  13. def x_min(self, T):
  14. return self.log_expectation(T) - 4.0 * self.sigma * math.sqrt(T)

描述我们这里设计到的产品的类:CallOption

  1. class CallOption:
  2. def __init__(self, strike):
  3. self.k = strike
  4. def ic(self, spot):
  5. return max(spot - self.k, 0.0)
  6. def bcl(self, spot, tau, model):
  7. return 0.0
  8. def bcr(self, spot, tau, model):
  9. return spot - math.exp(-model.r*tau) * self.k

完整的隐式格式:BSMScheme

  1. class BSMScheme:
  2. def __init__(self, model, payoff, T, M, N):
  3. self.model = model
  4. self.T = T
  5. self.M = M
  6. self.N = N
  7. self.dt = self.T / self.M
  8. self.payoff = payoff
  9. self.x_min = model.x_min(self.T)
  10. self.x_max = model.x_max(self.T)
  11. self.dx = (self.x_max - self.x_min) / self.N
  12. self.C = np.zeros((self.N+1, self.M+1)) # 全部网格
  13. self.xArray = np.linspace(self.x_min, self.x_max, self.N+1)
  14. self.C[:,0] = map(self.payoff.ic, np.exp(self.xArray))
  15. sigma_square = self.model.sigma*self.model.sigma
  16. r = self.model.r
  17. self.l_j = -(0.5*sigma_square*self.dt/self.dx/self.dx - 0.5 * (r - 0.5 * sigma_square)*self.dt/self.dx)
  18. self.c_j = 1.0 + sigma_square*self.dt/self.dx/self.dx + r*self.dt
  19. self.u_j = -(0.5*sigma_square*self.dt/self.dx/self.dx + 0.5 * (r - 0.5 * sigma_square)*self.dt/self.dx)
  20. def roll_back(self):
  21. for k in range(0, self.M):
  22. udiag = np.ones(self.N-1) * self.u_j
  23. ldiag = np.ones(self.N-1) * self.l_j
  24. cdiag = np.ones(self.N-1) * self.c_j
  25. mat = np.zeros((3,self.N-1))
  26. mat[0,:] = udiag
  27. mat[1,:] = cdiag
  28. mat[2,:] = ldiag
  29. rhs = np.copy(self.C[1:self.N,k])
  30. # 应用左端边值条件
  31. v1 = self.payoff.bcl(math.exp(self.x_min), (k+1)*self.dt, self.model)
  32. rhs[0] -= self.l_j * v1
  33. # 应用右端边值条件
  34. v2 = self.payoff.bcr(math.exp(self.x_max), (k+1)*self.dt, self.model)
  35. rhs[-1] -= self.u_j * v2
  36. x = solve_banded((1,1), mat, rhs)
  37. self.C[1:self.N, k+1] = x
  38. self.C[0][k+1] = v1
  39. self.C[self.N][k+1] = v2
  40. def mesh_grids(self):
  41. tArray = np.linspace(0, self.T, self.M+1)
  42. tGrids, xGrids = np.meshgrid(tArray, self.xArray)
  43. return tGrids, xGrids

应用在一起:

  1. model = BSMModel(100.0, 0.05, 0.2)
  2. payoff = CallOption(105.0)
  3. scheme = BSMScheme(model, payoff, 5.0, 100, 300)
  1. scheme.roll_back()
  1. from matplotlib import pylab
  2. pylab.figure(figsize=(12,8))
  3. pylab.plot(np.exp(scheme.xArray)[50:170], scheme.C[50:170,-1])
  4. pylab.xlabel('$S$')
  5. pylab.ylabel('$C$')
  6. <matplotlib.text.Text at 0x76ea7d0>

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

4. 收敛性测试

首先使用BSM模型的解析解获得精确解:

  1. analyticPrice = BSMPrice(1, 105., 100., 0.05, 0.0, 0.2, 5.)
  2. analyticPrice
pricedeltagammavegarhotheta
126.7618440.7496940.0071171.10319241.037549-3.832439

我们固定时间方向网格数为3000,分别计算不同S网格数情形下的结果:

  1. xSteps = range(50,300,10)
  2. finiteResult = []
  3. for xStep in xSteps:
  4. model = BSMModel(100.0, 0.05, 0.2)
  5. payoff = CallOption(105.0)
  6. scheme = BSMScheme(model, payoff, 5.0, 3000, xStep)
  7. scheme.roll_back()
  8. interp = CubicNaturalSpline(np.exp(scheme.xArray), scheme.C[:,-1])
  9. price = interp(100.0)
  10. finiteResult.append(price)

我们可以画下收敛图:

  1. anyRes = [analyticPrice['price'][1]] * len(xSteps)
  2. pylab.figure(figsize = (16,8))
  3. pylab.plot(xSteps, finiteResult, '-.', marker = 'o', markersize = 10)
  4. pylab.plot(xSteps, anyRes, '--')
  5. pylab.legend([u'隐式差分格式', u'解析解(欧式)'], prop = font)
  6. pylab.xlabel(u'价格方向网格步数', fontproperties = font)
  7. pylab.title(u'Black - Scholes - Merton 有限差分法', fontproperties = font, fontsize = 20)
  8. <matplotlib.text.Text at 0x7857bd0>

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