第8章 高效能的Python

許多高效能庫可以用於加速Python程式碼的執行：
• Cython 用於合併Py由on和C語言靜態編譯範型。
• IPython.parallel 用於在本地或者在群集上並行執行程式碼／函式。
• numexpr 用於快速數值運算。
• multiprocessing Python內建的（本地）並行處理模組。
• Numba 用於為CPU動態編譯Python程式碼。
• NumbaPro 用於為多核CPU和GPU動態編譯Python程式碼。

定義一個方便的函式，可以系統性地比較在相同或者不同資料集上執行不同函式的效能：

def perf_comp_data
 
(func_list, data_list, rep=3, number=1):
    """
    Function to compare the performance of different function.
    :param func_list: list with function names as strings
    :param data_list: list with data set names as strings
    :param rep: number of repetitions of the whole comparison
    :param number: number of executions for every function
    :return:
    """
 

    from timeit import repeat
    res_list = {}
    for name in enumerate(func_list):
        stmt = name[1] + '(' + data_list[name[0]] + ')'
        setup = "from __main__ import " + name[1] + ', ' + data_list[name[0]]
        results = repeat(stmt=stmt, setup=setup, repeat=rep, number=number)
        res_list[name[1
 
]] = sum(results) / rep
    res_sort = sorted(res_list.items(), key=lambda item: item[1])
    for item in res_sort:
        rel = item[1] / res_sort[0][1]
        print('function:' + item[0][1] + ', av.item sec: %9.5f, ' % item[1] + 'relative: %6.1f' % rel)

8.1 Python範型與效能

在金融學中與其他科學及資料密集學科一樣， 大資料集上的數值計算相當費時。舉個例子， 我們想要在包含 50 萬個數值的陣列上求取某個複雜數學表示式的值。我們選擇公式中的表示式，它的每次計算都會帶來一定的計算負擔。除此之外，該公式沒有任何特殊的含義。
數學表示式示例

def perf_comp_data(func_list, data_list, rep=3, number=1):
    """
    Function to compare the performance of different function.
    :param func_list: list with function names as strings
    :param data_list: list with data set names as strings
    :param rep: number of repetitions of the whole comparison
    :param number: number of executions for every function
    :return:
    """
    from timeit import repeat
    res_list = {}
    for name in enumerate(func_list):
        stmt = name[1] + '(' + data_list[name[0]] + ')'
        setup = "from __main__ import " + name[1] + ', ' + data_list[name[0]]
        results = repeat(stmt=stmt, setup=setup, repeat=rep, number=number)
        res_list[name[1]] = sum(results) / rep
    res_sort = sorted(res_list.items(), key=lambda item: item[1])
    for item in res_sort:
        rel = item[1] / res_sort[0][1]
        print('function:' + item[0] + ', av.item sec: %9.5f, ' % item[1] + 'relative: %6.1f' % rel)


# 8.1 Python範型與效能

from math import *


# 很容易轉換為一個Python函式
def f(x):
    return abs(cos(x)) ** 0.5 + sin(2 + 3 * x)

# 使用range函式，我們可以高效地生成一個包含 50 萬個數值的列表物件
I = 500000
a_py = range(I)

# 包含顯式迴圈的標準Python函式
def f1(a):
    res = []
    for x in a:
        res.append(f(x))
    return res

# 包含隱含迴圈的迭代子方法
def f2(a):
    return [f(x) for x in a]

# 包含隱含迴圈、使用eval的選代子方法
def f3(a):
    ex = 'abs(cos(x))**0.5+sin(2+3*x)'
    return [eval(ex) for x in a]

# Numy向量化實現
import numpy as np
a_np = np.arange(I)
def f4(a):
    return (np.abs(np.cos(a)) ** 0.5 + np.sin(2 + 3 * a))

# 專用庫numexpr求數值表示式的值。 這個庫內建了多執行緒執行支援
# numexpr單執行緒實現
import numexpr as ne
def f5(a):
    ex='abs(cos(a))**0.5+sin(2+3*a)'
    ne.set_num_threads(1)
    return ne.evaluate(ex)

# nwexpr多執行緒實現
def f6(a):
    ex = 'abs(cos(a))**0.5+sin(2+3*a)'
    ne.set_num_threads(16)
    return ne.evaluate(ex)

