The Python Quants

Market-Based Valuation of Volatility Options

Square-Root Jump Diffusion Plus for All Maturities

Yves Hilpisch

yves@pythonquants.com

www.pythonquants.com

The Python Quants GmbH

Data Import and Selection

In [1]:
from dx import *
  # dx library
import numpy as np
import pandas as pd
import datetime as dt
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()
In [3]:
third_fridays = sorted(set(vstoxx_futures['MATURITY']))
In [4]:
tol = 0.20  # otm & itm options
def get_option_selection(pricing_date, tol=tol):
    futures = vstoxx_futures[vstoxx_futures.DATE == pricing_date]# ['PRICE'].values[0]
    option_selection = pd.DataFrame()
    for mat in third_fridays[3:10]:
        forward = futures[futures.MATURITY == mat]['PRICE'].values[0]
        option_selection = option_selection.append(
            vstoxx_options[(vstoxx_options.DATE == pricing_date)
                         & (vstoxx_options.MATURITY == mat)
                         & (vstoxx_options.TYPE == 'C')  # only calls
                         & (vstoxx_options.STRIKE > (1 - tol) * forward)
                         & (vstoxx_options.STRIKE < (1 + tol) * forward)])
    return option_selection, futures
In [5]:
def srd_forwards(v0, p0, t):
    ''' Function for forward vols/rates in SRD model.
    kappa: float
        mean-reversion factor
    theta:
        long-run mean
    sigma:
        volatility factor (vol-vol)
    '''
    kappa, theta, sigma = p0
    t = get_year_deltas(t)
    g = math.sqrt(kappa ** 2 + 2 * sigma ** 2)
    sum1 = ((kappa * theta * (np.exp(g * t) - 1)) /
          (2 * g + (kappa + g) * (np.exp(g * t) - 1)))
    sum2 = v0 * ((4 * g ** 2 * np.exp(g * t)) /
            (2 * g + (kappa + g) * (np.exp(g * t) - 1)) ** 2)
    return sum1 + sum2
In [6]:
def srd_forward_error(p0):
    global initial_value, f, t
    if p0[0] < 0 or p0[1] < 0 or p0[2] < 0:
        return 100
    f_model = srd_forwards(initial_value, p0, t)
    MSE = np.sum((f - f_model) ** 2) / len(f)
    return MSE
In [7]:
from scipy.optimize import fmin
def generate_shift_base(pricing_date, futures):
    global initial_value, f, t
    # futures price array
    f = list(futures['PRICE'].values)
    f.insert(0, initial_value)
    f = np.array(f)
    # date array
    t = [_.to_pydatetime() for _ in futures['MATURITY']]
    t.insert(0, pricing_date)
    t = np.array(t)
    # calibration
    opt = fmin(srd_forward_error, (2., 15., 2.))
    # shift_calculation
    f_model = srd_forwards(initial_value, opt, t)
    shifts = f - f_model
    shift_base = np.array((t, shifts)).T
    return shift_base

Example Term Structure Calibration

In [8]:
pdate = pd.Timestamp('2014-03-31')
os, futures = get_option_selection(pdate)
initial_value = vstoxx_index['V2TX'].loc[pdate]
# futures price array
f = list(futures['PRICE'].values)
f.insert(0, initial_value)
f = np.array(f)
# date array
t = [_.to_pydatetime() for _ in futures['MATURITY']]
t.insert(0, pdate)
t = np.array(t)
In [9]:
from scipy.optimize import fmin
opt = fmin(srd_forward_error, (2., 15., 2.))
Optimization terminated successfully.
         Current function value: 0.042650
         Iterations: 312
         Function evaluations: 567

In [10]:
import matplotlib.pyplot as plt
%matplotlib inline
In [11]:
f_model = srd_forwards(initial_value, opt, t)
fig = plt.figure(figsize=(9, 5))
plt.plot(t, f_model, 'ro', label='model')
plt.plot(t, f, label='market'); plt.legend(loc=0)
plt.grid(True); fig.autofmt_xdate()
In [12]:
generate_shift_base(pdate, futures)
  # base shift values = diffs between model and market
Optimization terminated successfully.
         Current function value: 0.042650
         Iterations: 312
         Function evaluations: 567

