In [1]:
import numpy as np
import matplotlib.pyplot as plt
In [2]:
import time
def new_timer(id='', disabled=False):
beg_ts = time.time ()
def finish(id2='', reset=False):
nonlocal beg_ts
if disabled: return
end_ts = time.time ()
timediff = end_ts - beg_ts
print("%s%s elapsed time: %f" % (id, id2, timediff))
if reset: beg_ts = time.time() # needs `nonlocal` defn from Py3
return timediff
return finish
We first define a number of different European option payoff structures (ie payoff structures that only depende on the value of the spot at maturity), notable
In [3]:
def payoff_call(k=100):
def payoff(spot):
if spot > k:
return spot - k
else:
return 0
return payoff
def payoff_put(k=100):
def payoff(spot):
if spot < k:
return k - spot
else:
return 0
return payoff
def payoff_fwd(k=100):
def payoff(spot):
return spot - k
return payoff
def payoff_straddle(k=100):
def payoff(spot):
if spot < k:
return k - spot
else:
return spot - k
return payoff
We define a factory function that then provides us with a calculation function for our specific model. The factory function most importantly takes the number of Monte Carlo samples to be drawn N. For convenience it also takes a default payoff function as often the payoff will not be varied over the course of one analysis.
Importantly, the vector of Standard Gaussian random variables z that will be the basis of the Monte Carlo simulation is drawn immediately, and it is remembered in the closure This allows us to run the calculation function (the one returned by the factory) with different parameters but against the same sampling set, thereby greatly reducing Monte Carlo errors.
In [4]:
def calc_factory (N=1000, default_payoff=None):
z = np.random.standard_normal((N))
def run (fwd, sig, mat, payoff=None):
read_timer = new_timer(disabled=True)
x = fwd * exp(- 0.5 * sig * sig * mat + sig * sqrt(mat) * z)
read_timer("a")
if payoff == None:
payoff = default_payoff
po = list(map(payoff,x))
read_timer("b")
fv = mean(po)
read_timer("c")
return fv
return run
In [5]:
calc = calc_factory(N=1000)
calc(fwd=100, sig=0.2, mat=1, payoff=payoff_straddle(100))
Out[5]:
We also implement the Black Schole formula for testing the code
In [6]:
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)
return fv
In [7]:
strike = 100
mat = 1
forward = 100
vol = 0.1
call = payoff_call (k=strike)
put = payoff_put (k=strike)
fwd = payoff_fwd (k=strike)
straddle = payoff_straddle (k=strike)
In [8]:
callfv = calc(fwd=forward, sig=vol, mat=mat, payoff=call)
callfvc = bscall(fwd=100,strike=100,sig=0.1,mat=1)
putfv = calc(fwd=forward, sig=vol, mat=mat, payoff=put)
fwdfv = calc(fwd=forward, sig=vol, mat=mat, payoff=fwd)
stradfv = calc(fwd=forward, sig=vol, mat=mat, payoff=straddle)
#print ("Call = %f", callfv)
In [9]:
"""
print "Call = %f" % callfv
print " check = %f" % callfvc
print "Put = %f" % putfv
print "Forward = %f" % fwdfv
print " check = %f" % (callfv-putfv)
print "Straddle = %f" % stradfv
print " check = %f" % (callfv+putfv)
"""
Out[9]:
We have the calc function and we can compute it on a range of values. For this we need to freeze all but one parameters what we do with an ad-hoc closure: the function freeze freezes all parameters of its first argument (which must be a function with the signature of calc) except the paramter fwd. It returns a calculation function that only takes one parameter. This function is then applies to forwards.
TODO Use functools.partial instead of freeze see here
In [10]:
forwards = np.linspace(75, 125, 25)
def freeze(thecalc, sig=0.20, mat=1.0, payoff=call):
def calc1(fwd):
return thecalc(fwd=fwd, sig=sig, mat=mat, payoff=payoff)
return calc1
#fvs = calc(fwd=forwards, sig=0.20, mat=1.0, payoff=call)
payoff1 = straddle
In [11]:
plt.plot(forwards, list(map (freeze(calc, sig=0.20, mat=2.0, payoff=payoff1), forwards)), label="2y")
plt.plot(forwards, list(map (freeze(calc, sig=0.20, mat=1.0, payoff=payoff1), forwards)), label="1y")
plt.plot(forwards, list(map (freeze(calc, sig=0.20, mat=0.5, payoff=payoff1), forwards)), label="0.5y")
plt.plot(forwards, list(map (freeze(calc, sig=0.20, mat=0.25, payoff=payoff1), forwards)), label="0.25y")
plt.plot(forwards, list(map (freeze(calc, sig=0.20, mat=0.1, payoff=payoff1), forwards)), label="0.1y")
plt.title('Premium of a European Option')
plt.ylabel('forward premium')
plt.xlabel('forward')
plt.legend()
plt.show()
(c) Stefan Loesch / oditorium 2014; all rights reserved (license)
In [12]:
import sys
print(sys.version)