Implied Volatility Surface

Overview

The following vignette demonstrates the steps needed to visualize implied volatility surface.

Option Data

Note that the quantmod library allows you to download the option chain at each point of time. This means that for each underlying, there is a chain of options that move accordingly, depending on liquidity, as the underlying moves. This covers different strike prices and maturities. The following works with the source code of getOptionChain for the SPY ETF:

library(quantmod)
library(lubridate)
library(plyr)
library(plotly)

rm(list = ls())

tic <- "SPY"

`getOptionChain` <-
  function(Symbols, Exp=NULL, src="yahoo", ...) {
    Call <- paste("getOptionChain",src,sep=".")
    if(missing(Exp)) {
      do.call(Call, list(Symbols=Symbols, ...))
    } else {
      do.call(Call, list(Symbols=Symbols, Exp=Exp, ...))
    }
  }

getOptionChain.yahoo <- function(Symbols, Exp, ...)
{
  if(!requireNamespace("jsonlite", quietly=TRUE))
    stop("package:",dQuote("jsonlite"),"cannot be loaded.")
  
  NewToOld <- function(x) {
    if(is.null(x) || length(x) < 1)
      return(NULL)
    # clean up colnames, in case there's weirdness in the JSON
    names(x) <- tolower(gsub("[[:space:]]", "", names(x)))
    # set cleaned up colnames to current output colnames
    d <- with(x, data.frame(Strike=strike, Last=lastprice, Chg=change,
                            Bid=bid, Ask=ask, Vol=volume, OI=openinterest,
                            row.names=contractsymbol, stringsAsFactors=FALSE))
    
    # remove commas from the numeric data
    d[] <- lapply(d, gsub, pattern=",", replacement="", fixed=TRUE)
    d[] <- lapply(d, type.convert, as.is=TRUE)
    d
  }
  
  # Don't check the expiry date if we're looping over dates we just scraped
  checkExp <- !hasArg(".expiry.known") || !match.call(expand.dots=TRUE)$.expiry.known
  # Construct URL
  urlExp <- paste0("https://query2.finance.yahoo.com/v7/finance/options/", Symbols[1])
  # Add expiry date to URL
  if(!checkExp)
    urlExp <- paste0(urlExp, "?&date=", Exp)
  
  # Fetch data (jsonlite::fromJSON will handle connection)
  tbl <- jsonlite::fromJSON(urlExp)
  
  # Only return nearest expiry (default served by Yahoo Finance), unless the user specified Exp
  if(!missing(Exp) && checkExp) {
    all.expiries <- tbl$optionChain$result$expirationDates[[1]]
    all.expiries.posix <- .POSIXct(as.numeric(all.expiries), tz="UTC")
    
    # this is a recursive command
    if(is.null(Exp)) {
      # Return all expires if Exp = NULL
      out <- lapply(all.expiries, getOptionChain.yahoo, Symbols=Symbols, .expiry.known=TRUE)
      # Expiry format was "%b %Y", but that's not unique with weeklies. Change
      # format to "%b.%d.%Y" ("%Y-%m-%d wouldn't be good, since names should
      # start with a letter or dot--naming things is hard).
      return(setNames(out, format(all.expiries.posix, "%b.%d.%Y")))
    }     
    
    else {
      # Ensure data exist for user-provided expiry date(s)
      if(inherits(Exp, "Date"))
        valid.expiries <- as.Date(all.expiries.posix) %in% Exp
      else if(inherits(Exp, "POSIXt"))
        valid.expiries <- all.expiries.posix %in% Exp
      else if(is.character(Exp)) {
        expiry.range <- range(unlist(lapply(Exp, .parseISO8601, tz="UTC")))
        valid.expiries <- all.expiries.posix >= expiry.range[1] &
          all.expiries.posix <= expiry.range[2]
      }
      if(all(!valid.expiries))
        stop("Provided expiry date(s) not found. Available dates are: ",
             paste(as.Date(all.expiries.posix), collapse=", "))
      
      expiry.subset <- all.expiries[valid.expiries]
      if(length(expiry.subset) == 1)
        return(getOptionChain.yahoo(Symbols, expiry.subset, .expiry.known=TRUE))
      else {
        out <- lapply(expiry.subset, getOptionChain.yahoo, Symbols=Symbols, .expiry.known=TRUE)
        # See comment above regarding the output names
        return(setNames(out, format(all.expiries.posix[valid.expiries], "%b.%d.%Y")))
      }
    }
  }
  
  dftables <- lapply(tbl$optionChain$result$options[[1]][,c("calls","puts")], `[[`, 1L)
  #dftables <- mapply(NewToOld, x=dftables, SIMPLIFY=FALSE)
  
  
  fix_date <- function(x) {
    if(class(x) == "list") 
      return(NULL)
    x$expiration <- .POSIXct(as.numeric(x$expiration), tz="UTC")
    x$lastTradeDate <- .POSIXct(as.numeric(x$lastTradeDate), tz="UTC")
    x <- x[,sort(names(x))]
    return(x)
  }
  
  dftables <- lapply(dftables,fix_date)
  dftables <- dftables[!sapply(dftables,is.null)]
  dftables
}

