{% raw %} # 量化分析師的Python日記【第7天：Q Quant 之初出江湖】 > 來源：https://uqer.io/community/share/5514fc98f9f06c8f33904449 > 通過前幾日的學習，我們已經熟悉了Python中一些常用數值計算庫的用法。本篇中，作為Quant中的Q宗（[P Quant 和 Q Quant 到底哪個是未來？](http://www.zhihu.com/question/24820388))，我們將嘗試把之前的介紹的工具串聯起來，小試牛刀。 > 您將可以體驗到： > 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] 公式：  其中`S`為標的價格，`K`為執行價格，`r`為無風險利率，`τ=T?t`為剩余到期時間。 `N(x)`為標準正態分布的累積概率密度函數。`Call(S,K,r,τ,σ)`為看漲期權的價格。 ```py # 參數 spot = 2.45 strike = 2.50 maturity = 0.25 r = 0.05 vol = 0.25 ``` 觀察上面的公式，需要使用一些數學函數，我們把它分為兩部分： + `log`,`sqrt`,`exp`，這三個函數我們可以從標準庫`math`中找到 + 標準正態分布的累計概率密度函數，我們使用`scipy`庫中的`stats.norm.cdf`函數 ```py # 基于Black - Scholes 公式的期權定價公式 from math import log, sqrt, exp from scipy.stats import norm def call_option_pricer(spot, strike, maturity, r, vol): d1 = (log(spot/strike) + (r + 0.5 * vol *vol) * maturity) / vol / sqrt(maturity) d2 = d1 - vol * sqrt(maturity) price = spot * norm.cdf(d1) - strike * exp(-r*maturity) * norm.cdf(d2) return price ``` 我們可以使用這個函數計算我們關注期權的結果： ```py print '期權價格 : %.4f' % call_option_pricer(spot, strike, maturity, r, vol) 期權價格 : 0.1133 ``` ## 2. 使用numpy加速批量計算 大部分的時候，我們不止關心一個期權的價格，而是關心一個組合（成千上萬）的期權。我們想知道， 隨著期權組合數量的增長，我們計算時間的增長會有多塊？ ### 2.1 使用循環的方式 ```py import time import numpy as np portfolioSize = range(1, 10000, 500) timeSpent = [] for size in portfolioSize: now = time.time() strikes = np.linspace(2.0,3.0,size) for i in range(size): res = call_option_pricer(spot, strikes[i], maturity, r, vol) timeSpent.append(time.time() - now) ``` 從下圖中可以看出，計算時間的增長可以說是隨著組合規模的增長線性上升。 ```py from matplotlib import pylab import seaborn as sns font.set_size(15) sns.set(style="ticks") pylab.figure(figsize = (12,8)) pylab.bar(portfolioSize, timeSpent, color = 'r', width =300) pylab.grid(True) pylab.title(u'期權計算時間耗時（單位：秒）', fontproperties = font, fontsize = 18) pylab.ylabel(u'時間（s)', fontproperties = font, fontsize = 15) pylab.xlabel(u'組合數量', fontproperties = font, fontsize = 15) <matplotlib.text.Text at 0xdbad950> ```  ### 2.2 使用`numpy`向量計算 `numpy`的內置數學函數可以天然的運用于向量： ```py sample = np.linspace(1.0,100.0,5) np.exp(sample) array([ 2.71828183e+00, 1.52434373e+11, 8.54813429e+21, 4.79357761e+32, 2.68811714e+43]) ``` 利用 `numpy` 的數學函數，我們可以重寫原先的計算公式 `call_option_pricer`，使得它接受向量參數。 ```py # 使用numpy的向量函數重寫Black - Scholes公式 def call_option_pricer_nunmpy(spot, strike, maturity, r, vol): d1 = (np.log(spot/strike) + (r + 0.5 * vol *vol) * maturity) / vol / np.sqrt(maturity) d2 = d1 - vol * np.sqrt(maturity) price = spot * norm.cdf(d1) - strike * np.exp(-r*maturity) * norm.cdf(d2) return price ``` ```py timeSpentNumpy = [] for size in portfolioSize: now = time.time() strikes = np.linspace(2.0,3.0, size) res = call_option_pricer_nunmpy(spot, strikes, maturity, r, vol) timeSpentNumpy.append(time.time() - now) ``` 再觀察一下計算耗時，雖然時間仍然是隨著規模的增長線性上升，但是增長的速度要慢許多： ```py pylab.figure(figsize = (12,8)) pylab.bar(portfolioSize, timeSpentNumpy, color = 'r', width = 300) pylab.grid(True) pylab.title(u'期權計算時間耗時（單位：秒）- numpy加速版', fontproperties = font, fontsize = 18) pylab.ylabel(u'時間（s)', fontproperties = font, fontsize = 15) pylab.xlabel(u'組合數量', fontproperties = font, fontsize = 15) <matplotlib.text.Text at 0xe0ba090> ```  讓我們把兩次計算時間進行比對，更清楚的了解 `numpy` 計算效率的提升！ ```py fig = pylab.figure(figsize = (12,8)) ax = fig.gca() pylab.plot(portfolioSize, np.log10(timeSpent), portfolioSize, np.log(timeSpentNumpy)) pylab.grid(True) from matplotlib.ticker import FuncFormatter def millions(x, pos): 'The two args are the value and tick position' return '$10^{%.0f}$' % (x) formatter = FuncFormatter(millions) ax.yaxis.set_major_formatter(formatter) pylab.title(u'期權計算時間耗時（單位：秒）', fontproperties = font, fontsize = 18) pylab.legend([u'循環計算', u'numpy向量加速'], prop = font, loc = 'upper center', ncol = 2) pylab.ylabel(u'時間（秒)', fontproperties = font, fontsize = 15) pylab.xlabel(u'組合數量', fontproperties = font, fontsize = 15) <matplotlib.text.Text at 0xe0b6390> ```  ## 3. 使用`scipy`做仿真計算 期權價格的計算方法中有一類稱為 蒙特卡洛 方法。這是利用隨機抽樣的方法，模擬標的股票價格隨機游走，計算期權價格（未來的期望）。假設股票價格滿足以下的隨機游走：  仿真的方法可以模擬到期日的股票價格  這里的`z`是一個符合標準正態分布的隨機數。這樣我們可以計算最后的期權價格：  標準正態分布的隨機數獲取，可以方便的求助于 `scipy` 庫： ```py import scipy scipy.random.randn(10) array([ 0.36802702, 1.09560268, -1.0235275 , 0.15722882, 0.83718188, -0.27193135, -0.03485659, 1.02705248, 0.69479874, -0.35967107]) ``` ```py pylab.figure(figsize = (12,8)) randomSeries = scipy.random.randn(1000) pylab.plot(randomSeries) print u'均 值：%.4f' % randomSeries.mean() print u'標準差：%.4f' % randomSeries.std() 均 值：0.0336 標準差：0.9689 ```  結合 `scipy` `numpy` 我們可以定義基于蒙特卡洛的期權定價算法。 ```py # 期權計算的蒙特卡洛方法 def call_option_pricer_monte_carlo(spot, strike, maturity, r, vol, numOfPath = 5000): randomSeries = scipy.random.randn(numOfPath) s_t = spot * np.exp((r - 0.5 * vol * vol) * maturity + randomSeries * vol * sqrt(maturity)) sumValue = np.maximum(s_t - strike, 0.0).sum() price = exp(-r*maturity) * sumValue / numOfPath return price ``` ```py print '期權價格（蒙特卡洛） : %.4f' % call_option_pricer_monte_carlo(spot, strike, maturity, r, vol) 期權價格（蒙特卡洛） : 0.1102 ``` 我們這里實驗從1000次模擬到50000次模擬的結果，每次同樣次數的模擬運行100遍。 ```py pathScenario = range(1000, 50000, 1000) numberOfTrials = 100 confidenceIntervalUpper = [] confidenceIntervalLower = [] means = [] for scenario in pathScenario: res = np.zeros(numberOfTrials) for i in range(numberOfTrials): res[i] = call_option_pricer_monte_carlo(spot, strike, maturity, r, vol, numOfPath = scenario) means.append(res.mean()) confidenceIntervalUpper.append(res.mean() + 1.96*res.std()) confidenceIntervalLower.append(res.mean() - 1.96*res.std()) ``` 蒙特卡洛方法會有收斂速度的考量。這里我們可以看到隨著模擬次數的上升，仿真結果的置信區間也在逐漸收斂。 ```py pylab.figure(figsize = (12,8)) tabel = np.array([means,confidenceIntervalUpper,confidenceIntervalLower]).T pylab.plot(pathScenario, tabel) pylab.title(u'期權計算蒙特卡洛模擬', fontproperties = font, fontsize = 18) pylab.legend([u'均值', u'95%置信區間上界', u'95%置信區間下界'], prop = font) pylab.ylabel(u'價格', fontproperties = font, fontsize = 15) pylab.xlabel(u'模擬次數', fontproperties = font, fontsize = 15) pylab.grid(True) ```  ## 4. 計算隱含波動率 作為BSM期權定價最重要的參數，波動率`σ`是標的資產本身的波動率。是我們更關心的是當時的報價所反映的市場對波動率的估計，這個估計的波動率稱為隱含波動率（Implied Volatility）。這里的過程實際上是在BSM公式中，假設另外4個參數確定，期權價格已知，反解`σ`:  由于對于歐式看漲期權而言，其價格為對應波動率的單調遞增函數，所以這個求解過程是穩定可行的。一般來說我們可以類似于試錯法來實現。在`scipy`中已經有很多高效的算法可以為我們所用，例如Brent算法: ```py # 目標函數，目標價格由target確定 class cost_function: def __init__(self, target): self.targetValue = target def __call__(self, x): return call_option_pricer(spot, strike, maturity, r, x) - self.targetValue # 假設我們使用vol初值作為目標 target = call_option_pricer(spot, strike, maturity, r, vol) cost_sampel = cost_function(target) # 使用Brent算法求解 impliedVol = brentq(cost_sampel, 0.01, 0.5) print u'真實波動率： %.2f' % (vol*100,) + '%' print u'隱含波動率： %.2f' % (impliedVol*100,) + '%' 真實波動率： 25.00% 隱含波動率： 25.00% ``` {% endraw %}