# -*- coding: utf-8 -*-
"""
Created on Fri Jan 27 17:08:49 2023

@author: turinici
"""

# -*- 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 bachelier_price_call(sigma, S, K, r, t):
    ''' price of a call under the Bachelier model   
    imported from
    https://github.com/kpmooney/numerical_methods_youtube/blob/master/bachelier_model/Bachelier%20Model%201.ipynb
    '''

    if t <= 1e-6: # the option already expired
        return  np.maximum(S-K,0)

    d = (S * np.exp(r*t) - K) / np.sqrt(sigma**2/(2 * r) * (np.exp(2*r*t)-1) )   

    C = np.exp(-r * t) * (S * np.exp(r * t) - K) * norm.cdf(d) + \
        np.exp(-r * t) * np.sqrt(sigma**2/(2*r) * (np.exp(2*r*t)-1) ) * norm.pdf(d)

    return C

def bachelier_delta_call(sigma, S, K, r, t):
    deltax = 1e-7
    return (bachelier_price_call(sigma,S+deltax,K,r,t)-bachelier_price_call(sigma, S-deltax,K,r,t))/(2.*deltax)


T=1.0
N=255
M=400 # 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_Black_Scholes=0.25
sigma=sigma_Black_Scholes*S0
print('setting sigma Bachelier coherent with Black-Scholes')

taux_r=0.05

def st_bachelier(mu,sigma,t,wt,S0=1):
    """
    Calcule S_t de modele Bachelier en utilisant la formule exacte a partir
    d'un brownien.
    
    t = timerange, type (N+1,)
    wt np array shape (N+1,M)
    """
    st=np.zeros_like(wt)
    st[0,:]=S0
    for ii in range(wt.shape[0]-1):
        st[ii+1,:]= st[ii,:]*(1+mu*dt)+sigma*(wt[ii+1,:]-wt[ii,:]) 

    return st

St = st_bachelier(mu,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')

K=110

#delta hedging
cash=np.zeros_like(St)
delta=np.zeros_like(St)

#t=0 : mise en place du portef de hedging
#vente call
#cash[0,:] = blsprice(S0,K,taux_r,T,sigma)[0]
#delta[0,:] = blsdelta(S0,K,taux_r,T,sigma)[0]
cash[0,:] = bachelier_price_call(sigma,S0,K,taux_r,T)
delta[0,:] = bachelier_delta_call(sigma,S0,K,taux_r,T)

cash[0,:] -= S0*delta[0,:]

#t entre 0 et T 
for ii in range(N):
    cash[ii+1,:]=np.exp(taux_r*dt)*cash[ii,:]#capitalisation du cash
    delta[ii+1,:]= bachelier_delta_call(sigma,St[ii+1,:],K,taux_r,T-(ii+1)*dt)#calcul nouveau delta    
    cash[ii+1,:] = cash[ii+1,:]- (delta[ii+1,:]-delta[ii,:])*St[ii+1,:]
    
#livraison en T et plot de hist finale
#paiement de (St-K)+
cash[-1,:] -= np.maximum(St[-1,:]-K,0)
cash[-1,:] += delta[-1,:]*St[-1,:]

plt.subplot(2,2,3)
plt.title('hedging_cash')
plt.hist(cash[-1,:])