In [1]:
import numpy as np
import matplotlib.pyplot as plt
Those are our option and market parameters: the exercise price of the option strike, the forward price of the underlying security forward and its volatility vol (as the model is lognormal, the volatility is a percentage number; eg 0.20 = 20%)
In [2]:
strike = 100
mat = 1
forward = 100
vol = 0.3
We now define our payoff function using a closure: the variable payoff represents a function with one parameter spot with the strike k being frozen at whatever value it had when the outer function call was called to set payoff
In [3]:
def call(k=100):
def payoff(spot):
if spot > k:
return spot - k
else:
return 0
return payoff
payoff = call(k=strike)
We also define an analytic function for calculation the price of a call using the Black Scholes formula, allowing us to benchmark our results
In [4]:
from scipy.stats import norm
def bscall(fwd=100,strike=100,sig=0.1,mat=1):
lnfs = log(fwd/strike)
sig2t = sig*sig*mat
sigsqrt = sig*sqrt(mat)
d1 = (lnfs + 0.5 * sig2t)/ sigsqrt
d2 = (lnfs - 0.5 * sig2t)/ sigsqrt
fv = fwd * norm.cdf (d1) - strike * norm.cdf (d2)
#print "d1 = %f (N = %f)" % (d1, norm.cdf (d1))
#print "d2 = %f (N = %f)" % (d2, norm.cdf (d2))
return fv
#bscall(fwd=100, strike=100, sig=0.1, mat=1)
We now generate a set of Standard Gaussian variables $z$ as a basis for our simulation...
In [5]:
N = 10000
z = np.random.standard_normal((N))
#z
...and transform it in a lognormal variable with the right mean and log standard deviation, ie a variable that is distributed according to $LN(F,\sigma\sqrt{T})$. Specifically, to transform a Standard Gaussian $Z$ into a lognormal $X$ with the above parameters we use the following formula $$ X = F \times \exp ( -0.5 \sigma^2 T + \sigma \sqrt{T} Z ) $$
In [6]:
x = forward * exp(- 0.5 * vol * vol * mat + vol * sqrt(mat) * z)
min(x), max(x), mean(x)
Out[6]:
We first look at the histogram of the spot prices $x$ (the function trim_vals simply deals with with the fact that histogram returns the starting and the ending point of the bin, ie overall one point too many)
In [7]:
def trim_xvals(a):
a1 = np.zeros(len(a)-1)
for idx in range(0,len(a)-1):
#a1[idx] = 0.5*(a[idx]+a[idx+1])
a1[idx] = a[idx]
return a1
hg0=np.histogram(x, bins=50)
xvals0 = trim_xvals(hg0[1])
In [8]:
fwd1 = mean(x)
print ("forward = %f" % (fwd1))
plt.bar(xvals0,hg0[0], width=0.5*(xvals0[1]-xvals0[0]))
plt.title('forward distribution')
plt.xlabel('forward')
plt.ylabel('occurrences')
plt.show()
We now determine the payoff values from our draws of the final spot price. Note that we need to use the map command rather than simply writing po = payoff(x). The reason for this is that this latter form is not compatible with the if statement in our payoff function. We also already compute the forward value of the option, which is simply the average payoff over all simulations.
In [9]:
po = list(map(payoff,x))
fv = mean(po)
#po
Now we produce the histogram of the payoffs
In [10]:
hg = np.histogram(po,bins=50)
xvals = trim_xvals(hg[1])
In [11]:
plt.bar(xvals,hg[0], width=0.9*(xvals[1]-xvals[0]))
plt.title('payout distribution')
plt.xlabel('payout')
plt.ylabel('occurrences')
plt.show()
In the next step we compute our "Greeks", ie a number of derivatives of the forward value with respect to the underlying parameters. What is crucial here is that those derivative are calculated on the same draw random numbers $z$, otherwise the Monte Carlo sampling error will dwarf the signal. The sensitivities we compute are to increase / decrease the forward by one currency unit (for Delta and Gamma), to increase the volatility by one currency unit (for Vega), and to decrease the time to maturity by 0.1y (for Theta)
In [12]:
x = (forward+1) * exp(- 0.5 * vol * vol * mat + vol * sqrt(mat) * z)
po = list(map(payoff,x))
fv_plus = mean(po)
x = (forward-1) * exp(- 0.5 * vol * vol * mat + vol * sqrt(mat) * z)
po = list(map(payoff,x))
fv_minus = mean(po)
x = forward * exp(- 0.5 * (vol+0.01) * (vol+0.01) * mat + (vol+0.01) * sqrt(mat) * z)
po = list(map(payoff,x))
fv_volp = mean(po)
x = forward * exp(- 0.5 * vol * vol * (mat-0.1) + vol * sqrt(mat-0.1) * z)
po = list(map(payoff,x))
fv_timep = mean(po)
In [13]:
print ("Strike = %f" % strike)
print ("Maturity = %f" % mat)
print ("Forward = %f" % forward)
print ("Volatility = %f" % vol)
print ("FV = %f" % fv)
print (" check = %f" % bscall(fwd=forward, strike=strike, sig=vol, mat=mat))
print ("Delta = %f" % ((fv_plus - fv_minus)/2))
print ("Gamma = %f" % ((fv_plus + fv_minus - 2 * fv)))
print ("Theta = %f" % ((fv_timep - fv)))
print ("Vega = %f" % ((fv_volp - fv)))
(c) Stefan Loesch / oditorium 2014; all rights reserved (license)
In [14]:
import sys
print(sys.version)