Binomial European Option Pricing in R

The code for this is available at linanqiu/binomial-european-option-r.

This post by Intel https://software.intel.com/en-us/articles/binomial-options-pricing-model-code-for-intel-xeon-phi-coprocessor is a surprisingly good explanation of the theoretical background of this subject and an insanely fast CPU level implementation of this.

My Financial Engineering class was working on Binomial European Option Pricing, and the prof insisted that we show the entire tree matrix for the intermediate steps (and not just the final price). Turns out the packages available on https://cran.r-project.org/ don’t seem to do that. Either way, good exercise to implement this from scratch and revise for the midterms.

First, let’s talk about our model.

Let’s say time is \(T\) where \(T=1\) represents a year, and \(T=0.25\) represents 3 months. The number of periods we simulate is \(N\). Then, the amount of time represented by each period is \(\Delta t = \frac{T}{N}\).

We assume that in each period, the stock can go up by

\[U = e^{\sigma * \sqrt{\Delta t}}\]

where \(\sigma\) is the volatility, and \(\Delta t\) is \(\frac{T}{N}\).

Then, if the stock value is \(S\) at period 1, it would go up to

\[S_U = US = e^{\sigma * \sqrt{\Delta t}}S\] in period 2.

It can also go down by

\[D = e^{-\sigma * \sqrt{\Delta t}}\].

Astute readers will recognize this as a Geometric Brownian Motion (I will probably make another post about this next time).

That means, say \(S=10\) and \(\sigma = 0.2\), \(T=1\) (1 year) and \(N=2\), then \(U = e^{\sigma * sqrt{\Delta t}} = 1.15191\) and \(D = 0.8681234\).

1
2
3
period 0: 10.0
period 1: 8.68 11.5
period 2: 7.54 10.0 13.3

The number directly below represents \(S_D = DS\) and the number directly above represents \(S_U = US\), and the second period is calculated recursively.

Then, we can build a stock tree!

1
2
3
4
5
6
7
8
9
10
11
12
13
build_stock_tree = function(S, sigma, delta_t, N) {
tree = matrix(0, nrow=N+1, ncol=N+1)
u = exp(sigma*sqrt(delta_t))
d = exp(-sigma*sqrt(delta_t))
for (i in 1:(N+1)) {
for (j in 1:i) {
tree[i,j] = S * u^(j-1) * d^((i-1)-(j-1))
}
}
return(tree)
}

To see what this gives us, let’s try it with the parameters above:

1
2
3
4
5
> build_stock_tree(S=10, sigma=0.2, delta_t=1/2, N=2)
[,1] [,2] [,3]
[1,] 10.000000 0.0000 0.00000
[2,] 8.681234 11.5191 0.00000
[3,] 7.536383 10.0000 13.26896

Cool! We just did what we wanted to do but way faster. In fact, we can expand the number of levels really easily. Note that when N changes, delta_t should change too. This is not too much of a problem – we’ll be calculating delta_t programmatically soon.

1
2
3
4
5
6
7
8
> build_stock_tree(S=10, sigma=0.2, delta_t=1/5, N=5)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 10.000000 0.000000 0.000000 0.00000 0.00000 0.00000
[2,] 9.144406 10.935647 0.000000 0.00000 0.00000 0.00000
[3,] 8.362017 10.000000 11.958837 0.00000 0.00000 0.00000
[4,] 7.646568 9.144406 10.935647 13.07776 0.00000 0.00000
[5,] 6.992333 8.362017 10.000000 11.95884 14.30138 0.00000
[6,] 6.394073 7.646568 9.144406 10.93565 13.07776 15.63948

Now how do we get an option’s price from the underlier’s (stock’s) price? Let’s say that we have a call option. Then, we know that at termination, the value of the call option is \(\max(S-X, 0)\) where \(S\) here is the price of the stock at termination, and \(X\) is the strike price. So say \(X=10\), then in our example above, the various values of the call option would be

1
0.0000000 0.0000000 0.0000000 0.9356469 3.0777623 5.6394832

But how do we advance up the option back to the price at the beginning? We call upon the powers of the q-measure where \(q\) is the risk-neutral no-arbitrage probability which is independent of the actual probability of the stock going up or down. Finance textbooks tell us that

\[q \equiv \frac{R-D}{U-D}\]

where \(R\) is \(e^{r*\Delta t}\) where \(r\) is the interest rate per annum continuously compounded. \(R\) is then the discount factor. The derivation of \[q\] will be left to another time.

Then, we can write a function to calculate \[q\].

