####################################################################################################################################
### Program: SimTrade_Implied_Volatility_Surface_Python_code_2026_06_28_SB                                                       ###
### Author: Saral Bindal                                                                                                         ###
### Contact: saralbindal.24@kgpian.iitkgp.ac.in                                                                                  ###
### Article on the SimTrade blog: https://www.simtrade.fr/blog_simtrade/implied-volatility-surface-smiles-smirks-term-structure/ ###                                                           ###
####################################################################################################################################


# Python program to process S&P 500 option chain data and generate the implied volatility surface.

# Objectives of the program:
# 1) Download, filter and clean options data
# 2) Implement the deterministic volatility function
# 3) Plot the Implied Volatility Surfaces, curves and term structures


#################################################################################################
### Organization of the program:                                                              ###
###   STEP 1: Environment setup and package checking                                          ###
###   STEP 2: Downloading the data from yFinance                                              ###
###   STEP 3: Cleaning and Filtering data                                                     ###
###   STEP 4: Implementing Implied Volatility Calculator                                      ###
###   STEP 5: Implementing the Deterministic Volatility Function                              ###
###   STEP 6: Plotting the implied volatility surface                                         ###
###   STEP 7: Plotting the implied volatility curve                                           ###
###   STEP 8: Plotting the implied volatility term structure                                  ###
###   STEP 9: Plotting the Volatility Surface in different stress conditions                  ###
#################################################################################################


#################
# Documentation #
#################

# Python packages
# https://pandas.pydata.org/docs/
# https://numpy.org/doc/stable/
# https://matplotlib.org/stable/
# https://docs.scipy.org/doc/scipy/
# https://ranaroussi.github.io/yfinance/

# Academic articles
# Black, F., & Scholes, M. (1973). The Pricing of Options and Corporate Liabilities. Journal of Political Economy, 81(3), 637-654.
# Merton, R. C. (1973). Theory of Rational Option Pricing. Bell Journal of Economics and Management Science, 4(1), 141-183.
# Dumas, B., Fleming, J., and Whaley, R.E. (1998). Implied volatility functions: Empirical tests. The Journal of Finance, 53(6), 2059-2106.

# Note:
# This code was used to generate the implied volatility surface (IVS) for the
# S&P 500 options data on 2026-06-18. To reproduce the results presented in
# the accompanying article, use the cleaned dataset
# (clean_SP500_options_data_20260618.csv) provided with the article and run
# the code directly from STEP 4 onwards. To generate an IVS for a different
# trading day, modify the dates and rerun the entire program.

# Important:
# If using the program for a different trading day, check and update:
# 1) Dates used throughout the program
# 2) Spot price (S), risk-free rate (r) and dividend yield (q), if using fixed values
# 3) Time-to-maturity (T) and log-moneyness (M) values used in STEPS 7 and 8


#################################################
# STEP 1: Environment setup and package loading #        
#################################################

import sys
import subprocess

# Automatically install missing Python packages
def install(package):
    subprocess.check_call([sys.executable, "-m", "pip", "install", package])


# Core package for working with directories and files
import os

# Package for numerical calculations
try:
    import numpy as np
except ImportError:
    install("numpy")
    import numpy as np


# Package for data manipulation and dataframes
try:
    import pandas as pd
except ImportError:
    install("pandas")
    import pandas as pd


# Package for creating charts and plots
try:
    import matplotlib.pyplot as plt
    from matplotlib import cm
except ImportError:
    install("matplotlib")
    import matplotlib.pyplot as plt
    from matplotlib import cm


# Package for downloading financial market data
try:
    import yfinance as yf
except ImportError:
    install("yfinance")
    import yfinance as yf


# Package for scientific computing and optimization
try:
    import scipy
    from scipy.stats import norm
    from scipy.optimize import brentq
except ImportError:
    install("scipy")
    import scipy
    from scipy.stats import norm
    from scipy.optimize import brentq


# Standard Python libraries for working with dates and mathematical functions
from datetime import datetime
import math


#################################################
# STEP 2: Downloading the data from yFinance    #        
#################################################

