Synopsis

Estimating market volatility is a fundamental task in financial engineering, particularly in the pricing of options and risk management. Since the early 2000s, volatility indices like the VIX and VXV have enabled market participants to gain direct exposure to market volatility. These indices are derived from option prices and represent the market’s expectations of future variance under the risk-neutral measure.

This case study aims to reconstruct the VIX and VXV indices from S&P 500 option data using the official CBOE methodology. By computing these indices ourselves, we deepen our understanding of implied volatility surfaces and how derivative pricing on volatility indices differs from standard equity-based instruments.

We will also compare our calculated indices to real-world data pulled from Yahoo Finance to validate our methodology.

Objectives

Load necessary packages and determined days per year and interest rate.

 library(data.table)
 library(tidyverse)
 library(tinytex)
 library(lubridate)
 library(tidyquant)

Data Import and Preprocessing

This code block performs the initial data setup for VIX computation. It: - Sets constants like the risk-free interest rate and trading days in a year. - Loads the raw SP500 options data from a CSV file using fread(). - Calculates the MID_PRICE as the average of the ask and bid prices. - Adjusts STRIKEPRICE by dividing it by 1000 to match standard strike formatting. - Parses Time and MATURITYDATE into proper Date formats. - Computes the time to expiry (years_to_expiry) for each option, which is crucial for forward pricing and volatility estimation.

in_interest_rate <- 0.01
 in_days_in_year <- 365
 SPY500 <- fread("C:/Users/Rithik/Desktop/Dominic/sp500_options_2010 (1).csv")
SPY500[, MID_PRICE := (ASK_PRICE + BID_PRICE) / 2]
SPY500[, STRIKEPRICE := STRIKEPRICE / 1000]
SPY500[, Time := as.Date(Time, format = "%d-%b-%y")]
SPY500[, MATURITYDATE := as.Date(as.character(MATURITYDATE), format = "%Y%m%d")]
SPY500[, years_to_expiry := as.numeric(MATURITYDATE - Time) / in_days_in_year]

Initial Option Chain Visualization

This section generates a scatter plot of SP500 option prices for a specific date (2010-01-04) and maturity (2010-02-20). - It filters the dataset for options with the selected Time and MATURITYDATE. - Then it plots STRIKEPRICE vs MID_PRICE, differentiating option types (calls vs puts) using color coding.

This visualization helps us understand the price distribution of call and put options around different strikes on a given date—an important diagnostic before filtering for valid options in VIX calculation.

Data_Time <- ymd("2010-01-04")
Data_MATURITYDATE <- ymd("2010-02-20")

ggplot(SPY500[Time == Data_Time & MATURITYDATE == Data_MATURITYDATE],
       aes(x = STRIKEPRICE, y = MID_PRICE)) +
  geom_point(aes(color = OPTIONTYPE))

### Forward Price Estimation using Put-Call Parity This block calculates the forward index level for each date-maturity pair using put-call parity. Here’s how it works: - The dataset is reshaped into a wide format where call (C) and put (P) prices are columns. - For each (Time, MATURITYDATE) pair, the strike price where call-put difference is smallest is identified (k_min). - Using this k_min and its corresponding price difference, the forward price is computed via the formula:
forward = k_min + exp(rT) * (C - P) - The forward prices are then merged back into the original SPY500 dataset.

Finally, the forward level is visualized with a vertical line on a plot of STRIKEPRICE vs MID_PRICE. This forward value is critical for identifying at-the-money (ATM) strikes used in the VIX formula.

SPY500_wide <- dcast(SPY500, Time + MATURITYDATE + STRIKEPRICE ~ OPTIONTYPE, value.var = "MID_PRICE")
SPY500_wide[, diff := abs(C - P)]
SPY500_wide[, k_min := .SD[which.min(diff), STRIKEPRICE], by = .(Time, MATURITYDATE)]
SPY500_wide <- SPY500_wide[k_min == STRIKEPRICE, .(Time, MATURITYDATE, k_min, diff)]

SPY500 <- merge(SPY500, SPY500_wide, by = c("Time", "MATURITYDATE"))
SPY500[, forward := k_min + exp(in_interest_rate * years_to_expiry) * diff]