Out[12]:
array([[Timestamp('2014-03-31 00:00:00'), 0.0],
       [datetime.datetime(2014, 4, 18, 0, 0), -0.5258174028562053],
       [datetime.datetime(2014, 5, 16, 0, 0), 0.3216763257413966],
       [datetime.datetime(2014, 6, 20, 0, 0), -0.02332543663536768],
       [datetime.datetime(2014, 7, 18, 0, 0), 0.011483931309037132],
       [datetime.datetime(2014, 8, 15, 0, 0), 0.009521601492053833],
       [datetime.datetime(2014, 9, 19, 0, 0), -0.005069145560565147],
       [datetime.datetime(2014, 10, 17, 0, 0), -0.05281655212520775],
       [datetime.datetime(2014, 11, 21, 0, 0), 0.017621432524546066]], dtype=object)

Option Modelling with DX Analytics

In [13]:
def get_option_models(pricing_date, option_selection, futures):
    global initial_value
    me_vstoxx = market_environment('me_vstoxx', pricing_date)
    me_vstoxx.add_constant('initial_value', initial_value)
    final_date = max(futures.MATURITY).to_pydatetime()
    me_vstoxx.add_constant('final_date', final_date)
    me_vstoxx.add_constant('currency', 'EUR')
    me_vstoxx.add_constant('frequency', 'W')
    me_vstoxx.add_constant('paths', 5000)
    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)
    
    me_vstoxx.add_constant('lambda', 0.5)
    me_vstoxx.add_constant('mu', -0.2)
    me_vstoxx.add_constant('delta', 0.1)
    
    me_vstoxx.add_curve('term_structure', futures)
    
    vstoxx_model = square_root_jump_diffusion_plus('vstoxx_model', me_vstoxx)
    
    vstoxx_model.shift_base = generate_shift_base(pricing_date, futures)
    vstoxx_model.update_shift_values()
    
    
    # option parameters and payoff
    payoff_func = 'np.maximum(maturity_value - strike, 0)'
    
    option_models = {}
    for option in option_selection.index:
        maturity = option_selection['MATURITY'].ix[option]
        me_vstoxx.add_constant('maturity', maturity)
        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
In [14]:
def calculate_model_values(p0, fixed_seed=True):
    ''' Returns all relevant option values.
    
    Parameters
    ===========
    p0 : tuple/list
        tuple of kappa, theta, volatility, lamb, mu, delt
    
    Returns
    =======
    model_values : dict
        dictionary with model values
    '''
    kappa, theta, volatility, lamb, mu, delt = p0
    vstoxx_model.update(kappa=kappa,
                        theta=theta,
                        volatility=volatility,
                        lamb=lamb,
                        mu=mu,
                        delt=delt)
    results = [option_models[option].present_value(
               fixed_seed=fixed_seed)
               for option in option_models]
    model_values = dict(zip(option_models, results))
    return model_values

Model Calibration

In [15]:
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.
        or p0[3] < 0 or p0[4] < 0 or p0[5] < 0):
        return 100
    global option_selection, vstoxx_model, option_models, first, last
    pdate = str(option_selection['DATE'].iloc[0])[:10]
    mat = str(option_selection['MATURITY'].iloc[0])[:10]
    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 first is True:
        penalty = 0.0
    else:
        penalty = (np.sum((p0 - last) ** 2)) / 100
    return MSE + penalty        
