# -*- coding: utf-8 -*-
"""
Created on Wed Oct 28 14:17:38 2020

@author: Gabriel

Cours M2 ISF  P21: Gestion de risques et construction de portefeuille - 
seance stop-loss / CPPI

A faire: simulation Stop-Loss.
"""

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,3,1)
plt.plot(Wt)
plt.title('Realisations du brownien')

S0=100.0
mu=0.1
sigma=0.15

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,timerange,Wt)
plt.subplot(1,3,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 delta_stop_loss(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
    """
    
    
    call =  np.heaviside(Price-Strike,0)
    return call,1.0-call

blsdelta=delta_stop_loss 


K=110 # niveau stop-loss
cash = np.zeros_like(St)
deltaancien,_=blsdelta(S0,K,taux_r,T,sigma) #calcul de delta
cash[0,:]=S0 - deltaancien*S0# cash initial est 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)
  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,:] += deltanouveau*St[-1,:]#cash obtenu par la vente du portef de couverture 

# histogramme du cash
plt.figure
plt.subplot(1,3,3)
plt.hist(cash[-1,:] - np.maximum(St[-1,:],K))

# gain espere est max(ST,K)-S0= max(ST-K,0)+K-S0 = (ST-K)_+ - (S0-K)

# si croisement St <K on perd (K-St)