In this tutorial, I wanted to demonstrate how to set up the binomial tree valuation in Python and explain the mathematical calculation happening in the background. In this example, I will use a simple European option contract.
To code the binomial tree in Python, we can refer to the following three steps:
- Define the valuation parameters
- Initialize the matrix with the prices of the underlying asset
- Derive the matrix of the option based on the underlying asset matrix
Let’s start with the first step!
Step 1: Define the Valuation Parameters
The binomial tree calculation is based on the following parameters:
- S0 – the price of the underlying asset at period t=0;
- u & d – the factor of the asset price increasing (u) or decreasing (d) at each period; in most practice exercises, the u and d are defined as 0.1 and -0.1, respectively;
- p – the risk-neutral probability of the asset price increasing at each node; consequently, (1-p) signifies the risk-neutral probability at each node;
- r – the risk-free rate that is used in the calculation of p;
- T – time to expiry
- N – the number of periods that we observe the price movements of the underlying asset; (sum of the factors u and d)
- sigma – the volatility of the underlying asset
- K – strike price
Wrapping the parameters in the formula, we get:
u = e^{\sigma \times \sqrt{\Delta t}}[\katex]
We start the code by defining a function which takes the same parameters as the model we defined in the formulae.
As you can see in the first code snippet, the first step of parameter definition is rather straightforward. We can use the standard Python operators in combination with the numpy library which offers more functionality such as np.sqrt() function for square-root and np.exp() for raising the exponent (e) to a specified factor.
def calltree_numpy(S0, K, T, r, sigma, q, N):
"""European call price based on an N-step binomial tree"""
# setting the parameters
deltaT = T / float(N)
u = np.exp(sigma*np.sqrt(deltaT))
d = 1.0 / u
p = (np.exp((r-q)*deltaT)-d) / (u-d)
# probabilities * discounting factors
piu = np.exp(-(r-q)*deltaT) * p
pid = np.exp(-(r-q)*deltaT) * (1-p)Step 2: Initialize the Matrix S
Let’s continue with looking at the binomial tree from the underlying asset perspective. What we will do is set up the binomial tree scheme in the Python environment.
It can be done in three steps:
- Initialize a zero matrix with the (N+1) * (N+1) dimensions
- Populate the matrix with the values based on the valuation parameters S0, u and d
- Exclude the lower triangle from the matrix
The code looks as follows:
S = np.zeros((N+1, N+1))
for i in range(N+1):
for j in range(i, N+1):
S[i, j] = S0 * u**(j-i) * d**(i)
S = np.triu(S)
Let’s go line-by-line to understand the code better!
We start at initializing the zero-matrix of certain dimensions, i.e., we want to define a matrix S which is populated with zeros and has N+1 rows and N+1 columns.
Now, let’s see what happens at each node of the tree!
First, let’s agree that each node can be written as S[i, j], where S refers to the matrix, i is the index location of the row and j is the index location of the column. For instance, the node Su^2d is located at [1, 3].
Let us look separately at one node to see what is going on there. We can say that the value at each node is calculated as S[i, j] = S0*u^(j-i)*d^(i). Taking our previous example, S[1, 3] = S0*u^(3-1)*d^(1) = Su^2d.
Finally, we use a numpy function np.triu() to select only the upper triangle of the matrix. If you draw the scheme for yourself, you can see that now the matrix resembles the flipped scheme as outlined in Figure 1.
Step 3: Derive the Option Matrix
The option matrix, and, consequently the option price can be obtained in the following four steps:
- Initialize a zero-matrix with the same dimensions as the previously defined matrix S
- Look at the final matrix column where the asset price at option expiry is calculated. At this step, we want to see if the value of (S-K) is greater than zero. This step stems from the logic that at each node, the buyer will only execute the option contract if she can gain a profit from buying S for the price of K.
- Using backward induction, populate the option matrix with the values derived from the matrix S
- Return the value of the option at the node [0,0]
The Python code looks as follows:
#setting matrix
C = np.zeros((N+1, N+1))
C[:, N] = np.maximum(0, S[:, N]-K)
for j in range(N-1, -1, -1):
C[:j+1, j] = piu * C[:j+1, j+1] + pid * C[1:j+2, j+1]
return C[0, 0]
Let’s look at each step separately. After initializing the zero-matrix, we select the excerpt of the matrix – the last column. In Python, you can specify this excerpt as C[:, N].
At each node in the last column, we want to select the maximum between 0 and S-K. If S is higher than the strike price K, the profit from the option contract is realized and the investor will choose to execute the contract. In the opposite scenario, the investor abstains from realizing the contract and chooses to receive nothing (i.e., zero) rather than incur a loss.
If you print your option matrix at this point, you will see that the last column is populated whereas the other matrix nodes are zero. Hence, in the next step, we want to use backward induction starting from the last column to populate column N-1, and consequently, columns N-2, N-3, etc.
We are now looking only at columns, therefore, we take the previously defined j as an indication of a column. As we are moving backwards in the matrix, we start at the second last column N-1 and move until the first node at index 0 (specified as until -1 non-inclusive). And each time we move one step back (specified as -1). You can see that in Python code, the range is defined as range(N-1, -1, -1).
Next, we will be taking a snippet of matrix C and populating it based on the other matrix snippets. At this point, we can refer to the principles of matrix manipulation. Let’s visualize this step!
Suppose we are looking at j = 2, i.e., our current matrix snippet is C[:3, 2] which stands for the first three rows of column 2. In the picture, you can see that it is a 3*1 matrix, we can mark it as matrix A. To obtain the matrix of 3*1 we use the other two matrices, i.e., C[:3, 3] and C[1:4, 3]. You can see these matrices on the picture framed as B and C, respectively.
Now, matrices B and C are multiplied by piu and pid parameters, respectively. Both matrices have 3*1 dimensions. To obtain matrix A, we simply add up B and C. If you write out the calculation, it looks as follows:
In Python, we use the for-loop to perform the same calculation on all columns thus populating all matrix snippets until we arrive at the first node C[0, 0] which is a 1*1 matrix or our option price!
The complete code looks as follows:
def calltree_numpy(S0, K, T, r, sigma, q, N):
"""European call price based on an N-step binomial tree"""
#setting the parameters
deltaT = T / float(N)
u = np.exp(sigma*np.sqrt(deltaT))
d = 1.0 / u
p = (np.exp((r-q)*deltaT)-d) / (u-d)
# probabilities * discounting factors
piu = np.exp(-(r-q)*deltaT) * p
pid = np.exp(-(r-q)*deltaT) * (1-p)
#setting matrix
C = np.zeros((N+1, N+1))
S = np.zeros((N+1, N+1))
for i in range(N+1):
for j in range(i, N+1):
S[i, j] = S0 * u**(j-i) * d**(i)
S = np.triu(S)
C[:, N] = np.maximum(0, S[:, N]-K)
for j in range(N-1, -1, -1):
C[:j+1, j] = piu * C[:j+1, j+1] + pid * C[1:j+2, j+1] #backward induction
return[0,0]