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()
alternative link download
