0
$\begingroup$

Context

I'm a beginner quant and I'm trying to calibrate an vol surface using SPX Implied Vol data. The model is from Jim Gatheral and Antoine Jacquier's paper https://www.tandfonline.com/doi/full/10.1080/14697688.2013.819986. I'm testing the model for just one expiry date just for proof of concept. An area where I'm not so confident is the actual implementation of the model. Below is the code I used for the project.

Libraries Used

#numerical libraries
import pandas as pd
from scipy.stats import norm
from scipy.integrate import quad
import scipy.optimize as opt
import numpy as np
import time

#visualisation libraries
import matplotlib.pyplot as plt
import plotly.graph_objects as go 
import seaborn as sb
sb.set()

I use np arrays as much as I can just for speed purposes.

defining functions to use

#ATM Total Variance 
def Theta(ATM_IV,T):
    return ((ATM_IV)**2)*T

#Heston-like smoothing function
def Phi(theta, gamma):
    theta = Theta(ATM_IV, T)
    result = 1/(gamma*theta)*(1-(1-np.exp(-gamma*theta))/(gamma*theta))
    return result

#SSVI
def SSVI(k, ATM_IV, T, rho, gamma):
    theta = Theta(ATM_IV/1, T)
    phi = Phi(theta, gamma)
    result = (0.5 * theta) * (1 + rho * phi * k + np.sqrt((phi * k + rho)**2 + 1 - rho**2))
    return result

def black_scholes_call_price(T, k, S, IV):
    
    # T, k, S, IV = parameters
    r = 0.05  # Fixed interest rate

    #black-scholes call option pricing
    d1 = (np.log(S / k) + (r + ((IV/100)**2) / 2) * T) / ((IV/100) * np.sqrt(T))
    d2 = d1 - (IV/100) * np.sqrt(T)
    call_price = S * norm.cdf(d1) - k * np.exp(-r * T) * norm.cdf(d2)

    return call_price

#minimizes sum of square between implied vols
def objective_function_IV(params, k, ATM_IV, T):
        rho, gamma = params
        Est_IV = SSVI(k, ATM_IV, T, rho, gamma)
        return np.sum((IV - Est_IV)**2)
    
    #minimizes sum of sqaures between option prices
def objective_function_prices(params, k, ATM_IV, T):
    rho, gamma = params
    Est_Total_ImpVol = SSVI(k, ATM_IV, T, rho, gamma)
    Est_IV = np.sqrt(Est_Total_ImpVol/T)
    Est_Prices = black_scholes_call_price(T, S*(M/100), S, Est_IV)
    return np.sum((Benchmark_Prices - Est_Prices)**2)

def increasing_deriv(ATM_IV):
    if np.sum(np.diff(np.unique(ATM_IV))<0) == 0:
        return True
    else:
        return False
#no arbitrage condtion
def Heston_condition(params):
    rho, gamma = params
    return gamma - 0.25*(1.+np.abs(rho))

bnds = [(-1 + 1e-6, 1 - 1e-6), (0.1, 1)]
cons2 = [ {'type': 'ineq', 'fun': Heston_condition} ]

Data

Expiry Moneyness Stock Price IV ATM IV
    0.5 80.0    1202.08 14.8029 13.9493
    0.5 90.0    1202.08 14.8031 13.9493
    0.5 95.0    1202.08 14.4389 13.9493
    0.5 97.5    1202.08 14.1962 13.9493
    0.5 100.0   1202.08 13.9493 13.9493
    0.5 102.5   1202.08 13.7337 13.9493
    0.5 105.0   1202.08 13.5304 13.9493
    0.5 110.0   1202.08 13.1142 13.9493
    0.5 120.0   1202.08 12.9992 13.9493

Implementation

I take the above data frame and unzip each column into an array. For ease of implementation I just pasted the data I was using. I computed the strike and log moneyness to use for two separate calibrations.

T, M, S, IV, ATM_IV = data_transformation(data[[1202.08]]).query("Expiry == 0.5").to_numpy(dtype = 'float').T
k = S*(M/100) #strike prices for black-scholes
lnM = np.log(M/100) #log moneyness
theta = Theta(ATM_IV, T) #total atm variance
Benchmark_Prices = black_scholes_call_price(T, k, S, IV) # BS Option prices using market

Results

Here, when I calibrate using log moneyness and compute the SSVI vols using log moneyness the numbers aren't close to my market implied vols.

#implied vols using log moneyness
results = opt.minimize(objective_function_prices,[0.4 ,0.1], args = (lnM, ATM_IV, T), method = 'COBYLA', bounds = bnds, constraints = cons2 , options={'disp' : True})
print(SSVI(lnM, ATM_IV, T, results.x[0], results.x[1]))
array([97.72859768, 97.49787432, 97.39196288, 97.34107991, 97.29148524,
   97.24311527, 97.19591095, 97.10478361, 96.93433838])

