Search This Blog

Monday, October 24, 2016

Creating volatility smile in R

In this post, I try to use R to draw volatility curves. This is a draft version that I will be updating on a regular basis. I am using the Russell 2000 trade and option chain data for plotting these graphs.

Here are the outputs (and the PDF file)



 Here is my code


#!/usr/bin/Rscript
#!/usr/bin/Rscript
# russ2k_analysis.r
# source("russ2k_analysis.r")

library(RCurl)
library(e1071)
library(RQuantLib)
library(RcppBDT)
library(ggplot2)

home.folder <- '~/Work/finance/RUSS2K/'
#Load historic closing prices
today <- Sys.Date()
day <- as.numeric(format(today,'%d'))
month <- as.numeric(format(today,'%m'))
year <- as.numeric(format(today,'%Y'))

#Author: Arthgallo Wachs (arthgallo.wachs@gmail.com)
URL <- paste('http://real-chart.finance.yahoo.com/table.csv?s=%5ERUT&a=',(month-1),'&b=',(day),'&c=',(year-1),'&d=',(month-1),'&e=',day,'&f=',year,'&g=d&ignore=.csv',sep="")
URL
targetFileName <- paste(home.folder,"russ2k",today,".csv",sep="")
download.file(URL,targetFileName,quiet=FALSE,mode="w",cacheOK=TRUE)
russ2k <- read.csv(targetFileName, head=TRUE, sep=",")
russ2k$Date <- as.Date(russ2k$Date)
names(russ2k)
hist(russ2k$Close,breaks=50)
mean(russ2k$Close)
median(russ2k$Close)
mode(russ2k$Close)
kurtosis(russ2k$Close)
skewness(russ2k$Close)
moment(russ2k$Close, order=3, center=TRUE) 
underlying <- as.numeric(russ2k$Close[1])

source(paste(home.folder,"cSplit.r",sep=""))

# Load option chain
# http://www.cboe.com/delayedquote/QuoteTableDownload.aspx
option_chain_file <- paste(home.folder,"russ2k_optionchain.dat",sep="")

# russ2k.opt <- read.csv(option_chain_file, head=TRUE, sep=",")
# Read the last closing from option chain

russ2k.opt <- read.table(option_chain_file, 
                header=FALSE, 
                sep=",", 
                skip = "3", 
                col.names=c("Calls","Last.Sale","Net","Bid","Ask","Vol",
                "Open.Int","Puts","Last Sale.1","Net.1","Bid.1","Ask.1",
                "Vol.1","Open.Int.1","X"))

names(russ2k.opt)
russ2k.opt <- cSplit(russ2k.opt, c("Calls","Puts"), sep=" ")
names(russ2k.opt)
names(russ2k.opt)[names(russ2k.opt)=="Calls_1"] <- "year"
names(russ2k.opt)[names(russ2k.opt)=="Calls_2"] <- "month"
names(russ2k.opt)[names(russ2k.opt)=="Calls_3"] <- "strike"
names(russ2k.opt)[names(russ2k.opt)=="Calls_4"] <- "call_name"
names(russ2k.opt)[names(russ2k.opt)=="Puts_4"] <- "put_name"
names(russ2k.opt)[names(russ2k.opt)=="Bid"] <- "CBid"
names(russ2k.opt)[names(russ2k.opt)=="Ask"] <- "CAsk"
names(russ2k.opt)[names(russ2k.opt)=="Bid.1"] <- "PBid"
names(russ2k.opt)[names(russ2k.opt)=="Ask.1"] <- "PAsk"

names(russ2k.opt)

today <- Sys.Date()
russ2k.opt$year <- paste("20",russ2k.opt$year,sep="")

# Read the option chain and convert each expiry mm/YYYY to first of that month.
# This is needed since source month is "Jun" which needs to be translated to 06
russ2k.opt$expdate <- as.Date(paste(russ2k.opt$year,russ2k.opt$month,"01", sep="-"),format="%Y-%b-%d")
russ2k.opt$month <- format(russ2k.opt$expdate,"%m")

# Convert expiration date to third Friday for that month
russ2k.opt$expdate <- apply(russ2k.opt, 1, function(x) {getNthDayOfWeek(third,Fri,as.numeric(x[["month"]]),as.numeric(x[["year"]]))} )
russ2k.opt$expdate <- as.Date(russ2k.opt$expdate,origin="1970-01-01")

# Now convert expiration date to expiration days and expiration fractional years
russ2k.opt$expdays <- as.numeric(russ2k.opt$expdate - today,units="days")
russ2k.opt$expyrs <- russ2k.opt$expdays/365

# Now calculate implied Volatility for each strike price based on value
russ2k.opt$call_iv <- apply(russ2k.opt, 1, function(x) 
{
  tryCatch(
    EuropeanOptionImpliedVolatility(
      type  = "call",
      value = as.numeric(x[["CBid"]])+(as.numeric(x[["CAsk"]])-as.numeric(x[["CBid"]])/2),
      underlying = underlying,
      strike     = as.numeric(x[["strike"]]),
      dividendYield = 0,
      riskFreeRate  = 0.0003,
      maturity      = as.numeric(x[["expyrs"]]),
      volatility    = .3
  )[[1]],
  error = function(cond) {return(NA)},
  warning = function(cond) {return(NA)},
  finally = {
      #Do Nothing
  })
})

# Plot the volatility smile
#Drop column X 
russ2k.opt <- within(russ2k.opt, rm(X))