%%time
r1=f1(a_py)
r2=f2(a_py)
r3=f3(a_py)
r4=f4(a_np)
r5=f5(a_np)
r6=f6(a_np)
# Wall time: 35.1 s

# NumPy函式alJclose可以輕鬆地檢查兩個（類） ndarray物件是否包含相同資料
np.allclose(r1,r2)
# True
np.allclose(r1,r3)
# True
np.allclose(r1,r4)
# True
np.allclose(r1,r5)
# True
np.allclose(r1,r6)
# True

# 使用perf_comp_data函式
func_list=['f1','f2','f3','f4','f5','f6']
data_list=['a_py','a_py','a_py','a_np','a_np','a_np']
perf_comp_data(func_list,data_list)
# function:f6, av.item sec:   0.01623, relative:    1.0
# function:f5, av.item sec:   0.04650, relative:    2.9
# function:f4, av.item sec:   0.07293, relative:    4.5
# function:f2, av.item sec:   1.17137, relative:   72.2
# function:f1, av.item sec:   1.33291, relative:   82.1
# function:f3, av.item sec:  33.47790, relative: 2062.2

8.2 記憶體佈局與效能

import numpy as np
np.zeros((3,3),dtype=np.float64,order='C')
# array([[ 0.,  0.,  0.],
#        [ 0.,  0.,  0.],
#        [ 0.,  0.,  0.]])

# 元素在記憶體中儲存的順序：C表示類似C（行優先）
c=np.array([[1.,1.,1.],
            [2.,2.,2.],
            [3.,3.,3.]],order='C')
# F表示類似Fortran （列優先）
f=np.array([[1.,1.,1.],
            [2.,2.,2.],
            [3.,3.,3.]],order='F')

x = np.random.standard_normal((3, 1500000))
C = np.array(x, order='C')
F = np.array(x, order='F')
x = 0.0

%timeit C.sum(axis=0)  # 10 loops, best of 3: 19.3 ms per loop
%timeit C.sum(axis=1)  # 100 loops, best of 3: 10.3 ms per loop
# 第一個軸上計算總和比第二個軸慢了將近一倍
%timeit C.std(axis=0)  # 10 loops, best of 3: 112 ms per loop
%timeit C.std(axis=1)  # 10 loops, best of 3: 57.6 ms per loop

%timeit F.sum(axis=0)  # 10 loops, best of 3: 70.7 ms per loop
%timeit F.sum(axis=1)  # 10 loops, best of 3: 84.2 ms per loop
# 兩個軸的相對差值並不太大
%timeit F.std(axis=0)  # 1 loop, best of 3: 253 ms per loop
%timeit F.std(axis=1)  # 1 loop, best of 3: 227 ms per loop

# 與類似C的佈局相比， 類似F這種佈局的效能更差

8.3 平行計算

8.3.1 蒙特卡洛演算法

期權的蒙特卡洛估值是導致高計算負擔的金融演算法之一。作為特例，我們選擇Black-Scholes-Meron設定下的歐式看漲期權價值蒙特卡洛估值函式。在這種設定下，所要估值的期權標的遵循隨機微分方程式（SDE），如下公式。St是時間t的標的價值；r是一個常數——無風險短期利率；σ是恆定瞬時波動率；Z是布朗運動。
Black-Scholes-Metron SDE

歐式看漲期權的蒙特卡洛估算函式

def bsm_mcs_valuation(strike):
    """
    Dynamic Black-Scholes-Merton Monte Carlo estimator for European calls.
    :param strike:
    :return:
    """
    import numpy as np
    S0=100.;T=1.0;r=0.05;vola=0.2
    M=50;I=2000
    dt=T/M
    rand=np.random.standard_normal((M+1,I))
    S=np.zeros((M+1,I));S[0]=S0
    for t in range(1,M+1):
        S[t]=S[t-1]*np.exp((r-0.5*vola**2)*dt+vola*np.sqrt(dt)*rand[t])
    value=(np.exp(-r*T)*np.sum(np.maximum(S[-1]-strike,0))/I)
    return value

8.3.2 順序化計算

作為基準用例， 我們對不同行權價的100種期權進行估值。seq_value函式計算蒙特卡洛估算函式。返回包含行權價和估值結果的列表物件：

