BSM Model with GPD as the tail distribution in option pricing

Black Scholes Merton Model with Generalized Pareto Distribution as tail distribution in option pricing.

# USD/CAD Options Pricing with Black-Scholes and GPD Tails
# Install required packages if not already installed
# install.packages(c("quantmod", "evir", "MASS"))

library(quantmod)
library(evir)
library(MASS)
library(tidyverse)


# Step 1: Fetch USD/CAD data
cat("Fetching USD/CAD data...\n")

## Fetching USD/CAD data...

getSymbols("CADUSD=X", src = "yahoo", from = "2020-01-01", auto.assign = TRUE)

## [1] "CADUSD=X"

usdcad <- `CADUSD=X`
prices <- na.omit(Cl(usdcad))

# Calculate log returns
returns <- diff(log(prices))
returns <- na.omit(returns)

# Current spot price
S0 <- as.numeric(tail(prices, 1))
cat(sprintf("Current USD/CAD spot price: %.4f\n", S0))

## Current USD/CAD spot price: 0.7252

# Step 2: Fit Generalized Pareto Distribution to tails
cat("\nFitting GPD to tail distributions...\n")

## 
## Fitting GPD to tail distributions...

# Define threshold (95th percentile for right tail, 5th for left tail)
threshold_upper <- quantile(returns, 0.95)
threshold_lower <- quantile(returns, 0.05)

# Fit GPD to upper tail (for calls - right tail risk)
upper_exceedances <- returns[returns > threshold_upper] - threshold_upper
gpd_upper <- gpd(as.numeric(returns), threshold = as.numeric(threshold_upper))

# Fit GPD to lower tail (for puts - left tail risk)
lower_exceedances <- abs(returns[returns < threshold_lower] - threshold_lower)
gpd_lower <- gpd(-as.numeric(returns), threshold = -as.numeric(threshold_lower))

cat(sprintf("Upper tail threshold: %.6f\n", threshold_upper))

## Upper tail threshold: 0.006802

cat(sprintf(
  "GPD upper tail shape parameter (xi): %.4f\n",
  gpd_upper$par.ests[1]
))

## GPD upper tail shape parameter (xi): 0.2664

cat(sprintf(
  "GPD upper tail scale parameter (beta): %.6f\n",
  gpd_upper$par.ests[2]
))

## GPD upper tail scale parameter (beta): 0.002055

# Step 3: Calculate volatility incorporating tail behavior
# Standard volatility
vol_standard <- sd(returns) * sqrt(252)

# Adjusted volatility considering GPD tails
vol_adj <- sd(returns) * sqrt(252) * (1 + abs(gpd_upper$par.ests[1]) * 0.5)

cat(sprintf(
  "\nStandard volatility: %.4f (%.2f%%)\n",
  vol_standard,
  vol_standard * 100
))

## 
## Standard volatility: 0.0688 (6.88%)

cat(sprintf("GPD-adjusted volatility: %.4f (%.2f%%)\n", vol_adj, vol_adj * 100))

## GPD-adjusted volatility: 0.0779 (7.79%)

# Step 4: Black-Scholes with GPD adjustment
black_scholes_gpd <- function(S, K, T, r, sigma, type = "call", gpd_shape = 0) {
  # Adjust for tail risk
  tail_adj <- 1 + abs(gpd_shape) * sqrt(T)
  sigma_adj <- sigma * tail_adj

  d1 <- (log(S / K) + (r + 0.5 * sigma_adj^2) * T) / (sigma_adj * sqrt(T))
  d2 <- d1 - sigma_adj * sqrt(T)

  if (type == "call") {
    price <- S * pnorm(d1) - K * exp(-r * T) * pnorm(d2)
  } else {
    price <- K * exp(-r * T) * pnorm(-d2) - S * pnorm(-d1)
  }

  return(price)
}

# Step 5: Set parameters and calculate option prices
r <- 0.05 # Risk-free rate (5%)
T_values <- c(1 / 12, 3 / 12, 6 / 12, 1, 1.5, 2) # Time to maturity in years
K_values <- seq(S0 * 0.95, S0 * 1.05, length.out = 6) # Strike prices around spot

# Create results dataframe
results <- data.frame(
  Strike = numeric(),
  Maturity_Months = numeric(),
  Call_Price = numeric(),
  Put_Price = numeric(),
  Call_Standard = numeric(),
  Put_Standard = numeric()
)

