The Binomial Model

You have discussed the Binomial Tree Model in class and can even compute simple examples with pen and paper. This lecture will generalize the binomial method to arbitrary lengths of time period, etc. I have provided a simple example in an Excel Spreadsheet to get the gist of what we are doing. This works well for relatively simple models, but for more complicated examples, switching to a more sophisticated language is useful. However, in the beginning it is useful to have an example we know the answer to for computation. Let's start with our import statements.

In [1]:
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline

A Simple Example

The example that I compute in the Excel Spreadsheet is a 4 step, multiplicative binomial tree for a European Call Option. However, since there are no dividends, the same model could apply to an American Call as these will never be exercised early.

Exogenous Parameters:

  • The length of time to expiration: $T = 1$
  • The initial stock price: $S = \$100$
  • The strike price: $K = \$100$
  • The risk free rate: $r = 0.06$
  • The number of time steps: $N = 3$
  • The up proportion: $u = 1.1$
  • The down proportion: $d = 0.9091$

Endogenous Parameters:

  • The time step: $\Delta t = \frac{T}{N}$
  • The risk neutral probability: $p = \frac{e^{r \Delta t} - d}{u-d}$
  • The discount factor: $e^{-r\Delta t}$

The first step is to enter the computational parameters and compute the endogenous ones before proceeding.

In [4]:
#Exogenous Parameters - note how I am careful to include decimal places for real numbers.
T, S, K, r, N, u, d = 1.0, 100.0, 100.0, .06, 3, 1.1, .9091

dt   = T / N
disc = np.exp(-r * dt)
p    = (np.exp(r * dt) - d) / (u - d)

print 'The exogenous parameters: dt, disc, p'
print dt, disc, p
The exogenous parameters: dt, disc, p
0.333333333333 0.980198673307 0.581987113812

The computational parameters match those in our Excel example, which is a good first step. Now solving for the $t = 0$ price involves solving for the asset price and call price at all nodes on the tree, subject to the boundary condition. Lets first calculate the stock value at the $4$ grid points at the end of time.

In [10]:
St      = np.zeros(N + 1)
St[0]   = S * d ** N

for j in range(1, N+1):
    St[j] = St[j-1] * u / d

# Check that the end values for the stock price are correct
print St
[  75.13373406   90.9109091   110.0011      133.1       ]

Next we want to calculate the value of the option at the end of time.

In [20]:
C  = np.zeros(N + 1)

for j in range(N + 1):
    C[j] = max(0,St[j] - K)

#Check that the end values for the call price value are correct
print C
[  0.       0.      10.0011  33.1   ]

Now we want to step back through the tree to find the value at time $0$. This is an inherently recursive procedure. Therefore, we don't store all of the grid points, just those from the first period in the future.

In [34]:
C  = np.zeros(N + 1)

for j in range(N + 1):
    C[j] = max(0,St[j] - K)

for i in range(N):
    k = N - i - 1
    for j in range(k+1):
        C[j] = disc * (p * C[j+1] + (1- p) * C[j])
print 'The price of the option is:', C[0]
The price of the option is: 10.1454621546

Now what if we wanted to change some of the key parameters? This can be made easier by setting up our Call Option Pricing Example as a function. This can be done as follows:

In [2]:
def CallPrice(T = 1.0, S = 100.0, K = 100.0, r = 0.06, N = 3, u = 1.1, d =.9091):
    dt      = T / N
    disc    = np.exp(-r * dt)
    p       = (np.exp(r * dt) - d) / (u - d)
    St      = np.zeros(N + 1)
    C       = np.zeros(N + 1)
    St[0]   = S * d ** N
    
    for j in range(1, N+1):
        St[j] = St[j-1] * u / d
    
    for j in range(N + 1):
        C[j] = max(0,St[j] - K)
    
    for i in range(N):
        k = N - i - 1
        for j in range(k+1):
            C[j] = disc * (p * C[j+1] + (1- p) * C[j])
     
    return C[0], dt, p, disc
In [5]:
C = CallPrice()[0]
C
Out[5]:
10.145462154644788
In [54]:
C0, dt, p, disc = CallPrice()

print 'The Call Price is:              ', C0
print 'The risk neutral probability is:', p
The Call Price is:               10.1454621546
The risk neutral probability is: 0.581987113812

Notice how I have initialised the parameters for our example. This means that I don't need to provide arguments to the function, it will just use the defaults. However, if we want to check that the example is stable as the grid size is increased we can set $N = 100$

