Dr. Yves J. Hilpisch
The Python Quants GmbH
NYCPython Meetup "For Python Quants" – 16 January 2014
This talk focuses on
A brief bio:
Visit www.hilpisch.com.
You can find these Slides here:
Important Milestones (I), from a rather personal perspective.
Important Milestones (II), from a rather personal perspective.
Black Scholes Merton (1973) consider a financial model with a geometric Brownian motion representing the value of the only traded, risky asset:
\(dS_t = r S_t dt + \sigma S_t dZ_t\)
with \(S\) the stock price, \(r\) the risk-free rate, \(\sigma\) the volatility and \(Z\) a Brownian motion.
The other traded asset is a risk-less bond yielding the risk-free interest rate:
\(dB_t = r B_t\)
In this setting, Black Scholes Merton (1973) derive a closed-form solution for the price of a European call option as follows:
\[ \begin{eqnarray*} C(S_t, K, T, r, \sigma^{imp}) &=& S_{t} \cdot \mathbf{N}(d_{1})-e^{-r(T-t)} \cdot K \cdot \mathbf{N}(d_{2}) \\ \mathbf{N}(d)&=&\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{d}e^{-\frac{1}{2}x^{2}}dx \\ d_{1}&=&\frac{\log \frac{S_{t}}{K}+(r+\frac{\sigma^{2}}{2})(T-t)}{\sigma \sqrt{T-t}} \\ d_{2}&=&\frac{\log \frac{S_{t}}{K}+(r-\frac{\sigma^{2}}{2})(T-t)}{\sigma \sqrt{T-t}} \end{eqnarray*} \]
Here, \(\mathbf{N}\) stands for the cdf of the normal distribution.
class call_option_analytical(object):
''' Class for European call option valuation in BSM model.
Analytical fomula.'''
from math import log, sqrt, exp
from scipy import stats
global log, sqrt, exp, stats
def __init__(self, S0, K, T, r, sigma):
self.S0 = float(S0)
self.K = K
self.T = T
self.r = r
self.sigma = sigma
def present_value(self):
d1 = ((log(self.S0 / self.K) + (self.r + 0.5 * self.sigma ** 2) * T)
/ (self.sigma * sqrt(T)))
d2 = ((log(self.S0 / self.K) + (self.r - 0.5 * self.sigma ** 2) * T)
/ (self.sigma * sqrt(T)))
value = (self.S0 * stats.norm.cdf(d1, 0.0, 1.0)
- self.K * exp(-self.r * self.T) * stats.norm.cdf(d2, 0.0, 1.0))
return value
Let us generate a benchmark example for a European call option valuation.
S0 = 100
K = 100
T = 1.0
r = 0.05
sigma = 0.2
call = call_option_analytical(S0, K, T, r, sigma)
call.present_value()
10.450583572185565
One of the undeniable advantages of the analytical pricing formula is the speed when it comes to the numerical evaluation.
%%time
values = []
for K in np.linspace(80, 120, 5000):
call.K = K
values.append(call.present_value())
K = 100.
CPU times: user 973 ms, sys: 4.63 ms, total: 978 ms Wall time: 982 ms
The Monte Carlo estimator for the European call option's value with strike \(K\) is given by
\(C_0 \approx e^{-rT} \frac{1}{I} \sum_I \max[S_T(i)-K, 0]\)
where we can simulate the values of the risky asset at maturity using the log Euler discretization scheme
\(\log S_t = \log S_{t - \Delta t} + (r - 0.5 \sigma^2 ) \Delta t + \sigma \sqrt{\Delta t} z_t\)
class call_option_monte_carlo(object):
''' Class for European call option valuation in BSM model.
Monte Carlo simulation'''
import numpy as np
def __init__(self, S0, K, T, r, sigma, M, I):
self.S0 = float(S0)
self.K = K; self.T = T
self.r = r; self.sigma = sigma
self.M = M; self.I = I
def generate_paths(self):
dt = self.T / self.M
np.random.seed(0)
ran = np.random.standard_normal((self.M + 1, self.I))
ran = (ran - ran.mean()) / ran.std()
paths = S0 * np.exp(np.cumsum((self.r - 0.5 * self.sigma ** 2) * dt
+ self.sigma * np.sqrt(dt) * ran, axis=0))
paths[0] = S0
return paths
def present_value(self):
paths = self.generate_paths()
present_value = (np.exp(-self.r * self.T)
* np.sum(np.maximum(paths[-1] - self.K, 0)) / self.I)
return present_value
A first trial with the Monte Carlo based valuation approach.
call_mcs = call_option_monte_carlo(S0, K, T, r, sigma, M=50, I=10000)
%%time
print call_mcs.present_value()
10.4610933796 CPU times: user 48.7 ms, sys: 9.69 ms, total: 58.4 ms Wall time: 57.7 ms
Already quite good.
Adding more paths only improves the convergence of the MCS estimator in the limit.
%%time
call_mcs.I = 50000 # more paths
print call_mcs.present_value()
10.6068704869 CPU times: user 249 ms, sys: 41.5 ms, total: 291 ms Wall time: 290 ms
However, we also have the number/length of time intervals as a second parameter to improve convergence and precision; although at the cost of a longer computing time.
%%time
call_mcs.M = 500 # more steps
print call_mcs.present_value() # more time steps
10.4487236591 CPU times: user 2.47 s, sys: 679 ms, total: 3.14 s Wall time: 3.17 s
Potential problems with analytical formulae for pricing/valuation:
Sometimes significant advantages of analytical formulae for pricing/valuation:
Alleged drawbacks of alternative numerical methods (such as simulation):
On the other hand, many benefits from numerical methods (such as simulation):
Our class for the Monte Carlo valuation allows for valuable insights "behing the scenes".
plt.figure(figsize=(8, 5))
plt.hist(call_mcs.generate_paths()[-1], bins=30)
plt.grid(True); plt.xlabel('value at maturity'); plt.ylabel('frequency')
<matplotlib.text.Text at 0x106c8c0d0>
Not too long ago ...
Today, we are a bit more lucky ...
General ideas and approaches:
Global Valuation at its core:
Python specifics:
Lines of code (as of 12. Jan 2014): ca. 2,000
from DX_Lib_Valuation import *
One of the most simple classes is the one for constant short rate discounting.
r = constant_short_rate('short_rate', 0.05)
time_list = pd.date_range(start='2014/6/1', end='2014/12/31', freq='M').to_pydatetime()
time_list
array([datetime.datetime(2014, 6, 30, 0, 0), datetime.datetime(2014, 7, 31, 0, 0), datetime.datetime(2014, 8, 31, 0, 0), datetime.datetime(2014, 9, 30, 0, 0), datetime.datetime(2014, 10, 31, 0, 0), datetime.datetime(2014, 11, 30, 0, 0), datetime.datetime(2014, 12, 31, 0, 0)], dtype=object)
The class generates discount factors for a given list of datetime objects.
r.get_discount_factors(time_list)
array([[datetime.datetime(2014, 6, 30, 0, 0), 0.97510952647028781], [datetime.datetime(2014, 7, 31, 0, 0), 0.97925920727161508], [datetime.datetime(2014, 8, 31, 0, 0), 0.98329184073585585], [datetime.datetime(2014, 9, 30, 0, 0), 0.98747634223321412], [datetime.datetime(2014, 10, 31, 0, 0), 0.99154281422884716], [datetime.datetime(2014, 11, 30, 0, 0), 0.99576242860877573], [datetime.datetime(2014, 12, 31, 0, 0), 1.0]], dtype=object)
Alternatively, you can provide relative points in time as fractions of years.
time_list_delta = get_year_deltas(time_list)
r.get_discount_factors(time_list_delta, dtobjects=False)
array([[ 0. , 0.97510953], [ 0.08493151, 0.97925921], [ 0.16986301, 0.98329184], [ 0.25205479, 0.98747634], [ 0.3369863 , 0.99154281], [ 0.41917808, 0.99576243], [ 0.50410959, 1. ]])
For almost all (model & valuation) classes, we need to define a market environment.
me_jd = market_environment('jd', dt.datetime(2014, 7, 1))
type(me_jd)
DX_Lib_Frame.market_environment
A market environment mainly collects necessary constants and curves.
me_jd.add_constant('initial_value', 100.)
me_jd.add_constant('volatility', 0.2)
me_jd.add_constant('lambda', 0.7) # jump probability
me_jd.add_constant('mu', -0.8) # average jump size
me_jd.add_constant('delta', 0.2) # volatility of jump
me_jd.add_constant('currency', 'USD') # currency
me_jd.add_constant('frequency', 'B') # pandas date_range frequency string
me_jd.add_constant('paths', 10000) # number of paths
me_jd.add_constant('final_date', dt.datetime(2014, 12, 31))
me_jd.add_curve('discount_curve', r) # instance of discount class
Equipped with an instance of a market environment, you can initialize a model class, e.g. for a jump diffusion.
jd = jump_diffusion('jd_obj', me_jd)
paths = jd.get_instrument_values(fixed_seed=True)
paths[-1] # all values at final date
array([ 117.52748895, 105.12635975, 88.36302764, ..., 99.32911215, 92.31849354, 60.24410991])
Since you have access to all the data, detailed plots are only a couple of lines of code away.
plt.figure(figsize=(8, 5))
plt.plot(jd.time_grid, paths[:, 10:20])
plt.grid(True); plt.xlabel('date'); plt.ylabel('index level')
<matplotlib.text.Text at 0x10d2da910>
For valuation purposes, there are a number of valuation classes available, e.g. for European options with "arbitrary" payoffs. Again a market environment is needed.
me_call = market_environment('call', dt.datetime(2014, 7, 1))
me_call.add_constant('maturity', dt.datetime(2014, 12, 31))
me_call.add_constant('currency', 'USD')
call = valuation_mcs_european('call', me_call, jd, payoff_func='np.maximum(maturity_value - 100., 0)')
The valuation provides a number of methods, e.g. for valuation and Greeks estimation.
call.present_value()
4.808783
call.delta()
0.41847999999999885
call.vega()
19.11670000000001
Let us compile a couple of more statistics for the call option.
%%time
s_list = np.arange(80., 120.1, 5.)
pv = []; de = []; ve = []
for s in s_list:
call.update(initial_value=s)
pv.append(call.present_value())
de.append(call.delta())
ve.append(call.vega())
call.update(initial_value=100.)
CPU times: user 7.59 s, sys: 527 ms, total: 8.11 s Wall time: 8.14 s
The DX Analytics library provides a couple of standardized plotting functions.
plot_option_stats(s_list, pv, de, ve)
We now want to implement a more complex case with multiple, correlated risk factors and multiple derivatives.
The discounting class is central for risk-neutral pricing.
# constant short rate
r = constant_short_rate('r', 0.06)
Our market consists of three assets.
# market environments
me_gbm = market_environment('gbm', dt.datetime(2014, 1, 1))
me_jd = market_environment('jd', dt.datetime(2014, 1, 1))
me_sv = market_environment('sv', dt.datetime(2014, 1, 1))
Market environment for geometric Brownian motion.
# geometric Brownian motion
me_gbm.add_constant('initial_value', 36.)
me_gbm.add_constant('volatility', 0.2)
me_gbm.add_constant('currency', 'EUR')
me_gbm.add_constant('model', 'gbm')
Market environment for jump diffusion.
# jump diffusion
me_jd.add_constant('initial_value', 36.)
me_jd.add_constant('volatility', 0.2)
me_jd.add_constant('lambda', 0.5)
# probability for jump p.a.
me_jd.add_constant('mu', -0.75)
# expected jump size [%]
me_jd.add_constant('delta', 0.1)
# volatility of jump
me_jd.add_constant('currency', 'EUR')
me_jd.add_constant('model', 'jd')
Market environment for stochastic volatility model.
# stochastic volatility model
me_sv.add_constant('initial_value', 36.)
me_sv.add_constant('volatility', 0.2)
me_sv.add_constant('vol_vol', 0.1)
me_sv.add_constant('kappa', 2.5)
me_sv.add_constant('theta', 0.4)
me_sv.add_constant('rho', -0.5)
me_sv.add_constant('lambda', 0.0)
# probability for jump p.a.
me_sv.add_constant('mu', 0.0)
# expected jump size [%]
me_sv.add_constant('delta', 0.0)
# volatility of jump
me_sv.add_constant('currency', 'EUR')
me_sv.add_constant('model', 'svjd')
You can share common constants and curves by adding whole market environment objects to another one.
# valuation environment
val_env = market_environment('val_env', dt.datetime(2014, 1, 1))
val_env.add_constant('paths', 5000)
val_env.add_constant('frequency', 'W')
val_env.add_curve('discount_curve', r)
val_env.add_constant('starting_date', dt.datetime(2014, 1, 1))
val_env.add_constant('final_date', dt.datetime(2014, 12, 31))
# add valuation environment to market environments
me_gbm.add_environment(val_env)
me_jd.add_environment(val_env)
me_sv.add_environment(val_env)
We have now our market together. What is still missing are correlations.
underlyings = {'gbm' : me_gbm, 'jd' : me_jd, 'sv' : me_sv}
correlations = [['gbm', 'jd', 0.5], ['jd', 'sv', -0.75]]
First, the market model for the underlying.
gbm = geometric_brownian_motion('gbm_obj', me_gbm)
me_put = market_environment('put', dt.datetime(2014, 1, 1))
me_put.add_constant('maturity', dt.datetime(2014, 12, 31))
me_put.add_constant('strike', 40.)
me_put.add_constant('currency', 'EUR')
me_put.add_environment(val_env)
Then the model for the option.
am_put = valuation_mcs_american('am_put', mar_env=me_put, underlying=gbm,
payoff_func='np.maximum(strike - instrument_values, 0)')
am_put.present_value(fixed_seed=False, bf=5)
4.495371
We need a second underlying for the Maximum Call option.
jd = jump_diffusion('jd_obj', me_jd)
me_max_call = market_environment('put', dt.datetime(2014, 1, 1))
me_max_call.add_constant('maturity', dt.datetime(2014, 9, 15))
me_max_call.add_constant('currency', 'EUR')
me_max_call.add_environment(val_env)
Payoff definitions in this case must specify the underlying.
payoff_call = "np.maximum(np.maximum(maturity_value['gbm'], maturity_value['jd']) - 34., 0)"
We use a sub-set of the whole market for the valuation.
assets = {'gbm' : me_gbm, 'jd' : me_jd}
asset_corr = [correlations[0]]
max_call = valuation_mcs_multi_european('max_call', me_max_call, assets, asset_corr,
payoff_func=payoff_call)
max_call.present_value(fixed_seed=False)
5.69469
Of course, you now have two different Delta values.
max_call.delta('gbm')
0.5525599999999997
max_call.delta('jd')
0.32917000000000307
Let us generate a somehow broader overview of the option statistics. We start with the present value of the option given different initial values for the two underlyings.
asset_1 = np.arange(28., 46.1, 2.)
asset_2 = asset_1
a_1, a_2 = np.meshgrid(asset_1, asset_2)
value = np.zeros_like(a_1)
%%time
for i in range(np.shape(value)[0]):
for j in range(np.shape(value)[1]):
max_call.update('gbm', initial_value=a_1[i, j])
max_call.update('jd', initial_value=a_2[i, j])
value[i, j] = max_call.present_value()
CPU times: user 4.15 s, sys: 194 ms, total: 4.35 s Wall time: 4.35 s
And the graphical representation of the option present values.
plot_greeks_3d([a_1, a_2, value], ['gbm', 'jd', 'present value'])
Now the surfaces for the Greeks. Delta values first.
delta_1 = np.zeros_like(a_1)
delta_2 = np.zeros_like(a_1)
%%time
for i in range(np.shape(delta_1)[0]):
for j in range(np.shape(delta_1)[1]):
max_call.update('gbm', initial_value=a_1[i, j])
max_call.update('jd', initial_value=a_2[i, j])
delta_1[i, j] = max_call.delta('gbm')
delta_2[i, j] = max_call.delta('jd')
CPU times: user 13.5 s, sys: 636 ms, total: 14.1 s Wall time: 14.2 s
The surface for the gbm Delta.
plot_greeks_3d([a_1, a_2, delta_1], ['gbm', 'jd', 'delta gbm'])
The surface for the jd Delta.
plot_greeks_3d([a_1, a_2, delta_2], ['gbm', 'jd', 'delta jd'])
Finally, the values for the two Vega dimensions.
vega_1 = np.zeros_like(a_1)
vega_2 = np.zeros_like(a_1)
%%time
for i in range(np.shape(vega_1)[0]):
for j in range(np.shape(vega_1)[1]):
max_call.update('gbm', initial_value=a_1[i, j])
max_call.update('jd', initial_value=a_2[i, j])
vega_1[i, j] = max_call.vega('gbm')
vega_2[i, j] = max_call.vega('jd')
CPU times: user 12.1 s, sys: 564 ms, total: 12.7 s Wall time: 12.7 s
The surface for the gbm Vega.
plot_greeks_3d([a_1, a_2, vega_1], ['gbm', 'jd', 'vega gbm'])
The surface for the jd Vega.
plot_greeks_3d([a_1, a_2, vega_2], ['gbm', 'jd', 'vega jd'])
Finally, the third underlying needed.
sv = jump_diffusion('sv_obj', me_sv)
me_min_put = market_environment('min_put', dt.datetime(2014, 1, 1))
me_min_put.add_constant('maturity', dt.datetime(2014, 6, 17))
me_min_put.add_constant('currency', 'EUR')
me_min_put.add_environment(val_env)
This time, the derivative is a Minimum put option with American exercise.
payoff_put = "np.maximum(32. - np.minimum(instrument_values['jd'], instrument_values['sv']), 0)"
assets = {'jd' : me_jd, 'sv' : me_sv}
asset_corr = [correlations[1]]
asset_corr
[['jd', 'sv', -0.75]]
min_put = valuation_mcs_multi_american('min_put', val_env=me_min_put, assets=assets,
correlations=asset_corr, payoff_func=payoff_put)
min_put.present_value(fixed_seed=True)
4.620787
Again, two different values for the Greeks.
min_put.vega('jd')
4.802200000000045
min_put.vega('sv')
2.227300000000021
The derivative
class specifies a whole derivatives position.
am_put_pos = derivative('am_put_pos', 2, ['gbm'], me_put, 'American',
'np.maximum(instrument_values - 36., 0)')
max_call_pos = derivative('max_call_pos', 3, ['gbm', 'jd'], me_max_call,
'Multi European', payoff_call)
min_put_pos = derivative('min_put_pos', 5, ['sv', 'jd'], me_min_put,
'Multi American', payoff_put)
All is together to compose a portfolio with the derivatives positions as defined before.
positions = {'am_put_pos' : am_put_pos, 'max_call_pos' : max_call_pos,
'min_put_pos' : min_put_pos}
port = portfolio('portfolio', positions, val_env, underlyings, correlations)
DX Analytics allows for an integrated, consistent valuation of the portfolio given e.g. correlations between underlyings.
%%time
port.get_statistics()
CPU times: user 4.1 s, sys: 160 ms, total: 4.26 s Wall time: 4.24 s
position | name | quantity | value | currency | pos_value | pos_delta | pos_vega | |
---|---|---|---|---|---|---|---|---|
0 | max_call_pos | max_call_pos | 3 | 5.677053 | EUR | 17.031159 | {'jd': 0.9873, 'gbm': 1.68567} | {'jd': 18.3546, 'gbm': 25.8762} |
1 | am_put_pos | am_put_pos | 2 | 3.770114 | EUR | 7.540228 | 1.30004 | 26.019 |
2 | min_put_pos | min_put_pos | 5 | 4.617262 | EUR | 23.086310 | {'jd': -0.88965, 'sv': -0.84235} | {'jd': 25.251, 'sv': 9.617} |
Let's have a look at two paths of two correlated underlyings.
path_no = 20
paths1 = port.underlying_objects['gbm'].get_instrument_values()[:, path_no]
paths2 = port.underlying_objects['jd'].get_instrument_values()[:, path_no]
plt.figure(figsize=(9, 5))
plt.plot(port.time_grid, paths1, 'r', label='gbm')
plt.plot(port.time_grid, paths2, 'b', label='jd')
plt.gcf().autofmt_xdate()
plt.legend(loc=0); plt.grid(True)
# positive correlation case
The architecture and approach lends itself pretty well for risk management, scenario analysis and stress testing. For example, you could consider an alternative set of correlations.
correlations = [['gbm', 'jd', -0.7], ['jd', 'sv', 0.7]]
port = portfolio('portfolio', positions, val_env, underlyings, correlations)
port.get_statistics()
position | name | quantity | value | currency | pos_value | pos_delta | pos_vega | |
---|---|---|---|---|---|---|---|---|
0 | max_call_pos | max_call_pos | 3 | 6.833435 | EUR | 20.500305 | {'jd': 1.19421, 'gbm': 1.84566} | {'jd': 24.5532, 'gbm': 31.0476} |
1 | am_put_pos | am_put_pos | 2 | 3.780468 | EUR | 7.560936 | 1.27656 | 27.9452 |
2 | min_put_pos | min_put_pos | 5 | 4.182950 | EUR | 20.914750 | {'jd': -0.53865, 'sv': -1.0651} | {'jd': -0.392, 'sv': 4.795} |
Let's have a look at the same two paths of two correlated underlyings.
path_no = 20
paths1 = port.underlying_objects['gbm'].get_instrument_values()[:, path_no]
paths2 = port.underlying_objects['jd'].get_instrument_values()[:, path_no]
plt.figure(figsize=(9, 5))
plt.plot(port.time_grid, paths1, 'r', label='gbm')
plt.plot(port.time_grid, paths2, 'b', label='jd')
plt.gcf().autofmt_xdate()
plt.legend(loc=0); plt.grid(True)
# negative correlation case
DX Analytics leverages the experience from using Python for derivatives analytics of more than 10 years.
You can sign up for a DEXISION trial here www.dexision.com. All you have seen today and more in a standardized, Web-based fashion with a GUI.
DX Analytics is ...
10 years ago, Python was considered exotic in Finance – at best. Today, Python has become a major force in Finance due to a number of characteristics:
The Python Quants GmbH – the company Web site
Dr. Yves J. Hilpisch – my personal Web site
Derivatives Analytics with Python – my new book
Read an Excerpt and Order the Book
Contact Us
# Autosave Function for IPyton Notebook
def autosave(interval=3):
"""Autosave the notebook every interval (in minutes)"""
from IPython.core.display import Javascript
interval *= 60 * 1000 # JS wants intervals in miliseconds
tpl = 'setInterval ( "IPython.notebook.save_notebook()", %i );'
return Javascript(tpl % interval)
autosave()
<IPython.core.display.Javascript at 0x10d2df050>