# -*- coding: utf-8 -*-
"""trading-de-volatilite-M2ISFP21.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1WhXjKfB9XZnXd5IQldgi_qOij6hfyEVX
"""

# A faire: simulation avec volatilite realisee differente de la volatilite supposee (prevue)

"""Cours M2 ISF  : Gestion de risques et construction de portefeuille - séance valuation produits dérivés : Monte Carlo et delta hedging"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

T=1 # temps final
N=255*8*4
N=255
dt=T/N
timerange = np.linspace(0,T,N+1,endpoint=True)
timerange=timerange[:,None]
print(timerange.shape)

M=150 #no de realisations
#construction de Wt:
Wt=np.zeros((N+1,M))
W0=0
dW=np.sqrt(dt)*np.random.randn(N,M)
Wt[0,:]=W0
Wt[1:,:]=W0+np.cumsum(dW,0)

plt.figure()
plt.subplot(1,2,1)
plt.plot(Wt)
plt.title('Realisations du brownien')

S0=100.0
mu=0.1
sigma=0.15
sigma_i=0.15 # volatilite prevue, supposee
sigma_r=0.05 # volatilite reelle

taux_r=0.05

def sousjacent(mu,sigma,t,Wt):
  """ Rend les realisations de St une fois donne le tmps et Wt
  avec la formule St = S0 * exp( (mu-sigma*sigma/2.0)*t + sigma*Wt)  """ 
  return S0*np.exp( (mu-sigma**2.0/2.0)*t+ sigma*Wt)

St= sousjacent(mu,sigma_r,timerange,Wt)#volatilite reelle, car futur
plt.subplot(1,2,2)
plt.plot(St)
plt.title('Realisations du sous-jacent')
plt.savefig("realisations_Wt_St.jpg")

import numpy as np
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

K=110 # strike
#delta hedging
# on calcule / simule la situation d'un vendeur de call qui se couvre de maniere dynamique en achetant deltat parts de sous-jacent
# deltat = + d C_t / d S
# t=0 : vente option, cash = cash + prix option (celui de blsprice), mise en place du portefeuille de replication
cash = np.zeros_like(St)
prix_call,_=blsprice(S0,K,taux_r,T,sigma_i)
deltaancien,_=blsdelta(S0,K,taux_r,T,sigma_i) #calcul de delta
cash[0,:]=prix_call - deltaancien*S0

for jj in range(N):
  cash[jj+1,:]= cash[jj,:]*np.exp(taux_r*dt)# capitalisation
  #calcul du nouveau delta
  deltanouveau,_ = blsdelta(St[jj+1,:],K,taux_r,T-timerange[jj+1,0],sigma_i)
  cash[jj+1,:]= cash[jj+1,:]- (deltanouveau-deltaancien)*St[jj+1,:] #prix de la prise en compte du delta hedging
# t -> t + dt : ajustement de porfef. de replication, capitalisation cash
  deltaancien=deltanouveau

cash[-1,:] -=np.maximum(St[-1,:]-K,0)#perte due a la livraison 
cash[-1,:] += deltanouveau*St[-1,:]#cash obtenu par la vente du portef de couverture 

# histogramme du cash
plt.figure
plt.hist(cash[-1,:]*100/prix_call)