####################################################################################################################################
### Program: SimTrade_Implied_Volatility_Surface_R_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/ ###                                                                                            ###
####################################################################################################################################


# R 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 #
#################

# R packages
# https://www.rdocumentation.org/packages/quantmod/versions/0.4.26
# https://www.rdocumentation.org/packages/tidyverse/versions/2.0.0
# https://www.rdocumentation.org/packages/plot3D/versions/1.1

# 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 #        
#################################################

# Modify the path below to match the folder where your script and CSV files are stored
setwd("C:/Users/YourName/Documents/MyTradingProject")

install_and_load <- function(package) {
  if (!require(package, character.only = TRUE)) {
    install.packages(package, dependencies = TRUE, repos = "https://cloud.r-project.org")
    library(package, character.only = TRUE)
  }
}

# Quantitative financial modelling framework for loading equity quotes and option chain structures
install_and_load("quantmod")

# Standard toolkit library for data manipulation, cleaning, filtering and structural transformations
install_and_load("tidyverse")

# Specialized multi-dimensional visualization extension for rendering complex volumetric mesh projections
install_and_load("plot3D")


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

# Recording the current spot price (2026-06-18)
sp500_data <- getSymbols("^GSPC", src = "yahoo", from = "2026-06-17", to = "2026-06-19", auto.assign = FALSE)
#S = as.numeric(Cl(tail(sp500_data, 1)))
S <- 7496.0400390625
print(S)

