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:
Post a Comment