Source for this study
+ http://www.theresearchkitchen.com/archives/738 + http://www.frolovs.me/2016/05/26/binomial-pricing-trees-in-R.html.
The intuition behing the binomial tree simulation is to have at each point in time \((t+\Delta t)\) a probability p of prices going up, with the consequent probability \((1 - p)\) of the price going down.
Probability of going up or down remains constant during the whole process. A balanced binomial tree with a height \(h\) will have \(2^{h+1}-1\) nodes, while a binomial lattice with the same height will have \(\sum_{i=1} ^{h}i\). Height here refers to the number of steps the tree has.
To create the binomial lattice we will follow the variables below N <- number of steps of the tree, or simply the height X0 <- starting value u <- size of up price movement d <- size of down price movement
genlattice <- function(X0 = 100, u = 1.1, d = .75, N = 5) {
X <- c()
X[1] <- X0
count <- 2
for (i in 1:N) {
for (j in 0:i) {
X[count] <- X0 * u^j * d^(i-j)
count <- count + 1
}
}
return(X)
}
genlattice(N= 5, u=1.1, d=.9)## [1] 100.000 90.000 110.000 81.000 99.000 121.000 72.900 89.100
## [9] 108.900 133.100 65.610 80.190 98.010 119.790 146.410 59.049
## [17] 72.171 88.209 107.811 131.769 161.051Initially, the genlattice function simply simulates the asset price evolution from node to node based on the up and down probabilities. As a next step, it should also value the option and store its price, so asset and option prices pairs will be displayed in the lattice. Also, it should accept parameters for option valuation scenarios, such as volatility, interest rates, payoff function and so on.
Here is how the code will follow:
1. A public function (genlattice.vanilla.european.call) is defined to work as interface for the whole code. this function refers to 2 other functions (private)
+ payoff.vanilla.call: where the payoff for a call type option is returned. The values feed the option pricing calculation of genlattice.european.
+ genlattice.european: where asset price evolution and option price calculation is performed
genlattice.vanilla.european.call <- function(Asset, Volatility, IntRate, DividendRate, Strike, Expiry, NoSteps) {
return (genlattice.european(Asset=Asset, Volatility=Volatility, IntRate=IntRate, DividendRate=DividendRate, Strike=Strike, Expiry=Expiry, NoSteps=NoSteps, Payoff=payoff.vanilla.call))
}As mentioned before, the genlattice function needs now to accomodate the option price together with its corresponding asset price for each node. Results can no longer be stored in a simple vector as before (X0), so the final set will be placed on a dataframe. For that, one must first compute the amount of nodes in the tree and pre-populate it with NA values.
The code is based on a Cox-Ross-Rubenstein (CRR, 1979) method for pricing the option. In that the up tick movement is adjusted with volatility \(\sigma\) through time.
\(u = e^{\sigma \sqrt{\delta t}}\)
Down tick is the inverse of the up tick
\(d = \frac{1}{u} = e^{-\sigma \sqrt{\delta t}}\)
The synthetic probability (Q-measure) is calculated as
\(p = \frac{a- d}{u - d}\)
being
\(a = e^{-r \delta t}\)
Notice that r here represents the risk free interest rate for a Stock option, and for that the dividend rate must be considered, so that
\(r = IntRate - DividentRate\)
The first part of the code payoff functions for call/put options will be defined to get the boundaries conditions for the equations which will follow to perform the option price calculation.
# Generate a binomial lattice for an European option
payoff.vanilla.call <- function(Asset, Strike){
return( max(0, Asset - Strike))
}
payoff.vanilla.put <- function(Asset, Strike){
return( max(Strike - Asset, 0))
}
NoSteps <- 10
Asset <- 810
Expiry <- 0.5
IntRate <- 0.05
DividendRate <- 0.02
Volatility <- .2
Strike <- 800
genlattice.european <- function(Asset, Volatility, IntRate, DividendRate, Strike, Expiry, NoSteps, Payoff, Type){
# Count the number of nodes on the tree.
# For every step there are NoSteps + 1 node, so
# the total number of nodes would be summation of the
# sequence from 1 until NoSteps +1
count <- sum(1: (NoSteps + 1))
# The objective is to have Asset and Option prices stored in a dataframe. (where are payoffs stored ?)
# The mapping for the tree node (i,j) to a linear index inside the dataframe will need to be defined
X <- data.frame(matrix(NA, nrow = count, ncol = 2))
names(X) <- c("asset", "option")
# Time between price movements
dt = Expiry / NoSteps
# Option price (and Asset) discount factor
DiscountFactor <- exp(-IntRate * dt)
# The up and down ticks with corresponding (synthetic) probabilities (Cox, Ross, Rubenstein Method)
u = exp(Volatility * sqrt(dt))
d = 1/u
a = exp((IntRate - DividendRate) * dt)
p = (a - d) / (u - d)
# Up to now it was simply defining the variables and the rules for the option pricing calculation.
# From now on, for every node asset price and the corresponding option price is calculated.
# Starting from the last node of the tree (bottom right corner) and backwards
for (i in NoSteps:0) {
for (j in i:0) {
X$asset[count] <- Asset * u^(i-j) * d^j
# Compute payoff directly for last step's nodes,
# otherwise use a formula.
# More in detail, what happes is:
# 1) the last column's payoff derives from the function which finds the max between asset and strike prices difference or zero
# These last nodes are identified as those elements on column i from the binominal tree where i = NoSteps
# 2) walking backwards, the value of the option at a given node is the expected value of this option, meaning, the probability of this option going up times the value it will have when it goes up, plus the probability of this option going down times the value it will have when it comes down. Both values are then discounted for the delta time dt between the nodes.
# Since we have identified the option values at the last nodes, for each node position back we need to identify what are the future up and down movements.
if (i == NoSteps) {
X$option[count] <- Payoff(X$asset[count], Strike)
} else {
up <- X$option[sum(1:(i+1), j, 1)]
down <- X$option[sum(1:(i+1), j+1, 1)]
X$option[count] <- DiscountFactor * (p * up + (1-p) * down)
}
count <- count - 1
}
}
return(X)
}The last step is to display the data in the binomial tree The dotlattice function needs to render both asset and option prices to a desired level of precision. It must work with both European and American trees (the latter is different from the former by having an early exercise flag for each node).
# Generates a graph specification that can be fed into graphviz.
# Input: the binomial lattice produced by one of genlattice family functions.
dotlattice <- function(S, digits=2) {
shape <- "plaintext"
cat("digraph G {", "\n", sep="")
cat("node[shape=",shape,"];","\n", sep="")
cat("rankdir=LR;","\n")
cat("edge[arrowhead=none];","\n")
# Create a dot node for each element in the lattice
for (i in 1:nrow(S)) {
x <- round(S$asset[i], digits=digits)
y <- round(S$option[i], digits=digits)
# Detect the American tree and draw accordingly
early.exercise <- ""
if (("exercise" %in% colnames(S)) && S$exercise[i]) {
early.exercise <- "shape=oval,"
}
cat("node", i, "[", early.exercise, "label=\"", x, ", ", y, "\"];", "\n", sep="")
}
# The number of levels in a binomial lattice of length N
# is `$\frac{\sqrt{8N+1}-1}{2}$`
L <- ((sqrt(8*nrow(S)+1)-1)/2 - 1)
k<-1
for (i in 1:L) {
tabs <- rep("\t",i-1)
j <- i
while(j>0) {
cat("node",k,"->","node",(k+i),";\n",sep="")
cat("node",k,"->","node",(k+i+1),";\n",sep="")
k <- k + 1
j <- j - 1
}
}
cat("}", sep="")
}To finally plot the tree, first one must pass on the variables to genlattice function. Then data resulting from genlattice will be passed on to dotlattice function with a precision definition. Output is visualized through Graphiviz.
x <- genlattice.vanilla.european.call(Asset=100, Volatility=0.2, IntRate=0.05, DividendRate=0.02, Strike=90, Expiry=0.5, NoSteps=10)
y <- capture.output(dotlattice(x, digits=4))
cat(y, file="/Volumes/Vol1/Cloud_Juun/git/clone_repos/ARIMA1/efd/Binomial Trees/lattice.dot")