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

来源:https://uqer.io/community/share/5514fc98f9f06c8f33904449

通过前几日的学习,我们已经熟悉了Python中一些常用数值计算库的用法。本篇中,作为Quant中的Q宗(P Quant 和 Q Quant 到底哪个是未来?),我们将尝试把之前的介绍的工具串联起来,小试牛刀。

您将可以体验到:

  1. 如何使用python内置的数学函数计算期权的价格;
  2. 利用 numpy 加速数值计算;
  3. 利用 scipy 进行仿真模拟;
  4. 使用 scipy 求解器计算隐含波动率;

穿插着,我们也会使用matplotlib绘制精美的图标。

1. 关心的问题

我们想知道下面的一只期权的价格:

  • 当前价 spot : 2.45
  • 行权价 strike: 2.50
  • 到期期限 maturity : 0.25
  • 无风险利率 r : 0.05
  • 波动率 vol : 0.25 关于这样的简单欧式期权的定价,有经典的Black - Scholes [1] 公式:

量化分析师的Python日记【第7天:Q Quant 之初出江湖】 - 图1

其中S为标的价格,K为执行价格,r为无风险利率,τ=T−t为剩余到期时间。 N(x)为标准正态分布的累积概率密度函数。Call(S,K,r,τ,σ)为看涨期权的价格。

  1. # 参数
  2. spot = 2.45
  3. strike = 2.50
  4. maturity = 0.25
  5. r = 0.05
  6. vol = 0.25

观察上面的公式,需要使用一些数学函数,我们把它分为两部分:

  • log,sqrt,exp,这三个函数我们可以从标准库math中找到
  • 标准正态分布的累计概率密度函数,我们使用scipy库中的stats.norm.cdf函数
  1. # 基于Black - Scholes 公式的期权定价公式
  2. from math import log, sqrt, exp
  3. from scipy.stats import norm
  4. def call_option_pricer(spot, strike, maturity, r, vol):
  5. d1 = (log(spot/strike) + (r + 0.5 * vol *vol) * maturity) / vol / sqrt(maturity)
  6. d2 = d1 - vol * sqrt(maturity)
  7. price = spot * norm.cdf(d1) - strike * exp(-r*maturity) * norm.cdf(d2)
  8. return price

我们可以使用这个函数计算我们关注期权的结果:

  1. print '期权价格 : %.4f' % call_option_pricer(spot, strike, maturity, r, vol)
  2. 期权价格 : 0.1133

2. 使用numpy加速批量计算

大部分的时候,我们不止关心一个期权的价格,而是关心一个组合(成千上万)的期权。我们想知道, 随着期权组合数量的增长,我们计算时间的增长会有多块?

2.1 使用循环的方式

  1. import time
  2. import numpy as np
  3. portfolioSize = range(1, 10000, 500)
  4. timeSpent = []
  5. for size in portfolioSize:
  6. now = time.time()
  7. strikes = np.linspace(2.0,3.0,size)
  8. for i in range(size):
  9. res = call_option_pricer(spot, strikes[i], maturity, r, vol)
  10. timeSpent.append(time.time() - now)

从下图中可以看出,计算时间的增长可以说是随着组合规模的增长线性上升。

  1. from matplotlib import pylab
  2. import seaborn as sns
  3. font.set_size(15)
  4. sns.set(style="ticks")
  5. pylab.figure(figsize = (12,8))
  6. pylab.bar(portfolioSize, timeSpent, color = 'r', width =300)
  7. pylab.grid(True)
  8. pylab.title(u'期权计算时间耗时(单位:秒)', fontproperties = font, fontsize = 18)
  9. pylab.ylabel(u'时间(s)', fontproperties = font, fontsize = 15)
  10. pylab.xlabel(u'组合数量', fontproperties = font, fontsize = 15)
  11. <matplotlib.text.Text at 0xdbad950>

量化分析师的Python日记【第7天:Q Quant 之初出江湖】 - 图2

2.2 使用numpy向量计算

numpy的内置数学函数可以天然的运用于向量:

  1. sample = np.linspace(1.0,100.0,5)
  2. np.exp(sample)
  3. array([ 2.71828183e+00, 1.52434373e+11, 8.54813429e+21,
  4. 4.79357761e+32, 2.68811714e+43])

