Python

Calibration of Heston93 Stochastic Volatility Model

# Calibration of H93 Stochastic Volatility Model

# Calibration of Stoch Vol Jump Model to EURO STOXX Option Quotes via Numerical Integration
#Bakshi, Cao and Chen (1997)
# Data Source: www.eurexchange.com 

import sys
sys.path.append('09_gmm')
import math
import numpy as np
np.set_printoptions(suppress=True,
formatter={'all': lambda x: '%5.3f' % x})
import pandas as pd
from scipy.optimize import brute, fmin, minimize
import matplotlib as mpl
mpl.rcParams['font.family'] = 'serif'
from BCC_option_valuation import H93_call_value
from CIR_calibration import CIR_calibration, r_list
from CIR_zcb_valuation import B
#
# Calibrate Short Rate Model
#
kappa_r, theta_r, sigma_r = CIR_calibration()
#
# Market Data from www.eurexchange.com
# as of 30. September 2014
#
h5 = pd.HDFStore('08_m76/option_data.h5', 'r')
data = h5['data'] # European call & put option data (3 maturities)
h5.close()
S0 = 3225.93 # EURO STOXX 50 level 30.09.2014
r0 = r_list[0] # initial short rate

#
# Option Selection
#
tol = 0.02 # percent ITM/OTM options
options = data[(np.abs(data['Strike'] - S0) / S0) < tol]
# options = data[data['Strike'].isin([3100, 3150, 3225, 3300, 3350])]
#
# Adding Time-to-Maturity and Short Rates
#
for row, option in options.iterrows():
    T = (option['Maturity'] - option['Date']).days / 365.
    options.ix[row, 'T'] = T
    B0T = B([kappa_r, theta_r, sigma_r, r0, T])
    options.ix[row, 'r'] = -math.log(B0T) / T
#
# Calibration Functions
#
i = 0
min_MSE = 500
def H93_error_function(p0):
    ''' Error function for parameter calibration in BCC97 model via
    Lewis (2001) Fourier approach.
    Parameters
    ==========
    kappa_v: float
    mean-reversion factor
    theta_v: float
    long-run mean of variance
    sigma_v: float
    volatility of variance
    rho: float
    correlation between variance and stock/index level
    v0: float
    initial, instantaneous variance
    Returns
    =======
    MSE: float
    mean squared error
    '''
    global i, min_MSE
    kappa_v, theta_v, sigma_v, rho, v0 = p0
    
    if kappa_v < 0.0 or theta_v < 0.005 or sigma_v < 0.0 or \
                        rho < -1.0 or rho > 1.0:
        return 500.0
    if 2 * kappa_v * theta_v < sigma_v ** 2:
        return 500.0
    se = []
    for row, option in options.iterrows():
        model_value = H93_call_value(S0, option['Strike'], option['T'],
                                     option['r'], kappa_v, theta_v, sigma_v, rho, v0)
        se.append((model_value - option['Call']) ** 2)
    MSE = sum(se) / len(se)
    min_MSE = min(min_MSE, MSE)
    if i % 25 == 0:
        print '%4d |' % i, np.array(p0), '| %7.3f | %7.3f' % (MSE, min_MSE)
    i += 1
    return MSE

def H93_calibration_full():
    ''' Calibrates H93 stochastic volatility model to market quotes. '''
    # first run with brute force
    # (scan sensible regions)
    p0 = brute(H93_error_function,
                ((2.5, 10.6, 5.0), # kappa_v
                (0.01, 0.041, 0.01), # theta_v
                (0.05, 0.251, 0.1), # sigma_v
                (-0.75, 0.01, 0.25), # rho
                (0.01, 0.031, 0.01)), # v0
                finish=None)
    # second run with local, convex minimization
    # (dig deeper where promising)
    opt = fmin(H93_error_function, p0,
               xtol=0.000001, ftol=0.000001,
               maxiter=750, maxfun=900)
    np.save('11_cal/opt_sv', np.array(opt))
    return opt

def H93_calculate_model_values(p0):
    ''' Calculates all model values given parameter vector p0. '''
    kappa_v, theta_v, sigma_v, rho, v0 = p0
    values = []
    for row, option in options.iterrows():
        model_value = H93_call_value(S0, option['Strike'], option['T'],
                                     option['r'], kappa_v, theta_v, sigma_v, rho, v0)
        values.append(model_value)
    return np.array(values)