# Recording the current spot price (2026-06-18)
sp500 = yf.Ticker("^GSPC")
sp500_today = sp500.history(period="1d")
#S = sp500_today['Close'].iloc[-1]
S = 7496.0400390625
print(S)

# Recording the current risk-free rate (2026-06-18)
treasury_3m = yf.Ticker("^IRX")
treasury_today = treasury_3m.history(period="1d")
interest_rate_raw = treasury_today['Close'].iloc[-1]
#r = interest_rate_raw / 100
r = 0.036580002307891844
print(r)

# Recording the current dividend yield rate (2026-06-18)
q = 0.0104
print(q)

# Importing all the options data from yfinance for all the expiry dates available (starting from 2026-06-18)
tk = yf.Ticker('^SPX')

expirations = tk.options

data_import = pd.DataFrame()

for expiry in expirations:
    options = tk.option_chain(expiry)
    calls = options.calls
    puts = options.puts
    calls['optionType'] = 'C'
    puts['optionType'] = 'P'
    exp_data = pd.concat(objs=[calls, puts], ignore_index=True)
    exp_data['expiryDate'] = pd.to_datetime(expiry)
    data_import = pd.concat(objs=[data_import, exp_data], ignore_index=True)

data_import.drop(columns=['contractSymbol', 'contractSize', 'currency', 'lastTradeDate', 'lastPrice', 'change', 'percentChange'], inplace=True)
data_import.to_csv('SP500_options_data_20260618.csv', index = False)


#################################################
# STEP 3: Cleaning and Filtering data           #        
#################################################

data = pd.read_csv('SP500_options_data_20260618.csv', parse_dates=['expiryDate'])

# Filtering for OTM call and put options
data = data[((data['optionType'] == 'P') & (data['strike'] <= S)) | ((data['optionType'] == 'C') & (data['strike'] >= S))]

# Filtering out illiquid options with volume = 0 and open interest = 0
data = data[(data['volume'] != 0) | (data['openInterest'] != 0)]

# Creating the time-to-maturity column
today = pd.Timestamp('2026-06-18')
data['timeToMaturity'] = (data['expiryDate'] - today).dt.days/365.25
data = data[data['timeToMaturity']>0.01]

# Calculating the mid price
data['mid'] = (data['bid'] + data['ask'])/2

# Creating the risk-free discount factor
data['r_df'] = np.exp(-r*data['timeToMaturity'])

# Creating the dividend yield discount factor
data['q_df'] = np.exp(-q*data['timeToMaturity'])

# Applying the upper and lower bounds for the call and put options prices
data['upper_bound'] = 0.0
data['lower_bound'] = 0.0

data.loc[data['optionType'] == 'C', 'upper_bound'] = S * data['q_df']
data.loc[data['optionType'] == 'C', 'lower_bound'] = np.maximum(0, (S * data['q_df']) - (data['strike'] * data['r_df']))

data.loc[data['optionType'] == 'P', 'upper_bound'] = data['strike'] * data['r_df']
data.loc[data['optionType'] == 'P', 'lower_bound'] = np.maximum(0, (data['strike'] * data['r_df']) - (S * data['q_df']))

data = data[(data['mid'] < data['upper_bound']) & (data['mid'] > data['lower_bound'])]

# Sorting the values by time to maturity and within it by the strik price in ascending order
data = data.sort_values(by=['timeToMaturity', 'strike'], ascending=True)

# Filtering for option contracts not following calendar arbitrage for every expiration date
data['price_change'] = data.groupby('timeToMaturity')['mid'].diff()
data = data[((data['price_change'] > 0) & (data['optionType'] == 'P')) | 
            ((data['price_change'] < 0) & (data['optionType'] == 'C')) | 
            (data['price_change'].isna())]

data.to_csv('clean_SP500_options_data_20260618.csv', index = False)


######################################################
# STEP 4: Implementing Implied Volatility Calculator #        
######################################################

data_clean = pd.read_csv('clean_SP500_options_data_20260618.csv', parse_dates=['expiryDate'])

# Calculating the forward price
data_clean['forwardPrice'] = S * np.exp(data_clean['timeToMaturity']*(r-q))

# Calculating the log-moneyness
data_clean['log_moneyness'] = np.log(data_clean['strike']/data_clean['forwardPrice'])