利用 numpy 的数学函数,我们可以重写原先的计算公式 call_option_pricer,使得它接受向量参数。

  1. # 使用numpy的向量函数重写Black - Scholes公式
  2. def call_option_pricer_nunmpy(spot, strike, maturity, r, vol):
  3. d1 = (np.log(spot/strike) + (r + 0.5 * vol *vol) * maturity) / vol / np.sqrt(maturity)
  4. d2 = d1 - vol * np.sqrt(maturity)
  5. price = spot * norm.cdf(d1) - strike * np.exp(-r*maturity) * norm.cdf(d2)
  6. return price
  1. timeSpentNumpy = []
  2. for size in portfolioSize:
  3. now = time.time()
  4. strikes = np.linspace(2.0,3.0, size)
  5. res = call_option_pricer_nunmpy(spot, strikes, maturity, r, vol)
  6. timeSpentNumpy.append(time.time() - now)

再观察一下计算耗时,虽然时间仍然是随着规模的增长线性上升,但是增长的速度要慢许多:

  1. pylab.figure(figsize = (12,8))
  2. pylab.bar(portfolioSize, timeSpentNumpy, color = 'r', width = 300)
  3. pylab.grid(True)
  4. pylab.title(u'期权计算时间耗时(单位:秒)- numpy加速版', fontproperties = font, fontsize = 18)
  5. pylab.ylabel(u'时间(s)', fontproperties = font, fontsize = 15)
  6. pylab.xlabel(u'组合数量', fontproperties = font, fontsize = 15)
  7. <matplotlib.text.Text at 0xe0ba090>

量化分析师的Python日记【第7天:Q Quant 之初出江湖】 - 图3

让我们把两次计算时间进行比对,更清楚的了解 numpy 计算效率的提升!

  1. fig = pylab.figure(figsize = (12,8))
  2. ax = fig.gca()
  3. pylab.plot(portfolioSize, np.log10(timeSpent), portfolioSize, np.log(timeSpentNumpy))
  4. pylab.grid(True)
  5. from matplotlib.ticker import FuncFormatter
  6. def millions(x, pos):
  7. 'The two args are the value and tick position'
  8. return '$10^{%.0f}$' % (x)
  9. formatter = FuncFormatter(millions)
  10. ax.yaxis.set_major_formatter(formatter)
  11. pylab.title(u'期权计算时间耗时(单位:秒)', fontproperties = font, fontsize = 18)
  12. pylab.legend([u'循环计算', u'numpy向量加速'], prop = font, loc = 'upper center', ncol = 2)
  13. pylab.ylabel(u'时间(秒)', fontproperties = font, fontsize = 15)
  14. pylab.xlabel(u'组合数量', fontproperties = font, fontsize = 15)
  15. <matplotlib.text.Text at 0xe0b6390>

量化分析师的Python日记【第7天:Q Quant 之初出江湖】 - 图4

3. 使用scipy做仿真计算

期权价格的计算方法中有一类称为 蒙特卡洛 方法。这是利用随机抽样的方法,模拟标的股票价格随机游走,计算期权价格(未来的期望)。假设股票价格满足以下的随机游走:

量化分析师的Python日记【第7天:Q Quant 之初出江湖】 - 图5

仿真的方法可以模拟到期日的股票价格

量化分析师的Python日记【第7天:Q Quant 之初出江湖】 - 图6

这里的z是一个符合标准正态分布的随机数。这样我们可以计算最后的期权价格:

量化分析师的Python日记【第7天:Q Quant 之初出江湖】 - 图7

标准正态分布的随机数获取,可以方便的求助于 scipy 库:

  1. import scipy
  2. scipy.random.randn(10)
  3. array([ 0.36802702, 1.09560268, -1.0235275 , 0.15722882, 0.83718188,
  4. -0.27193135, -0.03485659, 1.02705248, 0.69479874, -0.35967107])
  1. pylab.figure(figsize = (12,8))
  2. randomSeries = scipy.random.randn(1000)
  3. pylab.plot(randomSeries)
  4. print u'均 值:%.4f' % randomSeries.mean()
  5. print u'标准差:%.4f' % randomSeries.std()
  6. 值:0.0336
  7. 标准差:0.9689