def seq_value(n):
    """
    Sequential option valuation
    :param n: number of option valuations/strikes
    :return:
    """
    strikes=np.linspace(80,120,n)
    option_values=[]
    for strike in strikes:
        option_values.append(bsm_mcs_valuation(strike))
    return strikes,option_values

n=100 # number of options to be valued
%time strikes,option_values_seq=seq_value(n)
# Wall time: 1.39 s

import matplotlib.pyplot as plt
plt.figure(figsize=(8,4))
plt.plot(strikes,option_values_seq,'b')
plt.plot(strikes,option_values_seq,'r.')
plt.grid(True)
plt.xlabel('strikes')
plt.ylabel('European call option values')

8.4 多處理

有時候在本地並行執行程式碼是很有益的。 這就是 “標準” Pthon 模
塊 multiprocessing 的用武之地：

import numpy as np
import multiprocessing as mp
import math
import matplotlib.pyplot as plt
def simulate_geometric_brownian_motion(p):
    M,I=p
    # time steps,paths
    S0=100;r=0.05;sigma=0.2;T=1.0
    # model parameters
    dt=T/M
    paths=np.zeros((M+1,I))
    paths[0]=S0
    for t in range(1,M+1):
        paths[t]=paths[t-1]*np.exp((r-0.5*sigma**2)*dt+sigma*math.sqrt(dt)*np.random.standard_normal(I))
    return paths

paths=simulate_geometric_brownian_motion((5,2))
paths
# array([[ 100.        ,  100.        ],
#        [  98.75585496,   86.36316092],
#        [ 109.5045796 ,   82.00664539],
#        [  92.85348223,   81.23649105],
#        [  73.79002067,   81.99661207],
#        [  67.4225339 ,   89.39928928]])

if __name__ == '__main__':
    I=10000  # number of paths
    M=100  # number of time steps
    t=100  # number of tasks/simulations
    # running on server with 8 cores/16 threads
    from time import time
    times=[]
    for w in range(1,17):
        t0=time()
        pool=mp.Pool(processes=w)
        # the pool of workers
        result = pool.map(simulate_geometric_brownian_motion,t*[(M,I),])
        # the mapping of the function to the list of parameter tuples
        times.append(time()-t0)

    plt.plot(range(1, 17), times)
    plt.plot(range(1, 17), times, 'ro')
    plt.grid(True)
    plt.xlabel('number of processes')
    plt.ylabel('time in seconds')
    plt.title('%d Monte Carlo simulations' % t)

效能和可用核心數量成正比。 不過， 超執行緒在本例中不能帶來太多好處（甚至更糟）

簡易的並行化
金融學中的許多問題可以應用簡單的並行化技術， 例如， 在演算法的不同例項之間沒有共享資料時。Python的multiprocesing模組可以高效地利用現代硬體架構的能力， 一般不需要改變基本演算法或者並行執行的Python函式．

8.5 動態編譯

Numba (http://numba.pydata.org）是開源 、 NumPy 感知的優化Python 程式碼編譯器。它使用 LLVM 編譯器基礎架構，將Python 位元組程式碼編譯專門用於 NumPy執行時和SciPy模組的機器程式碼。

8.5.1 介紹性示例

from math import cos,log
def f_py(I,J):
    res=0
    for i in range(I):
        for j in range(J):
            res+=int(cos(log(1)))
    return res

I,J=5000,5000
%time f_py(I,J)
# Wall time: 30 s
# 25000000

def f_np(I, J):
    a = np.ones((I, J), dtype=np.float64)
    return int(np.sum(np.cos(np.log(a)))), a

%time res, a = f_np(I, J)
# Wall time: 1.34 s

a.nbytes
# 200000000

import numba as nb
f_nb=nb.jit(f_py)

%time f_nb(I,J)
# Wall time: 741 ms
# 25000000

func_list=['f_py','f_np','f_nb']
data_list=3*['I,J']
perf_comp_data(func_list,data_list)

# function:f_nb, av.item sec:   0.00001, relative:    1.0
# function:f_np, av.item sec:   1.36470, relative: 156714.8
# function:f_py, av.item sec:  29.53817, relative: 3391993.0

速效方法
改善（數值演算法）效能的許多方法都需要花費可觀的精力。利用Python和Numba， 就有了需要最少精力的 一種方法一一般來說， 只需要匯入境庫和一行附加程式碼 ． 它不能用於所有型別演算法，但是往往值得（簡單地） 一試， 有時候確實能夠快速取得效果。

8.5.2 二項式期權定價方法

前面使用蒙特卡洛模擬方法、 利用平行計算估計歐式看漲期權的價值。估算期權價值的另一種流行數值方法是二項式期權定價模型。 這種模型和Black-Scholes-Meron設定一樣有風險資產（指數或者股票）以及無風險資產（債券）。 和蒙特卡洛方法一樣，從當天到期權到期日的時間間隔被分為通常等距的子間隔Δt， 如果時間s的指數水平為Ss， 則 t = s+Δt 時的指數水平為St = Ss·m ， 其中m是從｛u， d } 中隨機選取( )。r是 一個常數一一無風險利率。 風險中立的上漲概率為。