# Removing contracts that are long-dated (> 1 year) as they are higly illiquid also called as LEAPs
data_clean = data_clean[data_clean['timeToMaturity'] <= 1]

# Removing deep OTM contracts as they have very less trading volume
data_clean = data_clean[(data_clean['log_moneyness'] >= -0.45) & (data_clean['log_moneyness'] <= 0.45)]

# Computing our own implied volatility using the BSM option pricer and the Brent's root finding method

def BSM_price_discounted(otype: str, S: float, K: float, T: float, r_df: float, q_df: float, vol: float):
    if vol <= 0 or T <= 0:
        return 0.0
    forward_moneyness = (S * q_df) / (K * r_df)
    d1 = (math.log(forward_moneyness) + 0.5 * vol**2 * T) / (vol * math.sqrt(T))
    d2 = d1 - (vol * math.sqrt(T))
    
    if otype == 'C':
        return (S * q_df * norm.cdf(d1)) - (K * r_df * norm.cdf(d2))
    elif otype == 'P':
        return (K * r_df * norm.cdf(-d2)) - (S * q_df * norm.cdf(-d1))
    return 0.0

def IV_calculator_brentq(otype: str, S: float, K: float, T: float, r_df: float, q_df: float, mid: float):
    # Set target market price 
    target_price = mid
    
    # Establish lower and upper search brackets
    vol_lower = 0.001
    vol_upper = 5.0
    
    try:
        # One-line lambda execution: (Theoretical Price) - (Target Market Price)
        custom_iv = brentq(
            lambda vol: BSM_price_discounted(otype, S, K, T, r_df, q_df, vol) - target_price, 
            vol_lower, 
            vol_upper
        )
        return custom_iv
        
    except ValueError:
        # Handle cases where boundary signs do not flip (arbitrage violations)
        return np.nan
    
# Execute the IV calculator row-by-row
data_clean['custom_iv'] = data_clean.apply(
    lambda row: IV_calculator_brentq(
        otype=row['optionType'],       
        S=S, 
        K=row['strike'], 
        T=row['timeToMaturity'],       
        r_df=row['r_df'], 
        q_df=row['q_df'], 
        mid=row['mid']
    ), 
    axis=1
)


###############################################################
# STEP 5: Implementing the Deterministic Volatitliy Function  #        
###############################################################

# Generate parametric polynomial features for Implied Volatility Surface (IVS) modeling
design_matrix = pd.DataFrame({
    'intercept': 1.0,
    'M': data_clean['log_moneyness'],
    'M2': data_clean['log_moneyness'] ** 2,
    'τ': data_clean['timeToMaturity'],
    'τ2': data_clean['timeToMaturity'] ** 2,
    'M_τ': data_clean['log_moneyness'] * data_clean['timeToMaturity'],
    'target_iv': data_clean['custom_iv']
})

# Calculating the parameters
X = design_matrix[['intercept', 'M', 'M2', 'τ', 'τ2', 'M_τ']]
Y = design_matrix['target_iv']

alpha_0, alpha_1, alpha_2, alpha_3, alpha_4, alpha_5 = np.linalg.lstsq(X, Y)[0]

print('beta_0:', alpha_0)
print('beta_1:', alpha_1)
print('beta_2:', alpha_2)
print('beta_3:', alpha_3)
print('beta_4:', alpha_4)
print('beta_5:', alpha_5)


###############################################################
# STEP 6: Plotting the implied volatltiy surface (2026-06-18) #        
###############################################################

m_min, m_max = design_matrix['M'].min(), design_matrix['M'].max()
t_min, t_max = design_matrix['τ'].min(), design_matrix['τ'].max()

M_grid = np.linspace(m_min, m_max, 10)
T_grid = np.linspace(t_min, t_max, 10)

M_mesh, T_mesh = np.meshgrid(M_grid, T_grid)

IV_surface = (
    alpha_0 + 
    alpha_1 * M_mesh + 
    alpha_2 * (M_mesh**2) + 
    alpha_3 * T_mesh + 
    alpha_4 * (T_mesh**2) + 
    alpha_5 * (M_mesh * T_mesh)
)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

