# -*- coding: utf-8 -*-
"""
Created on Thu Apr  8 14:13:24 2021

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


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=100 # nombre de scenarios
dt = T/N
W0=0

timerange = np.linspace(0,T,N+1,endpoint=True)
timerange=timerange[:,None]
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)

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

def st(mu,sigma,t,wt):
    """
    Calcule S_t de modele Black-Scholes en utilisant la formule exacte a partir
    d'un brownien.
    """

    return np.exp((mu-sigma**2/2.)*t + sigma*wt)

St = st(taux_r,sigma,timerange,W)*S0

plt.subplot(2,2,1)
plt.plot(timerange,W)
plt.title('Brownien')

plt.subplot(2,2,2)
plt.plot(timerange,St)
plt.title('St risque neutre')

plt.subplot(2,2,3)
plt.hist(St[-1,:],50)
plt.title('hist de St')

K=110
prixcall,_ = blsprice(S0,K,taux_r,T,sigma)
prixMC=np.exp(-taux_r*T)*np.mean(np.maximum(St[-1,:]-K,0))

print("prix Monte Carlo=",prixMC)

plt.subplot(2,2,4)
#version qui affiche toutes les valeurs, avec beaucoup de valeurs nulles
#plt.hist(np.exp(-taux_r*T) * 
#         np.maximum(St[-1,:]-K,0)-prixcall ,50)
#version qui affiche seulement les valeurs positives
plt.hist(np.exp(-taux_r*T) * 
         np.maximum(St[-1,:]-K,0)[np.maximum(St[-1,:]-K,0)>0]-prixcall ,50)
plt.title('hist du prix Monte Carlo')




erreur_MC =prixMC-prixcall
savefig("monte_carlo.jpg")
#exercice: modifier le code pour faire la meme chose pour le prix du put