對該模型進行引數化：

# model & option parameters
S0=100.  # initial index level
T=1.  # call option maturity
r=0.05  # constant short rate
vola=0.20  # constant volatility factor of diffusion

# time parameters
M=1000 #time steps
dt=T/M # length of time interval
df= exp(-r*dt)  # discount factor per time interval

# binomial parameters
u=exp(vola*sqrt(dt))  # up-movement
d=1/u # down-movement
q=(exp(r*dt)-d)/(u-d) # martingale probability

歐式期權二項式演算法的實現主要包含如下部分：
指數水平模擬
連步模擬指數水平。
內在價值計算
計算到期日和每個時間步的內在價值。
風險中性折算
逐步折算（預期）內在價值直到達到現值。

在Python中，這可能採取函式binomial_py中的形式。該面數使用NumPy ndarray物件作為基本資料結構， 並實現3個不同的巢狀迴圈，以實現上述的3個步驟：

from math import *
# model & option parameters
S0=100.  # initial index level
T=1.  # call option maturity
r=0.05  # constant short rate
vola=0.20  # constant volatility factor of diffusion

# time parameters
M=1000 #time steps
dt=T/M # length of time interval
df= exp(-r*dt)  # discount factor per time interval

# binomial parameters
u=exp(vola*sqrt(dt))  # up-movement
d=1/u # down-movement
q=(exp(r*dt)-d)/(u-d) # martingale probability

import numpy as np
def binomial_py(strike):
    """
    Binomial option pricing via looping
    :param strike:float
        strike price of the European call option
    :return:
    """
    # LOOP 1 - Index Levels
    S=np.zeros((M+1,M+1),dtype=np.float64)
    # index level array
    S[0,0]=S0
    z1=0
    for j in range(1,M+1,1):
        z1=z1+1
        for i in range(z1+1):
            S[i,j]=S[0,0]*(u**j)*(d**(i*2))

    # LOOP 2 - Inner Values
    iv = np.zeros((M+1,M+1),dtype=np.float64)
    # inner value array
    z2=0
    for j in range(0,M+1,1):
        for i in range(z2+1):
            iv[i,j]=max(S[i,j]-strike,0)
        z2=z2+1

    # LOOP 3 - Valuation
    pv = np.zeros((M+1,M+1),dtype=np.float64)
    # present value array
    pv[:,M]=iv[:,M] # initialize last time point
    z3=M+1
    for j in range(M-1,-1,-1):
        z3=z3-1
        for i in range(z3):
            pv[i,j]=(q*pv[i,j+1]+(1-q)*pv[i+1,j+1])*df
    return pv[0,0]

上函式使用前面指定的引數． 返回歐式看漲期權的現值：

 %time round(binomial_py(100),3)
# Wall time: 4.23 s
# 10.449

將這個結果與蒙特卡洛函式bsm_mcs_valuation返回的估算結果比較

%time round(bsm_mcs_valuation(100),3)
# Wall time: 15 ms
# 10.183

兩個值很類似，它們只是“相似” 而不是相同。是因為蒙特卡洛估值和bsm_mcs_valuaton所實現的演算法都不是很精確， 不同的隨機數會導致（稍有）不同的估算結果， 對於健全的蒙特卡洛估算來說， 每次模擬使用2000Q條路徑也可能略少一些（但是可以得到較高的估值速度）。 相比之下， 本例中的二項式期權估價使用 1000 個時間步已經相當精確，但是花費的時間也長得多。