cat("\n", strrep("=", 70), "\n")

## 
##  ======================================================================

cat("OPTIONS PRICING RESULTS\n")

## OPTIONS PRICING RESULTS

cat("\n", strrep("=", 70), "\n")

## 
##  ======================================================================

# Calculate first 6 call and put prices
for (i in 1:6) {
  K <- K_values[i]
  T <- T_values[min(i, length(T_values))]

  # GPD-adjusted prices
  call_gpd <- black_scholes_gpd(
    S0,
    K,
    T,
    r,
    vol_adj,
    "call",
    gpd_upper$par.ests[1]
  )
  put_gpd <- black_scholes_gpd(
    S0,
    K,
    T,
    r,
    vol_adj,
    "put",
    gpd_lower$par.ests[1]
  )

  # Standard Black-Scholes for comparison
  call_std <- black_scholes_gpd(S0, K, T, r, vol_standard, "call", 0)
  put_std <- black_scholes_gpd(S0, K, T, r, vol_standard, "put", 0)

  results <- rbind(
    results,
    data.frame(
      Strike = K,
      Maturity_Months = T * 12,
      Call_Price = call_gpd,
      Put_Price = put_gpd,
      Call_Standard = call_std,
      Put_Standard = put_std
    )
  )
}

# Display results
print(results[, 1:4], row.names = FALSE)

##     Strike Maturity_Months Call_Price    Put_Price
##  0.6889350               1 0.03918897 5.039225e-05
##  0.7034389               3 0.03324351 2.359645e-03
##  0.7179427               6 0.03371090 7.628362e-03
##  0.7324466              12 0.04449009 1.367433e-02
##  0.7469505              18 0.05417635 1.838785e-02
##  0.7614544              24 0.06327560 2.229187e-02

cat("\n", strrep("=", 70), "\n")

## 
##  ======================================================================

cat("COMPARISON: GPD-Adjusted vs Standard Black-Scholes\n")

## COMPARISON: GPD-Adjusted vs Standard Black-Scholes

cat("\n", strrep("=", 70), "\n")

## 
##  ======================================================================

results$Call_Diff_Pct <- ((results$Call_Price - results$Call_Standard) /
  results$Call_Standard) *
  100
results$Put_Diff_Pct <- ((results$Put_Price - results$Put_Standard) /
  results$Put_Standard) *
  100

print(
  results[, c("Strike", "Maturity_Months", "Call_Diff_Pct", "Put_Diff_Pct")],
  row.names = FALSE
)

##     Strike Maturity_Months Call_Diff_Pct Put_Diff_Pct
##  0.6889350               1     0.1373236    362.23009
##  0.7034389               3     4.7671427     90.79342
##  0.7179427               6    13.2845361     59.59883
##  0.7324466              12    20.3946395     61.18461
##  0.7469505              18    25.0833565     65.71610
##  0.7614544              24    28.5029061     70.97133

cat("\n\nInterpretation:\n")

## 
## 
## Interpretation:

cat(
  "- Positive difference % means GPD model prices options higher (accounts for tail risk)\n"
)

## - Positive difference % means GPD model prices options higher (accounts for tail risk)

cat("- The GPD adjustment captures fat-tail behavior in currency markets\n")

## - The GPD adjustment captures fat-tail behavior in currency markets

cat("- This is particularly important for out-of-the-money options\n")

## - This is particularly important for out-of-the-money options

## Histogram
returns_df <- data.frame(
  Date = index(returns),
  Returns = as.numeric(returns)
)

# Create histogram with density overlay
p <- ggplot(returns_df, aes(x = Returns)) +
  geom_histogram(
    aes(y = after_stat(density)),
    bins = 50,
    fill = "steelblue",
    alpha = 0.7,
    color = "white"
  ) +
  geom_density(color = "darkred", linewidth = 1.2) +
  geom_vline(
    xintercept = threshold_upper,
    linetype = "dashed",
    color = "darkgreen",
    linewidth = 0.8
  ) +
  geom_vline(
    xintercept = threshold_lower,
    linetype = "dashed",
    color = "darkgreen",
    linewidth = 0.8
  ) +
  labs(
    title = "USD/CAD Log Returns Distribution",
    subtitle = "Histogram with Kernel Density Estimate and GPD Thresholds",
    x = "Log Returns",
    y = "Density"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    plot.subtitle = element_text(hjust = 0.5, size = 10),
    panel.grid.minor = element_blank()
  )

print(p)