# Black-Scholes formula for European call option

from math import *
def normcdf(x):
    #'Cumulative distribution function for the standard normal distribution'
    return (1.0 + erf(x / sqrt(2.0))) / 2.0

def normpdf(x):
    return 1/sqrt(2*pi) * exp(- x**2/2)

def callBS(S,K,T, sigma, r):
    d_plus = ((r + 1/2 * sigma**2) * T - log(K/S))/ (sigma * sqrt(T) )
    d_minus = ((r - 1/2 * sigma**2) * T - log(K/S))/ (sigma * sqrt(T) )
    return S * normcdf(d_plus) - K * exp(-r*T) * normcdf(d_minus)

# Set values for parameters

S = 100;
K = 110;
T = 1;
sigma = 0.1;
r = 0.02;

print("Black-Scholes price is " + str(callBS(S,K,T,sigma,r)))

Black-Scholes price is 1.3537383655181436

# Numerical integration of Black-Scholes price using rectangular rule

import numpy as np

def callBS_Rect(S,K,T, sigma, r):
    L = 100 # define upper limit of the integration
    d_minus = ((r - 1/2 * sigma**2) * T - log(K/S))/ (sigma * sqrt(T) )
    N = 10000 # number of subintervals
    delta = (L + d_minus) / N # length of a subinterval
    X = np.linspace(-d_minus,L,N) # the integrand is 0 for x less than - d_minus
    
    BS_Rect = 0
    for x in X:
            BS_Rect += exp(-r*T)/sqrt(2*pi) * (S*exp((r-0.5*sigma**2)*T + sigma * sqrt(T) * x) - K) * exp(-0.5 *x**2)
    return BS_Rect*delta

print("Black-Scholes price using rectangular rule is " + str(callBS_Rect(S,K,T,sigma,r)))

Black-Scholes price using rectangular rule is 1.3535774398429419

# Numerical integration of Black-Scholes price using trapezoidal rule

def callBS_Trap(S,K,T, sigma, r):
    L = 100 # define upper limit of the integration
    d_minus = ((r - 1/2 * sigma**2) * T - log(K/S))/ (sigma * sqrt(T) )
    N = 10000 # number of subintervals
    delta = (L + d_minus) / N # length of a subinterval
    X = np.linspace(-d_minus,L,N) # the integrand is 0 for x less than - d_minus
    
    BS_Trap = 0
    for i in range(0,N-1):
            BS_Trap +=  0.5 * ( exp(-r*T)/sqrt(2*pi) * (S*exp((r-0.5*sigma**2)*T + sigma * sqrt(T) * X[i]) - K) * exp(-0.5 *X[i]**2) 
                               + exp(-r*T)/sqrt(2*pi) * (S*exp((r-0.5*sigma**2)*T + sigma * sqrt(T) * X[i+1]) - K) * exp(-0.5 *X[i+1]**2) )
                               
    return BS_Trap*delta

print("Black-Scholes price using trapezoidal rule is " + str(callBS_Trap(S,K,T,sigma,r)))

Black-Scholes price using rectangular rule is 1.3535774398429439

# Numerical integration of Black-Scholes price using Simpson's rule

def callBS_Simp(S,K,T, sigma, r):
    L = 100 # define upper limit of the integration
    d_minus = ((r - 1/2 * sigma**2) * T - log(K/S))/ (sigma * sqrt(T) )
    N = 10000 # number of subintervals
    delta = (L + d_minus) / N # length of a subinterval
    X = np.linspace(-d_minus,L,N) # the integrand is 0 for x less than - d_minus
    
    BS_Simp = 0
    for i in range(2,N,2):
        BS_Simp +=  4 *  exp(-r*T)/sqrt(2*pi) * (S*exp((r-0.5*sigma**2)*T + sigma * sqrt(T) * X[i]) - K) * exp(-0.5 *X[i]**2)
    for i in range(1, N-1, 2):
        BS_Simp +=  2 *  exp(-r*T)/sqrt(2*pi) * (S*exp((r-0.5*sigma**2)*T + sigma * sqrt(T) * X[i]) - K) * exp(-0.5 *X[i]**2)
    BS_Simp += exp(-r*T)/sqrt(2*pi) * (S*exp((r-0.5*sigma**2)*T + sigma * sqrt(T) * X[0]) - K) * exp(-0.5 *X[0]**2)
    BS_Simp += exp(-r*T)/sqrt(2*pi) * (S*exp((r-0.5*sigma**2)*T + sigma * sqrt(T) * X[N-1]) - K) * exp(-0.5 *X[N-1]**2)                                                                                                           
    return BS_Simp*delta/3

print("Black-Scholes price using Simpson rule is " + str(callBS_Simp(S,K,T,sigma,r)))

Black-Scholes price using Simpson rule is 1.3535518877870085