Optimize option pricing code

I changed the main part of the script to:

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)

The point is that Python will print how long the calculation took to run. Before modifications, the script takes 3.05 seconds to run on my computer.

Now I can modify lsmput (I renamed LSMPut and some other functions and variables because they names shouldn't contain capital letters according to the Python style guide) and see which changes makes it run faster.

Here follows the changes that impacted the runtime the most. Changing

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

to

z = np.random.standard_normal(paths)

brings the runtime down to 2.43 seconds. Changing

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

to

ones = np.ones(len(X))

reduces the runtime further to 2.21 seconds. Precalculating the powers of X in basis_funct like this

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), ...

saves an additional 100 milliseconds. Changing

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

to

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

brings the runtime down to about 1.60 seconds. I think you see the pattern now; anytime you loop over a numpy array there's a numpy builtin function to do the job better and faster.

Then in the calculation of rest_paths:

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

Is the paths - 1 intentional here? Looks like a bug to me. Assuming it is a bug, updating P can be done more efficiently, like this:

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

Final code

The modified code runs in about 420 milliseconds on my computer. Note that I have removed all branches except for the k == 5 one from the basis_funct function. It should be straightforward for you to add it back.

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)