量化分析师的Python日记【Q Quant兵器谱之偏微分方程3的具体金融学运用】
来源:互联网 发布:虚拟网络映射 编辑:程序博客网 时间:2024/05/18 18:47
欢迎来到Black - Scholes — Merton的世界!本篇中我们将把第11天学习到的知识应用到这个金融学的具体方程之上!
import numpy as npimport mathimport seaborn as snsfrom matplotlib import pylabfrom CAL.PyCAL import *font.set_size(15)
1. 问题的提出
BSM模型可以设置为如下的偏微分方差初值问题:
做变量替换
2. 算法
按照之前介绍的隐式差分格式的方法,用离散差分格式代替连续微分:
\begin{align}&\frac{C_{j,k+1} - C_{j,k}}{\Delta \tau} = (r - \frac{1}{2}\sigma^2)\frac{C_{j+1,k+1} - C_{j-1,k+1}}{2\Delta x} + \frac{1}{2}\sigma^2 \frac{C_{j+1,k+1} - 2C_{j,k+1} + C_{j-1,k+1}}{\Delta x^2} - rC_{j,k+1}, \\\[5pt]\Rightarrow& \quad C_{j,k+1} - C_{j,k} - (r - \frac{1}{2}\sigma^2)\frac{\Delta \tau}{2\Delta x}(C_{j+1,k+1} - C_{j-1,k+1}) - \frac{1}{2}\sigma^2 \frac{\Delta \tau}{\Delta x^2}(C_{j+1,k+1} - 2C_{j,k+1} + C_{j-1,k+1}) + r\Delta \tau C_{j,k+1} = 0, \\\[5pt]\Rightarrow& \quad - (\frac{1}{2}(r - \frac{1}{2}\sigma^2)\frac{\Delta \tau}{\Delta x} + \frac{1}{2}\sigma^2\frac{\Delta \tau}{\Delta x^2})C_{j+1,k+1} + (1 + \sigma^2\frac{\Delta \tau}{\Delta x^2} + r\Delta \tau)C_{j,k+1} - (\frac{1}{2}\sigma^2\frac{\Delta \tau}{\Delta x^2} - \frac{1}{2}(r - \frac{1}{2}\sigma^2 )\frac{\Delta \tau}{\Delta x})C_{j-1,k+1} = C_{j,k}, \\\[5pt]\Rightarrow& \quad l_j C_{j-1,k+1} + c_j C_{j,k+1} + u_j C_{j+1,k+1} = C_{j,k}\end{align}
其中\begin{align}l_j &= - (\frac{1}{2}\sigma^2\frac{\Delta \tau}{\Delta x^2} - \frac{1}{2}(r - \frac{1}{2}\sigma^2 )\frac{\Delta \tau}{\Delta x}), \\\[5pt]c_j &= 1 + \sigma^2\frac{\Delta \tau}{\Delta x^2} + r\Delta \tau, \\\[5pt]u_j &= - (\frac{1}{2}(r - \frac{1}{2}\sigma^2)\frac{\Delta \tau}{\Delta x} + \frac{1}{2}\sigma^2\frac{\Delta \tau}{\Delta x^2})\end{align}
以上即为差分方程组。
这里有些细节需要处理,就是左右边界条件,我们这里使用Dirichlet边界条件,则:
\begin{align*}C_{0,k} = C(x_{min},\tau), \\\[5pt]C_{N,k} = C(x_{max},\tau)
\end{align*}
3.实现
import scipy as spfrom scipy.linalg import solve_banded
描述BSM方程结构的类:BSModel
class BSMModel: def __init__(self, s0, r, sigma): self.s0 = s0 self.x0 = math.log(s0) self.r = r self.sigma = sigma def log_expectation(self, T): return self.x0 + (self.r - 0.5 * self.sigma * self.sigma) * T def expectation(self, T): return math.exp(self.log_expectation(T)) def x_max(self, T): return self.log_expectation(T) + 4.0 * self.sigma * math.sqrt(T) def x_min(self, T): return self.log_expectation(T) - 4.0 * self.sigma * math.sqrt(T)
描述我们这里设计到的产品的类:
CallOption
class CallOption: def __init__(self, strike): self.k = strike def ic(self, spot): return max(spot - self.k, 0.0) def bcl(self, spot, tau, model): return 0.0 def bcr(self, spot, tau, model): return spot - math.exp(-model.r*tau) * self.k完整的隐式格式:
BSMScheme
class BSMScheme:2 def __init__(self, model, payoff, T, M, N):3 self.model = model4 self.T = T5 self.M = M6 self.N = N7 self.dt = self.T / self.M8 self.payoff = payoff9 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.N12 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 16 sigma_square = self.model.sigma*self.model.sigma17 r = self.model.r18 19 self.l_j = -(0.5*sigma_square*self.dt/self.dx/self.dx - 0.5 * (r - 0.5 * sigma_square)*self.dt/self.dx)20 self.c_j = 1.0 + sigma_square*self.dt/self.dx/self.dx + r*self.dt21 self.u_j = -(0.5*sigma_square*self.dt/self.dx/self.dx + 0.5 * (r - 0.5 * sigma_square)*self.dt/self.dx)22 23 def roll_back(self):24 25 for k in range(0, self.M):26 udiag = np.ones(self.N-1) * self.u_j27 ldiag = np.ones(self.N-1) * self.l_j28 cdiag = np.ones(self.N-1) * self.c_j29 30 mat = np.zeros((3,self.N-1))31 mat[0,:] = udiag32 mat[1,:] = cdiag33 mat[2,:] = ldiag34 rhs = np.copy(self.C[1:self.N,k])35 36 # 应用左端边值条件37 v1 = self.payoff.bcl(math.exp(self.x_min), (k+1)*self.dt, self.model)38 rhs[0] -= self.l_j * v139 40 # 应用右端边值条件41 v2 = self.payoff.bcr(math.exp(self.x_max), (k+1)*self.dt, self.model)42 rhs[-1] -= self.u_j * v243 44 x = solve_banded((1,1), mat, rhs)45 self.C[1:self.N, k+1] = x46 self.C[0][k+1] = v147 self.C[self.N][k+1] = v248 49 def mesh_grids(self):50 tArray = np.linspace(0, self.T, self.M+1)51 tGrids, xGrids = np.meshgrid(tArray, self.xArray)52 return tGrids, xGrids
应用在一起:
model = BSMModel(100.0, 0.05, 0.2)payoff = CallOption(105.0)scheme = BSMScheme(model, payoff, 5.0, 100, 300)
scheme.roll_back()
from matplotlib import pylabpylab.figure(figsize=(12,8))pylab.plot(np.exp(scheme.xArray)[50:170], scheme.C[50:170,-1])pylab.xlabel('$S$')pylab.ylabel('$C$')
4. 收敛性测试
首先使用BSM模型的解析解获得精确解:analyticPrice = BSMPrice(1, 105., 100., 0.05, 0.0, 0.2, 5.)analyticPrice
我们固定时间方向网格数为3000,分别计算不同
xSteps = range(50,300,10)finiteResult = []for xStep in xSteps: model = BSMModel(100.0, 0.05, 0.2) payoff = CallOption(105.0) scheme = BSMScheme(model, payoff, 5.0, 3000, xStep) scheme.roll_back() interp = CubicNaturalSpline(np.exp(scheme.xArray), scheme.C[:,-1]) price = interp(100.0) finiteResult.append(price)
我们可以画下收敛图:
anyRes = [analyticPrice['price'][1]] * len(xSteps)pylab.figure(figsize = (16,8))pylab.plot(xSteps, finiteResult, '-.', marker = 'o', markersize = 10)pylab.plot(xSteps, anyRes, '--')pylab.legend([u'隐式差分格式', u'解析解(欧式)'], prop = font)pylab.xlabel(u'价格方向网格步数', fontproperties = font)pylab.title(u'Black - Scholes - Merton 有限差分法', fontproperties = font, fontsize = 20)
- 量化分析师的Python日记【Q Quant兵器谱之偏微分方程3的具体金融学运用】
- 量化分析师的Python日记【Q Quant兵器谱 -之偏微分方程1】
- 量化分析师的Python日记【Q Quant兵器谱之偏微分方程2】
- 量化分析师的Python日记【Q Quant兵器谱之函数插值】
- 量化分析师的Python日记【Q Quant兵器谱之二叉树】
- Python--量化分析师的Python日记【第7天:Q Quant 之初出江湖】
- 量化分析师的Python日记【Q Quant 之初出江湖】
- 量化分析师的Python日记
- 量化分析师的Python日记
- 量化分析师的Python日记【第3天:一大波金融Library来袭之numpy篇】
- 量化分析师的Python日记【第3天:一大波金融Library来袭之numpy篇】
- 量化投资领域--关于Quant的FAQ
- 量化分析师的Python日记【第4天:一大波金融Library来袭之scipy篇】
- 量化分析师的Python日记【第4天:一大波金融Library来袭之scipy篇】
- 量化分析师的Python日记【第1天:谁来给我讲讲Python?】
- 量化分析师的Python日记【第1天:谁来给我讲讲Python?】
- 量化分析师的Python日记【第2天:再接着介绍一下Python呗】
- 量化分析师的Python日记【第1天:谁来给我讲讲Python?】
- R语言相关分析
- Java出现No enclosing instance of type E is accessible. Must qualify the allocation with an enclosing
- Android N 来电流程(MT)
- 放苹果
- SystemUI下的快速设置面板显示异常
- 量化分析师的Python日记【Q Quant兵器谱之偏微分方程3的具体金融学运用】
- git
- HTM-16.2代码(1)——编码端一些函数的说明
- 确定主题域
- 汽车权威网站
- Webapp中1px边框在retina屏中变粗的问题
- EXCEL疑难记录
- Java反射机制——学习总结
- Drupal 模块的 Hooks(钩子)。