Question1
Saeed Rahman
March 7, 2017
A)
Construct code to calculate option values using an additive binomial tree.For this part you need for versions European and American as well as Call and Put. You may use the same tree construction for all options.
BinomialTree = function(isCall, isAmerican=FALSE, K=100, Tm=1,
S0=100, r=0.06, sig=0.2, N=3,div=0,show=FALSE)
{
# Precompute constants ----
dt = Tm/N
nu=r-div-0.5*sig*sig
dxu=sqrt(sig*sig*dt+((nu*dt)^2))
dxd=-dxu
pu=0.5+0.5*(nu*dt/dxu)
pd=1-pu
disc=exp(-r*dt)
nRows = 2*N+1 #number of rows for the matrix
nCols = N+1 #number of columns
cp = ifelse(isCall, 1, -1) # to check ifts a call or a put
# Intialize asset prices ----
# Creating a matrix of nRows*nColumns with zeros and headings
V = S = matrix(0, nrow=nRows, ncol=nCols, dimnames=list(
paste("NumUps", N:-N, sep="="), paste("T", 0:N, sep="=")))
S[nCols, 1] = S0 # initial stock price
# iterating the elements of the matrix in a conical manner starting
# from the position of initial stock price
# For n=3, S[i,j]= S0, then update the forward diagonal elements
# Code is similar to the one used for trinomial tree
for (j in 1:N) {
for(i in (nCols-j+1):(nCols+j-1)) {
S[i-1, j+1] = S[i, j]*exp(dxu)
S[i+1, j+1] = S[i, j] *exp(dxd)
}
}
for (i in 1:nRows) {
V[i, nCols] = max( 0, cp * (S[i, nCols]-K))
}
# Step backwards through the tree ----
for (j in (nCols-1):1) {
for(i in (nCols-j+1):(nCols+j-1)) {
# V[i, j] = disc * (p*V[i-1,j+1] + (1-p)*V[i+1,j+1])
V[i, j] = disc * (pu*V[i-1,j+1] + pd*V[i+1,j+1])
if(isAmerican) {
# if american option, then take the Value at each node as the max of the
# value of option or the payoff at that period
V[i, j] = max(V[i, j], cp * (S[i, j] - K))
}
}
}
if(show)
{
print("Stock Tree")
print(S)
print("Option Value Tree")
print(V)
}
else
return(V[nCols,1])
}
#To display the matrix of Stock process and Option Value
BinomialTree(isCall=TRUE,K=100,Tm =1 ,S0 =100 ,sig = .2,N =4,r=.06,show = TRUE)## [1] "Stock Tree"
## T=0 T=1 T=2 T=3 T=4
## NumUps=4 0 0.00000 0.00000 0.00000 149.48039
## NumUps=3 0 0.00000 0.00000 135.18801 0.00000
## NumUps=2 0 0.00000 122.26217 0.00000 122.26217
## NumUps=1 0 110.57223 0.00000 110.57223 0.00000
## NumUps=0 100 0.00000 100.00000 0.00000 100.00000
## NumUps=-1 0 90.43862 0.00000 90.43862 0.00000
## NumUps=-2 0 0.00000 81.79145 0.00000 81.79145
## NumUps=-3 0 0.00000 0.00000 73.97106 0.00000
## NumUps=-4 0 0.00000 0.00000 0.00000 66.89841
## [1] "Option Value Tree"
## T=0 T=1 T=2 T=3 T=4
## NumUps=4 0.00000 0.000000 0.000000 0.00000 49.48039
## NumUps=3 0.00000 0.000000 0.000000 36.67122 0.00000
## NumUps=2 0.00000 0.000000 25.207510 0.00000 22.26217
## NumUps=1 0.00000 16.547632 0.000000 12.05646 0.00000
## NumUps=0 10.53007 0.000000 6.529383 0.00000 0.00000
## NumUps=-1 0.00000 3.536099 0.000000 0.00000 0.00000
## NumUps=-2 0.00000 0.000000 0.000000 0.00000 0.00000
## NumUps=-3 0.00000 0.000000 0.000000 0.00000 0.00000
## NumUps=-4 0.00000 0.000000 0.000000 0.00000 0.00000print(paste("European Call Price",round(BinomialTree(isCall=TRUE,K=100,Tm =1 ,S0 =100 ,sig = .2,N=200,r=.06),4)))## [1] "European Call Price 10.98"print(paste("European Put Price",round(BinomialTree(isCall=FALSE,K=100,Tm =1 ,S0 =100 ,sig = .2,N =200,r=.06),4)))## [1] "European Put Price 5.1568"print(paste("American Call Price",round(BinomialTree(isCall=TRUE,isAmerican = TRUE,K=100,Tm =1 ,S0 =100 ,sig = .2,N=200,r=.06),4)))## [1] "American Call Price 10.98"print(paste("American Put Price",round(BinomialTree(isCall=FALSE,isAmerican = TRUE,K=100,Tm =1 ,S0 =100 ,sig = .2,N =200,r=.06),4)))## [1] "American Put Price 5.796"B)
Download Option prices (you can use the Bloomberg Terminal, Yahoo! Fi- nance, etc.) for an equity, for 3 different maturities (1 month, 2 months, and 3 months) and 20 strike prices close to the value at the money. If 3 months does not exist use next one available. Please download the data DURING THE TRADING DAY (9:00am to 4:30pm ET). Otherwise your values will be way off. Do not forget to download the value of the underlying. For each strike price in the data, use the implied vol values in Homework 1 (see Problem 1c) and the current short-term interest rate (today the Feds Fund rate is r = 0.75%). Calculate the option price (European Calls and Puts) using the binomial tree, and compare the results with the Black Scholes price. Use at least 200 steps in your tree construction. Treat the options as American as well and plot these values side by side with the European and Black Scholes values. When you create the plot do not forget to plot the bid-ask values as well
Calculating Implied Volatility
Importing the data from the CSV files and calculating the IV
option_chain_call_csv <- read.csv(file="call.csv",header=TRUE, sep=",")
option_chain_put_csv <-read.csv(file="put.csv",header=TRUE, sep=",")
#call.csv and put.csv are available in the project folder
#Call
option_chain_call_csv$days_till_expiry <- as.Date(option_chain_call_csv$Expiry,"%m/%d/%Y")-as.Date("2017-02-19")
#Put
option_chain_put_csv$days_till_expiry <- as.Date(option_chain_put_csv$Expiry,"%Y/%m/%d")-as.Date("2017-02-19")
#Calculating the days till expiry
option_chain_call_csv$premium<-(option_chain_call_csv$Bid+option_chain_call_csv$Ask)/2
option_chain_put_csv$premium<-(option_chain_put_csv$Bid+option_chain_put_csv$Ask)/2
df=ImpliedVol_Secant(option_chain = option_chain_call_csv)
df=as.data.frame(df[1])
options.df_call=df[complete.cases(df$Implied_Volatility),]
df=ImpliedVol_Secant(option_chain = option_chain_put_csv)
df=as.data.frame(df[1])
options.df_put=df[complete.cases(df$Implied_Volatility),]Calculating the Option Prices Using Binomial Tree and Black Sholes Merton Formula
PriceCalculator<-function(options.df,show=FALSE){
eur_binomial <-{}
eur_bsm <-{}
ame_binomial <-{}
stock_df<-as.data.frame(getSymbols("AAPL",from = as.Date("2017-01-01"),
to=as.Date("2017-02-19"), env = NULL))
for (i in 1:nrow(options.df)){
eur_binomial <- append(eur_binomial,BinomialTree(
isCall = as.logical(options.df[i,"Type"]=="Call"),
K = as.numeric(options.df[i,"Strike"]),
Tm = as.numeric(options.df[i,"Days_till_Expiry"])/252,
S0 = as.numeric(tail(stock_df,1)[6]),
sig = as.numeric(options.df[i,"Implied_Volatility"]),
r = .75/100,
N = 200))
ame_binomial <- append(ame_binomial,BinomialTree(
isAmerican = TRUE,
isCall = as.logical(options.df[i,"Type"]=="Call"),
K = as.numeric(options.df[i,"Strike"]),
Tm = as.numeric(options.df[i,"Days_till_Expiry"])/252,
S0 = as.numeric(tail(stock_df,1)[6]),
sig = as.numeric(options.df[i,"Implied_Volatility"]),
r = .75/100,
N = 200))
eur_bsm <- append(eur_bsm,BSM(
type = ifelse(options.df[i,"Type"]=="Call",'c','p'),
K = as.numeric(options.df[i,"Strike"]),
t = as.numeric(options.df[i,"Days_till_Expiry"])/252,
S = as.numeric(tail(stock_df,1)[6]),
sigma = as.numeric(options.df[i,"Implied_Volatility"]),
r = .75/100))
}
options.df$BinomialEuropean <- eur_binomial
options.df$BinomialAmerican <- ame_binomial
options.df$BlackSholesMertonEuropean <- eur_bsm
if(show)
head(options.df)
else
return(options.df)
}
PriceCalculator(options.df_call,show=TRUE)## Days_till_Expiry Type Specification
## 8 5 Call 130 - Call Expiring On: 2/24/2017
## 9 5 Call 131 - Call Expiring On: 2/24/2017
## 10 5 Call 132 - Call Expiring On: 2/24/2017
## 11 5 Call 133 - Call Expiring On: 2/24/2017
## 12 5 Call 134 - Call Expiring On: 2/24/2017
## 13 5 Call 135 - Call Expiring On: 2/24/2017
## Implied_Volatility Strike Bid Ask BinomialEuropean BinomialAmerican
## 8 0.1360770 130 5.55 5.95 5.749815 5.749815
## 9 0.1146620 131 4.60 4.90 4.749936 4.749936
## 10 0.1267717 132 3.65 3.95 3.799360 3.799360
## 11 0.1373564 133 2.87 3.00 2.935358 2.935358
## 12 0.1350209 134 2.09 2.15 2.120356 2.120356
## 13 0.1337578 135 1.41 1.45 1.428714 1.428714
## BlackSholesMertonEuropean
## 8 5.749994
## 9 4.749993
## 10 3.799978
## 11 2.934957
## 12 2.119940
## 13 1.429927option_chain_call <-PriceCalculator(options.df_call)
option_chain_put <-PriceCalculator(options.df_put)Plotting the Call Option Prices that are expiring in 5 days
Zoomin to see the difference in prices. Also select legend to deactivate and activate individual series of option prices
# {r,echo=FALSE,results='asis',comment=NA, echo=FALSE, message=FALSE, warning=FALSE, comment=NA}
option_chain_call$Implied_Volatility<-NULL
df_split <-split(option_chain_call,option_chain_call$Days_till_Expiry)
first <-df_split[1]
first <- first$`5`
first[1:2] <- list(NULL)
rownames(first) <- first$Strike
first$Strike <- NULL
a <- Highcharts$new()
a$chart(type = "line")
a$chart(zoomType="xy")
a$title(text = "Apple Call Expiring in 5 Days")
a$xAxis(categories = rownames(first),title = list(text = "Strike"),replace = T)
a$yAxis(title = list(text = "Option Price"))
# a$xAxis(title = "Strikes")
a$data(first)
a$print(include_assets = TRUE)Plotting the Call Option Prices Expiring in 26 days
# {r,echo=FALSE,results='asis',comment=NA, echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results='asis'}
df_split <-split(option_chain_call,option_chain_call$Days_till_Expiry)
first <-df_split[2]
first <- first$`26`
first[1:2] <- list(NULL)
rownames(first) <- first$Strike
first$Strike <- NULL
a <- Highcharts$new()
a$chart(type = "line")
a$chart(zoomType="xy")
a$title(text = "Apple Call Expiring in 26 Days")
a$xAxis(categories = rownames(first),title = list(text = "Strike"),replace = T)
a$yAxis(title = list(text = "Option Price"))
# a$xAxis(title = "Strikes")
a$data(first)
a$print(include_assets = TRUE)Plotting the Option Prices for the Call Expiring in 152 days
df_split <-split(option_chain_call,option_chain_call$Days_till_Expiry)
first <-df_split[3]
first <- first$`152`
first[1:2] <- list(NULL)
rownames(first) <- first$Strike
first$Strike <- NULL
a <- Highcharts$new()
a$chart(type = "line")
a$chart(zoomType="xy")
a$title(text = "Apple Call Expiring in 152 Days")
a$xAxis(categories = rownames(first),title = list(text = "Strike"),replace = T)
a$yAxis(title = list(text = "Option Price"))
# a$xAxis(title = "Strikes")
a$data(first)
a$print(include_assets = TRUE)Plotting the Put Options Prices that are expiring in 5 days
Zoomin to see the difference in prices. Also select legend to deactivate and activate individual series of option prices
option_chain_put$Implied_Volatility<-NULL
df_split <-split(option_chain_put,option_chain_put$Days_till_Expiry)
first <-df_split[1]
first <- first$`5`
first[1:2] <- list(NULL)
rownames(first) <- first$Strike
first$Strike <- NULL
a <- Highcharts$new()
a$chart(type = "line")
a$chart(zoomType="xy")
a$title(text = "Apple Put Options Expiring in 5 Days")
a$xAxis(categories = rownames(first),title = list(text = "Strike"),replace = T)
a$yAxis(title = list(text = "Option Price"))
a$data(first)
a$print(include_assets = TRUE)Plotting the Option Prices for the Puts Expiring in 26 days
df_split <-split(option_chain_put,option_chain_put$Days_till_Expiry)
first <-df_split[2]
first <- first$`26`
first[1:2] <- list(NULL)
rownames(first) <- first$Strike
first$Strike <- NULL
a <- Highcharts$new()
a$chart(type = "line")
a$chart(zoomType="xy")
a$title(text = "Apple Put Options Expiring in 26 Days")
a$xAxis(categories = rownames(first),title = list(text = "Strike"),replace = T)
a$yAxis(title = list(text = "Option Price"))
a$data(first)
a$print(include_assets = TRUE)Plotting the Option Prices for the Puts Expiring in 152 days
df_split <-split(option_chain_put,option_chain_put$Days_till_Expiry)
first <-df_split[3]
first <- first$`152`
first[1:2] <- list(NULL)
rownames(first) <- first$Strike
first$Strike <- NULL
a <- Highcharts$new()
a$chart(type = "line")
a$chart(zoomType="xy")
a$title(text = "Apple Put Options Expiring in 152 Days")
a$xAxis(categories = rownames(first),title = list(text = "Strike"),replace = T)
a$yAxis(title = list(text = "Option Price"))
a$data(first)
a$print(include_assets = TRUE)C)
The table obtained above in the previous question clearly shows that even though there exist a small error in the option prices obtained, the binomial tree(simulation method) closely predicts the price as per the Analytical equation, that is the Black sholes equation. We can also validate the accuracy of the Binomial tree and Blacksholes model used in this program by observing the small error in price difference between the calculated and the original bid and ask values.
D)
QuestionD<-function(){
bsm_price <- {}
binomial_price <-{}
error <-{}
iter <-c(10, 20, 30, 40, 50, 100, 150, 200, 250, 300, 350, 400)
for(i in iter){
binomial_temp <-BinomialTree(isCall=FALSE,K=100,Tm =1 ,S0 =100 ,sig = .2,N =i,r=.06)
bsm_temp <-BSM(S=100,K=100,t=1,r=.06,sigma=.2,type="p")
binomial_price <- append(binomial_price,binomial_temp)
bsm_price <- append(bsm_price,bsm_temp)
error <- append(error,bsm_temp-binomial_temp)
}
print(paste("Error=",tail(error,n=1)))
plot(iter,error,xlab = "Number of Iterations",ylab = "Error", type="o",col="blue")
}
QuestionD()## [1] "Error= 0.00460598068031715"
It is very clear that that as the number of steps in the tree increases, the error comes down. And the error is lowest at the highest number of steps which is 400
Bonus Question
Using the binomial tree for American Calls and Puts, calculate the implied volatility corresponding to the data you have downloaded in part (b). You will need to use the bisection or Newton/secant method of finding roots with the respective binomial trees. Compare these values of the implied volatility with the volatilities calculated using the usual Black Scholes formula (as in Homework 1, Problem 1c). Write detailed observations.
Secant_Binomial <- function(S, K, t, r, type, option_price
, x0=0.1, x1=3, tolerance=1e-07, max.iter=10000){
x1=3
theta=.00001
fun.x1=BinomialTree(isCall=is.logical(type=="Call"),K=K,Tm =t ,S0 =S ,sig = x1,N =200,r=r)-option_price
count=1
start.time <- Sys.time()
while(abs(fun.x1) > tolerance && count<max.iter) {
x2=x1-theta
fun.x1=BinomialTree(isCall=is.logical(type=="Call"),K=K,Tm =t ,S0 =S ,sig = x1,N =200,r=r)-option_price
fun.x2=BinomialTree(isCall=is.logical(type=="Call"),K=K,Tm =t ,S0 =S ,sig = x2,N =200,r=r)-option_price
x1 <- x1- fun.x1/((fun.x1-fun.x2)/theta)
count <-count+1
print(count)
}
end.time <- Sys.time()
time.taken <- end.time - start.time
if(x2<0 || count>=max.iter)
return(list(NA,time.taken,count))
else
return(list(x2,time.taken,count))
}
ImpliedVol_Secant_Binomial<-function(symbol="AAPL",option_chain,rate=.75/100){
symbol="AAPL"
stock_df<-as.data.frame(getSymbols(symbol,from = as.Date("2017-01-01"),to=as.Date("2017-02-19"), env = NULL))
iv <- {}
original_iv <-{}
optionName <-{}
strike <-{}
days_till_expiry <-{}
time.taken <- 0
iterations <- 0
type <-{}
bid<-{}
ask<-{}
for (i in 1:nrow(option_chain))
{
try({
#SecantMethodUsingBinomial---
secant <- Secant_Binomial(
S = as.numeric(tail(stock_df,1)[6]),
K = as.numeric(option_chain[i,"Strike"]),
t = as.numeric(option_chain[i,"days_till_expiry"])/252,
r = rate,
type = ifelse((option_chain[i,"Type"]=="Call"), "c", "p"),
option_price = as.numeric(option_chain[i,"premium"]))
iv <- append(iv,as.numeric(secant[1]))
if(!is.na(secant[1])){
time.taken <- as.numeric(secant[2])+time.taken
iterations <- as.numeric(secant[3])+iterations
}
type <- append(type,as.character(option_chain[i,"Type"]))
strike<-append(strike,as.numeric(option_chain[i,"Strike"]))
optionName <- append(optionName,paste(option_chain[i,"Strike"],"-",
option_chain[i,"Type"],"Expiring On:",
option_chain[i,"Expiry"]))
days_till_expiry <- append(days_till_expiry,as.numeric(option_chain[i,"days_till_expiry"]))
bid<-append(bid,as.numeric(option_chain[i,"Bid"]))
ask<-append(ask,as.numeric(option_chain[i,"Ask"]))
# original_iv <- append(original_iv,as.numeric(option_chain[i,"Implied.Volatility"]))
})
}
option_chain_df <- data.frame(days_till_expiry,type,optionName,iv,strike,bid,ask)
names(option_chain_df)<-c("Days_till_Expiry","Type","Specification","Implied_Volatility","Strike","Bid","Ask")
time.taken <- time.taken/as.numeric(colSums(!is.na(option_chain_df))[3])
iterations <- iterations/as.numeric(colSums(!is.na(option_chain_df))[3])
list(option_chain_df,time.taken,iterations)
}
df=ImpliedVol_Secant_Binomial(option_chain = option_chain_csv)
df=as.data.frame(df[1])
options.df=df[complete.cases(df$Implied_Volatility),]