ggplot(SPY500[Time == Data_Time & MATURITYDATE == Data_MATURITYDATE],
       aes(x = STRIKEPRICE, y = MID_PRICE)) +
  geom_point(aes(color = OPTIONTYPE)) +
  geom_vline(aes(xintercept = forward))

### Identifying ATM Strike (k_zero) This section determines the at-the-money (ATM) strike price, denoted as k_zero, which is closest to the computed forward price: - Calculates the difference between each strike price and the forward price. - Filters out any negative differences (strike above forward). - For each (Time, MATURITYDATE) pair, the strike with the minimum non-negative difference is selected as k_zero.

A plot is generated showing STRIKEPRICE vs MID_PRICE with a vertical line indicating the selected k_zero.
This ATM strike serves as a reference point to filter valid options for VIX computation.

SPY500[, diff := forward - STRIKEPRICE]
SPY500[diff < 0, diff := NA]
SPY500[, k_zero := .SD[which.min(diff), STRIKEPRICE], by = .(Time, MATURITYDATE)]

ggplot(SPY500[Time == Data_Time & MATURITYDATE == Data_MATURITYDATE],
       aes(x = STRIKEPRICE, y = MID_PRICE)) +
  geom_point(aes(color = OPTIONTYPE)) +
  geom_vline(aes(xintercept = k_zero))

### Filtering Out Illiquid Options (Double-Zero Bids) This code identifies and flags illiquid options with zero bid prices that have not changed from the previous quote—referred to as “double-zero bids”. These are often deep out-of-the-money contracts with no trading activity.

SPY500 <- SPY500[order(Time, MATURITYDATE, STRIKEPRICE, OPTIONTYPE)]

SPY500[, last_BID_PRICE := shift(BID_PRICE, type = "lag"), 
       by = .(Time, MATURITYDATE, OPTIONTYPE)]
SPY500[, is_double_zero_bid := fifelse(
  !is.na(last_BID_PRICE) & last_BID_PRICE == BID_PRICE & BID_PRICE <= 1e-31,
  TRUE, FALSE)]


SPY500_filtered <- SPY500[Time == Data_Time & MATURITYDATE == Data_MATURITYDATE]

ggplot(SPY500_filtered, aes(x = STRIKEPRICE, y = MID_PRICE)) +
  geom_point(aes(color = is_double_zero_bid)) +
  geom_vline(xintercept = unique(SPY500_filtered$k_zero), 
             linetype = "dashed", color = "red") +
  labs(title = "Visualizing Double-Zero Bids",
       subtitle = paste("Date:", Data_Time))

### Filtering Options Around ATM Strike (k_zero) This code block focuses on selecting out-of-the-money options around the at-the-money strike (k_zero), a standard step in VIX construction The plot visually confirms this filtering with color-coded option types and a vertical black line marking the k_zero threshold. This selection is crucial for computing reliable implied volatility using only liquid, relevant options.

 SPY500 <- SPY500 %>%
  filter((OPTIONTYPE == "C" & STRIKEPRICE >= k_zero) |
         (OPTIONTYPE == "P" & STRIKEPRICE <= k_zero))
SPY500_filtered <- SPY500 %>%
  filter(Time == Data_Time, MATURITYDATE == Data_MATURITYDATE)

ggplot(SPY500_filtered, aes(x = STRIKEPRICE, y = MID_PRICE)) +
  geom_point(aes(color = OPTIONTYPE)) +
  geom_vline(xintercept = unique(SPY500_filtered$k_zero), color = "black") +
  labs(title = "Filtered Options Around K_zero",
       subtitle = paste("Date:", Data_Time))

### Filtering Near-Term Options for VIX Construction This section selects and cleans near-term options for calculating implied volatility: The final cleaned dataset (FilteredNearOptions) contains only relevant, liquid options that meaningfully contribute to the VIX formula. A scatter plot helps visually verify the result, showing option prices by strike with a dashed red line marking the k_near ATM level.

MaturityFrame <- SPY500 %>%
  mutate(MaturityGap = abs(MATURITYDATE - Time - 30)) %>%
  group_by(Time) %>%
  arrange(MaturityGap, .by_group = TRUE) %>%
  distinct(MATURITYDATE, .keep_all = TRUE) %>%
  slice_head(n = 2) %>%
  ungroup()

