full code example for implied volatility surface

ivFullCode.R

#load packages
library(quantmod)
library(rvest)
library(reshape2)
library(plotly)
library(akima)

#get underlying stock info and last trade date
symbol <- "AAPL"
priceInfo <- getQuote(symbol)
lastPrice <- priceInfo$Last
divYield <- getQuote(symbol, what = yahooQF("Dividend Yield"))$`Dividend Yield`
if(is.na(divYield)){divYield <- 0}
date <- as.Date(priceInfo$`Trade Time`)

#settings for moneyness and time to maturity
moneynessBoundaries <- c(0.85,1.15)
ttmBoundaries <- c(7, 180)

#scrape full site
baseUrl <- paste0("https://finance.yahoo.com/quote/",symbol,"/options")
baseHTML <- read_html(baseUrl)

#get available expiries and convert to time to maturity
expiriesUNIX <- baseHTML %>% html_nodes("option") %>% html_attr("value")
expiries <- as.Date((baseHTML %>% html_nodes("option") %>% html_text()), format = "%b %d, %Y")
timeToMats <- as.numeric(expiries - date)

#select applicable expiries
sel <- timeToMats >= ttmBoundaries[1] & timeToMats <= ttmBoundaries[2]

expiriesUNIX <- expiriesUNIX[sel]
expiries <- expiries[sel]
timeToMats <- timeToMats[sel]

#loop over expiries to get calls and puts
calls <- NULL
puts <- NULL

for(i in 1:length(expiriesUNIX)){
  expiryUrl <- paste0(baseUrl,"?date=",expiriesUNIX[i])
  expiryHTML <- read_html(expiryUrl)
  tmpCalls <- expiryHTML %>% html_nodes(".calls") %>% html_table()
  if(length(tmpCalls) > 0){
    tmpCalls <- tmpCalls[[1]]
    #sometimes column names are in uppercase, sometimes not
    colnames(tmpCalls) <- tolower(colnames(tmpCalls))
    #remove thousand separator and convert to numeric if applicable
    tmpCalls$strike <- as.numeric(gsub(",","",tmpCalls$strike))
    #add time to maturity
    tmpCalls$ttm <- timeToMats[i]
    #calculate moneyness
    tmpCalls$moneyness <- lastPrice/tmpCalls$strike
    #convert yahoo finance IV to numeric
    tmpCalls$ivOrig <- as.numeric(gsub("%","",gsub(",","",tmpCalls$`implied volatility`)))/100
    calls <- rbind(calls, tmpCalls)
  }
  
  tmpPuts <- expiryHTML %>% html_nodes(".puts") %>% html_table()
  if(length(tmpPuts) > 0){
    tmpPuts <- tmpPuts[[1]]
    #sometimes column names are in uppercase, sometimes not
    colnames(tmpPuts) <- tolower(colnames(tmpPuts))
    #remove thousand separator and convert to numeric if applicable
    tmpPuts$strike <- as.numeric(gsub(",","",tmpPuts$strike))
    #add time to maturity
    tmpPuts$ttm <- timeToMats[i]
    #calculate moneyness
    tmpPuts$moneyness <- tmpPuts$strike/lastPrice
    #convert yahoo finance IV to numeric
    tmpPuts$ivOrig <- as.numeric(gsub("%","",gsub(",","",tmpPuts$`implied volatility`)))/100
    puts <- rbind(puts, tmpPuts)
  }
}

#select only calls within the applicable moneyness boundaries
calls <- calls[calls$moneyness >= moneynessBoundaries[1] & calls$moneyness <= moneynessBoundaries[2],]
#select only calls that have traded during the last 5 minutes of the stocks last trading day
calls <- calls[strptime(calls$`last trade date`,format = "%Y-%m-%d %I:%M%p") >= strptime(paste0(date," 3:55PM EDT"), format = "%Y-%m-%d %I:%M%p"),]

#select only puts within the applicable moneyness boundaries
puts <- puts[puts$moneyness >= moneynessBoundaries[1] & puts$moneyness <= moneynessBoundaries[2],]
#select only puts that have traded during the last 5 minutes of the stocks last trading day
puts <- puts[strptime(puts$`last trade date`,format = "%Y-%m-%d %I:%M%p") >= strptime(paste0(date," 3:55PM EDT"), format = "%Y-%m-%d %I:%M%p"),]


#generalized black scholes merton model
gBSM <- function(S, X, sigma, r, q, ttm, type){
  #S = stock price
  #X = strike price
  #sigma = volatility
  #r = risk free interest rate
  #q = dividend yield
  #ttm = time to maturity in days
  #type = option type
  
  b <- r - q
  t <- ttm/365.25
  
  d1 <- (log(S / X) + (b + sigma^2 / 2) * t) / (sigma * sqrt(t))
  d2 <- d1 - sigma * sqrt(t)
  
  if(type == "call"){
    price <- S * exp((b - r) * t) * pnorm(d1) - X * exp(-r * t) * pnorm(d2)
  }else if (type == "put"){
    price <-  (X * exp(-r * t) * pnorm(-d2) - S * exp((b - r) * t) * pnorm(-d1))
  }
  
  return(price)
}

