Calibrating an SSVI Volatility Surface with Consistent Volatility Units
Article Quant Q&A · Author: khubquant
Summary
The document describes a one-expiry SPX implied-volatility calibration using the SSVI parameterization. The initial implementation produced implausible estimates when fitting market option prices, and the author also found that the fitting objective had compared implied volatility with total implied variance. A later correction identified a further unit mismatch: Black–Scholes pricing used volatility as a decimal, while the SSVI input used percentage points.
The revised approach aligns volatility units before calibration, computes total variance from at-the-money volatility and expiry, and uses log moneyness as the SSVI input. The described fit minimizes squared differences between market and model option prices, with calls and puts selected according to moneyness. The example reports a fit to historical crisis-period data but notes shortcomings at the front of the term structure. It is a single-expiry implementation example, not evidence that the parameterization or calibration will fit other dates or maturities reliably.
Key ideas
- SSVI uses log moneyness and total at-the-money variance as inputs in the described setup.
- Volatility units must be consistent across the SSVI model and Black–Scholes pricing.
- Comparing implied volatility with total implied variance creates a mismatched calibration objective.
- The revised example fits model parameters by minimizing squared differences in option prices.
- The reported example covers one expiry and retains fit limitations at the front of the term structure.
Tags
Full text
# How to calibrate a volatility surface using SSVI with market data?
# How to calibrate a volatility surface using SSVI with market data?
## 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])
```
## Answer by khubquant (score 1)
https://quant.stackexchange.com/a/78892
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.Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.