NearExpiry <- MaturityFrame %>%
  filter(Time == Data_Time) %>%
  slice(1) %>%
  pull(MATURITYDATE)

NextExpiry <- MaturityFrame %>%
  filter(Time == Data_Time) %>%
  slice(2) %>%
  pull(MATURITYDATE)

k_near <- SPY500 %>%
  filter(Time == Data_Time & MATURITYDATE == NearExpiry) %>%
  slice(1) %>%
  pull(k_zero)

k_next <- SPY500 %>%
  filter(Time == Data_Time & MATURITYDATE == NextExpiry) %>%
  slice(1) %>%
  pull(k_zero)

PutDropZone <- SPY500 %>%
  filter(Time == Data_Time & MATURITYDATE == NearExpiry &
         OPTIONTYPE == "P" & STRIKEPRICE < k_near & is_double_zero_bid == TRUE) %>%
  arrange(desc(STRIKEPRICE)) %>%
  slice_head(n = 2) %>%
  pull(STRIKEPRICE)

CallDropZone <- SPY500 %>%
  filter(Time == Data_Time & MATURITYDATE == NearExpiry &
         OPTIONTYPE == "C" & STRIKEPRICE > k_near & is_double_zero_bid == TRUE) %>%
  arrange(STRIKEPRICE) %>%
  slice_head(n = 2) %>%
  pull(STRIKEPRICE)

CleanPuts <- SPY500 %>%
  filter(Time == Data_Time & MATURITYDATE == NearExpiry & OPTIONTYPE == "P" & is_double_zero_bid == FALSE)
if (length(PutDropZone) > 0) {
  CleanPuts <- CleanPuts %>% filter(STRIKEPRICE > min(PutDropZone))
}

CleanCalls <- SPY500 %>%
  filter(Time == Data_Time & MATURITYDATE == NearExpiry & OPTIONTYPE == "C" & is_double_zero_bid == FALSE)
if (length(CallDropZone) > 0) {
  CleanCalls <- CleanCalls %>% filter(STRIKEPRICE < max(CallDropZone))
}

FilteredNearOptions <- bind_rows(CleanPuts, CleanCalls) %>% arrange(STRIKEPRICE)

CleanedPlot <- FilteredNearOptions %>%
  filter(Time == Data_Time & MATURITYDATE == NearExpiry)

ggplot(CleanedPlot, aes(x = STRIKEPRICE, y = MID_PRICE)) +
  geom_point(aes(color = is_double_zero_bid)) +
  geom_vline(xintercept = k_near, linetype = "dashed", color = "red") +
  labs(title = "Filtered Near-Term Options for VIX Calculation",
       subtitle = paste("Date:", Data_Time))

### Filtering Next-Term Options for Volatility Computation This part prepares the next-term options data for use in VIX estimation by: This cleaned next-term set ensures only reliable, liquid options contribute to the next-expiry variance estimate, aligning with CBOE’s methodology. A scatter plot visually confirms the data, with a dashed red line marking k_next, the central strike used for further computation.

put_zero_strikes <- SPY500 %>%
  filter(Time == Data_Time, MATURITYDATE == NextExpiry,
         OPTIONTYPE == "P", STRIKEPRICE < k_next, is_double_zero_bid) %>%
  arrange(-STRIKEPRICE) %>%
  head(2) %>%
  pull(STRIKEPRICE)

call_zero_strikes <- SPY500 %>%
  filter(Time == Data_Time, MATURITYDATE == NextExpiry,
         OPTIONTYPE == "C", STRIKEPRICE > k_next, is_double_zero_bid) %>%
  arrange(STRIKEPRICE) %>%
  head(2) %>%
  pull(STRIKEPRICE)

valid_puts <- SPY500 %>%
  filter(Time == Data_Time, MATURITYDATE == NextExpiry,
         OPTIONTYPE == "P", is_double_zero_bid == FALSE) %>%
  filter(if (length(put_zero_strikes) == 0) TRUE else STRIKEPRICE > min(put_zero_strikes))

valid_calls <- SPY500 %>%
  filter(Time == Data_Time, MATURITYDATE == NextExpiry,
         OPTIONTYPE == "C", is_double_zero_bid == FALSE) %>%
  filter(if (length(call_zero_strikes) == 0) TRUE else STRIKEPRICE < max(call_zero_strikes))

SPY500_next_term_filtered <- bind_rows(valid_puts, valid_calls) %>%
  arrange(STRIKEPRICE)

ggplot(SPY500_next_term_filtered, aes(x = STRIKEPRICE, y = MID_PRICE)) +
  geom_point(aes(color = is_double_zero_bid)) +
  geom_vline(aes(xintercept = k_next), color = "red", linetype = "dashed") +
  labs(title = "Filtered Next-Term Options", subtitle = paste("Date:", Data_Time))

### Calculating Average Mid Prices at ATM Strikes This section computes the average mid-price for options at the at-the-money (ATM) strike k_near and k_next for near- and next-term maturities, respectively. - These averages will later be used to replace the raw call price at the ATM strike to ensure accurate and stable variance contribution in the VIX formula.

These will serve as the base for constructing final datasets used in the implied variance calculation for each term.

NearStrikeRows <- FilteredNearOptions %>%
  filter(MATURITYDATE == NearExpiry & STRIKEPRICE == k_near)

AvgPrice_Near <- NearStrikeRows %>%
  summarise(avg_mid = mean(MID_PRICE, na.rm = TRUE)) %>%
  pull(avg_mid)

NextStrikeRows <- SPY500_next_term_filtered %>%
  filter(MATURITYDATE == NextExpiry & STRIKEPRICE == k_next)

AvgPrice_Next <- NextStrikeRows %>%
  summarise(avg_mid = mean(MID_PRICE, na.rm = TRUE)) %>%
  pull(avg_mid)

NearTermOptionsSet <- FilteredNearOptions %>%
  filter(Time == Data_Time & MATURITYDATE == NearExpiry)

NextTermOptionsSet <- SPY500_next_term_filtered %>%
  filter(Time == Data_Time & MATURITYDATE == NextExpiry)

VIX Calculation Using Prepared Option Data

This section performs the core VIX computation using the cleaned near- and next-term option datasets. - A synthetic “P/C” row is created at the ATM strike (k_near / k_next) using the average mid-price of puts and calls. This replaces the actual call at-the-money to ensure symmetry.

We extract: - Time to maturity (tau_near, tau_next) - Forward price (fwd_near, fwd_next) - Variance is calculated per term, adjusted by the forward price.

PC_row_near <- NearTermOptionsSet %>%
  filter(OPTIONTYPE == "C", STRIKEPRICE == k_near) %>%
  mutate(OPTIONTYPE = "P/C", MID_PRICE = AvgPrice_Near)

NearOptionsPrepared <- NearTermOptionsSet %>%
  filter(STRIKEPRICE != k_near) %>%
  bind_rows(PC_row_near) %>%
  arrange(STRIKEPRICE)

NearOptionsPrepared <- NearOptionsPrepared %>%
  mutate(StrikeSpacing = lead(STRIKEPRICE) - STRIKEPRICE,
         Contribution = (StrikeSpacing / (STRIKEPRICE^2)) *
         exp(in_interest_rate * years_to_expiry) * MID_PRICE)

tau_near <- NearOptionsPrepared %>%
  slice(1) %>%
  pull(years_to_expiry)

fwd_near <- NearOptionsPrepared %>%
  slice(1) %>%
  pull(forward)

variance_near <- 2 / tau_near * sum(NearOptionsPrepared$Contribution, na.rm = TRUE) -
  1 / tau_near * ((fwd_near / k_near - 1)^2)

PC_row_next <- NextTermOptionsSet %>%
  filter(OPTIONTYPE == "C", STRIKEPRICE == k_next) %>%
  mutate(OPTIONTYPE = "P/C", MID_PRICE = AvgPrice_Next)

NextOptionsPrepared <- NextTermOptionsSet %>%
  filter(STRIKEPRICE != k_next) %>%
  bind_rows(PC_row_next) %>%
  arrange(STRIKEPRICE)

NextOptionsPrepared <- NextOptionsPrepared %>%
  mutate(StrikeSpacing = lead(STRIKEPRICE) - STRIKEPRICE,
         Contribution = (StrikeSpacing / (STRIKEPRICE^2)) *
         exp(in_interest_rate * years_to_expiry) * MID_PRICE)

