优化期权定价码

Question

我是Python新手，正在尝试优化以下代码：

import sys
import numpy as np

def LSMPut(T, r, sigma, K, S0, TimeSteps, Paths, k):
    dt = T/TimeSteps
    t = np.arange(0, T+dt, dt).tolist()

    z=[np.random.standard_normal() for _ in range(Paths)]
    w = (r-sigma**2/2)*T + sigma*np.sqrt(T)*np.array(z)

    S = S0*np.exp(np.array(w))
    P=np.maximum(K-np.array(S),0)

    for i in range(TimeSteps-1, -1, -1):
        z=[np.random.standard_normal() for _ in range(Paths)]
        w = t[i]*np.array(w)/t[i+1] + sigma*np.sqrt(dt*t[i]/t[i+1])*np.array(z)

        S = S0*np.exp(np.array(w))
        itmPaths = [index for index,value in enumerate(K-np.array(S)) if value > 0]

        itmS = S[itmPaths]
        Pt = K - np.array(itmS)

        itmDiscP = P[itmPaths]*np.exp(-r*dt)

        A = BasisFunct(itmS, k)
        beta = np.linalg.lstsq(A,itmDiscP)[0]
        C = np.dot(A,beta)

        exPaths = [itmPaths[i] for i, value in enumerate(zip(Pt, C)) if value[0] > value[1]]
        restPaths = np.setdiff1d(np.arange(0, Paths-1, 1).tolist(), exPaths) # Rest of the paths

        P[exPaths] = [Pt[i] for i, value in enumerate(zip(Pt, C)) if value[0] > value[1]]
        P[restPaths] = np.array(P[restPaths])*np.exp(-r*dt)

        u=np.mean(P*np.exp(-r*dt))

    return u

def BasisFunct(X, k):
    Ones=[1 for _ in range(len(X))]
    if k == 1:
        A = np.column_stack((Ones,1 - np.array(X)))
    elif k == 2:
        A = np.column_stack((Ones,1 - np.array(X),1/2*(2-4*np.array(X) + np.array(X)**2)))
    elif k == 3:
        A = np.column_stack((Ones,1 - np.array(X),1/2*(2-4*np.array(X) + np.array(X)**2), 1/6*(6-18*np.array(X) + 9*np.array(X)**2-np.array(X)**3)))
    elif k == 4:
        A = np.column_stack((Ones,1 - np.array(X),1/2*(2-4*np.array(X) + np.array(X)**2), 1/6*(6-18*np.array(X) + 9*np.array(X)**2-np.array(X)**3),1/24*(24 - 96*np.array(X) + 72*np.array(X)**2 - 16*np.array(X)**3 + np.array(X)**4)))
    elif k == 5:
        A = np.column_stack((Ones,1 - np.array(X),1/2*(2-4*np.array(X) + np.array(X)**2), 1/6*(6-18*np.array(X) + 9*np.array(X)**2-np.array(X)**3),1/24*(24 - 96*np.array(X) + 72*np.array(X)**2 - 16*np.array(X)**3 + np.array(X)**4),1/120*(120-600*np.array(X)+600*np.array(X)**2-200*np.array(X)**3+25*np.array(X)**4-np.array(X)**5)))
    else:
            sys.exit("Too many basis functions requested")
    return A

print(LSMPut(1, 0.06, 0.15, 100, 90, 20, 1000000, 5))

代码的目的是计算美国看跌期权的价格。如果Paths > 1,000,000需要很长时间来执行，特别是当我执行灵敏度分析时。我想知道是否有任何方法来优化代码或加快处理时间。谢谢。

Answer 1

我将脚本的主要部分更改为：

from time import time
start = time()
val = LSMPut(1, 0.06, 0.15, 100, 90, 20, 40_000, 5)
print(f'Done in {time() - start:.3} seconds')
print(val)

重点是Python将打印计算运行所需的时间。修改前，脚本需要3.05秒才能在我的计算机上运行。