1
2
3
4
5
6
q_prob = function(r, delta_t, sigma) {
u = exp(sigma*sqrt(delta_t))
d = exp(-sigma*sqrt(delta_t))
return((exp(r*delta_t) - d)/(u-d))
}

Then, we know that

\[C = \frac{qC_U + (1-q)C_D}{R}\]

where \(C\) is the price of the option in a previous time period, and \(C_U\) is the value of the option when the underlier went up, and \(C_D\) is the value of the option when the underlier went down.

That means we can derive the step above the last row of the call option. Let’s say that \(r = 0.1\) for 10 percent interest continuously compounded a year. First we find that

\[R = e^{r * \Delta t} = e^{\frac{0.1}{5}} = 1.020201\]

Then, \(q = \frac{R-D}{U-D} = 0.5904327\) using our q_prob function above with q_prob(0.1, 1/5, 0.2).

1
2
0.0000000 0.0000000 0.5414976 2.1568506 4.499392 0.000000
0.0000000 0.0000000 0.0000000 0.9356469 3.077762 5.639483

You should verify the first row using the q-measure method.

We can code this into a method (account for both “put” and “call”).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
value_binomial_option = function(tree, sigma, delta_t, r, X, type) {
q = q_prob(r, delta_t, sigma)
option_tree = matrix(0, nrow=nrow(tree), ncol=ncol(tree))
if(type == 'put') {
option_tree[nrow(option_tree),] = pmax(X - tree[nrow(tree),], 0)
} else {
option_tree[nrow(option_tree),] = pmax(tree[nrow(tree),] - X, 0)
}
for (i in (nrow(tree)-1):1) {
for(j in 1:i) {
option_tree[i, j] = ((1-q)*option_tree[i+1,j] + q*option_tree[i+1,j+1])/exp(r*delta_t)
}
}
return(option_tree)
}

This function takes in a stock tree generated in the build_stock_tree function, and returns an option tree. For the previous option, it would generate

1
2
3
4
5
6
7
8
9
> tree = build_stock_tree(S=10, sigma=0.2, delta_t=1/5, N=5)
> value_binomial_option(tree, sigma=0.2, delta_t=1/5, r=0.1, X=10, type='call')
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.3515415 0.0000000 0.0000000 0.0000000 0.000000 0.000000
[2,] 0.6365307 1.8937675 0.0000000 0.0000000 0.000000 0.000000
[3,] 0.1813700 0.9740419 2.5965507 0.0000000 0.000000 0.000000
[4,] 0.0000000 0.3133870 1.4656468 3.4698679 0.000000 0.000000
[5,] 0.0000000 0.0000000 0.5414976 2.1568506 4.499392 0.000000
[6,] 0.0000000 0.0000000 0.0000000 0.9356469 3.077762 5.639483

Now let’s put everything together and make the output a nice little list.

1
2
3
4
5
6
binomial_option = function(type, sigma, T, r, X, S, N) {
q = q_prob(r=r, delta_t=T/N, sigma=sigma)
tree = build_stock_tree(S=S, sigma=sigma, delta_t=T/N, N=N)
option = value_binomial_option(tree, sigma=sigma, delta_t=T/N, r=r, X=X, type=type)
return(list(q=q, stock=tree, option=option, price=option[1,1], delta=delta))
}

Usually, there’s also something that we want to compute: the \[\Delta\] of stocks – the amount of stocks that are required to replicate the portfolio.

\[\Delta = \frac{C_U - C_D}{S_U - S_D}\]

We code this as

1
2
3
4
5
delta = function(binomial_option, row, col) {
stock_tree = binomial_option$stock
option_tree = binomial_option$option
return((option_tree[row+1, col+1] - option_tree[row+1, col])/(stock_tree[row+1, col+1] - stock_tree[row+1, col]))
}

Then, the overall function becomes

1
2
3
4
5
6
7
binomial_option = function(type, sigma, T, r, X, S, N) {
q = q_prob(r=r, delta_t=T/N, sigma=sigma)
tree = build_stock_tree(S=S, sigma=sigma, delta_t=T/N, N=N)
option = value_binomial_option(tree, sigma=sigma, delta_t=T/N, r=r, X=X, type=type)
delta = (option[2,2]-option[2,1])/(tree[2,2]-tree[2,1])
return(list(q=q, stock=tree, option=option, price=option[1,1], delta=delta))
}