tau_next <- NextOptionsPrepared %>%
  slice(1) %>%
  pull(years_to_expiry)

fwd_next <- NextOptionsPrepared %>%
  slice(1) %>%
  pull(forward)

variance_next <- 2 / tau_next * sum(NextOptionsPrepared$Contribution, na.rm = TRUE) -
  1 / tau_next * ((fwd_next / k_next - 1)^2)

N_30 <- 30 * 1440
N_365 <- in_days_in_year * 1440
N_T1 <- tau_near * N_365
N_T2 <- tau_next * N_365

VIX_estimate <- sqrt((tau_near * variance_near * ((N_T2 - N_30) / (N_T2 - N_T1)) +
                      tau_next * variance_next * ((N_30 - N_T1) / (N_T2 - N_T1))) * N_365 / N_30) * 100

VIX_estimate
## [1] 16.28477

Reusable VIX Function for Any Date and Period

CalculateVIX <- function(Data_Time, PeriodLength) {
  Dates_Maturity <- SPY500 %>%
    mutate(DayDiff = abs(MATURITYDATE - Time - PeriodLength)) %>%
    group_by(Time) %>%
    arrange(DayDiff, .by_group = TRUE) %>%
    distinct(MATURITYDATE, .keep_all = TRUE) %>%
    slice_head(n = 2) %>%
    ungroup()

  NearExpiry <- Dates_Maturity %>% filter(Time == Data_Time) %>% slice(1) %>% pull(MATURITYDATE)
  NextExpiry <- Dates_Maturity %>% filter(Time == Data_Time) %>% slice(2) %>% pull(MATURITYDATE)

  k_near <- SPY500 %>% filter(Time == Data_Time & MATURITYDATE == NearExpiry) %>% slice(1) %>% pull(k_zero)
  k_next <- SPY500 %>% filter(Time == Data_Time & MATURITYDATE == NextExpiry) %>% slice(1) %>% pull(k_zero)

  # Filter and clean near-term options
  near_data <- SPY500 %>% filter(Time == Data_Time & MATURITYDATE == NearExpiry)
  put_bound_near <- near_data %>% filter(OPTIONTYPE == "P", STRIKEPRICE < k_near, is_double_zero_bid) %>% arrange(desc(STRIKEPRICE)) %>% head(2) %>% pull(STRIKEPRICE)
  call_bound_near <- near_data %>% filter(OPTIONTYPE == "C", STRIKEPRICE > k_near, is_double_zero_bid) %>% arrange(STRIKEPRICE) %>% head(2) %>% pull(STRIKEPRICE)

  CleanPutsNear <- near_data %>% filter(OPTIONTYPE == "P", is_double_zero_bid == FALSE)
  if (length(put_bound_near) > 0) CleanPutsNear <- CleanPutsNear %>% filter(STRIKEPRICE > min(put_bound_near))

  CleanCallsNear <- near_data %>% filter(OPTIONTYPE == "C", is_double_zero_bid == FALSE)
  if (length(call_bound_near) > 0) CleanCallsNear <- CleanCallsNear %>% filter(STRIKEPRICE < max(call_bound_near))

  NearOptions <- bind_rows(CleanPutsNear, CleanCallsNear) %>% arrange(STRIKEPRICE)

  # Filter and clean next-term options
  next_data <- SPY500 %>% filter(Time == Data_Time & MATURITYDATE == NextExpiry)
  put_bound_next <- next_data %>% filter(OPTIONTYPE == "P", STRIKEPRICE < k_next, is_double_zero_bid) %>% arrange(desc(STRIKEPRICE)) %>% head(2) %>% pull(STRIKEPRICE)
  call_bound_next <- next_data %>% filter(OPTIONTYPE == "C", STRIKEPRICE > k_next, is_double_zero_bid) %>% arrange(STRIKEPRICE) %>% head(2) %>% pull(STRIKEPRICE)

  CleanPutsNext <- next_data %>% filter(OPTIONTYPE == "P", is_double_zero_bid == FALSE)
  if (length(put_bound_next) > 0) CleanPutsNext <- CleanPutsNext %>% filter(STRIKEPRICE > min(put_bound_next))

  CleanCallsNext <- next_data %>% filter(OPTIONTYPE == "C", is_double_zero_bid == FALSE)
  if (length(call_bound_next) > 0) CleanCallsNext <- CleanCallsNext %>% filter(STRIKEPRICE < max(call_bound_next))

  NextOptions <- bind_rows(CleanPutsNext, CleanCallsNext) %>% arrange(STRIKEPRICE)

  # ATM price estimation
  NearATM <- NearOptions %>% filter(STRIKEPRICE == k_near) %>% summarise(avg = mean(MID_PRICE, na.rm = TRUE)) %>% pull(avg)
  NextATM <- NextOptions %>% filter(STRIKEPRICE == k_next) %>% summarise(avg = mean(MID_PRICE, na.rm = TRUE)) %>% pull(avg)

  # Final sets
  NearFull <- NearOptions %>% filter(Time == Data_Time, MATURITYDATE == NearExpiry)
  NextFull <- NextOptions %>% filter(Time == Data_Time, MATURITYDATE == NextExpiry)

  near_row <- NearFull %>% filter(OPTIONTYPE == "C", STRIKEPRICE == k_near) %>% mutate(OPTIONTYPE = "P/C", MID_PRICE = NearATM)
  NearSet <- NearFull %>% filter(STRIKEPRICE != k_near) %>% bind_rows(near_row) %>% arrange(STRIKEPRICE)

  NearSet <- NearSet %>% mutate(
    delta_K = lead(STRIKEPRICE) - STRIKEPRICE,
    contrib = (delta_K / (STRIKEPRICE^2)) * exp(in_interest_rate * years_to_expiry) * MID_PRICE
  )

  tau_near <- NearSet %>% slice(1) %>% pull(years_to_expiry)
  fwd_near <- NearSet %>% slice(1) %>% pull(forward)

  var_near <- 2 / tau_near * sum(NearSet$contrib, na.rm = TRUE) - 1 / tau_near * ((fwd_near / k_near - 1)^2)

  next_row <- NextFull %>% filter(OPTIONTYPE == "C", STRIKEPRICE == k_next) %>% mutate(OPTIONTYPE = "P/C", MID_PRICE = NextATM)
  NextSet <- NextFull %>% filter(STRIKEPRICE != k_next) %>% bind_rows(next_row) %>% arrange(STRIKEPRICE)

  NextSet <- NextSet %>% mutate(
    delta_K = lead(STRIKEPRICE) - STRIKEPRICE,
    contrib = (delta_K / (STRIKEPRICE^2)) * exp(in_interest_rate * years_to_expiry) * MID_PRICE
  )

  tau_next <- NextSet %>% slice(1) %>% pull(years_to_expiry)
  fwd_next <- NextSet %>% slice(1) %>% pull(forward)

  var_next <- 2 / tau_next * sum(NextSet$contrib, na.rm = TRUE) - 1 / tau_next * ((fwd_next / k_next - 1)^2)

  # Final VIX
  N_30 <- PeriodLength * 1440
  N_365 <- in_days_in_year * 1440
  N_T1 <- tau_near * N_365
  N_T2 <- tau_next * N_365

  VIX <- sqrt((tau_near * var_near * ((N_T2 - N_30) / (N_T2 - N_T1)) +
               tau_next * var_next * ((N_30 - N_T1) / (N_T2 - N_T1))) * N_365 / N_30) * 100
  return(VIX)
}