# Recording the current risk-free rate (2026-06-18)
treasury_3m <- getSymbols("^IRX", src = "yahoo", from = "2026-06-17", to = "2026-06-19", auto.assign = FALSE)
interest_rate_raw <- as.numeric(Cl(tail(treasury_3m, 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 Yahoo Finance for all the expiry dates available (starting from 2026-06-18)
options_raw <- getOptionChain("^SPX", Exp = NULL, src = "yahoo")

data_import <- data.frame()

for (expiry in names(options_raw)) {
  calls <- options_raw[[expiry]][["calls"]]
  puts <- options_raw[[expiry]][["puts"]]
  
  if (!is.null(calls) && nrow(calls) > 0) {
    calls$optionType <- 'C'
  }
  if (!is.null(puts) && nrow(puts) > 0) {
    puts$optionType <- 'P'
  }
  
  exp_data <- bind_rows(calls, puts)
  if (!is.null(exp_data) && nrow(exp_data) > 0) {
    # quantmod outputs expiry names in MM.DD.YYYY string formats
    exp_data$expiryDate <- as.Date(expiry, format = "%b.%d.%Y")
    data_import <- bind_rows(data_import, exp_data)
  }
}

# Standardizing column headers mapping quantmod names to match python schemas
data_import <- data_import %>%
  rename(strike = Strike, bid = Bid, ask = Ask, volume = Vol, openInterest = OI)

# Dropping unnecessary fields as requested by structural workflow
data_import <- data_import %>% 
  select(-one_of('Chg', 'PctChg', 'IV', 'Last', 'Symbol'))

write.csv(data_import, 'SP500_options_data_20260618.csv', row.names = FALSE)


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

data <- read.csv('SP500_options_data_20260618.csv')
data$expiryDate <- as.Date(data$expiryDate)

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

# Filtering out illiquid options with volume = 0 and open interest = 0
data <- data %>% filter((volume != 0) | (openInterest != 0))

# Creating the time-to-maturity column
today <- as.Date('2026-06-18')
data$timeToMaturity <- as.numeric(data$expiryDate - today) / 365.25
data <- data %>% filter(timeToMaturity > 0.01)

# Calculating the mid price
data$mid <- (data$bid + data$ask) / 2

# Creating the risk-free discount factor
data$r_df <- exp(-r * data$timeToMaturity)

# Creating the dividend yield discount factor
data$q_df <- 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 <- data %>% mutate(
  upper_bound = case_when(
    optionType == 'C' ~ S * q_df,
    optionType == 'P' ~ strike * r_df,
    TRUE ~ upper_bound
  ),
  lower_bound = case_when(
    optionType == 'C' ~ pmax(0, (S * q_df) - (strike * r_df)),
    optionType == 'P' ~ pmax(0, (strike * r_df) - (S * q_df)),
    TRUE ~ lower_bound
  )
)

data <- data %>% filter((mid < upper_bound) & (mid > lower_bound))

# Sorting the values by time to maturity and within it by the strik price in ascending order
data <- data %>% arrange(timeToMaturity, strike)

# Filtering for option contracts not following calendar arbitrage for every expiration date
data <- data %>% 
  group_by(timeToMaturity) %>% 
  mutate(price_change = c(NA, diff(mid))) %>% 
  ungroup()

data <- data %>% filter(
  ((price_change > 0) & (optionType == 'P')) | 
    ((price_change < 0) & (optionType == 'C')) | 
    is.na(price_change)
)

write.csv(data, 'clean_SP500_options_data_20260618.csv', row.names = FALSE)


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

data_clean <- read.csv('clean_SP500_options_data_20260618.csv')
data_clean$expiryDate <- as.Date(data_clean$expiryDate)

# Calculating the forward price
data_clean$forwardPrice <- S * exp(data_clean$timeToMaturity * (r - q))

# Calculating the log-moneyness
data_clean$log_moneyness <- 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 %>% filter(timeToMaturity <= 1)

# Removing deep OTM contracts as they have very less trading volume
data_clean <- data_clean %>% filter(log_moneyness >= -0.45, log_moneyness <= 0.45)

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

BSM_price_discounted <- function(otype, S, K, T, r_df, q_df, vol) {
  if (vol <= 0 || T <= 0) {
    return(0.0)
  }
  forward_moneyness <- (S * q_df) / (K * r_df)
  d1 <- (log(forward_moneyness) + 0.5 * vol^2 * T) / (vol * sqrt(T))
  d2 <- d1 - (vol * sqrt(T))
  
  if (otype == 'C') {
    return((S * q_df * pnorm(d1)) - (K * r_df * pnorm(d2)))
  } else if (otype == 'P') {
    return((K * r_df * pnorm(-d2)) - (S * q_df * pnorm(-d1)))
  }
  return(0.0)
}

IV_calculator_brentq <- function(otype, S, K, T, r_df, q_df, mid) {
  target_price <- mid
  vol_lower <- 0.001
  vol_upper <- 5.0
  
  obj_func <- function(vol) {
    BSM_price_discounted(otype, S, K, T, r_df, q_df, vol) - target_price
  }
  
  res <- tryCatch({
    uniroot(obj_func, lower = vol_lower, upper = vol_upper)$root
  }, error = function(e) {
    return(NA)
  })
  return(res)
}

# Execute the IV calculator row-by-row
data_clean$custom_iv <- sapply(1:nrow(data_clean), function(i) {
  IV_calculator_brentq(
    otype = data_clean$optionType[i],
    S = S,
    K = data_clean$strike[i],
    T = data_clean$timeToMaturity[i],
    r_df = data_clean$r_df[i],
    q_df = data_clean$q_df[i],
    mid = data_clean$mid[i]
  )
})


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

# Generate parametric polynomial features for Implied Volatility Surface (IVS) modeling
design_matrix <- data.frame(
  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 <- as.matrix(design_matrix[, c('intercept', 'M', 'M2', 'τ', 'τ2', 'M_τ')])
Y <- design_matrix$target_iv

# OLS execution via standard matrix projection mapping
coefficients_ols <- solve(t(X) %*% X) %*% t(X) %*% Y
alpha_0 <- coefficients_ols[1]
alpha_1 <- coefficients_ols[2]
alpha_2 <- coefficients_ols[3]
alpha_3 <- coefficients_ols[4]
alpha_4 <- coefficients_ols[5]
alpha_5 <- coefficients_ols[6]

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


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

m_min <- min(design_matrix$M); m_max <- max(design_matrix$M)
t_min <- min(design_matrix$τ); t_max <- max(design_matrix$τ)

M_grid <- seq(m_min, m_max, length.out = 10)
T_grid <- seq(t_min, t_max, length.out = 10)

M_mesh <- outer(M_grid, rep(1, length(T_grid)))
T_mesh <- outer(rep(1, length(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)
)

z_surface <- IV_surface

persp3D(x = M_grid, y = T_grid, z = z_surface,
        xlim = c(m_min, m_max), ylim = c(0, 1), zlim = c(0.0, max(z_surface, na.rm = TRUE) * 1.1),
        xlab = 'Log-Moneyness (ln(K/F))', ylab = 'Time to Maturity (T)', zlab = 'Implied Volatility (%)',
        main = 'Implied Volatility Surface of the S&P 500 Options (2026-06-18)',
        theta = -70, phi = 20, expand = 0.5, colkey = FALSE,
        bty = "b2", col = viridis(100), border = "black", lwd = 0.3)


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

print(unique(design_matrix$τ))

# Helper numeric precision function matching Python approach 
isclose <- function(a, b, atol = 1e-05) { abs(a - b) <= atol }

T_1 <- 0.02190281
m_min_1 <- min(design_matrix$M[isclose(design_matrix$τ, T_1)])
m_max_1 <- max(design_matrix$M[isclose(design_matrix$τ, T_1)])
M_grid_1 <- seq(m_min_1, 0.2, length.out = 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 <- min(design_matrix$M[isclose(design_matrix$τ, T_2)])
m_max_2 <- max(design_matrix$M[isclose(design_matrix$τ, T_2)])
M_grid_2 <- seq(m_min_2, m_max_2, length.out = 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 <- min(design_matrix$M[isclose(design_matrix$τ, T_3)])
m_max_3 <- max(design_matrix$M[isclose(design_matrix$τ, T_3)])
M_grid_3 <- seq(m_min_3, m_max_3, length.out = 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))

plot(M_grid_1, IV_curve_1 * 100, type = 'l', col = '#0A2540', lwd = 2, xlim = c(-0.4, 0.2),
     xlab = 'Log-Moneyness (ln(K/F))', ylab = 'Implied Volatility (%)',
     main = 'Implied Volatility Curves of the S&P 500 Options (2026-06-18)', panel.first = grid(lty = ":"))
lines(M_grid_2, IV_curve_2 * 100, col = '#20639B', lwd = 2)
lines(M_grid_3, IV_curve_3 * 100, col = '#6BAED6', lwd = 2)
legend("topright", legend = c('T = 0.02 Y', 'T = 0.42 Y', 'T = 0.83 Y'), col = c('#0A2540', '#20639B', '#6BAED6'), lty = 1, lwd = 2)


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

print(unique(design_matrix$M))

M_1 <- 0
t_min_1 <- min(design_matrix$τ[isclose(design_matrix$M, M_1, atol = 0.01)])
t_max_1 <- max(design_matrix$τ[isclose(design_matrix$M, M_1, atol = 0.01)])
t_grid_1 <- seq(t_min_1, t_max_1, length.out = 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 <- min(design_matrix$τ[isclose(design_matrix$M, M_2, atol = 0.01)])
t_max_2 <- max(design_matrix$τ[isclose(design_matrix$M, M_2, atol = 0.01)])
t_grid_2 <- seq(t_min_2, t_max_2, length.out = 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 <- min(design_matrix$τ[isclose(design_matrix$M, M_3, atol = 0.01)])
t_max_3 <- max(design_matrix$τ[isclose(design_matrix$M, M_3, atol = 0.01)])
t_grid_3 <- seq(t_min_3, t_max_3, length.out = 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))

plot(t_grid_1, IV_ts_1 * 100, type = 'l', col = '#0A2540', lwd = 2,
     xlab = 'Time to Maturity (T)', ylab = 'Implied Volatility (%)',
     main = 'Implied Volatility Term Structure of the S&P 500 Options (2026-06-18)', panel.first = grid(lty = ":"))
lines(t_grid_2, IV_ts_2 * 100, col = '#20639B', lwd = 2)
lines(t_grid_3, IV_ts_3 * 100, col = '#6BAED6', lwd = 2)
legend("topright", legend = c('M = 0.0 (ATM)', 'M = 0.2 (OTM Call)', 'M = -0.2 (OTM Put)'), col = c('#0A2540', '#20639B', '#6BAED6'), lty = 1, lwd = 2)


#########################################################################
# 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 <- c(alpha_0, alpha_1, alpha_2, alpha_3, alpha_4, alpha_5)

m_min <- min(design_matrix$M); m_max <- max(design_matrix$M)
t_min <- min(design_matrix$τ); t_max <- max(design_matrix$τ)

M_grid <- seq(m_min, m_max, length.out = 10)
T_grid <- seq(t_min, t_max, length.out = 10)

M_mesh <- outer(M_grid, rep(1, length(T_grid)))
T_mesh <- outer(rep(1, length(M_grid)), T_grid)

alphas_covid <- alphas_baseline
alphas_covid[1] <- alphas_covid[1] + 0.62   
alphas_covid[2] <- alphas_covid[2] * 0.90  
alphas_covid[3] <- alphas_covid[3] + 0.15   
alphas_covid[4] <- alphas_covid[4] - 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[1] + 
    alphas_covid[2] * M_mesh + 
    alphas_covid[3] * (M_mesh**2) + 
    alphas_covid[4] * T_mesh + 
    alphas_covid[5] * (T_mesh**2) + 
    alphas_covid[6] * (M_mesh * T_mesh)
)

persp3D(x = M_grid, y = T_grid, z = IV_covid,
        xlim = c(m_min, m_max), ylim = c(0, 1), zlim = c(0.0, max(IV_covid, na.rm = TRUE) * 1.1),
        xlab = 'Log-Moneyness (ln(K/F))', ylab = 'Time to Maturity (T)', zlab = 'Implied Volatility (%)',
        main = 'Implied Volatility Surface of the S&P 500 Options (2026-06-18)',
        theta = -70, phi = 20, expand = 0.5, colkey = FALSE,
        bty = "b2", col = viridis(100), border = "black", lwd = 0.3)


# Plot 2 (Escalation of US-Iran war)

alphas_baseline <- c(alpha_0, alpha_1, alpha_2, alpha_3, alpha_4, alpha_5)

m_min <- min(design_matrix$M); m_max <- max(design_matrix$M)
t_min <- min(design_matrix$τ); t_max <- max(design_matrix$τ)

M_grid <- seq(m_min, m_max, length.out = 10)
T_grid <- seq(t_min, t_max, length.out = 10)

M_mesh <- outer(M_grid, rep(1, length(T_grid)))
T_mesh <- outer(rep(1, length(M_grid)), T_grid)

alphas_future_escalation <- alphas_baseline
alphas_future_escalation[1] <- alphas_future_escalation[1] + 0.25   # Moderate-to-high persistent lift
alphas_future_escalation[2] <- alphas_future_escalation[2] * 1.20   # Steep downside skew (panic put-buying for earnings protection)
alphas_future_escalation[3] <- alphas_future_escalation[3] + 0.08   # Asymmetric smirk curvature
alphas_future_escalation[4] <- alphas_future_escalation[4] + 0.30   # Flat term structure (fear stays high across long horizons)
alphas_future_escalation[6] <- alphas_future_escalation[6] - 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[1] + 
    alphas_future_escalation[2] * M_mesh + 
    alphas_future_escalation[3] * (M_mesh**2) + 
    alphas_future_escalation[4] * T_mesh + 
    alphas_future_escalation[5] * (T_mesh**2) + 
    alphas_future_escalation[6] * (M_mesh * T_mesh)
)

persp3D(x = M_grid, y = T_grid, z = IV_future_escalation,
        xlim = c(m_min, m_max), ylim = c(0, 1), zlim = c(0.0, max(IV_covid, na.rm = TRUE) * 1.1),
        xlab = 'Log-Moneyness (ln(K/F))', ylab = 'Time to Maturity (T)', zlab = 'Implied Volatility (%)',
        main = 'Implied Volatility Surface of the S&P 500 Options (2026-06-18)',
        theta = -70, phi = 20, expand = 0.5, colkey = FALSE,
        bty = "b2", col = viridis(100), border = "black", lwd = 0.3)


# End of the code