Finding the market implied volatility

Recall the Black-Schole’s formula \[V = BS(S, K, r, T, \sigma)\] Now we have an equation where we only have one unknown variable, volatility. We could solve for the implied volatility and then find out what that would value the option at. We could do so using a root finding method like the Newton’s method.

Newton’s method is a method for finding increasingly improved approximations to the roots of a function.

If we have a function \(f(x)\) and its derivative \(f′(x)\), we can start with an initial guess and successively improve it with an updated guess by doing: \[x_{n+1} = x_n - \frac{f(x_n)}{f'(x_{n})}\]

Here, our \(f(x)\) function is the Black-Scholes equation for pricing an option: \[f(x) = BS(\sigma)\]

find_vol <- function(target_value, call_put, S, K, T, r) {
    MAX_ITERATIONS <- 100
    PRECISION <- 1.0e-5
    call_put <- tolower(call_put)
    
    sigma <-  0.5
    for (i in 1: MAX_ITERATIONS) {
        price <- bs_price(call_put, S, K, T, r, sigma)
        vega <- bs_vega(call_put, S, K, T, r, sigma)
        diff <- target_value - price  # our root
        if (abs(diff) < PRECISION | vega < PRECISION) {
            return(sigma)
        }
        sigma <- sigma + diff/vega # f(x) / f'(x)
    }
    return(sigma)
}    

n = dnorm
N = pnorm

bs_price <- function(cp_flag,S,K,T,r,v,q=0.0) {
    d1 <- (log(S/K)+(r+v*v/2.)*T)/(v*sqrt(T))
    d2 <- d1-v*sqrt(T)
    if (cp_flag == 'c') {
        price <- S*exp(-q*T)*N(d1)-K*exp(-r*T)*N(d2)
    }
    else {
        price <- K*exp(-r*T)*N(-d2)-S*exp(-q*T)*N(-d1)
    }
    return(price)
}
bs_vega <- function(cp_flag,S,K,T,r,v,q=0.0) {
    d1 <- (log(S/K)+(r+v*v/2.)*T)/(v*sqrt(T))
    return(S * sqrt(T)*n(d1))
}    

Volatility Smile for EURO STOXX 50 Index Dividend Options

First the data of EURO STOXX 50 Index Dividend Options and Futures are imported,

library(readr)
library(dplyr)
data <- read_csv("2016_FEXD Settlement Price by Product and Expiration Date%281%29.csv", skip = 9)
data <- data[1:9]
names_list = c('day', 'id', 'product', 'expiration', 'symbol', 'strike', 'exercise', 'metrics', 'price')
colnames(data) <- names_list
data <- data.frame(data)

Next we filter the data at specified date with fixed expiration month, then we select the data with specified option type, and compute its implied volatility using Newton’s method,

volatility_smile <- function(data, current_day, expiration_month) {
    data_by_day <- filter(data, day == current_day, 
                      expiration == expiration_month)
    r <- 0
    F <- filter(data_by_day, symbol == 'F')$price
    cp <- 'C'
    T <- as.numeric(as.Date(expiration_day, 
                "%m/%d/%y") - as.Date(current_day, "%m/%d/%y"), 
                units = "days")/365.
    data_by_day <-filter(data_by_day, symbol == cp)
    data_by_day$implied_volatility <- 0
    data_by_day$Moneyness <- 0
    for (i in 1:nrow(data_by_day)) {
            sigma <- find_vol(data_by_day$price[i], 
                              cp, S = F*exp(-r*T),
                              K = data_by_day$strike[i], T, r)
            data_by_day$implied_volatility[i] <- sigma
            data_by_day$Moneyness[i] <- data_by_day$strike[i]/F
    }
    return(data_by_day)
}

Finally, we plot the implied volatility,

expiration_month <- "201612"
expiration_day <- "12/23/2016"  # the third Friday of December

library(ggplot2)
g <- ggplot()
current_day <- "10/05/2016"
data_by_day1 <- volatility_smile(data = data, current_day, expiration_month)
g <- g + geom_point(aes(data_by_day1$Moneyness, 
                   data_by_day1$implied_volatility), colour = "red")
current_day <- "08/05/2016"
data_by_day2 <- volatility_smile(data = data, current_day, expiration_month)
g <- g + geom_point(aes(data_by_day2$Moneyness, 
                   data_by_day2$implied_volatility), colour = "blue")
current_day <- "06/06/2016"
data_by_day3 <- volatility_smile(data = data, current_day, expiration_month)
g <- g + geom_point(aes(data_by_day3$Moneyness, 
                   data_by_day3$implied_volatility), colour = "green")
current_day <- "04/05/2016"
data_by_day4 <- volatility_smile(data = data, current_day, expiration_month)
g <- g + geom_point(aes(data_by_day4$Moneyness, 
                   data_by_day4$implied_volatility), colour = "gold")
current_day <- "02/05/2016"
data_by_day5 <- volatility_smile(data = data, current_day, expiration_month)
g <- g + geom_point(aes(data_by_day5$Moneyness, 
                   data_by_day5$implied_volatility), colour = "violet")
g <- g + labs(x = 'Moneyness = K/F', y = 'Implied Volatility')
g