surf = ax.plot_surface(M_mesh, T_mesh, IV_surface, 
                        cmap=cm.viridis,        
                        edgecolor='black',  
                        linewidth=0.3,      
                        antialiased=True)

ax.set_xlim(m_min, m_max)
ax.set_ylim(0, 1)
ax.set_zlim(0.0, float(np.nanmax(IV_surface)) * 1.1)

ax.set_xlabel('Log-Moneyness (ln(K/F))', fontsize = 18, labelpad = 17)
ax.set_ylabel('Time to Maturity (T)', fontsize = 18, labelpad = 17)
ax.set_zlabel('Implied Volatility (%)', fontsize = 18, labelpad = 17)
ax.set_title('Implied Volatility Surface of the S&P 500 Options (2026-06-18)', fontsize = 20, pad=15)
ax.tick_params(labelsize=13)

ax.view_init(elev=20, azim=-70, roll=0) 

plt.tight_layout()
plt.show()


########################################################################################
# STEP 7: Plotting the implied volatltiy curve (2026-06-18) for 3 different maturities #        
########################################################################################

print(design_matrix['τ'].unique())

T_1 = 0.02190281
m_min_1, m_max_1 = design_matrix.loc[np.isclose(design_matrix['τ'], T_1 ),'M'].min(), design_matrix.loc[np.isclose(design_matrix['τ'], T_1 ),'M'].max()
M_grid_1 = np.linspace(m_min_1, 0.2, 40)
IV_curve_1 = (alpha_0 + alpha_1 * M_grid_1 + alpha_2 * (M_grid_1**2) + alpha_3 * T_1 + alpha_4 * (T_1**2) + alpha_5 * (M_grid_1 * T_1))

T_2 = 0.42436687
m_min_2, m_max_2 = design_matrix.loc[np.isclose(design_matrix['τ'], T_2 ),'M'].min(), design_matrix.loc[np.isclose(design_matrix['τ'], T_2 ),'M'].max()
M_grid_2 = np.linspace(m_min_2, m_max_2, 40)
IV_curve_2 = (alpha_0 + alpha_1 * M_grid_2 + alpha_2 * (M_grid_2**2) + alpha_3 * T_2 + alpha_4 * (T_2**2) + alpha_5 * (M_grid_2 * T_2))

T_3 = 0.82683094
m_min_3, m_max_3 = design_matrix.loc[np.isclose(design_matrix['τ'], T_3 ),'M'].min(), design_matrix.loc[np.isclose(design_matrix['τ'], T_3 ),'M'].max()
M_grid_3 = np.linspace(m_min_3, m_max_3, 40)
IV_curve_3 = (alpha_0 + alpha_1 * M_grid_3 + alpha_2 * (M_grid_3**2) + alpha_3 * T_3 + alpha_4 * (T_3**2) + alpha_5 * (M_grid_3 * T_3))


fig, ax = plt.subplots()
ax.plot(M_grid_1, IV_curve_1*100, color='#0A2540', lw=2)
ax.plot(M_grid_2, IV_curve_2*100, color='#20639B', lw=2)
ax.plot(M_grid_3, IV_curve_3*100, color='#6BAED6', lw=2)
ax.grid(True, linestyle=':', alpha=0.6)
ax.legend(['T = 0.02 Y', 'T = 0.42 Y', 'T = 0.83 Y'], frameon=True, fontsize = 18)
ax.set_xlabel('Log-Moneyness (ln(K/F))', fontsize = 18, labelpad = 12)
ax.set_ylabel('Implied Volatility (%)', fontsize = 18, labelpad = 12)
ax.set_title('Implied Volatility Curves of the S&P 500 Options (2026-06-18)', fontsize = 20, pad=15)
ax.tick_params(labelsize=15)
ax.set_xlim(-0.4,0.2)
plt.show()


###################################################################################################
# STEP 8: Plotting the implied volatltiy term structure (2026-06-18) for ATM, OTM and ITM options #        
###################################################################################################

print(design_matrix['M'].unique())

