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.
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.
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
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
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!
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.
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.