# -*- 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 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 # 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.figure('brownian',figsize=(6,3)) plt.subplot(1,2,1) plt.plot(W) plt.subplot(1,2,2) 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,:]-= delta[0,:]*S[0,:] #achat de Delta(0) parts #de sous jacent: mettre en memoire Delta(0) et mettre à jour le cash f=0.00#frais de transaction de 1/1000 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,:] = cash[jj+1,:]- S[jj+1,:]*(delta[jj+1,:]-delta[jj,:]) \ - f*np.abs(S[jj+1,:]*(delta[jj+1,:]-delta[jj,:])) #achat ou vente de sous-jacent et mise à jour cash #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.figure("hist") plt.hist(cash[-1,:]) #plt.hist(cash_old[-1,:])