# Omit all rows with implied volatility as NA
russ2k.opt <- russ2k.opt[complete.cases(russ2k.opt),]
russ2k.opt$exprank <- rank(russ2k.opt$expdays, ties.method="max")

nextexpdays <- min(russ2k.opt$expdays)
russ2k_nextexp <- russ2k.opt[russ2k.opt$expdays==nextexpdays]

#Plot the implied volatility smiles on the graph
pdf("russell2000.pdf",width=10,height=8,paper="USr",pointsize=12)
russ2k.opt$strike <- as.numeric(russ2k.opt$strike)
xmin <- min(russ2k.opt$strike)
xmax <- max(russ2k.opt$strike)
ymin <- min(russ2k.opt$call_iv)
ymax <- max(russ2k.opt$call_iv)
ymin <- 0.1;ymax <- 0.3
plot(russ2k_nextexp$strike, russ2k_nextexp$call_iv, typ="l", col="green", xlim=c(xmin, xmax), ylim= c(ymin,ymax),
     main=c("Russell 2000 Call", "Volatility Smile"), xlab="Strike", ylab="Implied Vol")
unk_expdays <- unique(russ2k.opt$expdays)
unk_expdates <- unique(as.Date(russ2k.opt$expdate))
numexp <- length(unk_expdays)
for(i in 2:numexp)
{
 nextexpdays <- unk_expdays[i] 
 russ2k_nextexp <- russ2k.opt[russ2k.opt$expdays==nextexpdays]
 lines(russ2k_nextexp$strike, russ2k_nextexp$call_iv, col=rainbow(16)[i])
}
legend("bottomleft", cex=0.9, legend=as.vector(as.character(unk_expdates)), lty=1, col=rainbow(16)[2:numexp])
abline(v=underlying, col="red",lty=2)
closest_strike <- russ2k.opt$strike[which.min(abs(russ2k.opt$strike-underlying))]
closest_exp <- min(russ2k.opt$expdays)
closest_iv <- russ2k.opt[russ2k.opt$strike == closest_strike & russ2k.opt$expdays == closest_exp][1]$call_iv

abline(h=closest_iv, col="black", lty=2)

#Function for calculating Greeks
FxCalcGreeks <- function (type, 
   value, 
   underlying, 
   strike, 
   dividendYield,
   riskFreeRate,
   maturity,
   implied_volatility)
{
  #Calculates the greeks when passed a set of arguments
 tryCatch({
  # Compute all the greeks
  eurOp<-EuropeanOption(
     type  = type, 
     underlying = underlying,
     strike = strike, 
     dividendYield = dividendYield,
     riskFreeRate  = riskFreeRate,
     maturity = maturity,
     volatility = implied_volatility
  )
  c(eurOp$bsmValue,eurOp$delta,eurOp$gamma,eurOp$vega,eurOp$theta,eurOp$rho,eurOp$divRho)
 }
 ,
 error=function(cond) {return(NA)},
 warning=function(cond) {return(NA)},
 finally =
 {
  #Do Nothing
 })
}

#Calculate Greeks for each row of data in russ2k.opt
  russ2k.opt[, c('c.bsm.value','c.delta','c.gamma','c.vega','c.theta','c.rho','c.divRho')] <- t(apply(russ2k.opt,1,
  function(x)
  FxCalcGreeks(
   type="call",
   underlying=underlying,
   strike=as.numeric(x[["strike"]]),
   dividendYield = 0,
   riskFreeRate=0.0003,
   maturity = as.numeric(x[["expyrs"]]),
     implied_volatility    = as.numeric(x[["call_iv"]]) 
  )
  ))


#Identify the 16,50 & 80 delta points
list.delta.16 <- c()
list.delta.50 <- c()
list.delta.84 <- c()
unk.expdays <- unique(russ2k.opt$expdays)
unk.expdates <- unique(as.Date(russ2k.opt$expdate))
numexp <- length(unk_expdays)
for(i in 1:numexp)
{
  nextexpdays <- unk.expdays[i] 
  russ2k.nextexp <- russ2k.opt[russ2k.opt$expdays == nextexpdays]
  strike.16 <- russ2k.nextexp[which.min(abs(russ2k.nextexp$c.delta - 0.16))]$strike
  strike.50 <- russ2k.nextexp[which.min(abs(russ2k.nextexp$c.delta - 0.5))]$strike
  strike.84 <- russ2k.nextexp[which.min(abs(russ2k.nextexp$c.delta - 0.84))]$strike
  list.delta.16 <- c(list.delta.16, list(strike.16))
  list.delta.50 <- c(list.delta.50, list(strike.50))
  list.delta.84 <- c(list.delta.84, list(strike.84))
}
  xmin <- min(unk.expdates)
  xmax <- max(unk.expdates)
  ymin <- min(unlist(list.delta.84))
  ymax <- max(unlist(list.delta.16))

  plot(unk.expdates, list.delta.16, typ="l", col="green", xlim=c(xmin, xmax), ylim= c(0,ymax),
     main=c("Russell 2000 Option Cone", "Option Cone"), xlab="Expiration Dates", ylab="Strike")
  xticks <- as.character(unk.expdates,format='%m/%y')
  
  tryCatch({ 
      axis(side=1,at=xticks[seq(1,numexp,2)],tcl=0.4,lwd.ticks=1,mgp=c(0,0.25,0), las=2, cex.axis=0.35)
    },
    error= function(cond) {return(NA)},
    warning= function(cond) {return(NA)}
  )
  
  lines(unk.expdates, list.delta.50, col="black")
  lines(unk.expdates, list.delta.84, col="red")
  dev.off()

No comments: