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

# MonteCarlo
Nsim = 100000
def payoff(S):
    return np.maximum(0, S-E)
num0 = np.random.normal(0,sigma,Nsim)
S1 = S0*np.exp(r-sigma**2/2+num0)
C0 = math.exp(-r)*np.mean(payoff(S1))
print(f'Call value computed by {Nsim} simulations:',C0)

# Black-Scholes
def d1(t,S):
    return (np.log(S/E) + (r+sigma**2/2)*t)/(sigma*np.sqrt(t))
def d2(t,S):
    return d1(t,S) - sigma*np.sqrt(t)
def call_value(t,S):
    return norm.cdf(d1(t,S))*S-norm.cdf(d2(t,S))*E*np.exp(-r*t)
 
t = 1
C0 = call_value(t,S0)
print('Call value computed by Black & Scholes:',C0)

space=np.linspace(0.001,2,100)
t = 0.001
plt.plot(space,call_value(t,space))
t = 0.5
plt.plot(space,call_value(t,space))
t = 1
plt.plot(space,call_value(t,space))
t = 1.5
plt.plot(space,call_value(t,space))
plt.show()