In [58]:
C0, dt, p, disc = CallPrice(N = 100)
print 'The Call price is:              ', C0
print 'The step size is:               ', dt
print 'The risk neutral probability is:', p
The Call price is:               38.453570645
The step size is:                0.01
The risk neutral probability is: 0.479309481592

Notice that the call price changes a lot as we change the grid size. This method turns out to be unstable. Therefore, we need to be a little more careful about how we alter the step size. It turns out that when we alter the step size we also must alter the up and down jumps and the risk neutral probability accordingly. There are four methods proposed in Clemlow and Strickland. Each of these methods takes as given the risk neutral mean of the process ($\nu = r - \frac{1}{2}\sigma^2$) and the risk neutral variance, $\sigma^2$. The problem arises from the fact that these two parameters are needed to tie down $u, d, p$ leading to overidentification. We therefore require additional assumptions to pin down the parameter values. The four methods below use different assumptions.

  • Cox, Ross and Rubenstein

In this method we set $p = \frac{1}{2}$ and then find: $$ u = e^{(r - \frac{1}{2}\sigma^2)\Delta t + \sigma \sqrt{\Delta t}} \qquad d = e^{(r - \frac{1}{2}\sigma^2)\Delta t - \sigma \sqrt{\Delta t}} $$

  • Jarrow and Rudd

In this method we set $u = -d = e^{\sigma \sqrt{\Delta t}}$ and $p = \frac{1}{2} + \frac{r - \frac{1}{2}\sigma^2}{2\sigma}\sqrt{\Delta t}$.

  • Equal Probabilities/Additive

The second two methods employ an additive binomial tree structure. Previously we have assumed a multiplicative one. With an additive structure the natural logarithm of $S$, which we denote $x$ changes to $x + \Delta x_u$ with probability $p_u$. This leads to the following two equations in three unknowns:

$$\mathbb{E}[\Delta x] = p \Delta x_u + (1-p)\Delta x_d = \nu \Delta t $$ $$\mathbb{E}[\Delta x^2] = p \Delta x_u^2 + (1-p)\Delta x_d ^2 = \sigma^2 \Delta t + \nu^2 \Delta t^2 $$

Where $\nu = r - \frac{1}{2} \sigma^2$ as before. We once again can set $p = \frac{1}{2}$ to find: $$\Delta x_u = \frac{1}{2}\nu\Delta t + \frac{1}{2}\sqrt{4\sigma^2\Delta t - 3 \nu^2 \Delta t^2}$$ $$\Delta x_d = \frac{3}{2}\nu\Delta t - \frac{1}{2}\sqrt{4\sigma^2\Delta t - 3 \nu^2 \Delta t^2}$$

  • Trigeorgis

In this method we use equal jump sizes like in the Jarrow and Rudd method. This leads to:

$$\Delta x = \sqrt{\sigma^2 \Delta t + \nu^2 \Delta t^2} \qquad p = \frac{1}{2} + \frac{\nu \Delta t}{2\Delta x} $$

I don't expect you guys to understand why these approximations make sense at this point. However, once you learn about Geometric Brownian Motion and its random walk approximation you may see where the math comes from. Note that the key problem with binomial trees is that we have to make additional assumptions to pin down either the probability or the step size. The finite difference methods we introduce later will solve this problem.

Example: Using the Trigeorgis Method to Price a European Call

  • $K = \$100$
  • $T = 1$
  • $S = \$100$
  • $\sigma = 0.2$
  • $r = 0.06$
  • $N = 3$

Notice that now we don't specify the up and down jumps. They are based on the parameter values. As well, we will start by computing a simple 3 period example and hope to scale it up to a finer grid. The idea is that this time the solution should be more stable for large $N$.

In [8]:
def EuroCallPriceTrig(K = 100.0, T = 1.0, S = 100.0, sigma = 0.2, r = 0.06, N = 3):
    #Precompute constants
    dt  = T/N
    nu  = r - 0.5 * sigma ** 2
    dxu = np.sqrt(dt * sigma ** 2 + (nu * dt) ** 2)
    dxd = - dxu #symmetric jump sizes
    p   = 0.5 + 0.5 * (nu * dt / dxu)
    disc = np.exp(-r * dt)
    
    St = np.zeros(N + 1)
    C  = np.zeros(N + 1)
    
    #Initialize end of time Stock price
    St[0] = S * np.exp(N * dxd)
    for j in range(1,N+1):
        St[j] = St[j-1] * np.exp(dxu - dxd)
    
    #Initialize end of time call price - boundary condition
    for j in range(N + 1):
        C[j] = max(0, St[j] - K)
    
    #Step back through the tree
    for i in range(N):
        k = N - i - 1
        for j in range(k + 1):
            C[j] = disc * (p * C[j+1] + (1- p) * C[j])
     
    return C[0]

