High Performance Python: Numpy and Numba

2024年3月11日 03:45

Pythonの高速化には、Numpyによるベクトル化と、llvmを使い動的にコンパイルするNumbaの二つがコードを大きく書き換える必要もなく手軽と言える。

素数判定

def isPrime(I):
    if I % 2 == 0: return False  
    for i in range(3, int(I ** 0.5) + 1, 2):  
        if I % i == 0: return False  
    return True

ループ関数の中にif文が入っていることから、ベクトル化には向いていないこの関数をNumbaで高速化する。

from numba import jit

isPrime_nb = jit(isPrime)
%time isPrime_nb(p2)
%timeit isPrime(p2)

Numbaの最初の呼び出しは、コンパイルが入るから時間がかかるが、2回目はすでにコンパイルされていることから、早くなる。
　Numba無しと有りとでの比較が下の結果である。

$${\pi}$$推定

モンテカルロ法によって、単位円領域と$${2\times2}$$正方形領域の面積比から$${\pi}$$を推定する方法をNumpyとNumbaで比較する。
単位円の面積は$${\pi r^2}$$より$${\pi}$$で、単位円が内接する$${2\times 2}$$の正方形の面積は$${4}$$である。中心を原点に置く正方形内にランダムで$${N}$$個の点をとり、同じく原点を中心とする単位円の中に入っている点の数$${N_c}$$との比が、$${lim_{N\rightarrow\infty} \frac{N_c}{N}\rightarrow \frac{\pi}{4}}$$であることを使う。

def mcs_pi_py(n):
    circle = 0
    for _ in range(n):
        x, y = random.random(), random.random()
        if (x ** 2 + y ** 2)  <= 1:
            circle += 1
    return (4 * circle) / n

これを同じく、Numba無しと有りで比較した結果が以下である。

二項モデル

リスク資産と無リスク資産からなるポートフォリオの価値を、リスク資産の上下モデルから計算する。
　金利を$${r}$$、リスク資産のボラティリティを$${\sigma}$$とする。各ノード間を$${\Delta t}$$とする。
　ある時点$${t}$$でリスク資産の価値が$${S_t}$$で与えられたとすれば、$${t+\Delta t}$$でのリスク資産の価値は、$${S_{t+\Delta t}=m \cdot S_t }$$で与えられる。ここで、$${m}$$は、$${0 < d < e^{r\Delta t} < u = e^{\sigma \sqrt{\Delta t}}}$$を満たす$${u}$$か$${d}$$のどちらかを取る。

S0 = 36.  
T = 1.0  
r = 0.06  
sigma = 0.2 

def simulate_tree(M):
    dt = T / M  
    u = math.exp(sigma * math.sqrt(dt))  
    d = 1 / u  
    S = np.zeros((M + 1, M + 1))
    S[0, 0] = S0
    z = 1
    for t in range(1, M + 1):
        for i in range(z):
            S[i, t] = S[i, t-1] * u
            S[i+1, t] = S[i, t-1] * d
        z += 1
    return S

　結果の行列$${S}$$の上三角にしか意味はないことに注意する。

np.set_printoptions(formatter={'float':
                               lambda x: '%6.2f' % x})  
simulate_tree(4)

これをNumpyを使ってベクトル化する。

def simulate_tree_np(M):
    dt = T / M
    up = np.arange(M + 1)
    up = np.resize(up, (M + 1, M + 1))
    down = up.transpose() * 2
    S = S0 * np.exp(sigma * math.sqrt(dt) * (up - down))
    return np.triu(S)

Numbaを使った二項ツリー

simulate_tree_nb = jit(simulate_tree)
simulate_tree_nb(4)

これらの時間測定の結果は以下の通り。

モンテカルロシミュレーション

ブラックショールズ微分方程式によるオプション評価式
$${dS=r S_t dt + \sigma S_t dZ_t}$$を離散しする。
$${S_t=\displaystyle{S_{t-\Delta t}\exp\left((r-\frac{\sigma^2}{2}\Delta t + \sigma\sqrt{\Delta t}z\right)}}$$
$${z\sim \mathcal{N}[0,1]}$$
　これをモンテカルロシミュレーションに実装する。

M = 100  
I = 50000
def mcs_simulation_py(p):
    M, I, S0 = p
    dt = T / M
    S = np.zeros((M + 1, I))
    S[0] = S0
    rn = np.random.standard_normal(S.shape)  
    for t in range(1, M + 1):  
        for i in range(I):  
            S[t, i] = S[t-1, i] * math.exp((r - sigma ** 2 / 2) * dt +
                                         sigma * math.sqrt(dt) * rn[t, i])  
    return S

Numpyを使って、$${I}$$をベクトル化し、時間ループのみにする。

def mcs_simulation_np(p):
    M, I, S0 = p
    dt = T / M
    S = np.zeros((M + 1, I))
    S[0] = S0
    rn = np.random.standard_normal(S.shape)
    for t in range(1, M + 1):  
        S[t] = S[t-1] * np.exp((r - sigma ** 2 / 2) * dt +
                               sigma * math.sqrt(dt) * rn[t]) 
    return S

EWMA

pandasのDataFrameから、指数荷重移動平均EWMAを計算する。
$${EWMA_0=S_0}$$
$${EWMA_t = \alpha S_t +(1-\alpha)EWMA{t-1}, \ t =[1,\cdots,T]}$$

def ewma_py(x, alpha):
    y = np.zeros_like(x)
    y[0] = x[0]
    for i in range(1, len(x)):
        y[i] = alpha * x[i] + (1-alpha) * y[i-1]
    return y

この計算には、一時間毎の暗号資産を使う。

import yfinance as yfin
from pandas_datareader import data as pdr

yfin.pdr_override()

sym='BTC-USD'
cryptocurrencies = ['BTC-USD','ETH-USD']

df_data= pdr.get_data_yahoo(cryptocurrencies, start='2024-01-01', end='2024-03-01',interval='1h')['Adj Close']
df_data.columns=cryptocurrencies
df_data

これを使ったPythonのループ計算とNumbaを通したパフォーマンスの向上は以下の通り。

この記事が気に入ったらサポートをしてみませんか？