Volatility smile

Following code downloads option price data for a day $d$ from NSE website and plot volatility smile  :)

	library(data.table)
	library(RQuantLib)
	#### For OS x libcurl may not work###
	# This code work without quantmode#
	NSEIndex = "^NSEI"
	# Risk free interest rate fixed at 7.5%
	# Source NSE Mibor
	riskfree = 0.075
	#Specify date in mm/dd/yyyy
	dates = c("01/06/2015")
	tdata <- as.Date(dates, format = "%m/%d/%Y")
	# Find URL link at NSE website for the bhavcopy of the day
	urllink = paste("http://www.nseindia.com/content/historical/DERIVATIVES/",(format(tdata, "%Y")),"/",toupper(format(tdata, "%b")),"/fo",toupper(format(tdata, "%d%b%Y")),"bhav.csv.zip", sep="")
	tempzip = paste(toupper(format(tdata, "%d%b%Y")),"bhav.csv.zip", sep="")
	tempfilename = paste("fo",toupper(format(tdata, "%d%b%Y")), "bhav.csv", sep="")
	# download the compressed file and extract
	download.file(urllink,tempzip, mode="wb" , method = "libcurl")
	unzip(tempzip, tempfilename, overwrite = T)
	optiondf <- read.table(tempfilename, sep=",", header=T)
	print("Downloaded derivative market data for the day")
	print(head(optiondf))
	# Download or read indices' close price
	indexlnk = paste("http://www.nseindia.com/content/indices/ind_close_all_", toupper(format(tdata, "%d%m%Y")),".csv", sep="")
	tempindexzip = paste(toupper(format(tdata, "%d%m%Y")),"index.csv", sep="")
	download.file(indexlnk,tempindexzip, mode="wb" , method = "libcurl")
	indexdf <- read.table(tempindexzip, sep=",", header=T)
	print("Downloaded market index data for the day")
	print(head(indexdf))
	# get NIFTY 50 Index price
	indprice = 0
	# Dataframe for put/call options
	optionindxdf = NULL
	otmput = NULL
	otmcall =NULL
	# find NIFTY index price
	if(length(indexdf)>0) {
	colmhead = "Nifty 50"
	if(colmhead %in% indexdf[, 1]){
	indprice = indexdf[indexdf$Index.Name =="Nifty 50","Closing.Index.Value"]
	}
	else{
	indprice = indexdf[indexdf$Index.Name =="CNX Nifty","Closing.Index.Value"]
	}
	}
	if(length(indexdf)>0) {
	# Near expirty date
	optiondf$EDT = as.Date(optiondf$EXPIRY_DT, "%d-%b-%Y")
	# Days to maturity
	optiondf$dtm = as.numeric(optiondf$EDT - tdata)
	nearexp = min(optiondf[(optiondf$dtm>20) & (optiondf$dtm<100), "EDT"])
	optselindx = (optiondf$INSTRUMENT=="OPTIDX") & (optiondf$SYMBOL=="NIFTY") & (optiondf$CONTRACTS>0) & (optiondf$EDT==nearexp)
	optionindxdf = data.table(optiondf[optselindx, ])
	# Out of money call and put options
	otmput = optionindxdf[(optionindxdf$OPTION_TYP=="PE") & (optionindxdf$STRIKE_PR<indprice), ][order(STRIKE_PR)]
	otmcall = optionindxdf[(optionindxdf$OPTION_TYP=="CE") & (optionindxdf$STRIKE_PR>indprice), ][order(STRIKE_PR)]
	otmput$X = NULL
	otmcall$X = NULL
	print("Call option details:")
	print(head(otmcall))
	print("Put option details:")
	print(head(otmput))
	# Calculate implied volatility
	idx = seq(1, dim(otmput)[1])
	otmput$IV = 0.01
	for(j in idx){
	oval = otmput[j, CLOSE]
	stk = otmput[j, STRIKE_PR]
	# Annualized days
	matrty = otmput[j, dtm]/360.0
	dvdyld = 0
	expvol = 0.1
	civ = EuropeanOptionImpliedVolatility(type="put", value=oval, underlying=indprice, strike=stk, dividendYield=dvdyld, riskFreeRate=riskfree, maturity=matrty, volatility=expvol)
	otmput$IV[j] = civ
	}
	idx = seq(1, dim(otmcall)[1])
	otmcall$IV = 0.01
	for(j in idx){
	oval = otmcall[j, CLOSE]
	stk = otmcall[j, STRIKE_PR]
	matrty = otmcall[j, dtm]/360.0
	dvdyld = 0
	expvol = 0.
	civ = EuropeanOptionImpliedVolatility(type="call", value=oval, underlying=indprice, strike=stk, dividendYield=dvdyld, riskFreeRate=riskfree, maturity=matrty, volatility=expvol)
	otmcall$IV[j] = civ
	}
	# Plot the volatility smile/smirk
	smirkval = c(otmput$IV, otmcall$IV)
	smirkstk = c(otmput$STRIKE_PR, otmcall$STRIKE_PR)
	plot(smirkstk, smirkval, xlab="Strike", ylab= "Implied volatility", main ="Volatility smile", type="l", col = "dark red")
	abline(v=indprice, col="blue")
	text(indprice, 0.2, "Stock Price")
	}

view raw VolatilitySmirk.R hosted with ❤ by GitHub