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 | period 0: 10.0 |

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 | build_stock_tree = function(S, sigma, delta_t, N) { |

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

1 | > build_stock_tree(S=10, sigma=0.2, delta_t=1/2, N=2) |

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 | > build_stock_tree(S=10, sigma=0.2, delta_t=1/5, N=5) |

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 | q_prob = function(r, delta_t, sigma) { |

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 | 0.0000000 0.0000000 0.5414976 2.1568506 4.499392 0.000000 |

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 | value_binomial_option = function(tree, sigma, delta_t, r, X, type) { |

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 | > tree = build_stock_tree(S=10, sigma=0.2, delta_t=1/5, N=5) |

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

1 | binomial_option = function(type, sigma, T, r, X, S, N) { |

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 | delta = function(binomial_option, row, col) { |

Then, the overall function becomes

1 | binomial_option = function(type, sigma, T, r, X, S, N) { |

Look how pretty the output is for our example!

1 | > binomial_option(type='call', sigma=0.2, T=1, r=0.1, X=10, S=10, N=5) |

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

1 | > binomial_option(type='put', sigma=0.2, T=1, r=0.1, X=10, S=10, N=5) |

As a sanity check, we ensure that

- 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!
- \(\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 | periods = seq(100, 120) |

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 | library(parallel) |

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.