量化分析师的Python日记【Q Quant兵器谱之偏微分方程3的具体金融学运用】

原创 2017年01月03日 11:05:05

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

import numpy as np
import math
import seaborn as sns
from matplotlib import pylab
from CAL.PyCAL import *
font.set_size(15)

1. 问题的提出


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

{C(S,t)t+rSC(S,t)S+12σ2S22C(S,t)S2rC(S,t)=0,0tT\[5pt]C(S,T)=max(SK,0)

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

C(x,τ)τ\[5pt]C(x,0)\[5pt]C(xmax,τ)\[5pt]C(xmin,τ)=(r12σ2)C(x,τ)x+12σ22C(x,τ)x2rC(x,τ)=0,0τT=max(exK,0),=exmaxKeτ,=0

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 sp
from 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 = 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
        
16
        sigma_square = self.model.sigma*self.model.sigma
17
        r = self.model.r
18
        
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.dt
21
        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_j
27
            ldiag =  np.ones(self.N-1) * self.l_j
28
            cdiag =  np.ones(self.N-1) * self.c_j
29
            
30
            mat = np.zeros((3,self.N-1))
31
            mat[0,:] = udiag
32
            mat[1,:] = cdiag
33
            mat[2,:] = ldiag
34
            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 * v1
39
            
40
            # 应用右端边值条件
41
            v2 = self.payoff.bcr(math.exp(self.x_max), (k+1)*self.dt, self.model)
42
            rhs[-1] -= self.u_j * v2
43
            
44
            x = solve_banded((1,1), mat, rhs)
45
            self.C[1:self.N, k+1] = x
46
            self.C[0][k+1] = v1
47
            self.C[self.N][k+1] = v2
48
            
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 pylab
pylab.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,分别计算不同S网格数情形下的结果:

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兵器谱之二叉树】

通过之前几天的学习,Q Quant们应该已经熟悉了Python的基本语法,也了解了Python中常用数值库的算法。到这里为止,小Q们也许早就对之前简单的例子不满意,希望能在Python里面玩票大的!O...

量化分析师的Python日记【第6天:数据处理的瑞士军刀pandas下篇

####第二篇:快速进阶 在上一篇中我们介绍了如何创建并访问pandas的Series和DataFrame类型的数据,本篇将介绍如何对pandas数据进行操作,掌握这些操...

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

从今天开始我们将进入一个系列 —— 偏微分方程。作为这一系列的开篇,我们以热传导方差为引子,引出: 如何提一个偏微分方程的初边值问题;利用差分格式将偏微分方程离散化;显示差分格式;显示差分格式的条件...

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

这是量化分析师的偏微分方程系列的第二篇,在这一篇中我们将解决上一篇显式格式留下的稳定性问题。本篇将引入隐式差分算法,读者可以学到: 隐式差分格式描述三对角矩阵求解如何使用scipy加速算法实现 ...

量化分析师的Python日记【Q Quant兵器谱之函数插值】

在本篇中,我们将介绍Q宽客常用工具之一:函数插值。接着将函数插值应用于一个实际的金融建模场景中:波动率曲面构造。...

量化分析师的Python日记【Q Quant 之初出江湖】

本篇中,作为Quant中的Q宗(P Quant 和 Q Quant 到底哪个是未来?),我们将尝试把之前的介绍的工具串联起来,小试牛刀。 您将可以体验到: 如何使用python内置的数学函...

量化分析师的Python日记【第4天:一大波金融Library来袭之scipy篇】

上一篇介绍了numpy,本篇中着重介绍一下另一个量化金融中常用的库 scipy   一、SciPy概述 前篇已经大致介绍了NumPy,接下来让我们看看SciPy能...

量化分析师的Python日记【第4天:一大波金融Library来袭之scipy篇】

####一、SciPy概述 前篇已经大致介绍了NumPy,接下来让我们看看SciPy能做些什么。NumPy替我们搞定了向量和矩阵的相关操作,基本上算是一个高级的科学计算器。SciPy基于NumP...

量化分析师的Python日记【第3天:一大波金融Library来袭之numpy篇】

接下来要给大家介绍的系列中包含了Python在量化金融中运用最广泛的几个Library: numpy scipy pandas matplotlib 会给初学者一一介绍 ...
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:量化分析师的Python日记【Q Quant兵器谱之偏微分方程3的具体金融学运用】
举报原因:
原因补充:

(最多只允许输入30个字)