However, when I use the strikes, I get something that looks better, but its going in the wrong direction !

#implied vols using strike
results = opt.minimize(objective_function_prices,[0.4 ,0.1], args = (k, ATM_IV, T), method = 'COBYLA', bounds = bnds, constraints = cons2 , options={'disp' : True})
print(SSVI(k, ATM_IV, T, results.x[0], results.x[1]))
array([11.266532  , 12.56532181, 13.21530702, 13.54041906, 13.86560229,
       14.19085117, 14.51616076, 15.16694458, 16.46906994])

I think the issue may not be with my code, but with my implementation. Though, I could be wrong. Is anyone familiar with the implementation or can see where I'm going wrong in my implementation ? I've been troubleshooting and comparing this to Antoine Jacquier's code https://github.com/JackJacquier/SSVI/blob/master/SSVILocalVol.ipynb but I can't figure out why the numbers don't make sense.

Edit: I found that I was not minimizing between the estimated IV, but the estimated TOTAL IMPLIED VARIANCE. I modified the objective function and got the following estimated implied vols, though they are still off by some.

results = opt.minimize(objective_function_prices,[0.4 ,0.1], args = (lnM, ATM_IV, T), method = 'COBYLA', bounds = bnds, constraints = cons2 , options={'disp' : True})
np.sqrt(SSVI(lnM, ATM_IV, T, results.x[0], results.x[1])/0.5)
array([13.98060068, 13.96408782, 13.9565012 , 13.9528549 , 13.9493    ,
       13.94583201, 13.94244677, 13.93590927, 13.92367325])
Improve this question
$\endgroup$
1

1 Answer 1

Reset to default
0
$\begingroup$

The issue to my problem was a misalignment between the volatility used in my Black Scholes pricing function and the one used for the ATM IV in the SSVI. My BS vol was a decimal (e.g. 0.1345) whereas the ATM IV feed to the SSVI was 13.45. This messed with the dimensions of the model calibration and was the reason my IV estimates were all over the place and even in the wrong direction. To resolve this, I removed the scaling inside the function and scaled the ATM IV before I entered it into the model. Below is a revised implementation of my code as well as an example of a fit.

#ATM Total Variance 
    def Theta(ATM_IV,T):
      return ((ATM_IV)**2)*T

#This parameterization is arbitrage free for the entire set of parameters.
    def SQRT_SSVI(log_k, ATM_IV, T, rho, gamma):
      theta = Theta(ATM_IV, T)
      phi = gamma/(theta**0.5) #sqrt smoothing function from power law
      result = (0.5 * theta) * (1 + rho * phi * log_k + np.sqrt( (phi * log_k + rho)**2 + 1 - rho**2) )
      return result

#prices implied volatilities greater than 0 as calls and puts otherwise
    def black_scholes_pricing_div(T, lnM, k, S, IV, q, r):
      d1 = (np.log(S / k) + (r - q + (0.5 * (IV)**2)) * T) / ((IV) * np.sqrt(T))
      d2 = d1 - (IV) * np.sqrt(T)
      #pricing anything less than strike as puts and calls otherwise
      put_id = np.flatnonzero(lnM <= 0)
      call_id = np.flatnonzero(lnM > 0)

      put_price = k[put_id] * np.exp(-r[put_id] * T[put_id]) * norm.cdf(-d2[put_id]) - S[put_id] * np.exp(-q[put_id] *T[put_id]) * norm.cdf(d1[put_id])
      call_price = S[call_id] * np.exp(-q[call_id] * T[call_id]) * norm.cdf(d1[call_id]) - k[call_id] * np.exp(-r[call_id] * T[call_id]) * norm.cdf(d2[call_id])

      combined_id = np.concatenate([put_id, call_id])
      combined_price = np.concatenate([put_price, call_price])

      # Get the indices that would sort the 'combined_id' array
      sort_indices = np.argsort(combined_id)

# Apply the sorting indices to both 'combined_id' and 'combined_price'
      sorted_combined_id = combined_id[sort_indices]
      sorted_combined_price = combined_price[sort_indices]

# Create a new 2D array where each row is [ID, Price]
      prices = np.column_stack((sorted_combined_id, sorted_combined_price)).round(6)[:,1]

      return prices

The model is then fit minimizing the SSR for market prices and no penalty term was used. The code is also vectorized by using numpy and employs broadcasting to maximize performance.

Attached is the model fit during the GFC. The fit is surprisingly good despite its shortcomings on the front end of the term structure. enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.