Remember that before when we changed the grid size we saw large fluctuations in the price of the option. Lets do a test to confirm that this is a stable method.

In [117]:
pricestrig = np.zeros(100)
nrange     = np.zeros(100)

for n in range(1,101):
    NN = 5 * n   
    nrange[n-1] = NN
    pricestrig[n-1] = EuroCallPriceTrig(N = NN)

We can also check the stability of the equal probabilities method.

In [118]:
def EuroCallPriceEqP(K = 100.0, T = 1.0, S = 100.0, sigma = 0.2, r = 0.06, N = 3):
    #Precompute constants
    dt  = T/N
    nu  = r - 0.5 * sigma ** 2
    dxu = 0.5 * nu * dt + 0.5 * np.sqrt(4 * dt * sigma ** 2 - 3 * (nu*dt) **2) 
    dxd = 1.5 * nu * dt - 0.5 * np.sqrt(4 * dt * sigma ** 2 - 3 * (nu*dt) **2) 
    p   = 0.5 
    disc = np.exp(-r * dt)
    
    St = np.zeros(N + 1)
    C  = np.zeros(N + 1)
    
    #Initialize end of time Stock price
    St[0] = S * np.exp(N * dxd)
    for j in range(1,N+1):
        St[j] = St[j-1] * np.exp(dxu - dxd)
    
    #Initialize end of time call price - boundary condition
    for j in range(N + 1):
        C[j] = max(0, St[j] - K)
    
    #Step back through the tree
    for i in range(N):
        k = N - i - 1
        for j in range(k + 1):
            C[j] = disc * (p * C[j+1] + (1- p) * C[j])
     
    return C[0]

priceseq = np.zeros(100)
nrange   = np.zeros(100)

for n in range(1,101):
    NN = 5 * n   
    nrange[n-1] = NN
    priceseq[n-1] = EuroCallPriceEqP(N = NN)
In [119]:
plt.plot(nrange, pricestrig, 'b-', linewidth = 3, label ='Trigeorgis Method')
plt.plot(nrange, priceseq, 'g-', linewidth = 3, label ='Equal Prob. Method')
plt.xlabel('Number of Grid Points')
plt.ylabel('Price of a European Call Option')
plt.title('Stability of the Binomial Tree Method')
plt.legend()
plt.show()

As you can see, the price of the call settles down around $\$11$ as the grid gets finer. As well, there are still minor fluctuations in the price, showing that the approximation is not exact. As I mentioned before, we will use finite difference methods to come up with a more stable solution. As well, note how handy it is to define a function for pricing this type of option. We can easily change the key parameters with minimal hassle. This is not the case in a program like Excel. Note how there seems to be some discrepancy between the two methods, even for fine grids. This is another reason to use more accurate methods as we should worry that the assumption we make regarding the free parameter is influencing our results. From now on we will use the Trigeorgis Method as it seems to settle down the quickest.

Computing the Greeks

We will now use our EuroCallTrig() function to determine the hedge sensitivities:

  • $$\delta = \frac{\partial C}{\partial S}$$

For a small enough time step, $\Delta t$ we can use the one step ahead finite difference approximation: $$\frac{\partial C}{\partial S} = \frac{C_{1,1} - C_{1,0}}{S_{1,1}- S_{1,0}} $$ This will require saving the one step ahead values for $C$ and $S$. We can easily add this to our pricing function.

  • $$\gamma = \frac{\partial^2 C}{\partial S^2}$$

For the second partial derivative we can use the following formula:

$$ \frac{\partial^2 C}{\partial S^2} = \frac{(C_{2,2} - C_{2,1})/(S_{2,2} - S_{2,1}) - (C_{2,1} - C_{2,0}) / (S_{2,1} - S_{2,0})}{\frac{1}{2}(S_{2,2} - S_{2,0})}$$ This requires two step ahead information from our binomial tree, however, for small $\Delta t$ this approximation should work well.

  • $$\nu = \frac{\partial C}{\partial \sigma}$$

For the last two Greeks we can just use small perturbations of the option pricing function to approximate the derivatives: $$ \frac{\partial C}{\partial \sigma} = \frac{C(\sigma + \Delta \sigma) - C(\sigma - \Delta \sigma)}{2 \Delta \sigma}$$

  • $$\rho = \frac{\partial C}{\partial r}$$

Likewise:

$$ \frac{\partial C}{\partial r} = \frac{C(r + \Delta r) - C(r - \Delta r)}{2 \Delta r}$$

These are just central finite difference approximations to the derivative.

Let's amend our pricing function to return $\delta$ and $\gamma$ as outputs as well.

In [30]:
def EuroCallPriceTrig(K = 100.0, T = 1.0, S = 100.0, sigma = 0.2, r = 0.06, N = 3):
    #Precompute constants
    dt  = T/N
    nu  = r - 0.5 * sigma ** 2
    dxu = np.sqrt(dt * sigma ** 2 + (nu * dt) ** 2)
    dxd = - dxu #symmetric jump sizes
    p   = 0.5 + 0.5 * (nu * dt / dxu)
    disc = np.exp(-r * dt)
    
    St = np.zeros(N + 1)
    C  = np.zeros(N + 1)
    Cd = np.zeros(5)
    Sd = np.zeros(5)
    s  = 0
    
    Sd[1] = np.exp(dxu) * S  
    Sd[0] = np.exp(dxd) * S
    Sd[4] = np.exp(2 * dxu) * S
    Sd[3] = S
    Sd[2] = np.exp(2 * dxd) * S
    #Initialize end of time Stock price
    St[0] = S * np.exp(N * dxd)
    for j in range(1,N+1):
        St[j] = St[j-1] * np.exp(dxu - dxd)
    
    #Initialize end of time call price - boundary condition
    for j in range(N + 1):
        C[j] = max(0, St[j] - K)
    
    #Step back through the tree
    for i in range(N):
        k = N - i - 1
        for j in range(k + 1):
            C[j] = disc * (p * C[j+1] + (1- p) * C[j])
        if k == 2:
            gamma1 =  (C[2] - C[1]) / (Sd[4] - Sd[3])
            gamma2 =  (C[1] - C[0])/(Sd[3] - Sd[2]) 
            gamma3 = .5 * (Sd[4] - Sd[2])
            gamma  = (gamma1 - gamma2) / gamma3
        if k == 1:
             delta = (C[1] - C[0]) / (Sd[1] - Sd[0])
        
    return C[0], delta, gamma
In [46]:
C, delta, gamma = EuroCallPriceTrig(N = 200)
print 'The Solution: Call Price, delta, gamma'
print C, delta, gamma

sig  = .2
dsig = .0001 * sig
sigmaplus = sig + dsig
sigmamin  = sig - dsig

r    = .06
dr   = r * .0001
rplus = r + dr
rmin  = r - dr

vega = (EuroCallPriceTrig(N=200, sigma = sigmaplus)[0] - EuroCallPriceTrig(N=200, sigma = sigmamin)[0]) / (dsig * 2.0)
rho  = (EuroCallPriceTrig(N=200, r = rplus)[0] - EuroCallPriceTrig(N=200, r = rmin)[0]) / (dr * 2.0)
print 'The Estimates of Vega and Rho are:'
print vega, rho
The Solution: Call Price, delta, gamma
10.9800066195 0.655230638912 0.0184934321947
The Estimates of Vega and Rho are:
36.7683026989 54.5754665916

An American Put Option: Adjusting for Early Exercise

I am not going to write the function for the put option as this will be left as an exercise. However, I will discuss the key changes that occur.

  • The time $T$ payoffs are different:

$$P_{N,j} = \max\{0, K - S_{N,j} \} $$

  • We have to account for early exercise. Thus at time node $i < N we have that:

$$ P_{n,j} = \max \{K - S_{n,j}, P_{n,j}\}$$

These two simple modifications should allow you to complete the exercise listed below.

Dealing with Dividends

It is very easy to incorporate a constant flow rate of dividends, $\delta$, into the binomial model. We simply have to adjust the drift of the process, $\nu = r - \delta - \frac{1}{2} \sigma^2$ to account for the dividend payout. Once we make this simple correction the rest of the results carry through. We won't discuss this simple modification as it is relatively straightforward. But here is a European Call Price function which incorporates this feature:

In [48]:
def EuroCallPrice(K = 100.0, T = 1.0, S = 100.0, sigma = 0.2, r = 0.06, N = 3, delta = .02):
    #Precompute constants
    dt  = T/N
    nu  = r - delta - 0.5 * sigma ** 2
    dxu = np.sqrt(dt * sigma ** 2 + (nu * dt) ** 2)
    dxd = - dxu #symmetric jump sizes
    p   = 0.5 + 0.5 * (nu * dt / dxu)
    disc = np.exp(-r * dt)
    
    St = np.zeros(N + 1)
    C  = np.zeros(N + 1)
    Cd = np.zeros(5)
    Sd = np.zeros(5)
    s  = 0
    
    Sd[1] = np.exp(dxu) * S  
    Sd[0] = np.exp(dxd) * S
    Sd[4] = np.exp(2 * dxu) * S
    Sd[3] = S
    Sd[2] = np.exp(2 * dxd) * S
    #Initialize end of time Stock price
    St[0] = S * np.exp(N * dxd)
    for j in range(1,N+1):
        St[j] = St[j-1] * np.exp(dxu - dxd)
    
    #Initialize end of time call price - boundary condition
    for j in range(N + 1):
        C[j] = max(0, St[j] - K)
    
    #Step back through the tree
    for i in range(N):
        k = N - i - 1
        for j in range(k + 1):
            C[j] = disc * (p * C[j+1] + (1- p) * C[j])
        if k == 2:
            gamma1 =  (C[2] - C[1]) / (Sd[4] - Sd[3])
            gamma2 =  (C[1] - C[0])/(Sd[3] - Sd[2]) 
            gamma3 = .5 * (Sd[4] - Sd[2])
            gamma  = (gamma1 - gamma2) / gamma3
        if k == 1:
             delta = (C[1] - C[0]) / (Sd[1] - Sd[0])
        
    return C[0], delta, gamma

C0 = EuroCallPrice(delta = 0.0, N = 200)[0]
C1 = EuroCallPrice(delta = .02, N = 200)[0]

print 'The Euro Call Price with No dividends: ', C0
print 'The Euro Call Price with delta = .02:  ', C1
The Euro Call Price with No dividends:  10.9800066195
The Euro Call Price with delta = .02:   9.71884531718

Adding dividends clearly reduces the price of the call. This is because a constant flow rate of dividends slows the growth of the stock price, lowering the likelihood of a high stock price at the exercise date. Thus, it is less likely the call will be exercised. Even if it is the payoff will be smaller on average, leading to a lower call price.

If the asset pays a known discrete proportional dividend at a known time in the future we can adjust the binomial tree to account for this. We will assume that the dividend date corresponds with a node on the tree. If the time is before the dividend date, $\tau$, then the grid remains unchanged. However, the value of the stock at all days after $\tau$ is $(1-\hat{\delta})S_{t,i}$, where $\hat{\delta}$ is the proportional dividend. Once we adjust the end of time stock price downwards by the dividend factor we can compute the tree as above. However, we had to make a lot of assumptions regarding the dividend (proportional, known date, known size). To accomodate more realistic dividend structures is complicated (and not discussed in Clemlow and Strickland. We will not further discuss dividends.

Time Varying Volatility

A common requirement in practice is the incorporation of time varying volatility into a derivative pricing model. In order to price securities consistently we must make adjustments to the binomial tree model to account for time varying volatilities. However, for time varying volatilities it is probably best to use a more advanced method, such as finite differences.

The simplest way to adjust the binomial model for time varying volatility is to fix the space step and adjust the probabilities and the length of the time step. This ensures that the tree recombines. If we fix the space step as $\Delta x$ and have a time varying volatility, $\sigma_i$ for time step $i$, this leads to time varying risk neutral probabilities, time step and drifts: $$p_i \Delta x - (1-p_i)\Delta x = \nu_i \Delta t_i $$ $$p_i \Delta x^2 + (1-p_i)\Delta x^2 = \sigma_i^2\Delta t_i + \nu_i^2\Delta t_i^2 $$ $$\nu_i = r - \frac{1}{2}\sigma_i^2$$

Solving the above system for $\nu_i, \Delta t_i$ and $p_i$ for each time and for an appropriately chosen $\Delta x$ leads to a time varying volatility version of the binomial method. A good choice of $\Delta x$ is: $$\Delta x = \sqrt{\bar{\sigma}^2 \Delta \bar{t} + \bar{\nu}^2 \Delta \bar{t}^2} \qquad\mbox{Where:} \;\; \bar{\sigma} = \frac{1}{N}\sum_{i=1}^{N}\sigma_i\qquad \bar{\nu} = \frac{1}{N}\sum_{i=1}^{N}\nu_i $$

Exercise: Pricing an American Put

  1. Write a function which uses the binomial method to price an American put option with the following parameters:
  • $K = 100$, $T = 1$, $N = 100$, $S = 100$, $r = .06$, $\sigma = 0.2$

You can use either the Equal Probability or the Trigeorgis method. Report the price for $N = \{100, 200, 300\}$. Is the method convergent? Use your function to calculate vega.