In [24]:
# 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
In [22]:
# Compare values of different Greeks using numerical differentiation versus theoretical results

d = 0.0001;

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

delta_numerical = 1/d * (callBS(S+d, K,T,sigma,r) - callBS(S, K,T,sigma,r))
delta_theoretical = normcdf(d_plus)

print("Difference in deltas " + str(abs(delta_numerical - delta_theoretical)))

vega_numerical = 1/d * (callBS(S, K,T,sigma+d,r) - callBS(S, K,T,sigma,r))
vega_theoretical1 = K*exp(-r*T)*normpdf(d_minus)*sqrt(T)
vega_theoretical2 = S*normpdf(d_plus)*sqrt(T)

print("Difference in vegas " + str(abs(vega_numerical - vega_theoretical1)))

rho_numerical = 1/d * (callBS(S, K,T,sigma,r + d) - callBS(S, K,T,sigma,r))
rho_theoretical = K*T*exp(-r*T)*normcdf(d_minus)

print("Difference in rhos " + str(abs(rho_numerical - rho_theoretical)))

gamma_numerical = 1/d**2*(callBS(S+d, K,T,sigma,r) -  2* callBS(S, K,T,sigma,r) + callBS(S-d, K,T,sigma,r) )
gamma_theoretical = normpdf(d_plus)/(sigma*sqrt(T)*S)

print("Difference in gammas " + str(abs(gamma_numerical - gamma_theoretical)))

Difference in deltas 1.5579801374010493e-06
Difference in vegas 0.008789550730725182
Difference in rhos 0.014444634497685627
Difference in gammas 2.267266009520208e-06
In [26]:
# Finding implied volatility using Bisection method

a = 0.05;
b = 0.15;
tol = 0.0005;

cMarket = 1.35;

while (b-a > tol):
    if callBS(S,K,T,(a+b)/2, r) - cMarket > 0:
        b = (a+b)/2
    else :
        a = (a+b)/2

print("The implied volatiltiy at market price " + str(cMarket) + " is "  + str(a))

The implied volatiltiy at market price 1.35 is 0.099609375
In [29]:
# Finding implied volatility using Newton method

vol = 0.2
n = 0

while (abs(callBS(S,K,T,vol, r) - cMarket) > tol):
    n = n+1;
    vega = K*exp(-r*T)*normpdf(d_minus)*sqrt(T)
    vol = vol - (callBS(S,K,T,vol, r) - cMarket)/vega # Newton step

print("The implied volatiltiy at market price " + str(cMarket) + " is "  + str(vol))

The implied volatiltiy at market price 1.35 is 0.09987782168826839
In [ ]: