This R script calculates and visualizes Gamma Exposure (GEX) for SPY options based on option chain data. It is primarily inspired by the methodologies described in:

Thomas Moran’s RPubs post on Gamma Exposure

Perfiliev’s blog on calculating Gamma Exposure and Zero Gamma level

The workflow of the script includes:

🔹 Data Preparation The script reads an option chain CSV file containing SPY options data and extracts the current spot price of SPY.

It cleans and formats the data, including parsing expiration dates and computing the days till expiration for each contract.

It also fetches the 10-Year US Treasury rate from FRED to be used as the risk-free rate in the Black-Scholes model.

🔹 GEX Calculation Gamma Exposure is computed for both calls and puts using the formula: GEX = Gamma × Open Interest × 100 × Spot^2 × 0.01 Call GEX is positive, while Put GEX is negative.

The script aggregates GEX per strike to visualize the Total Gamma Exposure distribution across strikes.

🔹 Gamma Surface The core feature of the script is the construction of a Gamma Surface, which simulates how total gamma exposure evolves as:

The spot price of SPY varies across a defined range (e.g., 500 to 560), and

The time to expiration decreases, modeled in reverse time from market open to close in 3-minute intervals (total 405 minutes).

For each spot level and time slice:

The script applies a vectorized Black-Scholes gamma formula for each strike using implied volatility (IV), time to expiry, and open interest.

It computes the sum of gamma exposure for all call and put contracts.

The results are stored in a matrix and converted into a long-format dataframe for plotting.

🔹 Visualization A heatmap is generated using ggplot2, where:

The x-axis represents minutes left until expiration

The y-axis shows the simulated SPY spot price

The color gradient indicates the net gamma exposure, from negative (red) to positive (green)

This surface plot helps visualize the dynamic behavior of gamma exposure across both time and price dimensions, offering traders and analysts a better understanding of dealer hedging pressures and potential gamma inflection zones.

# Set working directory and load libraries
##setwd("E:/GEX/GEX-main")

library(data.table)
library(timeDate)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:data.table':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(stringr)
library(formattable)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:formattable':
## 
##     style
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(tidyr)
library(purrr)
## 
## Attaching package: 'purrr'
## The following object is masked from 'package:data.table':
## 
##     transpose
library(writexl)
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:data.table':
## 
##     yearmon, yearqtr
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## The following objects are masked from 'package:data.table':
## 
##     first, last
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
## The following objects are masked from 'package:data.table':
## 
##     dcast, melt
library(ggplot2)

# Load custom libraries
source("00_scripts/Libraries.R")
libraries()

# === Load option chain and spot ===
file <- "250412spy_quotedata.csv"
option_chain <- fread(file)
spotLine <- fread(file, skip = 1, nrows = 1)
## Warning in fread(file, skip = 1, nrows = 1): Stopped early on line 3. Expected
## 3 fields but found 8. Consider fill=TRUE and comment.char=. First discarded
## non-empty line: <<"Date: April 12, 2025 at 9:10 AM EDT",Bid: 534.1,Ask:
## 534.54,Size: 2*1,"Volume: 97,798,087">>
spotLineData <- strsplit(as.character(colnames(spotLine[,2])), ":")
spotPrice <- as.numeric(str_trim(spotLineData[[1]][2]))

# === Get 10Y risk-free rate from FRED ===
risk_free_rate <- getSymbols(Symbols = "DGS10", src = "FRED", auto.assign = FALSE) %>%
  tail(1) %>% as.data.frame() %>% pull(1)

# (Optional) Get current SPY price
spy_spot <- getSymbols(Symbols = "SPY", src = "yahoo", auto.assign = FALSE) %>%
  tail(1) %>% as.data.frame() %>% pull(1)

# === Define strike range around spot ===
fromStrike <- 0.8 * spotPrice
toStrike <- 1.2 * spotPrice

# === Extract current date from file ===
dateLine <- fread(file, skip = 2, nrows = 1)
## Warning in fread(file, skip = 2, nrows = 1): Stopped early on line 4. Expected
## 5 fields but found 22. Consider fill=TRUE and comment.char=. First discarded
## non-empty line: <<Expiration Date,Calls,Last
## Sale,Net,Bid,Ask,Volume,IV,Delta,Gamma,Open Interest,Strike,Puts,Last
## Sale,Net,Bid,Ask,Volume,IV,Delta,Gamma,Open Interest>>
todayDateData <- strsplit(as.character(colnames(dateLine[,1])), " ")
todayYear <- todayDateData[[1]][4]
todayMonth <- todayDateData[[1]][2]
todayDay <- todayDateData[[1]][3]
todayDate <- as.Date(ymd(paste0(todayYear, '-', todayMonth, '-', todayDay)))

# === Rename option chain columns ===
colnames(option_chain) <- c('ExpirationDate','Calls','CallLastSale','CallNet','CallBid','CallAsk','CallVol',
                            'CallIV','CallDelta','CallGamma','CallOpenInt','StrikePrice','Puts','PutLastSale',
                            'PutNet','PutBid','PutAsk','PutVol','PutIV','PutDelta','PutGamma','PutOpenInt')

# === Parse expiration date ===
option_chain$ExpirationDate <- as.Date(mdy(substring(option_chain$ExpirationDate,5)))

# === Compute time to expiration in years ===
option_chain$daysTillExp <- as.numeric(option_chain$ExpirationDate - todayDate) / 365

# === Compute GEX for calls and puts ===
option_chain$CallGEX <- option_chain$CallGamma * option_chain$CallOpenInt * 100 * spotPrice^2 * 0.01
option_chain$PutGEX  <- option_chain$PutGamma  * option_chain$PutOpenInt  * 100 * spotPrice^2 * 0.01 * -1
option_chain$TotalGamma <- option_chain$CallGEX + option_chain$PutGEX

# === Aggregate GEX by strike price ===
dfAgg <- option_chain %>%
  filter(StrikePrice > fromStrike, StrikePrice < toStrike) %>%
  group_by(StrikePrice) %>%
  summarise(
    CallGEX = sum(CallGEX, na.rm = TRUE),
    PutGEX = sum(PutGEX, na.rm = TRUE),
    TotalGamma = sum(TotalGamma, na.rm = TRUE)
  )

# === Compute total gamma value ===
TotalGammaSum <- round(sum(dfAgg$TotalGamma), 3)

# === Chart: Total Gamma Exposure ===
chart1 <- plot_ly(data = dfAgg, x = ~StrikePrice, y = ~TotalGamma, type = 'bar', name = 'Total Gamma') %>%
  add_lines(x = spotPrice,
            line = list(color = rgb(3, 74, 23, maxColorValue = 255), dash = 'dot'),
            name = paste('Spot Price:', round(spotPrice))) %>%
  layout(
    title = paste("Total Gamma Exposure: $", TotalGammaSum, " Bn per 1% SPX Move -", todayDate),
    xaxis = list(title = "Strike Price"),
    yaxis = list(title = "Gamma Exposure ($ billions/1%)")
  )

# === Chart: Call vs Put Gamma Exposure ===
chart2 <- plot_ly(data = dfAgg, x = ~StrikePrice, y = ~CallGEX, type = 'bar', name = 'Call Gamma') %>%
  add_bars(y = ~PutGEX, marker = list(color = 'red'), name = 'Put Gamma') %>%
  add_lines(x = spotPrice,
            line = list(color = rgb(3, 74, 23, maxColorValue = 255), dash = 'dot'),
            name = paste('Spot Price:', round(spotPrice))) %>%
  layout(
    title = paste("Call vs Put Gamma Exposure -", todayDate),
    xaxis = list(title = "Strike Price"),
    yaxis = list(title = "Gamma Exposure ($ billions/1%)")
  )

chart1
chart2
# === Set simulation range and step ===
step <- 0.5
SPOT_base <- round(spotPrice)
range_low <- round(spotPrice*.9)
range_high <- round(spotPrice*1.1)
levels <- seq(range_high, range_low, by = -step)

# === Select nearest expiration date ===
min_exp <- min(option_chain$daysTillExp[option_chain$daysTillExp > 0], na.rm = TRUE)
filtered_chain <- subset(option_chain, daysTillExp == min_exp)
filtered_chain <- subset(filtered_chain, StrikePrice >= range_low & StrikePrice <= range_high)

# === Function to compute gamma exposure ===
calcGammaEx <- Vectorize(function(S, K, vol, T, r, q, OI) {
  if (T <= 0 || vol == 0) return(0)
  d1 <- (log(S/K) + T*(r - q + (vol^2)/2)) / (vol * sqrt(T))
  gamma <- exp(-q*T) * dnorm(d1) / (S * vol * sqrt(T))
  return(OI * 100 * S^2 * 0.01 * gamma)
})

# === Time intervals in minutes until expiration ===
minutes_remaining <- seq(405, 0, by = -3)

# === Prepare matrix for gamma surface ===
totalGammaMatrix <- matrix(NA, nrow = length(levels), ncol = length(minutes_remaining), dimnames = list(levels, minutes_remaining))

# === Calculate gamma surface ===
for (t in seq_along(minutes_remaining)) {
  T <- minutes_remaining[t]/(60*24)/265  # convert minutes to trading years
  for (i in seq_along(levels)) {
    level <- levels[i]
    callGEX <- calcGammaEx(level, filtered_chain$StrikePrice, filtered_chain$CallIV,
                           T, risk_free_rate, 0, filtered_chain$CallOpenInt)
    putGEX <- calcGammaEx(level, filtered_chain$StrikePrice, filtered_chain$PutIV,
                          T, risk_free_rate, 0, filtered_chain$PutOpenInt)
    totalGammaMatrix[i, t] <- sum(callGEX, na.rm = TRUE) - sum(putGEX, na.rm = TRUE)
  }
}

# === Convert matrix to long-format data frame ===
df <- melt(totalGammaMatrix)
colnames(df) <- c("Level", "MinutesRemaining", "Gamma")
df$Level <- as.numeric(as.character(df$Level))
df$MinutesRemaining <- as.numeric(as.character(df$MinutesRemaining))

# === Plot gamma surface ===
ggplot(df, aes(x = MinutesRemaining, y = Level, fill = Gamma)) +
  geom_tile() +
  scale_fill_gradient2(low = "red", mid = "white", high = "green", midpoint = 0) +
  scale_x_reverse() +
  labs(title = "Gamma Surface for the nearest expedition 0DTE",
       x = "Min Left until Expiry",
       y = "Simulated Spot Price") +
  theme_minimal(base_size = 14)