Trinomial Trees and Finite Difference Methods

We showed how to solve option pricing problems using binomial tree methods in an earlier note. However, we showed that the choice of either equal probabilities or equal step sizes impacted the price of the option. This is not a desirable property. As well, it was challenging to deal with time varying parameters since the tree was not always recombining. To overcome these difficulties we will introduce a trinomial tree method in this note. This imporoves the accuracy of the approximation.

Next we go on to solve the Black-Scholes Partial Differential Equation with finite difference methods. It turns out this is analgous to using a trinomial tree. However, we may also view this as solving the PDE:

$$-\frac{\partial C(x,t)}{\partial t} +rC(x,t)= \nu\frac{\partial C(x,t)}{\partial x} + \frac{1}{2}\sigma^2\frac{\partial^2 C(x,t)}{\partial x^2} $$

Subject to appropriate boundary conditions.

To understand how finite difference methods work recall the definition of a derivative:

$$\frac{d f(x)}{dx} = \lim_{\Delta x \to 0} \frac{f(x+\Delta x) - f(x)}{\Delta x} $$

If we take the domain of the function to be $[0,1]$ and divide into $N$ discrete grid points we can think of

Example

Let's solve the ODE:

$$\frac{df(x)}{dx} = e^x \qquad \mbox{s.t.} \;\; f(0) = 1 $$

We know the solution to the ODE is:

$$f(x) = e^x $$

Thus, we can see how well a finite difference scheme works.

We can write out the approximation as:

$$\frac{f_{i+1} - f_i}{\Delta x} = e^{x_i} $$

Or:

$$f_{i+1} = f_i + \Delta x e^{x_i} $$

Subject to $f_0 = 1$. We can then compare the approximation to the known solution. We would hope for arbitrarily small $\Delta x$ the solution would converge to the true one.

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

xmin = 0.0
xmax = 1.0
N    = 3

xgrid = np.linspace(xmin, xmax, N)
dx    = 1.0 * (xmax - xmin) / (N - 1.0)
f    = np.ones(N)

for i in range(1,N):
    f[i] = f[i-1] + dx * np.exp(xgrid[i-1])
    
ftrue = np.exp(xgrid)

plt.plot(xgrid, ftrue, 'r--', linewidth = 4, label = 'Actual Solution')
plt.plot(xgrid, f, 'b-', linewidth = 4, label = 'Approximate Solution')
plt.legend(loc = 'upper left')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Finite Difference Solution: N = 3')
plt.show()

xmin = 0.0
xmax = 1.0
N    = 30

xgrid = np.linspace(xmin, xmax, N)
dx    = 1.0 * (xmax - xmin) / (N - 1.0)
f    = np.ones(N)

for i in range(1,N):
    f[i] = f[i-1] + dx * np.exp(xgrid[i-1])
    
ftrue = np.exp(xgrid)

plt.plot(xgrid, ftrue, 'r--', linewidth = 4, label = 'Actual Solution')
plt.plot(xgrid, f, 'b-', linewidth = 4, label = 'Approximate Solution')
plt.legend(loc = 'upper left')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Finite Difference Solution: N = 30')
plt.show()

Applying finite differences to partial differential equations generalizes this procedure. However, first we analyze a trinomial tree.

In a trinomial tree we allow the stock price to rise or fall by an equal amount, $\Delta x$. However, it may also stay the same. We need to solve for the risk neutral probabilities but not the step size!

For the differential equation of $x = \ln S$ we have:

$$dx = \nu dt + \sigma dz $$

For a given time increment and $\Delta x$:

$$\mathbb{E}[\Delta x ] = p_u \Delta x + p_m \times 0 - p_d \Delta x = \nu \Delta t $$ $$\mathbb{E}[\Delta x^2 ] = p_u \Delta x^2 + p_m \times 0 + p_d \Delta x^2 = \sigma^2 \Delta t + \nu^2 \Delta t^2$$ $$p_u + p_m + p_d =1$$