量化分析师的Python日记【第7天:Q Quant 之初出江湖】 - 图8

结合 scipy numpy 我们可以定义基于蒙特卡洛的期权定价算法。

  1. # 期权计算的蒙特卡洛方法
  2. def call_option_pricer_monte_carlo(spot, strike, maturity, r, vol, numOfPath = 5000):
  3. randomSeries = scipy.random.randn(numOfPath)
  4. s_t = spot * np.exp((r - 0.5 * vol * vol) * maturity + randomSeries * vol * sqrt(maturity))
  5. sumValue = np.maximum(s_t - strike, 0.0).sum()
  6. price = exp(-r*maturity) * sumValue / numOfPath
  7. return price
  1. print '期权价格(蒙特卡洛) : %.4f' % call_option_pricer_monte_carlo(spot, strike, maturity, r, vol)
  2. 期权价格(蒙特卡洛) : 0.1102

我们这里实验从1000次模拟到50000次模拟的结果,每次同样次数的模拟运行100遍。

  1. pathScenario = range(1000, 50000, 1000)
  2. numberOfTrials = 100
  3. confidenceIntervalUpper = []
  4. confidenceIntervalLower = []
  5. means = []
  6. for scenario in pathScenario:
  7. res = np.zeros(numberOfTrials)
  8. for i in range(numberOfTrials):
  9. res[i] = call_option_pricer_monte_carlo(spot, strike, maturity, r, vol, numOfPath = scenario)
  10. means.append(res.mean())
  11. confidenceIntervalUpper.append(res.mean() + 1.96*res.std())
  12. confidenceIntervalLower.append(res.mean() - 1.96*res.std())

蒙特卡洛方法会有收敛速度的考量。这里我们可以看到随着模拟次数的上升,仿真结果的置信区间也在逐渐收敛。

  1. pylab.figure(figsize = (12,8))
  2. tabel = np.array([means,confidenceIntervalUpper,confidenceIntervalLower]).T
  3. pylab.plot(pathScenario, tabel)
  4. pylab.title(u'期权计算蒙特卡洛模拟', fontproperties = font, fontsize = 18)
  5. pylab.legend([u'均值', u'95%置信区间上界', u'95%置信区间下界'], prop = font)
  6. pylab.ylabel(u'价格', fontproperties = font, fontsize = 15)
  7. pylab.xlabel(u'模拟次数', fontproperties = font, fontsize = 15)
  8. pylab.grid(True)

量化分析师的Python日记【第7天:Q Quant 之初出江湖】 - 图9

4. 计算隐含波动率

作为BSM期权定价最重要的参数,波动率σ是标的资产本身的波动率。是我们更关心的是当时的报价所反映的市场对波动率的估计,这个估计的波动率称为隐含波动率(Implied Volatility)。这里的过程实际上是在BSM公式中,假设另外4个参数确定,期权价格已知,反解σ:

量化分析师的Python日记【第7天:Q Quant 之初出江湖】 - 图10

由于对于欧式看涨期权而言,其价格为对应波动率的单调递增函数,所以这个求解过程是稳定可行的。一般来说我们可以类似于试错法来实现。在scipy中已经有很多高效的算法可以为我们所用,例如Brent算法:

  1. # 目标函数,目标价格由target确定
  2. class cost_function:
  3. def __init__(self, target):
  4. self.targetValue = target
  5. def __call__(self, x):
  6. return call_option_pricer(spot, strike, maturity, r, x) - self.targetValue
  7. # 假设我们使用vol初值作为目标
  8. target = call_option_pricer(spot, strike, maturity, r, vol)
  9. cost_sampel = cost_function(target)
  10. # 使用Brent算法求解
  11. impliedVol = brentq(cost_sampel, 0.01, 0.5)
  12. print u'真实波动率: %.2f' % (vol*100,) + '%'
  13. print u'隐含波动率: %.2f' % (impliedVol*100,) + '%'
  14. 真实波动率: 25.00%
  15. 隐含波动率: 25.00%