M_1 = 0
t_min_1, t_max_1 = design_matrix.loc[np.isclose(design_matrix['M'], M_1, atol = 0.01 ),'τ'].min(), design_matrix.loc[np.isclose(design_matrix['M'], M_1 , atol = 0.01),'τ'].max()
t_grid_1 = np.linspace(t_min_1, t_max_1, 40)
IV_ts_1 = (alpha_0 + alpha_1 * M_1 + alpha_2 * (M_1**2) + alpha_3 * t_grid_1 + alpha_4 * (t_grid_1**2) + alpha_5 * (t_grid_1 * M_1))

M_2 = 0.2
t_min_2, t_max_2 = design_matrix.loc[np.isclose(design_matrix['M'], M_2, atol = 0.01 ),'τ'].min(), design_matrix.loc[np.isclose(design_matrix['M'], M_2 , atol = 0.01),'τ'].max()
t_grid_2 = np.linspace(t_min_2, t_max_2, 40)
IV_ts_2 = (alpha_0 + alpha_1 * M_2 + alpha_2 * (M_2**2) + alpha_3 * t_grid_2 + alpha_4 * (t_grid_2**2) + alpha_5 * (t_grid_2 * M_2))

M_3 = -0.2
t_min_3, t_max_3 = design_matrix.loc[np.isclose(design_matrix['M'], M_3, atol = 0.01 ),'τ'].min(), design_matrix.loc[np.isclose(design_matrix['M'], M_3 , atol = 0.01),'τ'].max()
t_grid_3 = np.linspace(t_min_3, t_max_3, 40)
IV_ts_3 = (alpha_0 + alpha_1 * M_3 + alpha_2 * (M_3**2) + alpha_3 * t_grid_3 + alpha_4 * (t_grid_3**2) + alpha_5 * (t_grid_3 * M_3))

fig, ax = plt.subplots()
ax.plot(t_grid_1, IV_ts_1*100, color='#0A2540', lw=2)
ax.plot(t_grid_2, IV_ts_2*100, color='#20639B', lw=2)
ax.plot(t_grid_3, IV_ts_3*100, color='#6BAED6', lw=2)
ax.grid(True, linestyle=':', alpha=0.6)
ax.legend(['M = 0.0 (ATM)', 'M = 0.2 (OTM Call)', 'M = -0.2 (OTM Put)'], frameon=True, fontsize = 18)
ax.set_xlabel('Time to Maturity (T)', fontsize = 18, labelpad = 12)
ax.set_ylabel('Implied Volatility (%)', fontsize = 18, labelpad = 12)
ax.set_title('Implied Volatility Term Structure of the S&P 500 Options (2026-06-18)', fontsize = 20, pad=15)
ax.tick_params(labelsize=15)
plt.show()


#########################################################################
# STEP 9: Plotting the Volatltiy Surface in different stress conditions #        
#########################################################################

# Note: the values are not based on real data, just an assumption of possible changes

# Plot 1 (Covid-19 Pandemic)

alphas_baseline = [alpha_0, alpha_1, alpha_2, alpha_3, alpha_4, alpha_5]

m_min, m_max = design_matrix['M'].min(), design_matrix['M'].max()
t_min, t_max = design_matrix['τ'].min(), design_matrix['τ'].max()

M_grid = np.linspace(m_min, m_max, 10)
T_grid = np.linspace(t_min, t_max, 10)

M_mesh, T_mesh = np.meshgrid(M_grid, T_grid)

alphas_covid = alphas_baseline.copy()
alphas_covid[0] += 0.62   
alphas_covid[1] *= 0.90  
alphas_covid[2] += 0.15   
alphas_covid[3] -= 0.53

IV_surface = (
    alpha_0 + 
    alpha_1 * M_mesh + 
    alpha_2 * (M_mesh**2) + 
    alpha_3 * T_mesh + 
    alpha_4 * (T_mesh**2) + 
    alpha_5 * (M_mesh * T_mesh)
)

IV_covid = (
    alphas_covid[0] + 
    alphas_covid[1] * M_mesh + 
    alphas_covid[2] * (M_mesh**2) + 
    alphas_covid[3] * T_mesh + 
    alphas_covid[4] * (T_mesh**2) + 
    alphas_covid[5] * (M_mesh * T_mesh)
)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(M_mesh, T_mesh, IV_covid, 
                        cmap=cm.viridis,        
                        edgecolor='black',  
                        linewidth=0.3,      
                        antialiased=True, label = 'During the Covid-19 crash (2020)')

