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