In [16]:
import scipy.optimize as spo
def get_parameter_series(pricing_date_list):
    ''' Returns parameter series for the calibrated model over time.
    
    Parameters
    ==========
    pricing_date_list: pd.DatetimeIndex
        object with relevant pricing dates
    maturity_list: list
        list with maturities to be calibrated
    
    Returns
    =======
    parameters: pd.DataFrame
        DataFrame object with parameter series
    '''
    global initial_value, option_selection, vstoxx_model, option_models, first, last
    parameters = pd.DataFrame()
    first = True
    for pricing_date in pricing_date_list:
        initial_value = vstoxx_index['V2TX'][pricing_date]
        option_selection, futures = get_option_selection(pricing_date)
        vstoxx_model, option_models = get_option_models(pricing_date,
                                                        option_selection,
                                                        futures)
        if first is True:
            # global optimization to start with
            opt = spo.brute(mean_squared_error,
                ((0.5, 3.51, 0.75),   # range for kappa
                 (10., 20.1, 2.5),   # range for theta
                 (0.5, 10.51, 5.),  # range for volatility
                 (0.1, 0.71, 0.3),  # range for lambda
                 (0.1, 0.71, 0.3),  # range for mu
                 (0.1, 0.21, 0.1)),  # range for delta
                 finish=None)
        # local optimization
        opt = spo.fmin(mean_squared_error, opt,
                       maxiter=550, maxfun=650,
                       xtol=0.0000001, ftol=0.0000001);
        MSE = mean_squared_error(opt)
        parameters = parameters.append(
                 pd.DataFrame(
                 {'date' : pricing_date,
                  'initial_value' : vstoxx_model.initial_value,
                  'kappa' : opt[0],
                  'theta' : opt[1],
                  'sigma' : opt[2],
                  'lambda' : opt[3],
                  'mu' : opt[4],
                  'delta' : opt[5],
                  'MSE' : MSE},
                  index=[0]), ignore_index=True)
        first = False
        last = opt
        print ("Pricing Date %s" % str(pricing_date)[:10]
               + " | MSE %6.5f" % MSE)
    return parameters
In [17]:
%%time
pricing_date_list = pd.date_range('2014/3/1', '2014/3/31', freq='B')
parameters = get_parameter_series(pricing_date_list)
Optimization terminated successfully.
         Current function value: 0.888388
         Iterations: 139
         Function evaluations: 250
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-03 | MSE 0.01315
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-04 | MSE 0.01742
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-05 | MSE 0.00839
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-06 | MSE 0.00432
Optimization terminated successfully.
         Current function value: 0.520270
         Iterations: 100
         Function evaluations: 196
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-07 | MSE 0.01112
Warning: Maximum number of function evaluations has been exceeded.
Optimization terminated successfully.
         Current function value: 0.003540
         Iterations: 296
         Function evaluations: 563
Pricing Date 2014-03-10 | MSE 0.00354
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-11 | MSE 0.00459
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-12 | MSE 0.00513
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-13 | MSE 0.03272
Optimization terminated successfully.
         Current function value: 0.930028
         Iterations: 126
         Function evaluations: 238
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-14 | MSE 0.00630
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-17 | MSE 0.00857
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-18 | MSE 0.00856
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-19 | MSE 0.00545
Optimization terminated successfully.
         Current function value: 0.020705
         Iterations: 210
         Function evaluations: 372
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-20 | MSE 0.00394
Optimization terminated successfully.
         Current function value: 0.004574
         Iterations: 245
         Function evaluations: 434
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-21 | MSE 0.01374
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-24 | MSE 0.00326
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-25 | MSE 0.00588
Optimization terminated successfully.
         Current function value: 0.009320
         Iterations: 281
         Function evaluations: 493
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-26 | MSE 0.00586
Optimization terminated successfully.
         Current function value: 0.022939
         Iterations: 212
         Function evaluations: 388
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-27 | MSE 0.00847
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-28 | MSE 0.00981
Optimization terminated successfully.
         Current function value: 0.042650
         Iterations: 312
         Function evaluations: 567
Warning: Maximum number of function evaluations has been exceeded.
Pricing Date 2014-03-31 | MSE 0.00711
CPU times: user 12min 4s, sys: 24.5 s, total: 12min 29s
Wall time: 13min 2s

Parameter Time Series

