import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt 
import math

S0 = 1 #initial value of underlying asset
E = 1 #exercise price
sigma = 0.2 #volatility
mu = 0.06 #drift: it is ininfluent
r = 0.1 #interest rate

# Black-Scholes
def d1(S,t):
    return (np.log(S/E) + (r+sigma**2/2)*t)/(sigma*np.sqrt(t))
def d2(S,t):
    return d1(S,t) - sigma*np.sqrt(t)
def call_value(S,t):
    return S*norm.cdf(d1(S,t))-E*np.exp(-r*t)*norm.cdf(d2(S,t))
def put_value(S,t):
    return E*np.exp(-r*t)*norm.cdf(-d2(S,t))-S*norm.cdf(-d1(S,t))
def mix_value(S,t):
    return call_value(S,t)+put_value(S,t)

plt.rcParams['figure.figsize'] = [8, 12]
figure, axis = plt.subplots(3, 1)

space=np.linspace(0.001,2,100)
for t in [0.001,0.5,1,1.5]:
    axis[0].plot(space,call_value(space,t))
axis[0].set_title("European Call")
axis[0].set_xlabel("S",loc="right")
#axis[0].set_ylabel("C",rotation="horizontal",loc="top")
axis[0].legend(["t=0","t=0.5","t=1","t=1.5"])

for t in [0.001,0.5,1,1.5]:
    axis[1].plot(space,put_value(space,t))
axis[1].set_title("European Put")
#axis[1].set_xlabel("S",loc="right")
#axis[1].set_ylabel("P",rotation="horizontal",loc="top")

for t in [0.001,0.5,1,1.5]:
    axis[2].plot(space,mix_value(space,t))
axis[2].set_title("European Mix")
#axis[2].set_xlabel("S",loc="right")
#axis[2].set_ylabel("P",rotation="horizontal",loc="top")

plt.show()