#objective function
volOptimFun <- function(sigma, price, S, K, r, q, ttm, type){
  abs(price - gBSM(S, K, sigma, r, q, ttm, type))
}

#wrapper for the optimization function so it can be used with apply
getIV <- function(x, S, r, q, type){
  res <- optimize(volOptimFun, interval = c(0, 2), price = as.numeric(x["ask"]), S = S, K = as.numeric(x["strike"]), r = r, q = q, ttm = as.numeric(x["ttm"]), type = type)
  return(res$minimum)
}

#calculating IV
calls$iv <- apply(calls, 1, getIV, S = lastPrice, r = 0.0011, q = divYield, type = "call")
puts$iv <- apply(puts, 1, getIV, S = lastPrice, r = 0.0011, q = divYield, type = "put")

#create grids
ivGridCalls <- acast(calls, ttm ~ moneyness, value.var = "iv")
ivGridPuts <- acast(puts, ttm ~ moneyness, value.var = "iv")

#get coordinates of NAs in grid
toInterpolate <- which(is.na(ivGridCalls))
coords <- cbind(toInterpolate%%dim(ivGridCalls)[1], toInterpolate%/%dim(ivGridCalls)[1] + 1)
coords[coords[,1] == 0, 2] <- coords[coords[,1] == 0, 2] - 1 
coords[coords[,1] == 0, 1] <- dim(ivGridCalls)[1]

#loop through NAs and interpolate
for(i in 1:nrow(coords)){
  #get the coordinates of a 10x10 area around the missing value
  x1 <- max(coords[i,1] - 10, 1)
  x2 <- min(coords[i,1] + 10, dim(ivGridCalls)[1])
  y1 <- max(coords[i,2] - 10, 1)
  y2 <- min(coords[i,2] + 10, dim(ivGridCalls)[2])
  
  #get the moneyness/time to mat combination of the missing value
  x0 <- as.numeric(rownames(ivGridCalls)[coords[i,1]])
  y0 <- as.numeric(colnames(ivGridCalls)[coords[i,2]])
  
  #get the part of the grid that is used to interpolate and remove all missing values that are present
  interpGrid <- ivGridCalls[x1:x2,y1:y2]
  interpGrid <- melt(interpGrid)
  interpGrid <- na.omit(interpGrid)
  
  #interpolate linearly
  interpVal <- interp(x = interpGrid$Var1, y = interpGrid$Var2, z = interpGrid$value,
                      xo = x0, yo = y0,
                      linear = TRUE, extrap = TRUE)$z[1,1]
  
  #if linear interpolation doesnt yield a result, use spline interpolation
  if(is.na(interpVal)){
    interpVal <- interp(x = interpGrid$Var1, y = interpGrid$Var2, z = interpGrid$value,
                        xo = x0, yo = y0,
                        linear = FALSE, extrap = TRUE)$z[1,1]
  }
  
  #if the resulting value is clearly wrong, e.g. negative or way outside the values that are used to interpolate,
  #leave it as NA
  if(interpVal < 0 | interpVal > max(interpGrid$value * 1.5)){
    interpVal <- NA
  }
  
  #replace the value with the result of the interpolation
  ivGridCalls[coords[i,1],coords[i,2]] <- interpVal
}

#plot the resulting implied volatility surface
xaxx <- list(
  gridcolor='rgb(255, 255, 255)',
  zerolinecolor='rgb(255, 255, 255)',
  showbackground=TRUE,
  backgroundcolor='rgb(230, 230,230)',
  title = "Moneyness"
)

yaxx <- list(
  gridcolor='rgb(255, 255, 255)',
  zerolinecolor='rgb(255, 255, 255)',
  showbackground=TRUE,
  backgroundcolor='rgb(230, 230,230)',
  title = "Time To Maturity"
)

zaxx <- list(
  gridcolor='rgb(255, 255, 255)',
  zerolinecolor='rgb(255, 255, 255)',
  showbackground=TRUE,
  backgroundcolor='rgb(230, 230,230)',
  tickformat = "%",
  title = "Implied Volatility"
)

fig <- plot_ly(x = colnames(ivGridCalls), y =  rownames(ivGridCalls), z = ivGridCalls)
fig <- fig %>% add_surface()
fig <- fig %>% layout(scene = list(xaxis=xaxx, yaxis=yaxx, zaxis = zaxx))
fig <- fig %>% plotly::colorbar(title = "", x = 0.9, y = 0.75, tickformat = "%")
fig