from dx import *
import numpy as np
import pandas as pd
from pylab import plt
plt.style.use('seaborn')

"""
Explanation: <img src="http://hilpisch.com/tpq_logo.png" alt="The Python Quants" width="45%" align="right" border="4">
Implied Volatilities and Model Calibration
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.
End of explanation
"""


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()

"""
Explanation: VSTOXX Futures & Options Data
We start by loading VSTOXX data from a pandas HDFStore into DataFrame objects (source: Eurex, cf. http://www.eurexchange.com/advanced-services/).
End of explanation
"""


%matplotlib inline
vstoxx_index['V2TX'].plot(figsize=(10, 6))

"""
Explanation: VSTOXX index for the first quarter of 2014.
End of explanation
"""


vstoxx_futures.info()

vstoxx_futures.tail()

"""
Explanation: The VSTOXX futures data (8 futures maturities/quotes per day).
End of explanation
"""


vstoxx_options.info()

vstoxx_options.tail()

"""
Explanation: The VSTOXX options data. This data set is quite large due to the large number of European put and call options on the VSTOXX.
End of explanation
"""


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)

third_fridays = {}
for month in set(vstoxx_futures['EXP_MONTH']):
    third_fridays[month] = third_friday(dt.datetime(2014, month, 1))

third_fridays

"""
Explanation: As a helper function we need a function to calculate all relevant third Fridays for all relevant maturity months of the data sets.
End of explanation
"""


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

"""
Explanation: Implied Volatilities from Market Quotes
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.
End of explanation
"""


%%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)

"""
Explanation: The following loop now calculates the implied volatilities for all those options whose strike lies within the defined tolerance level.
End of explanation
"""


options_data[60:70]

"""
Explanation: A selection of the results.
End of explanation
"""


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()

"""
Explanation: And the complete results visualized.
End of explanation
"""


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

"""
Explanation: Market Modeling
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:

for a single maturity only
call options only
for a certain moneyness of the options

Relevant Market Data
The following following returns the relevant market data per calibration date:
End of explanation
"""


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

"""
Explanation: Options Modeling
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.
End of explanation
"""


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

"""
Explanation: 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.
End of explanation
"""


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(list(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

"""
Explanation: Calibration Functions
Mean-Squared Error Calculation
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.
End of explanation
"""


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

"""
Explanation: Implementing the Calibration Procedure
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.
End of explanation
"""


%%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)

"""
Explanation: The Calibration Itself
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.
End of explanation
"""


paramet = parameters.set_index('pricing_date')
paramet.tail()

"""
Explanation: Calibration Results
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.
End of explanation
"""


%matplotlib inline
paramet[['kappa', 'theta', 'sigma', 'MSE']].plot(subplots=True, color='b', figsize=(10, 12))
plt.tight_layout()

"""
Explanation: 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.
End of explanation
"""


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(list(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')

"""
Explanation: 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.
End of explanation
"""


