# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""

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

def blsprice(Price,Strike,Rate,TimeToMaturity,Volatility,DividendRate=0):
    """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
    """
    
    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

#test: blsprice(100,110,0.05,1,0.2)
#blsprice(np.array([100,101]),np.array([110,111]),0.05,1,0.2)

def blsdelta(Price,Strike,Rate,TimeToMaturity,Volatility,DividendRate=0):
    """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
    """
    
    if TimeToMaturity <= 1e-6: # the option already expired
        call = (Price>=Strike).astype(np.float) # 1 if in the money, zero otherwise
        put = (Price<=Strike).astype(np.float) # cf above
        return call,put
        
    
    d1 = np.log(Price/Strike)+(Rate-DividendRate + Volatility**2/2.0)*TimeToMaturity;
    d1 = d1/(Volatility* np.sqrt(TimeToMaturity))
    
    call =  np.exp(-DividendRate*TimeToMaturity) * norm.cdf(d1)
    put  = -np.exp(-DividendRate*TimeToMaturity) * norm.cdf(-d1)
    return call,put
T=1
N=255*8
dt=T/N
M=1000
print("M=",M)

S0=100
K=110
mu=0.15
sigma=0.2
r=0.05

#construction mouvement brownien par somme cummulative 
dW= np.sqrt(dt)*np.random.randn(N,M)
W= np.zeros((N+1,M))
W[0,:]=0
W[1:,:]=np.cumsum(dW,axis=0)

t=np.linspace(0,T,N+1,endpoint=True)
#plot
plt.figure('brownian')
plt.plot(W)

# construction du processus de prix S: dS/S = mu*dt + sigma* dWt
S= S0*np.exp( (mu-sigma**2/2)*t[:,None] + sigma*W)
plt.plot(S)

#DELTA HEDGING : t=0 : mise en place du portefeuille de replication
cash=np.zeros_like(S)# matrice du meme format que S
#Calcul prix initial et vente de call
cash[0,:]=blsprice(S0,K,r,T,sigma)[0] #remplacer par le prix initial
delta=np.zeros((N+1,M))
delta[0,:]=blsdelta(S0,K,r,T,sigma)[0]#calcul de Delta(t=0)
cash[0,:]=cash[0,:]-delta[0,:]*S[0,:]#achat de Delta(0) parts 
#de sous jacent: mettre en memoire Delta(0) et mettre à jour le cash
f=0.001#frais de transaction

for jj in range(N):#boucle en temps: 
    cash[jj+1,:]=cash[jj,:]*np.exp(r*dt)#capitalisation du cash
    delta[jj+1,:]=blsdelta(S[jj+1,:],K,r,T-jj*dt,sigma)[0]#recalcul de Delta nouveau
    cash[jj+1,:] -= (delta[jj+1,:]-delta[jj,:])*S[jj+1,:]#achat ou vente de sous-jacent et mise à jour cash
    cash[jj+1,:] -= f*np.abs((delta[jj+1,:]-delta[jj,:])*S[jj+1,:])

#t=T: paiement de (S(T)-K)+ du cash; vente des ss-jacents
cash[-1,:] -= np.maximum(S[-1,:]-K,0)
cash[-1,:] += delta[-1,:]*S[-1,:]

# histograme de la valeur finale du portefeuille c'est à dire
# du cash
#cash_old=cash
plt.hist(cash[-1,:])
#plt.hist(cash_old[-1,:])