Look how pretty the output is for our example!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
> binomial_option(type='call', sigma=0.2, T=1, r=0.1, X=10, S=10, N=5)
$q
[1] 0.5904327
$stock
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 10.000000 0.000000 0.000000 0.00000 0.00000 0.00000
[2,] 9.144406 10.935647 0.000000 0.00000 0.00000 0.00000
[3,] 8.362017 10.000000 11.958837 0.00000 0.00000 0.00000
[4,] 7.646568 9.144406 10.935647 13.07776 0.00000 0.00000
[5,] 6.992333 8.362017 10.000000 11.95884 14.30138 0.00000
[6,] 6.394073 7.646568 9.144406 10.93565 13.07776 15.63948
$option
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.3515415 0.0000000 0.0000000 0.0000000 0.000000 0.000000
[2,] 0.6365307 1.8937675 0.0000000 0.0000000 0.000000 0.000000
[3,] 0.1813700 0.9740419 2.5965507 0.0000000 0.000000 0.000000
[4,] 0.0000000 0.3133870 1.4656468 3.4698679 0.000000 0.000000
[5,] 0.0000000 0.0000000 0.5414976 2.1568506 4.499392 0.000000
[6,] 0.0000000 0.0000000 0.0000000 0.9356469 3.077762 5.639483
$price
[1] 1.351541
$delta
[1] 0.7018805

Let’s try the model for a put with the same parameters.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
> binomial_option(type='put', sigma=0.2, T=1, r=0.1, X=10, S=10, N=5)
$q
[1] 0.5904327
$stock
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 10.000000 0.000000 0.000000 0.00000 0.00000 0.00000
[2,] 9.144406 10.935647 0.000000 0.00000 0.00000 0.00000
[3,] 8.362017 10.000000 11.958837 0.00000 0.00000 0.00000
[4,] 7.646568 9.144406 10.935647 13.07776 0.00000 0.00000
[5,] 6.992333 8.362017 10.000000 11.95884 14.30138 0.00000
[6,] 6.394073 7.646568 9.144406 10.93565 13.07776 15.63948
$option
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.3999157 0.0000000 0.00000000 0 0 0
[2,] 0.7232877 0.1892841 0.00000000 0 0 0
[3,] 1.2369985 0.3916873 0.05535867 0 0 0
[4,] 1.9613263 0.7768750 0.13789428 0 0 0
[5,] 2.8096541 1.4399698 0.34348430 0 0 0
[6,] 3.6059268 2.3534319 0.85559356 0 0 0
$price
[1] 0.3999157
$delta
[1] -0.2981195

As a sanity check, we ensure that

  1. Price of call minus price of put (long a call, short a put) equals a forward. In other words \(C - P = S + \frac{X}{e^{r}}\). In our case, \(1.351541 - 0.3999157 = 10 - \frac{10}{e^{0.1}} = 0.9516258\) so that’s good!
  2. \(\Delta_C - \Delta_P = 1\), which again is true!

Woohoo!

What if we want to find out how option price evolves as we increase the number of periods? Well we can do that!

1
2
3
4
5
6
7
8
9
10
11
periods = seq(100, 120)
option_price_vary_period = function(period) {
print(period)
option = binomial_option(type='call', sigma=0.2, T=1, r=0, X=100, S=100, N=period)
return(option$price)
}
values = sapply(periods, option_price_vary_period)
library(ggplot2)
data = as.data.frame(list(periods=periods, values=values))
plot = ggplot(data=data) + geom_line(aes(x=periods, y=values)) + labs(title="Call Value", x="Periods", y="Value")
plot

However, this gets really slow as we increase the number of periods. This is rather unfortunate, since we are going at \(O(N^2)\). However, we can use a slightly faster parallel implementation using the library parallel.

1
2
3
4
5
6
7
8
9
10
library(parallel)
cl = makeCluster(8)
clusterEvalQ(cl, source('binomial.R'))
periods = seq(100, 1000)
periods = sample(periods)
valuesPar = parSapply(cl=cl, periods, option_price_vary_period)
data = as.data.frame(list(periods=periods, values=valuesPar))
plot = ggplot(data=data) + geom_line(aes(x=periods, y=values, alpha=0.1)) + geom_point((aes(x=periods, y=values))) + labs(title="Call Value", x="Periods", y="Value")
plot
stopCluster(cl)

This assumes that binomial.R is in the same folder. This should speed things up A LOT. Reason why I randomized periods in the 5th line is because the larger periods take WAY longer, so you’ll want to distribute that among the cores rather evenly (since parSapply segments the input into equal segments increasingly). This doesn’t affect the graphing at all, and if you want you can always sort the result.