In [18]:
parameters
Out[18]:
MSE date delta initial_value kappa lambda mu sigma theta
0 0.013154 2014-03-03 0.465984 21.8609 2.695352 0.697442 0.439569 1.428634 17.608632
1 0.017423 2014-03-04 0.358975 18.6455 2.636244 0.728487 0.515683 1.469509 17.587893
2 0.008386 2014-03-05 0.384300 18.7004 2.263364 0.728821 0.450962 1.588534 17.575541
3 0.004323 2014-03-06 0.294376 17.8166 2.059080 0.736781 0.507651 1.666605 17.581056
4 0.011118 2014-03-07 0.317709 20.4730 2.384964 0.682157 0.481501 1.598423 17.117888
5 0.003540 2014-03-10 0.316727 20.1237 2.319733 0.718487 0.498948 1.668622 17.077978
6 0.004593 2014-03-11 0.347467 19.7688 2.196480 0.742221 0.451550 1.790583 17.088220
7 0.005125 2014-03-12 0.353123 20.2120 2.298805 0.734982 0.443342 1.917061 17.089660
8 0.032722 2014-03-13 0.363692 22.0896 2.390056 0.728509 0.433385 1.943588 16.468850
9 0.006296 2014-03-14 0.484433 23.2095 2.296335 0.700847 0.328635 2.051268 16.428798
10 0.008571 2014-03-17 0.509279 21.4617 2.271931 0.724860 0.330999 1.980083 16.418118
11 0.008558 2014-03-18 0.502846 20.1089 1.968023 0.742881 0.324256 2.176096 16.584652
12 0.005449 2014-03-19 0.387386 19.2930 1.854628 0.750082 0.383034 2.000749 16.571851
13 0.003936 2014-03-20 0.378975 18.1292 1.740213 0.768905 0.384222 2.132805 16.792499
14 0.013737 2014-03-21 0.389176 17.1916 1.532092 0.772228 0.370253 2.495453 17.053771
15 0.003264 2014-03-24 0.480760 19.4248 1.715106 0.688444 0.285955 2.562484 17.003097
16 0.005885 2014-03-25 0.486194 18.2637 1.705710 0.692390 0.279246 2.683363 17.257420
17 0.005863 2014-03-26 0.439389 17.5869 1.596382 0.700394 0.279168 2.774905 17.465881
18 0.008470 2014-03-27 0.409179 17.6397 1.620665 0.726145 0.265277 3.017647 17.658831
19 0.009806 2014-03-28 0.358049 17.0324 1.605664 0.721946 0.287684 3.331746 17.826602
20 0.007108 2014-03-31 0.342546 17.6639 1.705822 0.725757 0.291088 3.181513 17.844350

Mean-Squared Errors

In [19]:
%matplotlib inline
dat = parameters.MSE
dat.hist(bins=30)
plt.axvline(dat.mean(), color='r', ls='dashed',
                lw=1.5, label='mean = %5.4f' % dat.mean())
plt.legend()
plt.savefig('plots/SRJDplus_Hist_MSE_%s.png' % 
                str(max(parameters.date))[:10])

Model Fit on Selected Pricing Dates

In [21]:
pdate = max(parameters.date)
  # last pricing date
opt = np.array(parameters[parameters.date == pdate][[
    'kappa', 'theta', 'sigma', 'lambda', 'mu', 'delta']])[0]
  # optimal parameters for that date and the maturity
option_selection, futures = get_option_selection(pdate, tol=tol)
vstoxx_model, option_models = get_option_models(pdate, option_selection, futures)
model_values = calculate_model_values(opt, fixed_seed=False)
model_values = pd.DataFrame(model_values.values(), index=model_values.keys(),
                            columns=['MODEL'])
option_selection = option_selection.join(model_values)
print len(option_selection)
mats = set(option_selection.MATURITY.values)
fig, axarr = plt.subplots(len(mats), 2, sharex=True, figsize=(10, 16))
for z, mat in enumerate(mats):
    os = option_selection[option_selection.MATURITY == mat]
    strikes = os.STRIKE.values
    axarr[z, 0].set_ylabel('%s' % str(mat)[:10])
    axarr[z, 0].plot(strikes, os.PRICE.values, label='Market Quotes')
    axarr[z, 0].plot(strikes, os.MODEL.values, 'ro', label='Model Prices')
    axarr[z, 0].grid()
    wi = 0.3
    axarr[z, 1].bar(strikes - wi / 2, os.MODEL.values - os.PRICE.values, width=wi)
    axarr[z, 1].grid()
plt.savefig('plots/SRJDplus_Fit_Last_%s.png' % str(pdate)[:10])
57
Optimization terminated successfully.
         Current function value: 0.042650
         Iterations: 312
         Function evaluations: 567
57