ax.plot_surface(M_mesh, T_mesh, IV_surface, 
                        cmap=cm.viridis,        
                        edgecolor='black',  
                        linewidth=0.3,      
                        antialiased=True,
                        alpha = 0.1, label = 'Baseline (2026-06-18)')

ax.set_xlim(m_min, m_max)
ax.set_ylim(0, 1)
ax.set_zlim(0.0, float(np.nanmax(IV_covid)) * 1.1)

ax.legend(frameon=True, fontsize = 15)

ax.set_xlabel('Log-Moneyness (ln(K/F))', fontsize = 18, labelpad = 17)
ax.set_ylabel('Time to Maturity (T)', fontsize = 18, labelpad = 17)
ax.set_zlabel('Implied Volatility (%)', fontsize = 18, labelpad = 17)
ax.set_title('Implied Volatility Surface of the S&P 500 Options (2026-06-18)', fontsize = 20, pad=15)
ax.tick_params(labelsize=13)

ax.view_init(elev=20, azim=-70, roll=0) 

plt.tight_layout()
plt.show()


# Plot 2 (Escalation of US-Iran war)

alphas_baseline = [alpha_0, alpha_1, alpha_2, alpha_3, alpha_4, alpha_5]

m_min, m_max = design_matrix['M'].min(), design_matrix['M'].max()
t_min, t_max = design_matrix['τ'].min(), design_matrix['τ'].max()

M_grid = np.linspace(m_min, m_max, 10)
T_grid = np.linspace(t_min, t_max, 10)

M_mesh, T_mesh = np.meshgrid(M_grid, T_grid)

alphas_future_escalation = alphas_baseline.copy()
alphas_future_escalation[0] += 0.25   # Moderate-to-high persistent lift
alphas_future_escalation[1] *= 1.20   # Steep downside skew (panic put-buying for earnings protection)
alphas_future_escalation[2] += 0.08   # Asymmetric smirk curvature
alphas_future_escalation[3] += 0.30   # Flat term structure (fear stays high across long horizons)
alphas_future_escalation[5] -= 0.25

IV_surface = (
    alpha_0 + 
    alpha_1 * M_mesh + 
    alpha_2 * (M_mesh**2) + 
    alpha_3 * T_mesh + 
    alpha_4 * (T_mesh**2) + 
    alpha_5 * (M_mesh * T_mesh)
)

IV_future_escalation = (
    alphas_future_escalation[0] + 
    alphas_future_escalation[1] * M_mesh + 
    alphas_future_escalation[2] * (M_mesh**2) + 
    alphas_future_escalation[3] * T_mesh + 
    alphas_future_escalation[4] * (T_mesh**2) + 
    alphas_future_escalation[5] * (M_mesh * T_mesh)
)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')


ax.plot_surface(M_mesh, T_mesh, IV_future_escalation, 
                        cmap=cm.viridis,        
                        edgecolor='black',  
                        linewidth=0.3,      
                        antialiased=True, label = 'Future war escalation')

ax.plot_surface(M_mesh, T_mesh, IV_surface, 
                        cmap=cm.viridis,        
                        edgecolor='black',  
                        linewidth=0.3,      
                        antialiased=True,
                        alpha = 0.1, label = 'Baseline (2026-06-18)')

ax.set_xlim(m_min, m_max)
ax.set_ylim(0, 1)
ax.set_zlim(0.0, float(np.nanmax(IV_covid)) * 1.1)
ax.legend(frameon=True, fontsize = 15)

ax.set_xlabel('Log-Moneyness (ln(K/F))', fontsize = 18, labelpad = 17)
ax.set_ylabel('Time to Maturity (T)', fontsize = 18, labelpad = 17)
ax.set_zlabel('Implied Volatility (%)', fontsize = 18, labelpad = 17)
ax.set_title('Implied Volatility Surface of the S&P 500 Options (2026-06-18)', fontsize = 20, pad=15)
ax.tick_params(labelsize=13)

ax.view_init(elev=20, azim=-70, roll=0) 

plt.tight_layout()
plt.show()


# End of the code