# -*- coding: utf-8 -*-
"""
Created on Mon Apr 12 14:51:49 2021

@author: turinici
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm


#implementation of the Black-Scholes formula
def blsprice(Price,Strike,Rate,TimeToMaturity,Volatility,DividendRate=0):
    """
    Computes price of option with analytic formula.
    input:
       S:Price - Current price of the underlying asset.
    
       Strike:Strike - Strike (i.e., exercise) price of the option.
    
       Rate: Rate - Annualized continuously compounded risk-free rate of return over
         the life of the option, expressed as a positive decimal number.
    
       TimeToMaturity:Time - Time to expiration of the option, expressed in years.
    
      Volatility: volatility
      DividendRate = continuous dividend rate

    output: price of a call and of a put (tuple)
    """
    
    if TimeToMaturity <= 1e-6: # the option already expired
        call = np.max(Price-Strike,0)
        put = np.max(Strike-Price,0)
        return call,put
        
    
    d1 = np.log(Price/Strike)+(Rate-DividendRate + Volatility**2/2.0)*TimeToMaturity;
    d1 = d1/(Volatility* np.sqrt(TimeToMaturity))
    d2 = d1-(Volatility*np.sqrt(TimeToMaturity))
    
    call = Price * np.exp(-DividendRate*TimeToMaturity) * \
        norm.cdf(d1)-Strike* np.exp(-Rate*TimeToMaturity) * norm.cdf(d2)
    put  = Strike* np.exp(-Rate*TimeToMaturity) * norm.cdf(-d2)\
        -Price* np.exp(-DividendRate*TimeToMaturity) * norm.cdf(-d1)
    return call,put



    
T=1.0
N=255
M=300#number of scenarios
dt= T/N
W0=0#standard brownian motion
trange=np.linspace(0,T,N+1,endpoint=True)

dW=np.sqrt(dt)*np.random.randn(N,M)
W=np.zeros((N+1,M)) 
W[0,:]=W0
W[1:,:]=W0+np.cumsum(dW,0)

plt.figure(2)
plt.plot(trange,W)


S0=100.
mu=0.1
sigma=0.25
taux_r=0.05

#compute S_t with dS_t = mu S_t dt + sigma S_t d W_t
St=np.zeros_like(W)
St[0,:]=S0

for ii in range(N):
    St[ii+1,:] =St[ii,:]+ mu*St[ii,:]*dt + sigma*St[ii,:]*dW[ii,:]
    
plt.figure(3)
plt.plot(trange,St)
plt.title('$S_t$')






#compute the Monte Carlo price of an option
# solve St in risk-neutral probability, denote rn_St

#compute rn_St with 
# d rn_St = mu rn_St dt + sigma rn_St d W_t
rn_St=np.zeros_like(W)
rn_St[0,:]=S0

for ii in range(N):
    rn_St[ii+1,:] =rn_St[ii,:]+ taux_r*rn_St[ii,:]*dt \
        + sigma*rn_St[ii,:]*dW[ii,:]
    
#compute the price of the call
K=110
prixcall,_ = blsprice(S0,K,taux_r,T,sigma)
prixMC=np.exp(-taux_r*T)*np.mean(np.maximum(rn_St[-1,:]-K,0))


print("prix Monte Carlo=",prixMC)

plt.subplot(2,2,4)
plt.hist(np.exp(-taux_r*T) * 
         np.maximum(rn_St[-1,:]-K,0)-prixcall ,50)
plt.title('hist du prix Monte Carlo')

erreur_MC =prixMC-prixcall
#plt.savefig("euler_maruyama_monte_carlo.jpg")