# EXAMPLE to get all expiration in a single data.frame object
ds <- getOptionChain.yahoo(tic,NULL)
ds <- lapply(ds, function(ds_i) lapply(1:length(ds_i), 
                                       function(i) data.frame(Type = names(ds_i)[i], ds_i[[i]]))  )
ds <- lapply(ds, function(ds_i) do.call(plyr::rbind.fill,ds_i)  )
ds <- do.call(plyr::rbind.fill,ds)
ds$Date <- date(ds$lastTradeDate )
ds$Expiration <- date(ds$expiration)

We have about rncol(ds)` variables in the data with 7283 observations.

head(ds)

Surface

smile.T <- function(T.) {
  
  # find the closest expiration data around T.
  which.exp <- ds[which.min( abs(  (ds$Expiration - ds$Date)/250 - T.) ),"Expiration"]
  T. <- min(unique(which.exp - ds$Date)/250)
  
  ds.i <- ds[ds$Expiration %in% which.exp,]
  dim(ds.i)
  
  ds.i <- ds.i[with(ds.i,order(Type,strike)),]
  
  ds.l <- dlply(ds.i,"Type",function(x) x[,c("strike","lastPrice")]  )
  ds.i <- Reduce(function(...) merge(..., by = "strike"), ds.l  )
  dim(ds.i)
  
  names(ds.i) <- c("K","C","P")
  x <- lm(I(C-P) ~ I(-K), data = ds.i)[[1]]
  disc <- as.numeric(x[2])
  PVF <- as.numeric(x[1])
  r <- (-log(disc)/T.)
  
  # compare with BS formula
  bs <- function(PVF, K, disc, T., sigma) {
    values <- c(2)
    
    d1 <- log(PVF/(K*disc))/(sigma*sqrt(T.)) + 0.5*sigma*sqrt(T.)
    d2 <- d1 - sigma * sqrt(T.)
    
    values[1] <- PVF*pnorm(d1) - K*disc*pnorm(d2)
    values[2] <- K*disc*pnorm(-d2) - PVF*pnorm(-d1)
    
    return(values)
  }
  
  
  K <- ds.i[,"K"]
  C <- ds.i[,"C"]
  P <- ds.i[,"P"]
  
  Sigma <- numeric()
  Sigma2 <- numeric()
  
  for(i in 1:length(K)) {
    # cat(i,"\n")
    sigma.i <- try(uniroot(function(x) bs(PVF, K[i], disc, T.,x)[1]  - C[i], c(-1,2))$root, silent = T)
    if(inherits(sigma.i,"try-error")) 
      sigma.i <- NA
    sigma.j <- try(uniroot(function(x) bs(PVF, K[i], disc, T.,x)[2]  - P[i], c(-1,2))$root, silent = T)
    if(inherits(sigma.j,"try-error")) 
      sigma.j <- NA
    
    
    Sigma <- c(Sigma,sigma.i) # call prices
    Sigma2 <- c(Sigma2,sigma.j) # put pries
    
  }
  
  smile <- data.frame(Sigma = Sigma,K)
  smile <- smile[!smile$Sigma < 0.05,]
  smile$T. <- T.
  smile <- na.omit(smile)
  
  # {
  #   plot(Sigma~K, data =smile, pch = 20, cex = 0.5, main = "Call Options", ylab = expression(sigma) )
  #   lines(fitted(loess(Sigma~K,data = smile,span = 1))~smile$K, col = 2)
  #   grid(10)
  # }
  
  smile2 <- data.frame(Sigma = Sigma2,K)
  smile2 <- smile2[!smile2$Sigma < 0.05,]
  smile2$T. <- T.
  smile2 <- na.omit(smile2)
  # {
  #   plot(Sigma~K, data =smile2, pch = 20, cex = 0.5, main = "Put Options", ylab = expression(sigma))
  #   lines(fitted(loess(Sigma~K,data = smile2,span = 1))~smile2$K, col = 2)
  #   grid(10)
  # }
  
  list(smile,smile2,PVF,T.)
  
}

T.seq <- seq(3,36,by = 3)/12
smile.l <- lapply(T.seq, smile.T)

smile.call <- lapply(smile.l, function(x) x[[1]] ) 
smile.put <- lapply(smile.l, function(x) x[[2]] ) 
smile.pvf <- sapply(smile.l, function(x) x[[3]] ) 
smile.T <- sapply(smile.l, function(x) x[[4]] ) 

smile.call <- unique(Reduce(rbind,smile.call))
smile.put <- unique(Reduce(rbind,smile.put))

# add underlying price
S_0 <- as.numeric(last(get(getSymbols(tic))[,6]))
smile.call$INM <- log(S_0/smile.call$K)
smile.put$INM <- log(S_0/smile.put$K)
smile.call <- smile.call[smile.call$INM < 0,]
smile.put <- smile.put[smile.put$INM > 0,]

plot_ds1 <- smile.call[smile.call$T. %in% min(smile.call$T.),]
plot_ds2 <- smile.put[smile.put$T. %in% min(smile.put$T.),]
plot_ds <- rbind(plot_ds1,plot_ds2)
plot_ds <- plot_ds[order(plot_ds$K),]
plot(Sigma~K, data =plot_ds, pch = 20, cex = 0.5, main = "Volatility Smile for OTM Calls and Puts", 
     ylab = expression(sigma) )
lines(fitted(loess(Sigma~K,data = plot_ds,span = 1))~plot_ds$K, col = 2)
grid(10)

Finally, we use plotly to visualize the surface

model <- loess(Sigma~K+T.,data = smile.call, control = loess.control(surface = "interpolate"))

K.seq <- seq(range(smile.call$K)[1],range(smile.call$K)[2],by = 1)
T.seq <- seq(range(smile.call$T.)[1],range(smile.call$T.)[2],by = 1/12)
X <- expand.grid(K.seq,T.seq)
names(X) <- c("K","T.")
Z <- predict(model,X)


{
  f <- list(
    family = "Courier New, monospace",
    size = 11,
    color = "#7f7f7f"
  )
  
  axx <- list(
    title = "Exercise Price",
    titlefont = f
  )
  
  axy <- list(
    title = "Expiration (Years)",
    titlefont = f
    
  )
  
  axz <- list(
    title = "Implied Volatility",
    titlefont = f
    
  )
  
  p <- plot_ly(x = K.seq, y = T.seq, z = Z) %>% add_surface() 
  p <- layout(p,title = paste("Volatility Surface for", tic),scene = list(xaxis=axx,yaxis=axy,zaxis=axz))
  p
}

## This version of 'bslib' is designed to work with 'shiny' >= 1.6.0.
##     Please upgrade via install.packages('shiny').