AllDates <- SPY500 %>% select(Time) %>% distinct()

VIX_30_DAY <- AllDates %>% rowwise() %>% mutate(VIX = CalculateVIX(Time, 30))
options(tibble.print_max = Inf)

🔍 VIX Correlation Analysis with Yahoo Finance Data

We begin by importing the adjusted closing prices of the VIX index, which reflect market expectations of volatility. These are then merged with our internally computed VIX_30_DAY dataset using the Time variable to align dates. This step acts as a robust validation mechanism for our model, ensuring that it captures implied volatility with a high degree of precision. If the correlation is low, it would signal a misalignment between model assumptions and market pricing dynamics.

YahooVIX <- tq_get("^VIX",
                   get = "stock.prices",
                   from = "2010-01-04",
                   to = "2010-12-31") %>%
  transmute(Time = date,
            Source_VIX = adjusted)
MergedVIX <- inner_join(YahooVIX, VIX_30_DAY, by = "Time")
VIX_Correlation <- with(MergedVIX, cor(Source_VIX, VIX, use = "complete.obs"))
print(VIX_Correlation)
## [1] 0.9892307

Visual Comparison: Yahoo VIX vs Calculated VIX

This plot offers a side-by-side time series comparison of the official VIX values from Yahoo Finance and the VIX values we computed using SP500 options data.

