barrier.numpy.py

barrier_numpy.py

from random import random
from math import sqrt, log, exp
from numba import njit
import numpy
from numpy.random import normal

def path_final_min_max(initial, steps, sdt, volatility, drift):
    randoms = numpy.exp(normal(size=steps)*volatility*sdt) * drift
    randoms[0] = 1
    factors = numpy.cumprod(randoms) * initial
    return factors[-1], numpy.min(factors), numpy.max(factors)

@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))

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 range(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))