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()

Making space on boot for Ubuntu

Often we need to create space on the Ubuntu's boot partition. The offending files are older versions of Ubuntu, still sitting on the boot partition. Deleting old versions of Ubuntu is straightforward. Simply, run the following command and boot drive will be cleaned of the old versions.

sudo apt-get autoremove --purge

The above command will delete all, but the current and previous version in the boot partition, allowing you to load new updates of the OS.

However, sometimes deleting the old updates through the previous command is not enough. I ran across this myself when my software updater claimed that I needed to free up still more space. I found the following command on a blog that cleans out all but the current version. Use this only when the previous command is not adequate for your task..

dpkg -l 'linux-*' | sed '/^ii/!d;/'"$(uname -r | sed "s/\(.*\)-\([^0-9]\+\)/\1/")"'/d;s/^[^ ]* [^ ]* \([^ ]*\).*/\1/;/[0-9]/!d' | xargs sudo apt-get -y purge

This did the trick for me, and I was able to install the latest updates.

All the best!