import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from scipy.stats import norm
import matplotlib.dates as mdates
import matplotlib.cbook as cbook
import pandas as pd
import datetime as dt
def N(d):
# return the Normal function up to value d
return norm().cdf(d)
def F(t,K,T,r):
# return the function F from the theory
return K * np.exp(r*(t-T))
S = 51.50 #stockprice
K = 50 #maturity price
T = 72 #maturity time, in days
t = np.linspace(0,T-0.00000000001,T+1)
r = 0.05 / 365 #risk-free interest rate
sigma = np.array([0.2010,0.1600,0.12,0.24,0.28,0.4840]) / 16 #volatility as of 31 march 2021 (CBOE index)
C = []
# calculate the d values for in the normal distribution
print(t)
for index in range(len(sigma)):
d1 = (np.log(S / K) + (T - t) * (r + sigma[index] ** 2 / 2)) / (sigma[index] * (T - t) ** (1/2))
d2 = (np.log(S / K) + (T - t) * (r - sigma[index] ** 2 / 2)) / (sigma[index] * (T - t) ** (1/2))
# calculate values for C(S,t)
C.append( S * N(d1) - F(t,K,T,r) * N(d2) )
C = np.array(C)
plt.figure(figsize=(12,9), dpi=100)
# plot the values for C and t
for i in range(len(sigma)):
plt.plot(t,C[i], label='$\sigma$ = {}'.format(sigma[i]))
plt.xlabel('Time ($\it{days}$)')
plt.ylabel('C ($\it{€}$)')
plt.grid()
plt.legend()
plt.xlim(0,75)
plt.ylim(0,8)
plt.title('Price of a European call option')
plt.savefig('Black-Scholes Theoretical application')
plt.show()
print(C[0][0])
def N(d):
# return the Normal function up to value d
return norm().cdf(d)
def F(t,K,T,r):
# return the function F from the theory
return K * np.exp(r*(t-T))
#stockprice of Total, each day, starting at march 1st, closing price. Via Euronext
S =np.array([38.835, 38.915, 39.310, 40.560, 40.975, 40.975, 40.975,
40.920, 40.665, 41.310, 41.740, 42.190, 42.190, 42.190,
41.290, 40.765, 40.610, 40.435, 40.230, 40.230, 40.230,
39.730, 39.440, 40.065, 38.745, 39.160, 39.160, 39.160,
39.745, 40.250, 39.775,
39.100, 39.100, 39.100, 39.100,
39.100, 38.820, 38.915, 38.260, 37.845, 37.845, 37.845,
37.670, 37.875, 38.295, 38.000, 37.950])
K = 40 #strike price
T = 47 #maturity time, in days
t = np.linspace(0,T-0.00000001,T+1)
r = np.array([-0.611094, -0.605278, -0.605430, -0.617836, -0.613456, -0.613456, -0.613456,
-0.608869, -0.611540, -0.606683, -0.604049, -0.603366, -0.603366, -0.603366,
-0.598771, -0.603709, -0.604606, -0.607932, -0.613101, -0.613101, -0.613101,
-0.614519, -0.610431, -0.611952, -0.622158, -0.625250, -0.625250, -0.625250,
-0.627750, -0.617208, -0.636103,
-0.628903, -0.628903, -0.628903, -0.628903,
-0.628903, -0.626207, -0.631266, -0.638866, -0.631205, -0.631205, -0.631205,
-0.632801, -0.630718, -0.636989, -0.642934, -0.638678]) / 36500 #risk-free interest rate as per ECB yield curve 3 month.
sigma = np.array([23.350000, 24.100000, 26.670000, 28.570000, 24.660000, 24.660000, 24.660000,
25.469999, 24.030001, 22.559999, 21.910000, 20.690001, 20.690001, 20.690001,
20.030001, 19.790001, 19.230000, 21.580000, 20.950001, 20.950001, 20.950001,
18.879999, 20.299999, 21.200001, 19.809999, 18.860001, 18.860001, 18.860001,
20.740000, 19.610001, 19.400000,
17.330000, 17.330000, 17.330000, 17.330000,
17.910000, 18.120000, 17.160000, 16.950000, 16.690000, 16.690000, 16.690000,
16.910000, 16.650000, 16.990000, 16.570000, 16.250000]) / 1600 #volatility each day from march 1st till april 1st
C_check = [[0.92, 1.11, 0.83, 0.47, 0.26, 0.23, 0.10, 0.03, 0.01, 0.01, 0.01],
[ 28, 29, 30, 31, 36, 37, 38, 39, 42, 43, 44]]
# calculate the d values for in the normal distribution
d1 = (np.log(S / K) + (T - t[:-1]) * (r + sigma ** 2 / 2)) / (sigma * (T - t[:-1]) ** (1/2))
d2 = (np.log(S / K) + (T - t[:-1]) * (r - sigma ** 2 / 2)) / (sigma * (T - t[:-1]) ** (1/2))
# calculate values for C(S,t)
C = ( S * N(d1) - F(t[:-1],K,T,r) * N(d2) )
plt.figure(figsize=(12,9), dpi=100)
fig, ax1 = plt.subplots(figsize=(12,9))
color = 'tab:red'
ax1.set_xlabel('Time from 1-Mar-2021 $(days)$')
ax1.set_ylabel('Option Price $(€)$', color=color)
C_plot = ax1.plot(C, color=color, label='Theoretical TO2 option price')
C_real_plot = ax1.scatter(C_check[1],C_check[0], color=color, label='Real TO2 option price')
ax1.tick_params(axis='y', labelcolor=color)
ax1.legend(loc='upper left')
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.set_ylabel('Stockprice $(€)$', color=color) # we already handled the x-label with ax1
S_plot = ax2.plot(S, color=color, label ='TOTAL stock price')
ax2.tick_params(axis='y', labelcolor=color)
ax1.grid()
ax1.set_xlim([0,T+1])
ax1.set_ylim([-0.5,7.5])
ax2.set_ylim([35.5,43.5])
fig.tight_layout() # otherwise the right y-label is slightly clipped
ax2.legend(loc='upper right')
plt.title('Price of a European call option and its underlying stock')
plt.savefig('Black-Scholes Real Life application v2', bbox_inches='tight')
plt.show()
overview = []
overview.append(S)
overview.append(r)
overview.append(sigma)
df1 = pd.DataFrame(overview,
index=['Stockprice (€)', 'Risk-free interest rate', 'Volatility'])
df1 = df1.T
df1.to_excel("Black-Scholes real life values.xlsx")
def N(d):
# return the Normal function up to value d
return norm().cdf(d)
def F(t,K,T,r):
# return the function F from the theory
return K * np.exp(r*(t-T))
S = 51.50 #stockprice
K = 50 #maturity price
T = 72 #maturity time, in days
t = np.linspace(0,T-0.00000000001,T+1)
r = 0.05 / 365 #risk-free interest rate
sigma = np.array([0.2010,0.1600,0.12,0.24,0.28,0.4840]) / 16 #volatility as of 31 march 2021 (CBOE index)
C = []
# calculate the d values for in the normal distribution
print(t)
for index in range(len(sigma)):
d1 = (np.log(S / K) - (T - t) * (r + sigma[index] ** 2 / 2)) / (sigma[index] * (T - t) ** (1/2))
d2 = (np.log(S / K) + (T - t) * (r - sigma[index] ** 2 / 2)) / (sigma[index] * (T - t) ** (1/2))
# calculate values for C(S,t)
C.append( K * (S/K)**((-2*r)/(sigma[index]**2)) * N(d1) - F(t,K,T,r) * N(d2) )
C = np.array(C)
plt.figure(figsize=(12,9), dpi=100)
# plot the values for C and t
for i in range(len(sigma)):
plt.plot(t,C[i], label='$\sigma$ = {}'.format(sigma[i]))
plt.xlabel('Time ($\it{days}$)')
plt.ylabel('C ($\it{€}$)')
plt.grid()
plt.legend()
plt.xlim(0,75)
#plt.ylim(0,8)
plt.title('Price of a European call option')
plt.savefig('Black-Scholes Theoretical application')
plt.show()
print(C[0][0])