# -*- coding: utf-8 -*-
"""
Copyright : Gabriel Turinici 2022

Do not distribute.
"""

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


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

def sousjacent(S0,W,dt,mu,sigma):
    """
    Input : mu,sigma, dt : real scalars
    dt : time step  
    W.shape = N+1,M N=time steps, M=no. of scenarios
    
    Returns : solution of dS/S = mu dt + sigma dWt

    """
    N=W.shape[0]-1
    t = dt * np.linspace(0,N*dt,N+1,endpoint=True)[:,None]
    St = S0*np.exp( (mu-sigma**2/2)*t +sigma*W)
    return St

T=1
N=250*8*4
dt=T/N
M=200 #no. realisations
W=np.zeros((N+1,M))
dW=np.sqrt(dt)*np.random.randn(N,M) #increments Browniens
W[1:,:]=np.cumsum(dW,axis=0)
#plt.plot(W)

S0=100
mu=0.15
sigma=0.2
sigmaestime=0.25
K=100
r=0.05

St=sousjacent(S0,W,dt,mu,sigma)
#plt.figure('S')    
#plt.plot(St)

#delta hedging
#mise en place du portefeuille de replication en t=0
cash=np.zeros_like(St)
#calcul du prix de l'option a vendre et ajout du prix dans le cash
cash[0,:] = blsprice(S0,K,r,T,sigmaestime)[0]
#calcul de Delta_0 et achat de Delta_0 parts de s/ jacent
titres=np.zeros_like(St)#nombre de titres en portefeuille
titres[0,:] = blsdelta(S0,K,r,T,sigmaestime)[0]
#mise a jour du cash
cash[0,:] -= titres[0,:]*S0

for k in range(N):# boucle sur le temps
    #mise à jour du cash, delta,cash
    cash[k+1,:]= cash[k,:]* np.exp(r*dt)
    titres[k+1,:] = blsdelta(St[k+1,:],K,r,T-k*dt,sigmaestime)[0]
    cash[k+1,:] -= St[k+1,:]*(titres[k+1,:]-titres[k,:])
    
#livraison à la fin en t=T: payer (ST-K)+, vendre le portef s/jacent
cash[-1,:] -= np.maximum(0,St[-1,:]-K)
cash[-1,:] += titres[-1,:]*St[-1,:]

#histogramme du cash final
plt.hist(cash[-1,:])

#cj = cash[-1,:]
#ch = cash[-1,:]
#cq = cash[-1,:]
#plt.hist([cj,ch,cq])