This setion of the documentation illustrates how to calculate implied volatilities and how to calibrate a model to VSTOXX volatility index call option quotes. The example implements the calibration for a total of one month worth of data.
In [1]:
from dx import *
import numpy as np
import pandas as pd
import seaborn as sns; sns.set()
We start by loading VSTOXX data from a pandas HDFStore into DataFrame objects (source: Eurex, cf. http://www.eurexchange.com/advanced-services/).
In [2]:
h5 = pd.HDFStore('./data/vstoxx_march_2014.h5', 'r')
vstoxx_index = h5['vstoxx_index']
vstoxx_futures = h5['vstoxx_futures']
vstoxx_options = h5['vstoxx_options']
h5.close()
VSTOXX index for the first quarter of 2014.
In [3]:
%matplotlib inline
vstoxx_index['V2TX'].plot(figsize=(10, 6))
Out[3]:
The VSTOXX futures data (8 futures maturities/quotes per day).
In [4]:
vstoxx_futures.info()
In [5]:
vstoxx_futures.tail()
Out[5]:
The VSTOXX options data. This data set is quite large due to the large number of European put and call options on the VSTOXX.
In [6]:
vstoxx_options.info()
In [7]:
vstoxx_options.tail()
Out[7]:
As a helper function we need a function to calculate all relevant third Fridays for all relevant maturity months of the data sets.
In [8]:
import datetime as dt
import calendar
def third_friday(date):
day = 21 - (calendar.weekday(date.year, date.month, 1) + 2) % 7
return dt.datetime(date.year, date.month, day)
In [9]:
third_fridays = {}
for month in set(vstoxx_futures['EXP_MONTH']):
third_fridays[month] = third_friday(dt.datetime(2014, month, 1))
In [10]:
third_fridays
Out[10]:
Often calibration efforts are undertaken to replicate the market implied volatilities or the so-called volatility surface as good as possible. With DX Analytics and the BSM_european_option class, you can efficiently calculate (i.e. numerically estimate) implied volatilities. For the example, we use the VSTOXX futures and call options data from 31. March 2014.
Some definitions, the pre-selection of option data and the pre-definition of the market environment needed.
In [11]:
V0 = 17.6639 # VSTOXX level on 31.03.2014
futures_data = vstoxx_futures[vstoxx_futures.DATE == '2014/3/31'].copy()
options_data = vstoxx_options[(vstoxx_options.DATE == '2014/3/31')
& (vstoxx_options.TYPE == 'C')].copy()
me = market_environment('me', dt.datetime(2014, 3, 31))
me.add_constant('initial_value', 17.6639) # index on 31.03.2014
me.add_constant('volatility', 2.0) # for initialization
me.add_curve('discount_curve', constant_short_rate('r', 0.01)) # assumption
options_data['IMP_VOL'] = 0.0 # initialization new iv column
The following loop now calculates the implied volatilities for all those options whose strike lies within the defined tolerance level.
In [12]:
%%time
tol = 0.3 # tolerance level for moneyness
for option in options_data.index:
# iterating over all option quotes
forward = futures_data[futures_data['MATURITY'] == \
options_data.loc[option]['MATURITY']]['PRICE'].values
# picking the right futures value
if (forward * (1 - tol) < options_data.loc[option]['STRIKE']
< forward * (1 + tol)):
# only for options with moneyness within tolerance
call = options_data.loc[option]
me.add_constant('strike', call['STRIKE'])
me.add_constant('maturity', call['MATURITY'])
call_option = BSM_european_option('call', me)
options_data.loc[option, 'IMP_VOL'] = \
call_option.imp_vol(call['PRICE'], 'call', volatility_est=0.6)
A selection of the results.
In [13]:
options_data[60:70]
Out[13]:
And the complete results visualized.
In [14]:
import matplotlib.pyplot as plt
%matplotlib inline
plot_data = options_data[options_data.IMP_VOL > 0]
plt.figure(figsize=(10, 6))
for maturity in sorted(set(options_data['MATURITY'])):
data = plot_data.isin({'MATURITY': [maturity,]})
data = plot_data[plot_data.MATURITY == maturity]
# select data for this maturity
plt.plot(data['STRIKE'], data['IMP_VOL'],
label=maturity.date(), lw=1.5)
plt.plot(data['STRIKE'], data['IMP_VOL'], 'r.')
plt.xlabel('strike')
plt.ylabel('implied volatility of volatility')
plt.legend()
plt.show()
This sub-section now implements the model calibration based on selected options data. In particular, we choose, for a given pricing date, the following options data:
The following following returns the relevant market data per calibration date:
In [15]:
tol = 0.2
def get_option_selection(pricing_date, maturity, tol=tol):
''' Function selects relevant options data. '''
forward = vstoxx_futures[(vstoxx_futures.DATE == pricing_date)
& (vstoxx_futures.MATURITY == maturity)]['PRICE'].values[0]
option_selection = \
vstoxx_options[(vstoxx_options.DATE == pricing_date)
& (vstoxx_options.MATURITY == maturity)
& (vstoxx_options.TYPE == 'C')
& (vstoxx_options.STRIKE > (1 - tol) * forward)
& (vstoxx_options.STRIKE < (1 + tol) * forward)]
return option_selection, forward
Given the options and their respective quotes to which to calibrate the model, the function get_option_models returns the DX Analytics option models for all relevant options. As risk factor model we use the square_root_diffusion class.
In [16]:
def get_option_models(pricing_date, maturity, option_selection):
''' Models and returns traded options for given option_selection object. '''
me_vstoxx = market_environment('me_vstoxx', pricing_date)
initial_value = vstoxx_index['V2TX'][pricing_date]
me_vstoxx.add_constant('initial_value', initial_value)
me_vstoxx.add_constant('final_date', maturity)
me_vstoxx.add_constant('currency', 'EUR')
me_vstoxx.add_constant('frequency', 'W')
me_vstoxx.add_constant('paths', 10000)
csr = constant_short_rate('csr', 0.01)
# somewhat arbitrarily chosen here
me_vstoxx.add_curve('discount_curve', csr)
# parameters to be calibrated later
me_vstoxx.add_constant('kappa', 1.0)
me_vstoxx.add_constant('theta', 1.2 * initial_value)
me_vstoxx.add_constant('volatility', 1.0)
vstoxx_model = square_root_diffusion('vstoxx_model', me_vstoxx)
# square-root diffusion for volatility modeling
# mean-reverting, positive process
# option parameters and payoff
me_vstoxx.add_constant('maturity', maturity)
payoff_func = 'np.maximum(maturity_value - strike, 0)'
option_models = {}
for option in option_selection.index:
strike = option_selection['STRIKE'].ix[option]
me_vstoxx.add_constant('strike', strike)
option_models[option] = \
valuation_mcs_european_single(
'eur_call_%d' % strike,
vstoxx_model,
me_vstoxx,
payoff_func)
return vstoxx_model, option_models
The function calculate_model_values estimates and returns model value estimates for all relevant options given a parameter set for the square_root_diffusion risk factor model.
In [17]:
def calculate_model_values(p0):
''' Returns all relevant option values.
Parameters
===========
p0 : tuple/list
tuple of kappa, theta, volatility
Returns
=======
model_values : dict
dictionary with model values
'''
kappa, theta, volatility = p0
vstoxx_model.update(kappa=kappa,
theta=theta,
volatility=volatility)
model_values = {}
for option in option_models:
model_values[option] = \
option_models[option].present_value(fixed_seed=True)
return model_values
The calibration of the pricing model is based on the minimization of the mean-squared error (MSE) of the model values vs. the market quotes. The MSE calculation is implemented by the function mean_squared_error which also penalizes economically implausible parameter values.
In [18]:
i = 0
def mean_squared_error(p0):
''' Returns the mean-squared error given
the model and market values.
Parameters
===========
p0 : tuple/list
tuple of kappa, theta, volatility
Returns
=======
MSE : float
mean-squared error
'''
if p0[0] < 0 or p0[1] < 5. or p0[2] < 0 or p0[2] > 10.:
return 100
global i, option_selection, vstoxx_model, option_models, first, last
pd = dt.datetime.strftime(
option_selection['DATE'].iloc[0].to_pydatetime(),
'%d-%b-%Y')
mat = dt.datetime.strftime(
option_selection['MATURITY'].iloc[0].to_pydatetime(),
'%d-%b-%Y')
model_values = calculate_model_values(p0)
option_diffs = {}
for option in model_values:
option_diffs[option] = (model_values[option]
- option_selection['PRICE'].loc[option])
MSE = np.sum(np.array(option_diffs.values()) ** 2) / len(option_diffs)
if i % 150 == 0:
# output every 0th and 100th iteration
if i == 0:
print '%12s %13s %4s %6s %6s %6s --> %6s' % \
('pricing_date', 'maturity_date', 'i', 'kappa',
'theta', 'vola', 'MSE')
print '%12s %13s %4d %6.3f %6.3f %6.3f --> %6.3f' % \
(pd, mat, i, p0[0], p0[1], p0[2], MSE)
i += 1
return MSE
The function get_parameter_series calibrates the model to the market data for every date contained in the pricing_date_list object for all maturities contained in the maturity_list object.
In [19]:
import scipy.optimize as spo
def get_parameter_series(pricing_date_list, maturity_list):
global i, option_selection, vstoxx_model, option_models, first, last
# collects optimization results for later use (eg. visualization)
parameters = pd.DataFrame()
for maturity in maturity_list:
first = True
for pricing_date in pricing_date_list:
option_selection, forward = \
get_option_selection(pricing_date, maturity)
vstoxx_model, option_models = \
get_option_models(pricing_date, maturity, option_selection)
if first is True:
# use brute force for the first run
i = 0
opt = spo.brute(mean_squared_error,
((0.5, 2.51, 1.), # range for kappa
(10., 20.1, 5.), # range for theta
(0.5, 10.51, 5.0)), # range for volatility
finish=None)
i = 0
opt = spo.fmin(mean_squared_error, opt,
maxiter=200, maxfun=350, xtol=0.0000001, ftol=0.0000001)
parameters = parameters.append(
pd.DataFrame(
{'pricing_date' : pricing_date,
'maturity' : maturity,
'initial_value' : vstoxx_model.initial_value,
'kappa' : opt[0],
'theta' : opt[1],
'sigma' : opt[2],
'MSE' : mean_squared_error(opt)}, index=[0,]),
ignore_index=True)
first = False
last = opt
return parameters
This completes the set of necessary function to implement such a larger calibration effort. The following code defines the dates for which a calibration shall be conducted and for which maturities the calibration is required.
In [20]:
%%time
pricing_date_list = pd.date_range('2014/3/1', '2014/3/31', freq='B')
maturity_list = [third_fridays[7]]
parameters = get_parameter_series(pricing_date_list, maturity_list)
The results are now stored in the pandas DataFrame called parameters. We set a new index and inspect the last results. Throughout the MSE is pretty low indicated a good fit of the model to the market quotes.
In [21]:
paramet = parameters.set_index('pricing_date')
paramet.tail()
Out[21]:
This is also illustrated by the visualization of the time series data for the calibrated/optimal parameter values. The MSE is below 0.01 throughout.
In [22]:
%matplotlib inline
paramet[['kappa', 'theta', 'sigma', 'MSE']].plot(subplots=True, color='b', figsize=(10, 12))
plt.tight_layout()
The following generates a plot of the calibration results for the last calibration day. The absolute price differences are below 0.10 EUR for all options.
In [23]:
index = paramet.index[-1]
opt = np.array(paramet[['kappa', 'theta', 'sigma']].loc[index])
option_selection = get_option_selection(index, maturity_list[0], tol=tol)[0]
model_values = np.sort(np.array(calculate_model_values(opt).values()))[::-1]
import matplotlib.pyplot as plt
%matplotlib inline
fix, (ax1, ax2) = plt.subplots(2, sharex=True, figsize=(10, 8))
strikes = option_selection['STRIKE'].values
ax1.plot(strikes, option_selection['PRICE'], label='market quotes')
ax1.plot(strikes, model_values, 'ro', label='model values')
ax1.set_ylabel('option values')
ax1.grid(True)
ax1.legend(loc=0)
wi = 0.25
ax2.bar(strikes - wi / 2., model_values - option_selection['PRICE'],
label='market quotes', width=wi)
ax2.grid(True)
ax2.set_ylabel('differences')
ax2.set_xlabel('strikes')
Out[23]:
Copyright, License & Disclaimer
© Dr. Yves J. Hilpisch | The Python Quants GmbH
DX Analytics (the "dx library") is licensed under the GNU Affero General Public License version 3 or later (see http://www.gnu.org/licenses/).
DX Analytics comes with no representations or warranties, to the extent permitted by applicable law.
http://tpq.io | team@tpq.io | http://twitter.com/dyjh
Quant Platform | http://quant-platform.com
Derivatives Analytics with Python (Wiley Finance) | http://derivatives-analytics-with-python.com
Python for Finance (O'Reilly) | http://python-for-finance.com