可以嘗試 NumPy 向量化技術，從二項式方法中得到同樣精確、但是速度更快的結果。 binomial_np 函式初看有些神祕，但是， 當執行單獨的構建步驟並檢查結果， 後臺（ NumPy ）發生的操作就顯而易見了：

def binomial_np(strike):
    """
    Binomial option pricing with NumPy
    :param strike: float
        strike price of the European call option
    :return: 
    """
    # Index Levels with NumPy
    mu=np.array(M+1)
    mu=np.resize(mu,(M+1,M+1))
    md=np.transpose(mu)
    mu=u**(mu-md)
    S=S0*mu*md

    # Valuation Loop
    pv=np.maximum(S-strike,0)
    z=0
    for t in range(M-1,-1,-1): # backward iteration
        pv[0:M-z,t]=(q*pv[0:M-z,t+1]+(1-q)*pv[1:M-z+1,t+1])*df
        z+=1
    return pv[0,0]

下面我們簡單地看看後臺的情況。 為了簡潔和易於理解， 只考慮 M = 4 的時間步。第一步如下：

# 第一步
M = 4  # four time steps only
mu = np.arange(M + 1)
mu
# array([0, 1, 2, 3, 4])

# 第二步
mu = np.resize(mu, (M + 1, M + 1))
mu
# array([[0, 1, 2, 3, 4],
#        [0, 1, 2, 3, 4],
#        [0, 1, 2, 3, 4],
#        [0, 1, 2, 3, 4],
#        [0, 1, 2, 3, 4]])

# 第三步
md = np.transpose(mu)
md
# array([[0, 0, 0, 0, 0],
#        [1, 1, 1, 1, 1],
#        [2, 2, 2, 2, 2],
#        [3, 3, 3, 3, 3],
#        [4, 4, 4, 4, 4]])


# 第四步
mu = u ** (mu - md)
mu.round(3)
# array([[ 1.   ,  1.006,  1.013,  1.019,  1.026],
#        [ 0.994,  1.   ,  1.006,  1.013,  1.019],
#        [ 0.987,  0.994,  1.   ,  1.006,  1.013],
#        [ 0.981,  0.987,  0.994,  1.   ,  1.006],
#        [ 0.975,  0.981,  0.987,  0.994,  1.   ]])


# 第五步
md = d ** md
md.round(3)
# array([[ 1.   ,  1.   ,  1.   ,  1.   ,  1.   ],
#        [ 0.994,  0.994,  0.994,  0.994,  0.994],
#        [ 0.987,  0.987,  0.987,  0.987,  0.987],
#        [ 0.981,  0.981,  0.981,  0.981,  0.981],
#        [ 0.975,  0.975,  0.975,  0.975,  0.975]])

# 最後將所有步驟聚合起來
S = S0 * mu * md
S.round(3)
# array([[ 100.   ,  100.634,  101.273,  101.915,  102.562],
#        [  98.743,   99.37 ,  100.   ,  100.634,  101.273],
#        [  97.502,   98.121,   98.743,   99.37 ,  100.   ],
#        [  96.276,   96.887,   97.502,   98.121,   98.743],
#        [  95.066,   95.669,   96.276,   96.887,   97.502]])

在ndarray物件S中，只有上三角矩陣是重要的。雖然這種方法進行的計算多於原則上的需要，但是這種方法和預計的一樣，比嚴重依賴Python級別巢狀迴圈的第一個版本快得多：

M=1000 # reset number of time steps
%time round(binomial_np(100),3)
# Wall time: 1.03 s
# 0.0 這個結果不太對

Numba在金融演算法中也是很重要

binomial_nb=nb.jit(binomial_py)
%time round(binomial_nb(100),3)
# Wall time: 1.58 s
# 10.449

我們還沒有看出對NumPy向盤化版本有多少加速效果，因為第一次呼叫編譯後的函式涉及一些開銷。因此，使用 prf_comp_func函式應該更現實地揭示出不同實現的效能對比。顯然，Numba編譯版本確實明顯快於NumPy版本：

func_list=['binomial_py','binomial_np','binomial_nb']
K=100
data_list=3*['K']
perf_comp_data(func_list,data_list)
# function:binomial_nb, av.item sec:   0.12641, relative:    1.0
# function:binomial_np, av.item sec:   0.96281, relative:    7.6
# function:binomial_py, av.item sec:   4.26450, relative:   33.7