By overlaying these two series: - We can visually assess how closely our model tracks the actual market expectations of volatility. - A close alignment between the two lines would support the validity and robustness of our approach.

This comparison not only enhances transparency but also helps identify periods of divergence, which may suggest structural shifts, data limitations, or model calibration issues.

Ultimately, this visualization serves as an intuitive confirmation of whether our model-based VIX aligns with real market benchmarks.

library(ggplot2)
VIX_Merged <- left_join(YahooVIX, VIX_30_DAY, by = "Time")

ggplot(VIX_Merged, aes(x = Time)) +
  geom_line(aes(y = Source_VIX, color = "Yahoo VIX"), linetype = "dashed") +
  geom_line(aes(y = VIX, color = "Calculated VIX"), linetype = "solid") +
  scale_color_manual(values = c("Yahoo VIX" = "blue", "Calculated VIX" = "red")) +
  labs(title = "VIX: Yahoo vs. Calculated",
       y = "VIX Value", color = "Legend") +
  theme_minimal()

### VXV (93-Day Volatility) Summary This code calculates the 93-day VIX (VXV) for each date using the CalculateVIX() function with PeriodLength = 93. - Applies the function row-wise across all dates. - Uses summary() to return key stats: mean, min, max, and quartiles.

Helps evaluate the overall distribution and range of the calculated VXV values.

VXV_summary <- AllDates %>%
  rowwise() %>%
  mutate(VIX = CalculateVIX(Time, 93)) %>%
  ungroup()

summary(VXV_summary$VIX)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.67   20.97   23.97   24.66   27.83   40.15

Comparing Calculated VXV with Yahoo VXV

In this block, we fetch the official 93-day VIX (VXV) from Yahoo Finance and compare it to our calculated values: - ^VXV data is pulled and renamed for clarity. - A joined dataset combines our estimates with Yahoo’s official values based on the date. - The correlation between the two series is computed using cor() with complete.obs. - Shows how closely aligned our VXV model is with the market benchmark. - A high correlation validates the accuracy of our VXV construction approach.

Yahoo_XVIX <- tq_get("^VXV", get = "stock.prices",
                     from = "2010-01-04", to = "2010-12-31") %>%
  transmute(Time = date, Ref_VXV = adjusted)

X_VIX_Comparison <- inner_join(Yahoo_XVIX, VXV_summary, by = "Time")

Correlation_Result <- cor(X_VIX_Comparison$Ref_VXV, X_VIX_Comparison$VIX, use = "complete.obs")

print(Correlation_Result)
## [1] 0.9967151

VXV Comparison — Side-by-Side Visualization

This section plots the Yahoo VXV and our calculated VXV on two vertically stacked plots for a clear visual comparison. - facet_wrap() is used to separate Yahoo and calculated values into distinct panels. - Each panel uses the same x-axis (date) but allows independent y-axis scales for better visibility. > Highlights: - Visually compares accuracy and trend similarity between our model and the market. - The use of faceting provides an uncluttered view, making deviations or tracking patterns easy to detect.

YahooVXV <- Yahoo_XVIX %>%
  transmute(Time = Time, VIX = Ref_VXV, Source = "Yahoo")

MyVXV <- VXV_summary %>%
  mutate(Source = "Calculated") %>%
  select(Time, VIX, Source)

VXV_Combined <- bind_rows(MyVXV, YahooVXV)

ggplot(VXV_Combined, aes(x = Time, y = VIX)) +
  geom_line(color = "steelblue") +
  facet_wrap(~ Source, ncol = 1, scales = "free_y") +
  labs(title = "VXV - Yahoo vs Calculated (Side-by-Side)",
       y = "VXV Value", x = "Date") +
  theme_minimal()