It turns out that a good choice of $\Delta x$ is $\sigma \sqrt{3\Delta t}$. The reason why we cannot choose it independently of $\Delta t$ has to do with the discrete time approximation of a Brownian Motion. Solving the three equations above for the three risk neutral probabilities allows us to construct the grid and apply the finite difference method. Solving the trinomial tree is very similar to solving the binomial model and doing so will be left as an exercise.

Solving the PDE With Finite Differences

We can view the finite difference method as simply extending the Trinomial Tree to a grid of values. If we consider the solution to the PDE to be $C(x, t)$ then we can think of the discretized solution as a matrix of values. Each column represents the call price as a function of the underlying asset price at the time step in question. The important thing we have to account for is boundary conditions. As will be shown in tutorial, we need to make assumptions about the value of the call option for high and low values of $x$. This is because a forward difference approximation requires information from outside the grid.

To understand the problem a bit further we will use the following notation. $i$ will index the $x$ grid and $t$ the time grid. Therefore, $C_{i,t}$ is really just $C(x_i, t)$. The Black SCholes PDE becomes:

$$\frac{C_{i,t+1} - C_{i,t}}{\Delta t} = 0.5\sigma^2 \frac{C_{i+1,t+1} - 2 C_{i, t+1} + C_{i-1,t+1}}{\Delta x^2} + \nu \frac{C_{i+1,t+1} - C_{i-1,t+1}}{2\Delta x} - r C_{i,t+1} $$

It can be shown that this is equivalent to:

$$C_{i,t} = p_u C_{i+1,t+1} + p_m C_{i,t+1} + p_d C_{i-1,t+1} $$

For appropriate choices of the probabilities. As well, we use the following boundary conditions:

$$\frac{\partial C}{\partial S} =1 \;\; \mbox{For S large} \qquad \frac{\partial C}{\partial S} = 0 \;\; \mbox{For S small}$$

This means that $C_{N,t} = C_{N-1,t} + S_{N,t} - S_{N-1,t}$ and $C_{-N,t} = C_{-N+1, t}$. We may want to change the size of the grid to make sure that the boundary condition is not impacting the call price. We will solve the American Call Pricing problem in the assignment.

Stability and Convergence

It is important to ensure that the probabilities, $p_u,p_d,p_m$ are positive and that the solution converges for different grid sizes. Imagine that $C$ represents the exact, unknown, solution of the PDE and $O$ represents the exact solution of the finite difference equation. The difference, $C-O$ is called the discretization error and is a result of approximating partial derivatives by finite differences. This can be made arbitrarily small by making the grid finer in both time and space. $\tilde{O}$ represents the actual numerical solution. Therefore, $\tilde{O}-O$ is the error due to rounding. A finite difference method is convergent if the discretisation error tends to zero as the time and space steps get smaller and the round off error is small and bounded at all times.

It turns out that if the volatility is small relative to the interest rate, $p_d$ can become negative. Furthermore, if the volatility is too large, then $p_m$ can become negative. Therefore, we cannot choose our time and space step sizes independently of $\sigma$. This leads to the recommendation (which we will not dwell on) that:

$$\Delta x = \sigma \sqrt{3 \Delta t} $$

Assignment Problem

Choose one of the following two questions to do.

  1. For a given time increment and $\Delta x$ solve for $p_u,p_d,p_m$ in the trinomial model. Use these values to write a function which prices a European call option with K = 100, S = 100, r = .06 , $\delta$ = .03, $\sigma$ = .2. Assume that $\Delta x = \sigma \sqrt{3\Delta t}$ This should closely resemble your code for the binomial tree example.

  2. For a given time increment and $\Delta x$ solve for $p_u,p_d,p_m$ in the finite difference model. Use these values to write a function which prices an American Call Option with K = 100, S = 100, r = .06 , $\delta$ = .03, $\sigma$ = .2. Assume that $\Delta x = \sigma \sqrt{3\Delta t}$. Plot the call price as a function of $x$ for a few different time steps. You may choose to start with a coarse grid and make it finer. The output should be an $I\times N$ matrix where $I$ is the number of $x_i$ gridpoints and $N$ is the number of time steps.