我們可以得出如下結論：
效率：使用Nnmba只需要花費很少的額外精力。原始函式往往完全不需要改變；你所需要做的就是呼叫jlt函式。
加速：Numba往往帶來執行速度的顯著提高．不僅和純Py曲。n相比是如此．即使對向量化的NumPy實現也有明顯優勢。
記憶體：使用Numba不需要初始化大型陣列物件；編譯器專門為手上的問題生成機器程式碼（和Numpy的 “通用 ” 函式相比）並維持和純Python相同的記憶體效率。

8.6 用Cython進行靜態編譯

Numba的優勢是對任意函式應用該方法毫不費力。但是，Numba只能為某些問題 “毫不費力 ” 地產生顯著的效能改善。另一種方法更為靈活，但是也需要更多精力，這就是通過Cython的靜態編譯， Cython是Python和C語言的混血兒。從Python的角度看，需要注意的主要不同是靜態型別宣告（和 C語言相同）和一個單獨的編譯步驟（和任何編譯語言相同）。

# 正常的Python程式碼：
def f_py(I, J):
    res = 0.  # we work on a float object
    for i in range(I):
        for j in range(J * I):
            res += 1
    return res

I,J=500,500
%time f_py(I,J)
# Wall time: 9 s
# Out[4]: 125000000.0

使用Cython靜態型別宣告的巢狀迴圈
建立一個名為nested_loop.pyx的檔案

# Nested loop example with Cython
# nested_loop.pyx
def f_cy(int I,int J):
    cdef double res = 0
    # double float nuch slower than in or long
    for i in range(I):
        for j in range(J*I):
            res+=1
    return res

在這個簡單的例子中不需要任何特殊的C模組 ，導人模組有一種簡單的方法一一就是通過pyximport

import pyximport

pyximport.install()

# 直接從Cython模組中導人
import sys

sys.path.append('E:\python_for_finance\chapter08')
# path to the Cython script
# not needed if in same directory

from nested_loop import f_cy

# 如果報錯： Unable to find vcvarsall.bat
# 是因為找不到vcvarsall.bat，這個問題件是在VS目錄下的，
# 我的VS2015的目錄是 E:\Program Files (x86)\Microsoft Visual Studio 14.0
# 但是裡面沒有這個檔案。這時可以在VS中新增一個C++專案，VS會提醒安裝新的軟體包，繼續安裝就可以了
# 安裝好後就可以搜到vcvarsall.bat了
# 然後新增環境變數
# VS90COMNTOOLS : E:\Program Files (x86)\Microsoft Visual Studio 14.0\VC

I, J = 500, 500
%time res = f_cy(I, J)
# Wall time: 131 ms
res
# 125000000.0

在IPython Notebook中工作時,使用Cython有一個更便利的方法：

%load_ext Cython
%%cython
def f_cy(int I,int J):
    cdef double res = 0
    # double float nuch slower than in or long
    for i in range(I):
        for j in range(J*I):
            res+=1
    return res

I, J = 500, 500
%time res = f_cy(I, J)
# Wall time: 131 ms 效能結果相同
res
# 125000000.0

看看Numba在這種情況下能起什麼作用

import numba as nb
f_nb=nb.jit(f_py)

%time res=f_nb(I,J)
# Wall time: 947 ms
# 第一次呼叫函式時， 效能比 Cython 版本差（第一次呼叫 Numba 編譯函式總是有某些開銷）
%time res=f_nb(I,J)
# Wall time: 131 ms 再次呼叫，效能就相同了

8.7 在GPU上生成隨機數

最後 一個主題是使用裝置進行大規模的並行操作——也就是通用圖形處理單元(GPGPU 或者簡稱 GPU ）。 要使用 Nvidia GPU ， 就必須安裝 CUDA （統一計算裝置架構， https://developer.nvidia..com ）。利用 Nvidia GPU 的簡單方法之一是使用 NumbaPro,這個由 Contnuum Analytics 開發的高效能庫為 GPU（或者多核 CPU ）動態編譯 Python。
有一個金融領域可以從 GPU 的使用中得到很大的好處：蒙特卡洛模擬。 特別是（偽）隨機數生成。