# -*- coding: utf-8 -*- """ Created on Sat Mar 28 14:50:39 2020 @author: Gabriel """ # -*- coding: utf-8 -*- """ GT 29/10/19 """ # %matplotlib auto import numpy as np import matplotlib.pyplot as plt from scipy.stats import norm def blsprice(Price,Strike,Rate,TimeToMaturity,Volatility,DividendRate=0): """ Computes price of option with analytic formula. 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 output: price of a call and of a put (tuple) """ 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 T=1.0 N=255 M=100 # nombre de scenarios dt = T/N W0=0 timerange = np.linspace(0,T,N+1,endpoint=True) timerange=timerange[:,None] dW= np.sqrt(dt)*np.random.randn(N,M) W=np.zeros((N+1,M)) W[0,:]=W0 W[1:,:]=W0+np.cumsum(dW,0) S0=100. mu=0.1 sigma=0.25 taux_r=0.05 def st(mu,sigma,t,wt): """ Calcule S_t de modele Black-Scholes en utilisant la formule exacte a partir d'un brownien. """ return np.exp((mu-sigma**2/2.)*t + sigma*wt) St = st(taux_r,sigma,timerange,W)*S0 plt.subplot(2,2,1) plt.plot(timerange,W) plt.title('Brownien') plt.subplot(2,2,2) plt.plot(timerange,St) plt.title('St risque neutre') plt.subplot(2,2,3) plt.hist(St[-1,:],50) plt.title('hist de St') K=110 prixcall,_ = blsprice(S0,K,taux_r,T,sigma) prixMC=np.exp(-taux_r*T)*np.mean(np.maximum(St[-1,:]-K,0)) print("prix Monte Carlo=",prixMC) plt.subplot(2,2,4) plt.hist(np.exp(-taux_r*T) * np.maximum(St[-1,:]-K,0)-prixcall ,50) plt.title('hist du prix Monte Carlo') erreur_MC =prixMC-prixcall savefig("monte_carlo.jpg")