barrier_numba.py

barrier_numba.py

from random import random
from math import sqrt, log, exp
from numba import njit, jit, prange

@njit(fastmath=True)
def box_muller_rand():
    while True:
        x = random() * 2.0 - 1
        y = random() * 2.0 - 1
        d = x * x + y * y 
        if d < 1:
            return x * sqrt(-2 * log(d) / d)

@njit(fastmath=True)
def path_final_min_max(initial, steps, sdt, volatility, drift):
    current = initial
    minimum = current
    maximum = current
    for i in range(steps - 1):
        current = current * drift * exp(sdt * volatility * box_muller_rand())
        if minimum > current:
            minimum = current
        if maximum < current:
            maximum = current

    return current, minimum, maximum

@njit(fastmath=True)
def payoff(value, minimum, maximum, cp, knockin, knockout, strike):
    if knockin and knockin > spot and maximum < knockin: # Up and In
        return 0
    elif knockin and knockin < spot and minimum > knockin: # Down and In
        return 0
    elif knockout and knockout < spot and minimum < knockin: # Down and Out
        return 0
    elif knockout and knockout > spot and maximum > knockout: # Up and Out
        return 0
    return max(0, cp * (value - strike))

@njit(parallel=True, fastmath=True)
def price_option(strike, spot, time, volatility, risk_free, call_or_put='c', knockin=0, knockout=0, simulations=2000, steps_per_unit = 365):
    if knockin and knockout:
        raise Exception("Unable to cope with 2 barriers!")

    cp = 1 if call_or_put == 'c' else -1
    dt = 1 / steps_per_unit
    steps = int(time * steps_per_unit)
    sdt = sqrt(dt)
    drift = exp((risk_free - 0.5 * volatility * volatility) * dt)

    total_premium = 0
    for i in prange(simulations):
        value, minimum, maximum = path_final_min_max(spot, steps, sdt, volatility, drift)
        total_premium += payoff(value, minimum, maximum, cp, knockin, knockout, strike)

    return total_premium / simulations * exp(-time * risk_free)

spot=100
vol=0.2
risk_free=0.05
print(price_option(105, spot, 1, vol, risk_free, 'c', simulations=50000))