Dr. Yves J. Hilpisch
The Python Quants GmbH
PyData London – 23 February 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
import numpy as np
values = []
for K in np.linspace(80, 120, 5000):
call.K = K
values.append(call.present_value())
K = 100.
CPU times: user 1.21 s, sys: 2.97 ms, total: 1.21 s Wall time: 1.21 s
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 41.2 ms, sys: 9.44 ms, total: 50.6 ms Wall time: 49.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 261 ms, sys: 39.9 ms, total: 301 ms Wall time: 300 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.32 s, sys: 715 ms, total: 3.04 s Wall time: 3.11 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".
import matplotlib.pyplot as plt
%matplotlib inline
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 0x1177141d0>
Not too long ago ...
Today, we are a bit more lucky ...
General ideas and approaches:
Global Valuation at its core:
Python specifics:
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) # unit zero-coupon bond values
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.
import matplotlib.pyplot as plt
%matplotlib inline
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 0x10d71a450>
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.63 s, sys: 526 ms, total: 8.15 s Wall time: 8.17 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.2 s, sys: 185 ms, total: 4.38 s Wall time: 4.4 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 11.4 s, sys: 517 ms, total: 11.9 s Wall time: 12 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 10.7 s, sys: 510 ms, total: 11.2 s Wall time: 11.2 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 3.5 s, sys: 132 ms, total: 3.64 s Wall time: 3.61 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} |
3 rows × 8 columns
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} |
3 rows × 8 columns
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
The necessary imports.
import numpy as np
import pandas as pd
import datetime as dt
import matplotlib.pyplot as plt
%matplotlib inline
from DX_Lib_Calibration import *
To illustrate the calibration of the Square Root Jump Diffusion model to VSTOXX call option data, we use real data from 28. December 2012.
term_structure = np.load('vstoxx_term_structure_28122012.npy')
option_data = np.load('vstoxx_option_data_28122012.npy')
len(option_data) # number of call options
39
option_data[:5] # maturity, strike, quote
array([[datetime.datetime(2013, 1, 16, 0, 0), 14.0, 7.45], [datetime.datetime(2013, 1, 16, 0, 0), 18.0, 3.6], [datetime.datetime(2013, 1, 16, 0, 0), 22.0, 1.4], [datetime.datetime(2013, 1, 16, 0, 0), 25.0, 0.75], [datetime.datetime(2013, 1, 16, 0, 0), 28.0, 0.45]], dtype=object)
The initialization for the calibration procedure.
r = constant_short_rate('r', 0.02)
srjd_me = market_environment('srjd_env', dt.datetime(2012, 12, 28))
srjd_me.add_curve('discount_curve', r)
srjd_me.add_curve('term_structure', term_structure)
srjd_me.add_curve('option_data', option_data)
srjd_me.add_constant('initial_value', 21.3534)
# VSTOXX level on 28.12.2012
srjd_me.add_constant('kappa', 5.)
srjd_me.add_constant('theta', 20.)
srjd_me.add_constant('volatility', 3.)
srjd_me.add_constant('lambda', 0.4)
srjd_me.add_constant('mu', 0.4)
srjd_me.add_constant('delta', .1)
srjd_me.add_constant('final_date', dt.datetime(2014, 6, 30))
srjd_me.add_constant('maturity', dt.datetime(2014, 6, 30))
srjd_me.add_constant('model', 'srjd')
srjd_me.add_constant('frequency', 'B')
srjd_me.add_constant('paths', 1000)
srjd_me.add_constant('strike', 21.)
The instantiation of the calibration class.
cali = srjd_calibration('cali', srjd_me)
First the VSTOXX term structure based on the futures data.
# dummy parametrization
kappa, theta, volatility, lamb, mu, delt = 2.0, 25., 3., 0.4, 0.3, 0.2
parameters = kappa, theta, volatility, lamb, mu, delt
cali.term_calibration(parameters=parameters)
Warning: Maximum number of function evaluations has been exceeded.
Let us have a look at the results.
plt.figure(figsize=(8, 4))
plt.subplot(211)
plt.plot(cali.forward_rates[:, 0], cali.forward_rates[:, 1], 'ro', label='Model')
plt.plot(term_structure[:, 0], term_structure[:, 1], label='Market')
plt.grid(True); plt.legend(loc=0)
plt.subplot(212)
plt.bar(cali.forward_rates[:, 0],
term_structure[:, 1] - cali.forward_rates[:, 1], width=3, label='Differences')
plt.grid(True); plt.legend(loc=0)
<matplotlib.legend.Legend at 0x10d9e7790>
Setting parameter ranges for the global optimization.
parameter_ranges = ((.25, 1.251, 1.0), (15., 20.01, 5.), (1., 2.01, 1.),
(0.4, 0.71, 0.3), (0.4, 0.71, 0.3), (0.1, 0.3, 0.2))
# kappa, theta, volatility,
# lambda, mu, delta
Running the global optimization (brute force).
%time cali.global_optimization(parameter_ranges)
[ 0.25 15. 1. 0.4 0.4 0.1 ] 4.29 4.29 [ 0.25 15. 1. 0.4 0.4 0.1 ] 4.29 4.29 [ 0.25 15. 1. 0.4 0.7 0.1 ] 0.39 0.39 [ 0.25 15. 1. 0.7 0.4 0.1 ] 0.39 1.54 [ 0.25 15. 1. 0.7 0.7 0.1 ] 0.39 4.11 [ 0.25 15. 2. 0.4 0.4 0.1 ] 0.39 2.43 [ 0.25 15. 2. 0.4 0.7 0.1 ] 0.11 0.11 [ 0.25 15. 2. 0.7 0.4 0.1 ] 0.11 0.59 [ 0.25 15. 2. 0.7 0.7 0.1 ] 0.11 6.1 [ 0.25 20. 1. 0.4 0.4 0.1 ] 0.11 3.63 [ 0.25 20. 1. 0.4 0.7 0.1 ] 0.11 0.23 [ 0.25 20. 1. 0.7 0.4 0.1 ] 0.11 1.13 [ 0.25 20. 1. 0.7 0.7 0.1 ] 0.11 5.34 [ 0.25 20. 2. 0.4 0.4 0.1 ] 0.11 1.94 [ 0.25 20. 2. 0.4 0.7 0.1 ] 0.11 0.16 [ 0.25 20. 2. 0.7 0.4 0.1 ] 0.11 0.41 [ 0.25 20. 2. 0.7 0.7 0.1 ] 0.11 7.61 [ 1.25 15. 1. 0.4 0.4 0.1 ] 0.11 8.33 [ 1.25 15. 1. 0.4 0.7 0.1 ] 0.11 3.65 [ 1.25 15. 1. 0.7 0.4 0.1 ] 0.11 5.64 [ 1.25 15. 1. 0.7 0.7 0.1 ] 0.11 0.26 [ 1.25 15. 2. 0.4 0.4 0.1 ] 0.11 6.49 [ 1.25 15. 2. 0.4 0.7 0.1 ] 0.11 2.56 [ 1.25 15. 2. 0.7 0.4 0.1 ] 0.11 4.11 [ 1.25 15. 2. 0.7 0.7 0.1 ] 0.11 0.19 [ 1.25 20. 1. 0.4 0.4 0.1 ] 0.11 5.5 [ 1.25 20. 1. 0.4 0.7 0.1 ] 0.11 1.51 [ 1.25 20. 1. 0.7 0.4 0.1 ] 0.11 2.98 [ 1.25 20. 1. 0.7 0.7 0.1 ] 0.11 0.7 [ 1.25 20. 2. 0.4 0.4 0.1 ] 0.11 3.9 [ 1.25 20. 2. 0.4 0.7 0.1 ] 0.11 0.83 [ 1.25 20. 2. 0.7 0.4 0.1 ] 0.11 1.86 [ 1.25 20. 2. 0.7 0.7 0.1 ] 0.11 1.3 CPU times: user 24.4 s, sys: 1.1 s, total: 25.5 s Wall time: 25.5 s
The "global" results (I).
cali.opt_parameters
array([ 0.25, 15. , 2. , 0.4 , 0.7 , 0.1 ])
cali.min_err_value
0.11093010628186292
The "global" results (II).
plot_calibration_results(cali)
Now the local optimization, using the results from the global optimization.
%time cali.local_optimization(maxi=100)
[ 0.25 15. 2. 0.4 0.7 0.1 ] 0.11 0.11 [ 0.26 15. 2. 0.4 0.7 0.1 ] 0.11 0.11 [ 0.25 15.75 2. 0.4 0.7 0.1 ] 0.11 0.11 [ 0.25 15. 2.1 0.4 0.7 0.1 ] 0.11 0.11 [ 0.25 15. 2. 0.42 0.7 0.1 ] 0.1 0.1 [ 0.25 15. 2. 0.4 0.74 0.1 ] 0.1 0.13 [ 0.25 15. 2. 0.4 0.7 0.11] 0.1 0.11 [ 0.25 15.25 2.03 0.41 0.67 0.1 ] 0.1 0.12 [ 0.25 15.19 2.02 0.4 0.68 0.1 ] 0.09 0.09 [ 0.24 15.31 2.04 0.41 0.69 0.1 ] 0.09 0.09 [ 0.23 15.47 2.06 0.41 0.69 0.1 ] 0.09 0.11 [ 0.25 14.42 2.06 0.41 0.69 0.1 ] 0.09 0.12 [ 0.25 15.42 2.01 0.4 0.7 0.1 ] 0.09 0.1 [ 0.25 15.31 2.06 0.41 0.69 0.1 ] 0.09 0.11 [ 0.25 15.41 2.08 0.42 0.69 0.1 ] 0.09 0.11 [ 0.24 15.54 1.97 0.42 0.69 0.1 ] 0.09 0.1 [ 0.25 15.18 1.96 0.41 0.7 0.11] 0.08 0.08 [ 0.25 15.07 1.9 0.4 0.7 0.11] 0.08 0.11 [ 0.25 15.24 1.94 0.41 0.69 0.1 ] 0.08 0.11 [ 0.25 15.29 2.03 0.41 0.69 0.1 ] 0.08 0.12 [ 0.24 15.25 2. 0.41 0.69 0.1 ] 0.08 0.08 [ 0.25 15.18 1.99 0.41 0.69 0.1 ] 0.08 0.09 [ 0.25 15.09 1.98 0.41 0.7 0.1 ] 0.08 0.16 [ 0.25 15.3 1.99 0.4 0.7 0.1 ] 0.08 0.09 [ 0.25 15.36 1.97 0.41 0.69 0.1 ] 0.08 0.12 [ 0.25 15.24 2.01 0.41 0.69 0.1 ] 0.08 0.1 [ 0.25 15.41 1.99 0.4 0.69 0.11] 0.08 0.11 [ 0.25 15.16 2.01 0.4 0.7 0.1 ] 0.08 0.1 [ 0.25 15.02 1.99 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 14.83 2. 0.41 0.7 0.1 ] 0.08 0.11 [ 0.25 15.12 1.97 0.4 0.7 0.1 ] 0.08 0.12 [ 0.25 15.21 2. 0.41 0.69 0.1 ] 0.08 0.08 [ 0.25 15.22 1.96 0.42 0.69 0.1 ] 0.08 0.11 [ 0.25 15.18 2. 0.4 0.7 0.1 ] 0.08 0.1 [ 0.25 15.1 1.98 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 15.12 2. 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 15.14 2. 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 15.16 1.99 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 15.1 1.99 0.41 0.69 0.1 ] 0.08 0.08 [ 0.25 15.09 2. 0.4 0.7 0.1 ] 0.08 0.1 [ 0.25 15.12 1.98 0.41 0.7 0.1 ] 0.08 0.11 [ 0.25 15.1 2. 0.41 0.7 0.1 ] 0.08 0.09 [ 0.25 15.11 1.99 0.41 0.7 0.1 ] 0.08 0.11 [ 0.25 15.1 1.99 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 15.11 1.99 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 15.07 1.98 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 15.04 1.98 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 15.07 1.99 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 15.08 1.99 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 15.05 1.98 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 14.97 1.98 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 14.99 1.97 0.41 0.7 0.1 ] 0.08 0.08 [ 0.26 14.93 1.96 0.41 0.7 0.1 ] 0.08 0.08 [ 0.25 14.96 1.97 0.41 0.7 0.1 ] 0.08 0.08 [ 0.26 14.89 1.97 0.41 0.7 0.1 ] 0.08 0.08 [ 0.26 14.78 1.97 0.41 0.71 0.1 ] 0.08 0.08 [ 0.26 14.89 1.95 0.41 0.7 0.1 ] 0.08 0.08 [ 0.26 14.8 1.96 0.41 0.71 0.1 ] 0.08 0.08 [ 0.26 14.68 1.95 0.41 0.71 0.11] 0.08 0.08 [ 0.26 14.79 1.94 0.41 0.71 0.1 ] 0.08 0.08 [ 0.26 14.64 1.93 0.41 0.71 0.1 ] 0.08 0.08 [ 0.27 14.43 1.91 0.41 0.72 0.1 ] 0.08 0.1 [ 0.27 14.61 1.93 0.41 0.71 0.1 ] 0.08 0.08 [ 0.28 14.44 1.92 0.41 0.72 0.1 ] 0.08 0.08 [ 0.27 14.47 1.93 0.41 0.72 0.1 ] 0.08 0.08 [ 0.27 14.38 1.93 0.41 0.72 0.1 ] 0.08 0.1 [ 0.26 14.76 1.94 0.41 0.71 0.1 ] 0.08 0.08 [ 0.27 14.47 1.9 0.41 0.72 0.11] 0.08 0.08 [ 0.27 14.37 1.91 0.41 0.72 0.11] 0.08 0.08 [ 0.28 14.16 1.9 0.41 0.72 0.11] 0.08 0.1 [ 0.28 14.26 1.9 0.41 0.72 0.11] 0.08 0.1 [ 0.27 14.64 1.93 0.41 0.71 0.1 ] 0.08 0.08 [ 0.28 14.33 1.89 0.41 0.72 0.1 ] 0.07 0.07 [ 0.29 14.15 1.87 0.41 0.73 0.1 ] 0.07 0.07 [ 0.29 14.21 1.89 0.41 0.72 0.1 ] 0.07 0.07 [ 0.29 14.07 1.87 0.41 0.73 0.1 ] 0.07 0.07 [ 0.3 13.78 1.84 0.41 0.74 0.11] 0.07 0.07 [ 0.29 14. 1.84 0.41 0.73 0.11] 0.07 0.07 [ 0.3 13.84 1.86 0.41 0.73 0.1 ] 0.07 0.07 [ 0.31 13.68 1.82 0.41 0.74 0.11] 0.07 0.07 [ 0.32 13.31 1.77 0.41 0.75 0.11] 0.07 0.07 [ 0.32 13.4 1.78 0.41 0.75 0.1 ] 0.07 0.07 [ 0.32 13.28 1.77 0.41 0.75 0.1 ] 0.07 0.07 [ 0.33 13.05 1.75 0.41 0.75 0.11] 0.07 0.07 [ 0.36 12.5 1.69 0.41 0.77 0.11] 0.07 0.07 [ 0.34 12.89 1.74 0.41 0.76 0.1 ] 0.07 0.07 [ 0.35 12.73 1.69 0.41 0.77 0.11] 0.07 0.07 [ 0.36 12.43 1.66 0.41 0.77 0.11] 0.07 0.07 [ 0.32 13.45 1.8 0.41 0.75 0.11] 0.07 0.07 [ 0.34 12.84 1.73 0.41 0.76 0.11] 0.07 0.07 [ 0.34 12.56 1.71 0.41 0.77 0.11] 0.07 0.07 [ 0.32 13.33 1.76 0.41 0.75 0.11] 0.07 0.07 [ 0.32 13.22 1.76 0.41 0.75 0.11] 0.07 0.07 [ 0.3 13.65 1.83 0.41 0.74 0.11] 0.07 0.07 [ 0.34 12.96 1.73 0.41 0.76 0.11] 0.07 0.07 [ 0.33 13.08 1.76 0.41 0.76 0.11] 0.07 0.07 [ 0.33 13. 1.76 0.41 0.76 0.1 ] 0.07 0.07 [ 0.34 12.73 1.71 0.41 0.76 0.11] 0.07 0.07 [ 0.32 13.27 1.77 0.41 0.75 0.11] 0.07 0.07 [ 0.34 12.88 1.74 0.41 0.76 0.11] 0.07 0.07 [ 0.34 12.72 1.72 0.41 0.76 0.11] 0.07 0.07 [ 0.34 12.86 1.73 0.41 0.76 0.11] 0.07 0.07 [ 0.32 13.24 1.77 0.41 0.75 0.11] 0.07 0.07 [ 0.33 13.1 1.77 0.41 0.75 0.11] 0.07 0.07 [ 0.32 13.16 1.79 0.41 0.75 0.11] 0.07 0.07 [ 0.34 12.73 1.73 0.41 0.76 0.11] 0.07 0.07 [ 0.35 12.59 1.72 0.41 0.77 0.11] 0.07 0.07 [ 0.33 13.08 1.75 0.41 0.76 0.11] 0.07 0.07 [ 0.33 13.07 1.76 0.41 0.76 0.11] 0.07 0.07 [ 0.33 13.19 1.78 0.4 0.75 0.11] 0.07 0.07 [ 0.32 13.29 1.78 0.41 0.75 0.11] 0.07 0.07 [ 0.34 12.87 1.74 0.41 0.76 0.11] 0.07 0.07 [ 0.32 13.13 1.77 0.41 0.75 0.11] 0.07 0.07 [ 0.32 13.25 1.79 0.41 0.75 0.11] 0.07 0.07 [ 0.33 12.96 1.75 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.99 1.76 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.95 1.77 0.41 0.76 0.1 ] 0.07 0.07 [ 0.34 12.89 1.74 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.88 1.75 0.41 0.76 0.11] 0.07 0.07 [ 0.33 13.09 1.77 0.41 0.76 0.11] 0.07 0.07 [ 0.33 13.05 1.78 0.41 0.76 0.11] 0.07 0.07 [ 0.33 13.09 1.79 0.41 0.75 0.11] 0.07 0.07 [ 0.32 13.19 1.79 0.41 0.75 0.11] 0.07 0.07 [ 0.33 12.96 1.76 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.93 1.76 0.41 0.76 0.11] 0.07 0.07 [ 0.33 13.06 1.77 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.95 1.77 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.89 1.77 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.87 1.76 0.41 0.76 0.11] 0.07 0.07 [ 0.34 12.84 1.76 0.41 0.76 0.11] 0.07 0.07 [ 0.33 13. 1.79 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.9 1.78 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.85 1.79 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.89 1.8 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.85 1.82 0.41 0.75 0.11] 0.07 0.07 [ 0.34 12.68 1.77 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.83 1.81 0.41 0.76 0.11] 0.07 0.07 [ 0.33 12.86 1.82 0.41 0.75 0.11] 0.07 0.07 [ 0.34 12.65 1.81 0.41 0.76 0.11] 0.07 0.07 [ 0.34 12.48 1.82 0.41 0.76 0.11] 0.07 0.09 [ 0.34 12.68 1.83 0.41 0.75 0.11] 0.07 0.07 [ 0.34 12.57 1.86 0.41 0.75 0.11] 0.07 0.07 [ 0.33 12.86 1.86 0.41 0.75 0.11] 0.07 0.07 [ 0.33 12.7 1.87 0.41 0.75 0.11] 0.07 0.07 [ 0.33 12.66 1.87 0.41 0.75 0.11] 0.07 0.07 [ 0.33 12.57 1.88 0.41 0.75 0.11] 0.07 0.07 [ 0.34 12.43 1.91 0.41 0.75 0.11] 0.07 0.07 [ 0.34 12.44 1.91 0.41 0.75 0.11] 0.07 0.1 [ 0.33 12.75 1.84 0.41 0.75 0.11] 0.07 0.07 [ 0.33 12.67 1.93 0.41 0.75 0.11] 0.07 0.07 Warning: Maximum number of function evaluations has been exceeded. CPU times: user 1min 44s, sys: 4.54 s, total: 1min 49s Wall time: 1min 49s
The "local" results (I).
cali.opt_parameters
array([ 0.33544549, 12.4302789 , 1.90633636, 0.40903462, 0.75270964, 0.11147443])
cali.min_err_value
0.065486741801247944
The "local" results (II).
plot_calibration_results(cali)
The "local" results (III).
plot_calibration_results(cali, relative=True)
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.derivatives-analytics.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 0x10b886fd0>