[This article was first published on stotastic » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

With the stock market freaking out and all, I figured I should take a look at how volatility was being priced in the option market. The CBOE generously provides snapshots of market data for anyone interested to download. By using this data, we can calculate the markets ‘implied volatility’, or level of ‘freaking out’. For those not familiar with the concept of implied volatility, essentially we can take the prices of options in the market and back out the volatility implied by those prices using the Black-Scholes formula. Its been shown over and over again that the assumptions of the Black-Scholes model don’t hold up to empirical data; but its an easy calculation to perform, and so implied volatility is a widely used metric. Anyway, below is my Black-Scholes option pricing function and the function used to back out implied volatility (written in R of course). Since implied volatility can only be found numerically, I used the Bisection Method to calculate it since it was easy to implement, but there are faster methods out there.

```## Black-Scholes Function
BS <-
function(S, K, T, r, sig, type="C"){
d1 <- (log(S/K) + (r + sig^2/2)*T) / (sig*sqrt(T))
d2 <- d1 - sig*sqrt(T)
if(type=="C"){
value <- S*pnorm(d1) - K*exp(-r*T)*pnorm(d2)
}
if(type=="P"){
value <- K*exp(-r*T)*pnorm(-d2) - S*pnorm(-d1)
}
return(value)
}

## Function to find BS Implied Vol using Bisection Method
implied.vol <-
function(S, K, T, r, market, type){
sig <- 0.20
sig.up <- 1
sig.down <- 0.001
count <- 0
err <- BS(S, K, T, r, sig, type) - market

## repeat until error is sufficiently small or counter hits 1000
while(abs(err) > 0.00001 && count<1000){
if(err < 0){
sig.down <- sig
sig <- (sig.up + sig)/2
}else{
sig.up <- sig
sig <- (sig.down + sig)/2
}
err <- BS(S, K, T, r, sig, type) - market
count <- count + 1
}

## return NA if counter hit 1000
if(count==1000){
return(NA)
}else{
return(sig)
}
}```

To apply the calculations above, I downloaded option data on the S&P 500 (SPX) as of midday Friday, May 21st. Data is available for multiple expiration dates, but I narrowed my analysis to only the near dated contract, which was the June contract. As with all market data, I had to do a bit of cleaning first, but you can find my clean data here. Once cleaned, the data was read into R and cranked though the implied volatility function. Note that the S&P 500 was at 1082 at the time of these quotes.

```## read in data

## calculate implied vol for Call
S <- 1082.74
T <- 28/365
r <- 0.01

n <- dim(dat)[1]
c.vol.Bid <- rep(0,n)
p.vol.Bid <- rep(0,n)

for(i in 1:n){
c.vol.Ask[i] <- implied.vol(S, dat\$K[i], T, r, dat\$C.Ask[i], "C")
c.vol.Bid[i] <- implied.vol(S, dat\$K[i], T, r, dat\$C.Bid[i], "C")
p.vol.Ask[i] <- implied.vol(S, dat\$K[i], T, r, dat\$P.Ask[i], "P")
p.vol.Bid[i] <- implied.vol(S, dat\$K[i], T, r, dat\$P.Bid[i], "P")
}```

Finally the fun part, I plotted the implied volatility vs the strikes for both the calls and puts. I also plotted the open interest (number of open contracts held in the market) vs strike. For reference, I put the current S&P 500 index value (at the time of these quotes) as a vertical dashed red line.

```## plot Volatility Smile and Open Interest
par(mfrow=c(2,2))
plot(dat\$K, c.vol.Bid, typ="l", col="green", xlim=c(800, 1350), ylim= c(.2, .60),
main=c("S&P 500 June Call", "Volatility Smile"), xlab="Strike", ylab="Implied Vol")
abline(v=S, col="red", lty=2)
legend("bottomleft", cex=0.9, c("Ask", "Bid"), lty=1, col=c("blue", "green"))

plot(dat\$K, p.vol.Bid, typ="l", col="green", xlim=c(800, 1350), ylim= c(.2, .6),
main=c("S&P 500 June Put", "Volatility Smile"), xlab="Strike", ylab="Implied Vol")
abline(v=S, col="red", lty=2)
legend("bottomleft", cex=0.9, c("Ask", "Bid"), lty=1, col=c("blue", "green"))

plot(dat\$K, dat\$C.Open, xlim=c(800, 1350), ylim=c(0, 260000), type="h",
main="Open Interest", xlab="Strike", ylab="Open Interest")
abline(v=S, col="red", lty=2)

plot(dat\$K, dat\$P.Open, xlim=c(800, 1350), ylim=c(0, 260000), type="h",
main="Open Interest", xlab="Strike", ylab="Open Interest")
abline(v=S, col="red", lty=2)```

The most obvious thing you should notice is the phenomenon known as the ‘volatility smile’. The assumptions of the Black-Scholes model tell us that volatility should be constant for all strikes. Clearly the market doesn’t believe this assumption and there are numerous theories on why this happens. Something you may have missed, but is rather shocking, most of the put contracts held are at a strike of 800. These are essential insurance contracts for people who are worried about the S&P 500 dropping over 25% within the next month. Scary!