The statistical programming language R is a very powerful tool for statistical and graphical techniques to explore, illustrate, study, and model almost any imaginable mathematical scenario. It is also very popular as a tool for the study of quantitative finance. Here we will take a look at the use of R to explore the Black-Scholes option pricing model and the Cox-Ross-Rubenstein binomial option pricing models by web scrapping option prices and modeling on the data returned.
Options derive their value from an underlying asset. They are of a class of financial instruments that also include forward and futures contracts and swaps. The underlying asset of an option can include stocks, bonds, currencies, futures, and many other assets with some determinable value.
In all options transactions there is a contract where, under agreed upon terms, party A will sell something to party B and party B will buy something from party A if certain conditions are met or occur. There is no guarantee that the transaction to buy or sell the underlying will actually occur as it is based on a number of factors. This creates an asymmetric relationship where one party has the right to buy the option and the other party has the obligation to sell the option if the contractual conditions are met. This means there is no guarantee that the transaction will take place or not at the outset.
Option pricing requires all models to make assumptions regarding the future price movements of the underlying product. This is a nonlinear process. One of the parties can cap their loss or gain or do both and the price behavior of the option is based on the change in a number of input variables.
There are different kinds of models that have been designed to model the price behavior of options. This is necessary so that the buyer of the option contract has a way to know if the premium they are paying is appropriate and the seller of the option contract knows that the premium they are receiving is sufficient compensation for the risk of having to deliver the underlying assets. In this paper we restrict the underlying to common stocks that trade on existing exchanges, such as the NYSE and NASDAQ, which we will simply call stocks.
There are two primary models that are used, and we restrict our discussion to these, the Black-Scholes model or the Black-Scholes-Merton model, which assumes the underlying prices changes continuously and the Cox-Ross-Rubenstein model, aka the binomial model, which assumes that the stock prices change at discrete intervals of time.
We will use the statisical programming language R to examine these two models and illustrate the similarities and differences between them.
The Black-Scholes model operates under as set of assumptions that includes
This means that the change in price of the asset \(dS = S_t-S_{t-1}\) is modeled by some mean, \(\mu\) times the original price, \(S\), times the change in time, \(dt\) plus the volatility, \(\sigma\), of the asset, \(S\), times the asset, \(S\), times the change of the Weiner process, \(dW\). The Weiner process is the assumption that the changes in the value follow a standard normal distribution with mean \(0\) and variance \(t\). Here the mean, \(\mu\), is also called the drift.
The drift and volatility are constant parameters and \(W\) is a standard Wiener process.
We assume the market is efficient and therefore arbitrage free.
The underlying is a common stock that pays no dividends. There are modification to the Black-Scholes model that do allow dividends to be included that we will discuss lightly.
The ability to buy or short sell the underlying in any amounts is possible. This includes fractional transactions.
There are no transaction costs to buy or sell the option contract or the underlying asset.
We know the risk free interest rate and it is constant.
It is obvious that these assumptions have limitations in the real world that need to be considered in real world decisions. Still, this model and its modifications are widely used by industry investors, hedgers, and speculators as a means of valuation and decisioning.
An option is called a European option if it can only be exercised at the time of expiration of the contact while an American option can be exercised at any time up to the expiration data. The story of the naming convention as told by Robert Merton, says that Paul Samuelson coined the terms on the idea that the simpler option that can only be exercised at expiration was more suitable for the “European mind” while the more complex option that can be exercised at any time was more in keeping with the “American mind” that he based on a discussion about his ideas for modifying the Black-Scholes model to allow prior to expiration exercise with a Swiss put and call dealer in 1973, in the pre options exchange days.
A European call option is an option that give the buyer of the option the right to buy the underlying at the agreed upon strike price at expiration.
For an easy to understand example, if I purchase a call option to buy Microsoft (MSFT) for the strike price of \(240\) with an expiration of January, 15th, 2021, the premium as of today, Sunday, December 6th, 2020 has an ask of \(1.13\) for the premium. The actual price of one contract is for a round lot of \(100\) shares for a cost of \(1.13 \times 100 = 113.00\). In practical terms this means that the buyer of this call option, should the price of MSFT exceed \(240\) per share, can exercise his right to buy MSFT at \(240\) per share and the seller must fulfill their obligation to deliver the shares.
If we assume MSFT is trading at \(255\) per share on January 15th, then the call owner can make a profit of \(255-240+1.13 = 16.13\) per share for a total profit of \(16.13 \times 100 = 1613.00\) should the buyer exercise and immediately sell his newly owned share at the \(255\) per share price.
We have a formula for a Black-Scholes-Merton European call option
\[S e^{-\delta (T-t)}N(d_1) - Ke^{-r (T-t)}N(d_2)\]
where \(d_1\) and \(d_2\) are found by
\[d_1=\frac{\text{ ln }(S/K) + (r-\delta+0.5\sigma^2)(T-t)}{\sigma \sqrt{T-t}}\text{ and }d_2 = d_1-\sigma\sqrt{T-t}=\frac{\text{ ln }(S/K)+(r-\delta - 0.5\sigma^2)(T-t)}{\sigma\sqrt{T-t}}\]
Where \(\delta\) is the dividend yield, which for the Black-Scholes model is non-existent to that \(\delta=0\) and \((T-t)\) is the time to expiration and \(N\) denotes the cumulative distribution function (cdf) of the standard normal distribution.
A put option gives the buyer of the option the right to sell the shares of the underlying at the strike price for the premium value of the option. A put on the same underlying with the same strike and expiration date is essentially a mirror image of the corresponding call option in terms of profit or loss.
The Black-Scholes-Mertom put formula is
\[Ke^{-r(T-t)}N(-d_2)-Se^{-\delta (T-t)}N(-d_1)\]
where \(N(-d_1)\) and \((-d_2)\) where \(1-N(d_1)=N(-d_1)\) and \(1-N(d_2) = N(-d_2)\)
R is a statistical programming language with over \(17,000\) different packages for analytical modeling of every imaginal type. There are a several packages specifically for option price modeling and we will restrict ourselves to just a few here.
Consider a call and put option on a stock I have been trading recently, Square, Inc., ticker symbol, \(SQ\). The closing price for Friday, December 4, 2020 was \(208.15\). Let’s consider a call and a put option with an expiration of March 19, 2021 with a striker price of \(210.00\) dollars
We are going to need the volatility of the stock and a risk free rate. We will use the quantmod package to get the stock prices for the last year and to acquire the current 90-day US Treasury Bill yield as the risk free rate. The volatility will be represented by the standard deviation of \(SQ\) which is \(40.01209\) and the risk free rate is \(.08\) percent. The R code used is in the notated code blocks below.
This is a very comprehensive package that allows for very easy “scraping” of pricing data, financial data, economic data, and more, from a variety of sources. Here we will pull in the prices for Square, Inc. and put them in a dataframe where we can perform more detailed analysis of the pricing information.
We will also pull in the 90-day US Treasury bill rates as a proxy for the risk-free rate from the Federal Reserve Economic Database, commonly known as FRED using the same quantmod function, getSymbols.
Using the __getSymbols() function, acquire the price data for \(SQ\) for one year, since December 4, 2019.
getSymbols("SQ",src="yahoo", from = "2019-12-04") # get SQ prices from yahoo finance
## [1] "SQ"
getSymbols("DGS3MO",src = "FRED") #get 90-day T-Bill yield from FRED
## [1] "DGS3MO"
tail(DGS3MO)
## DGS3MO
## 2021-05-27 0.02
## 2021-05-28 0.01
## 2021-05-31 NA
## 2021-06-01 0.02
## 2021-06-02 0.02
## 2021-06-03 0.02
A look at the summary(), head(), tail(), and sd() functions gives us more insight into the price data. The standard deviation of the adjusted closing prices is found by the sd() function. In this case the standard deviation for the past year is \(40.01209\).
head(SQ)
## SQ.Open SQ.High SQ.Low SQ.Close SQ.Volume SQ.Adjusted
## 2019-12-04 67.47 68.430 67.190 67.78 3715100 67.78
## 2019-12-05 67.59 67.890 66.660 67.14 4265000 67.14
## 2019-12-06 67.67 68.323 67.610 67.98 4058400 67.98
## 2019-12-09 67.78 68.060 66.720 67.06 4335100 67.06
## 2019-12-10 67.30 67.430 65.420 65.84 6047200 65.84
## 2019-12-11 65.90 66.070 64.875 65.75 6518000 65.75
tail(SQ)
## SQ.Open SQ.High SQ.Low SQ.Close SQ.Volume SQ.Adjusted
## 2021-05-26 218.04 223.830 215.800 222.34 11914400 222.34
## 2021-05-27 222.38 223.390 217.660 220.90 8555700 220.90
## 2021-05-28 221.88 225.900 221.650 222.52 6779800 222.52
## 2021-06-01 223.54 227.065 218.577 221.95 6902900 221.95
## 2021-06-02 221.33 224.880 219.490 220.41 7147500 220.41
## 2021-06-03 218.87 219.940 211.250 211.43 7231500 211.43
sd(SQ$SQ.Adjusted)
## [1] 70.44053
Next we will find the Black-Scholes model price for the put and/or call in question looking at three different approaches:
Price the call option for \(SQ\) with expiration date “Mar.19.2021” and strike price \(210\) with closing price \(208.15\) as the input for \(S\). The time until expiration from today, Sunday, December 6, 2020 is 103 days. Therefore \(T=103/365 = 0.2821917808\).
We build a function in R called BSMCall that takes in as arguments the variable components of the Black-Scholess-Merton model and then run it to generate a value for the option. A caveat is in order: The Yahoo Finance Implied Volatility (IV) is much higher than the computed annualized volatility (57.23% versus 40.01%) which would be an item for consideration if using this to model real cash money decisions.
S = 208.15 # Stock price
K = 210 # Strike price
Q = 1 # Use Q for terminal time, T is special character?
t = 0.2821917808 # Time until expiration
r = 0.008 # Risk free rate
b = 0 # dividend yield
sigma = .4001209
BSMCall <- function(S, K, r, sigma, delta, Q, t)
{
d1 <- (log(S/K) + (r - delta + sigma^2/2)*(Q-t))/(sigma*sqrt(Q-t))
d2 <- d1 - sigma*sqrt(Q-t)
CallValue <- S*exp(-delta*(Q-t))*pnorm(d1) - K*exp(-r*(Q-t))*pnorm(d2)
}
Plug in the values for the function BSMCall and run the model code to get the price of the call.
Call1 <- BSMCall(208.15, 210, 0.008, 0.4001209, 0.0, 1, 103/365)
Call1
## [1] 27.73699
Plug in the values for the function BSMPut
BSMPut <- function(S, K, r, sigma, delta, T, t)
{
d1 <- (log(S/K) +(r - delta + sigma^2/2)*(T-t))/(sigma*sqrt(T-t))
d2 <- d1 - sigma*sqrt(T-t)
PutValue <- K*exp(-r*(T-t))*pnorm(-d2) - S*exp(-delta*(T-t))*pnorm(-d1)
}
Put1 <- BSMPut(208.15, 210, 0.008, 0.4001209, 0.0, 1, 103/365)
Put1
## [1] 28.38453
The arguments for the plain vanilla (aka, ordinary) Black-Scholes option pricing model in fOptions are
TypeFlag = “c” for call or “p” for put
S = current or starting stock price
X = strike price. The function uses X instead of K.
Time = time to expiration in years
r = risk free interest rate
sigma = volatility of the stock
b = cost of carry rate. Generally simply use risk free rate here
Plug in the argument values used above into the function to output the computed value of the call option for the given arguments.
library(fOptions)
GBSOption(TypeFlag = "c", S = 208.15, X = 210, Time = 103/365, r = 0.008, sigma = 0.4001209, b = 0.008)
##
## Title:
## Black Scholes Option Valuation
##
## Call:
## GBSOption(TypeFlag = "c", S = 208.15, X = 210, Time = 103/365,
## r = 0.008, b = 0.008, sigma = 0.4001209)
##
## Parameters:
## Value:
## TypeFlag c
## S 208.15
## X 210
## Time 0.282191780821918
## r 0.008
## b 0.008
## sigma 0.4001209
##
## Option Price:
## 16.99555
##
## Description:
## Fri Jun 04 17:03:58 2021
Using the same arguments but changing the TypeFlag to “P” we can generate the put value for the arguments.
library(fOptions)
GBSOption(TypeFlag = "p", S = 208.15, X = 210, Time = 103/365, r = 0.008, sigma = 0.4001209, b = 0.008)
##
## Title:
## Black Scholes Option Valuation
##
## Call:
## GBSOption(TypeFlag = "p", S = 208.15, X = 210, Time = 103/365,
## r = 0.008, b = 0.008, sigma = 0.4001209)
##
## Parameters:
## Value:
## TypeFlag p
## S 208.15
## X 210
## Time 0.282191780821918
## r 0.008
## b 0.008
## sigma 0.4001209
##
## Option Price:
## 18.372
##
## Description:
## Fri Jun 04 17:03:58 2021
Using getOptionChain() from quantmod package requires the ticker symbol as the main input. This does not return historical price data and the specific expiration data for an option chain is needed.
This can be looked up manually on Yahoo Finance or by exploring expiration.
optChain <- getOptionChain("SQ", Exp = "2021")
names(optChain)
## [1] "Jun.04.2021" "Jun.11.2021" "Jun.18.2021" "Jun.25.2021" "Jul.02.2021"
## [6] "Jul.09.2021" "Jul.16.2021" "Jul.23.2021" "Aug.20.2021" "Sep.17.2021"
## [11] "Dec.17.2021"
The names in the list are assigned with respect to the expiration date going forward.
optChain1 <- optChain[["Aug.20.2021"]]
class(optChain1)
## [1] "list"
names(optChain1)
## [1] "calls" "puts"
Note that opt1 is also a list, where each item corresponds to the call and put options data, such that calls and puts are in a data.frame
sapply(optChain1, class)
## calls puts
## "data.frame" "data.frame"
To put this information together so that we have the data.frame for the option chain, we use the lapply() function.
# adds a column for option type and another for expiry
lapply(1:length(optChain1), function(i) data.frame(optChain1
[[i]],Type = names(optChain1)[i],
Expiration = "August.20.2021" ,
stringsAsFactors = F) )
optChain1 <- Reduce(rbind,optChain1) #combines the two lists by rows as one data.frame
class(optChain1)
Pull the rows in the data.frame containing the $210 strike price to view the call and put option values.
class(optChain1)
## [1] "data.frame"
call <- optChain1[26:31, ]
put <- optChain1[65:70, ]
call
## Strike Last Chg Bid Ask Vol OI LastTradeTime
## SQ210820C00290000 290 2.07 0.23999989 2.06 2.20 26 120 2021-06-04 12:46:43
## SQ210820C00300000 300 1.60 0.18000007 1.54 1.68 35 245 2021-06-04 15:39:47
## SQ210820C00310000 310 1.21 0.09000003 1.16 1.28 9 101 2021-06-04 15:49:28
## SQ210820C00320000 320 0.96 0.00000000 0.88 0.99 4 107 2021-06-04 13:14:00
## SQ210820C00330000 330 0.82 -0.20999998 0.59 0.82 2 41 2021-06-04 12:24:22
## SQ210820P00100000 100 0.34 0.00000000 0.04 0.25 1 84 2021-05-24 09:32:00
## IV ITM
## SQ210820C00290000 0.4841360 FALSE
## SQ210820C00300000 0.4898733 FALSE
## SQ210820C00310000 0.4948781 FALSE
## SQ210820C00320000 0.5009816 FALSE
## SQ210820C00330000 0.5125781 FALSE
## SQ210820P00100000 0.6972687 FALSE
put
## Strike Last Chg Bid Ask Vol OI LastTradeTime IV ITM
## NA NA NA NA NA NA NA NA <NA> NA NA
## NA.1 NA NA NA NA NA NA NA <NA> NA NA
## NA.2 NA NA NA NA NA NA NA <NA> NA NA
## NA.3 NA NA NA NA NA NA NA <NA> NA NA
## NA.4 NA NA NA NA NA NA NA <NA> NA NA
## NA.5 NA NA NA NA NA NA NA <NA> NA NA
The option prices found by these three approaches all had differing results. The Yahoo Finance prices in the table are the ask prices.
GUIDE is a gui derivatives package in R that provides an interactive interface when run in RStudio or some other IDE that runs R.
By simply running GUIDE in a code block a pop-up window, separate from RStudio, with various parameter input options appears. This window has a variety of input options for quantitative finance related instruments including the Black-Scholes model. I am not sure that it is functioning properly as it gave me very odd prices but it is an interesting package I came across that is worth taking a look at.
GUIDE()
We called this the binomial model in our class and I will take the liberty of using that terminology moving forward.
The binomial model assumes the price of the underlying asset follows a discrete binomial process. Prices either rise in a period, which we denote as \(u\), or they fall in a period which we denote as \(d\). Both \(u\) and \(d\) are fixed multipliers in the model that measure the underlying price changes as it goes either up or down. Note that \(u = 1/d\) and \(d=1/u\) in a reciprocal relationship. There is also an assumption that the tree is recombining, i.e. \(Sud = Sdu\).
The model requires a decision as to how many steps or periods, \(n\) are being used. This determines how many steps the time to maturity will be divided into and the length of one time step, \(h\). This manifests at \(\Delta t = \frac{T-t}{n} = h\) where \(T\) is total time, usually one year, \(t\) is the time to maturity, and \(n\) is the number of steps.
We compute \(u\) by \(u=e^{r(T-t)+\sigma\sqrt{(T-t)}}\) and \(d=e^{r(T-t)-\sigma\sqrt{(T-t)}}\). The tree is built by applying the up factor and down factor to find the next node prices of the underlying for as many periods as are required.
At maturity, we compute the option prices at the final nodes using Call \(=max(S_T-K, 0)\) and Put \(=max(K-S_T, 0)\). Working backwards to the starting node, we compute the option prices using the risk-neutral probability that allows us to discount with the risk free rate using \(p^*=\frac{e^-r(T-t)-d}{u-d}\).
To find the option value at the prior node, and letting C = P (calls and puts) in the formula being interchangeable, we use the formula
\[C=e^{-rh}\big[p^*C_u+(1-p^*)C_d\big]\]
Where \(C\) is the call (or put) option in the node that feeds into the call values, \(C_u\) and \(C_d\). at each node, it is the two further nodes option values that provide the information for the prior node. The process repeats until the value of the option at the original node is computed.
The “brute force” method in R involves a great deal of code so we will bypass that process for this model and use an established R package.
The R package fOptions has a function, CRRBinomialTreeOption() that will price a put or call option and build the binomial tree for as many time steps as specified.
The arguments for the binomial option pricing model in fOptions, CRRBinomialTreeOption() are
TypeFlag = “ce” for European call or “pe” for European put, “ca” for American call, “pa” for American put
S = current or starting stock price
X = strike price. The function uses X instead of K.
Time = time to expiration in years
r = risk free interest rate
b = cost of carry rate. Generally simply use risk free rate here
sigma = volatility of the stock
n = number of periods
We will use the model to find the price for the \(SQ\) March 19, 2021 options with strike price \(210\) as we did above.
Using just \(n=5\) periods, initially, we will first price the call.
CRRBinomialTreeOption(TypeFlag = "ce", S = 208.15, X = 210, Time = 103/365, r = 0.008, b = 0.008, sigma = .400121, n =5)
##
## Title:
## CRR Binomial Tree Option
##
## Call:
## CRRBinomialTreeOption(TypeFlag = "ce", S = 208.15, X = 210, Time = 103/365,
## r = 0.008, b = 0.008, sigma = 0.400121, n = 5)
##
## Parameters:
## Value:
## TypeFlag ce
## S 208.15
## X 210
## Time 0.282191780821918
## r 0.008
## b 0.008
## sigma 0.400121
## n 5
##
## Option Price:
## 17.87306
##
## Description:
## Fri Jun 04 17:04:03 2021
To price a European Put we only need to change the TypeFlag to “pe”.
CRRBinomialTreeOption(TypeFlag = "pe", S = 208.15, X = 210, Time = 103/365, r = 0.008, b = 0.008, sigma = .400121, n =5)
##
## Title:
## CRR Binomial Tree Option
##
## Call:
## CRRBinomialTreeOption(TypeFlag = "pe", S = 208.15, X = 210, Time = 103/365,
## r = 0.008, b = 0.008, sigma = 0.400121, n = 5)
##
## Parameters:
## Value:
## TypeFlag pe
## S 208.15
## X 210
## Time 0.282191780821918
## r 0.008
## b 0.008
## sigma 0.400121
## n 5
##
## Option Price:
## 19.24951
##
## Description:
## Fri Jun 04 17:04:03 2021
As the number of periods, \(n\), increases, the model will converge to the Black-Scholes model which we can illustrate by increasing \(n\) and running the .
First, we will look at the plot of the two trees with \(n=5\) that we just generated.
BinTreeEuroCall <- BinomialTreeOption(TypeFlag = "ce", S = 208.15, X = 210, Time = 103/365, r = 0.008, b = 0.008, n =5, sigma = .400121)
#?BinomialTreePlot
BinomialTreePlot(BinTreeEuroCall, dy = 1.0, xlab = "Time Periods", ylab = "Number of Steps", main = "European Call Option Tree")
We can illustrate a tree for the European Put option for our inputs for \(SQ\).
BinTreeEuroPut <- BinomialTreeOption(TypeFlag = "pe", S = 208.15, X = 210, Time = 103/365, r = 0.008, b = 0.008, n =5, sigma = .400121)
#?BinomialTreePlot
BinomialTreePlot(BinTreeEuroPut, dy = 1.0, xlab = "Time Periods", ylab = "Number of Steps", main = "European Put Option Tree")
To show the convergence to the Black-Scholes model, recall that the price generated from the GBSOption() function in the fOptions package was \(16.995\) for the call. If we increase the number of periods to \(n=50\) we see the call option price has essentially converged.
CRRBinomialTreeOption(TypeFlag = "ce", S = 208.15, X = 210, Time = 103/365, r = 0.008, b = 0.008, sigma = .400121, n =50)
##
## Title:
## CRR Binomial Tree Option
##
## Call:
## CRRBinomialTreeOption(TypeFlag = "ce", S = 208.15, X = 210, Time = 103/365,
## r = 0.008, b = 0.008, sigma = 0.400121, n = 50)
##
## Parameters:
## Value:
## TypeFlag ce
## S 208.15
## X 210
## Time 0.282191780821918
## r 0.008
## b 0.008
## sigma 0.400121
## n 50
##
## Option Price:
## 16.99531
##
## Description:
## Fri Jun 04 17:04:04 2021
Both models have a lot of assumptions that need to hold and both are similar in many ways as well. If we can hedge the option perfectly by holding the appropriate amount of the underlying asset we can create a risk-free portfolio. The yield of the risk-free portfolio must equal the risk-free rate.
The correct hedging ratio is the delta of the option which is the partial derivative of the option with respect to the underlying price.
The value of the derivative is the risk-neutral expectation of the value of future possible values, discounted back by the risk-free rate.
As the number of steps in the binomial model increases, the difference between the two models goes away and the option values converge. It is worth noting that the binomial model prices changes are discrete while price changes in the Black-Scholes model are continuous.
We can illustrate the convergence of the two models for a European call using the following plot:
crrprices <- sapply(1:200, function(n){
CRRBinomialTreeOption(TypeFlag = "ce",
S = 208.15, X = 210,
Time = 103/365,
r = 0.008,
b = 0.008,
sigma = 0.400121,
n = n)@price
})
bsprices <- GBSOption(TypeFlag = "c", S = 208.15, X = 210, Time = 103/365, r = 0.008, sigma = 0.400121, b = 0.008)@price
plot(1:200, crrprices, type = 'l', xlab = "Number of Steps", ylab = "Prices")
abline(h = bsprices, col = "red")
legend("bottomright", legend = c('CRR-price', 'BS-price'), col = c("black", "red"), pch = 19)
We can see from the plot that the CRR Binomial pricing model converges to the Black-Scholes pricing model for this call option by about 120 steps.
Running the same for the put results in what appears to be an indentical convergence rate.
crrprices <- sapply(1:200, function(n){
CRRBinomialTreeOption(TypeFlag = "pe",
S = 208.15, X = 210,
Time = 103/365,
r = 0.008,
b = 0.008,
sigma = 0.400121,
n = n)@price
})
bsprices <- GBSOption(TypeFlag = "p", S = 208.15, X = 210, Time = 103/365, r = 0.008, sigma = 0.400121, b = 0.008)@price
plot(1:200, crrprices, type = 'l', xlab = "Number of Steps", ylab = "Prices")
abline(h = bsprices, col = "red")
legend("bottomright", legend = c('CRR-price', 'BS-price'), col = c("black", "red"), pch = 19)
Website for R package quantmod. https://www.quantmod.com/
CRAN Repository .pdf for quantmod package: https://cran.r-project.org/web/packages/quantmod/quantmod.pdf
Relevant YouTube Video Series on Using Quantmod, https://www.youtube.com/playlist?list=PL_-iSWCyeCajX6FNjeZ0W2LgN61KNKDqC
CRAN Repository package guide for GUIDE package in R, https://cran.r-project.org/web/packages/GUIDE/index.html
Website with R tutorial tips on using quantmod functions in R, https://ntguardian.wordpress.com/2017/03/27/introduction-stock-market-data-r-1/
A. Penati, G. Pennacchi, The Cox-Ross_rubenstein Option Pricing Model, .pdf file of paper hosted at http://home.cerge-ei.cz/petrz/FM/f400n10.pdf
A. M. Belachew, Dr. M.B.Martin, Dr. V. R. Kota Convergence of the multi-step binomial market model to the Black-Scholes financial model International Journal of Advance Research in Science and Engineering, V7, Issue 2, February 2018
Trenca, Ioan & Pochea, Miruna & FILIP, Angela. (2010). Options evaluation - Black-Scholes model vs. binomial options pricing model. Finance : Challenges of the Future.
Makamo, Event. (2018). Binomial Models And Option Pricing valuation in R. University of South Africa, ResearchGate
Derivative Markets, 3rd Edition, McDonald, Addison Wesly, 2012, ISBN-10: 0321543084