现在，我可以修改lsmput (我重新命名了LSMPut和其他一些函数和变量，因为它们的名称不应该包含大写字母，根据Python指南)，并查看哪些更改使其运行得更快。

下面是对运行时影响最大的更改。改变

z=np.array([np.random.standard_normal() for _ in range(paths)])

至

z = np.random.standard_normal(paths)

使运行时缩短到2.43秒。改变

ones=[1 for _ in range(len(X))]

至

ones = np.ones(len(X))

将运行时进一步减少到2.21秒。像这样重新计算X在basis_funct中的幂

X2 = X**2
X3 = X**3
X4 = X**4
X5 = X**5
...
A = np.column_stack((ones, 1 - X, 1/2 * (2 - 4*X + X2), ...

节省额外100毫秒。改变

itmPaths = [index for index,value in enumerate(K - S) if value > 0]

至

itmPaths = np.nonzero((K - S) > 0)[0]

使运行时降到大约1.60秒。我认为您现在看到了模式；每当您在numpy数组上循环时，就会有一个numpy内置函数来更好、更快地完成这项工作。

然后在计算rest_paths时：

rest_paths = np.setdiff1d(np.arange(0, paths - 1, 1), exPaths)

这里的paths - 1是故意的吗？在我看来是个虫子。假设它是一个bug，更新P可以更有效，如下所示：

mask = np.zeros(P.shape, dtype = bool)
...
mask.fill(False)
mask[itmPaths[Pt > C]] = True

P[mask] = Pt[Pt > C]
P[~mask] *= e_r_dt

最终代码

修改后的代码在我的计算机上运行约420毫秒。注意，我已经从k == 5函数中删除了除basis_funct 1之外的所有分支。对您来说，将其添加回应该很简单。

import sys
import numpy as np

# Uncomment this line during development to always get the same random
# numbers.
# np.random.seed(1234)

def lsmput(T, r, sigma, K, S0, TimeSteps, paths, k):
    dt = T/TimeSteps
    t = np.arange(0, T+dt, dt).tolist()

    z = np.random.standard_normal(paths)
    w = (r-sigma**2/2)*T + sigma*np.sqrt(T)*z

    S = S0*np.exp(w)
    P = np.maximum(K - S,0)

    e_r_dt = np.exp(-r * dt)
    mask = np.zeros(P.shape, dtype = bool)
    for i in range(TimeSteps-1, -1, -1):
        z = np.random.standard_normal(paths)
        w = t[i]*w/t[i+1] + sigma*np.sqrt(dt*t[i]/t[i+1])*z

        S = S0 * np.exp(w)
        itmPaths = np.nonzero(K > S)[0]
        itmS = S[itmPaths]
        Pt = K - itmS

        itmDiscP = P[itmPaths] * e_r_dt

        A = basis_funct(itmS, k)
        beta = np.linalg.lstsq(A, itmDiscP)[0]
        C = np.dot(A, beta)

        mask.fill(False)
        mask[itmPaths[Pt > C]] = True

        P[mask] = Pt[Pt > C]
        P[~mask] *= e_r_dt
    return np.mean(P * e_r_dt)

def basis_funct(X, k):
    ones = np.ones(len(X))
    assert k == 5
    X2 = X**2
    X3 = X**3
    X4 = X**4
    X5 = X**5
    A = np.column_stack((ones, 1 - X,
                         1/2 * (2 - 4*X + X2),
                         1/6 * (6 - 18*X + 9*X2 - X3),
                         1/24 * (24 - 96*X + 72*X2 - 16*X3 + X4),
                         1/120 * (120 - 600*X + 600*X2 - 200*X3 + 25*X4 - X5)))
    return A

from time import time
start = time()
val = lsmput(1, 0.06, 0.15, 100, 90, 20, 40_000, 5)
print(f'Done in {time() - start:.3} seconds')
print(val)