Solution

#Simple formula for computing a binomial coeffecient
def binomial_coefficient(n, k):
    return math.factorial(n) / (math.factorial(n-k) * math.factorial(k))
print(binomial_coefficient(3, 2))

#Probability density function for a binomial distribution
def binomial_distribution_pdf(p, n, k):
    q = 1-p
    return (p ** k) * (q ** (n-k)) * binomial_coefficient(n, k)
print(binomial_distribution_pdf(p_up, 3, 2))

We could also make a function to compute for each level all the probabilities based on this model.

#Build a function to find associated probabilities for each node
def create_node_probabilities(rf, T, sigma, T_increment):
    node_levels = []

    p_up = compute_up_probability(rf, T_increment, sigma)

    for t in range(T+1):
        nodes = [binomial_distribution_pdf(p_up, t, up) for up in range(t+1)]
        nodes = nodes[::-1]
        node_levels.append(nodes)
    return node_levels
probabilities = create_node_probabilities(.04, 12, .08, 1/12)
print(probabilities)

With all of this information in mind, we can mutliply all probabilities of nodes in levels by their values and then sum to get the expected value of each level. Using the numpy function dot will multiply numbers in the same position of two lists and then sum all these values. If you look at the expected value, it is clear that it is going up as a linear function related to the risk-free growth rate.

node_levels = create_nodes(100, 12, .08, 1/12)
probabilities = create_node_probabilities(.04, 12, .08, 1/12)

#Compute expected value of the stock in each state
expected_value = []
for prob, vals in zip(probabilities, node_levels):
    #Converting to numpy arrays and doing dot applies the dot product to vectors
    expected_value.append(np.array(prob).dot(np.array(vals)))
expected_value = pd.Series(expected_value)
expected_value.plot(kind='line')
plt.show()

We can confirm this by discounting back all expected values by the discount factor applicable given the risk-free rate and the total time.

The code below will find the value at each node if a $100 strike price call were to be exercised at that time.

#Notice discounting each shows risk-neural probabilities mean no gain or loss
discount_factors = [np.exp(.04 * 1/12*t) for t in range(13)]
discount_factors = pd.Series(discount_factors)

(expected_value / discount_factors).plot(kind='line')
plt.ylim([99,101])
plt.show()

With these trees set up, we begin to explore how they are used for options. Given an option we can value each node according to the value it would have if exercised at that point in time. To aid in the computation, let’s use np.where. This functions takes as a first argument the list of true/false values for a list. The second argument is what the value should be for each element of the list if the value is True in the position, and the third argument is what the value should be if that boolean value is instead False.

probabilities = create_node_probabilities(.04, 12, .08, 1/12)
node_levels = create_nodes(100, 12, .08, 1/12)

#Iterate through what each node would be valued at as a call option
for prob, vals in zip(probabilities, node_levels):
    prob = np.array(prob)
    vals = np.array(vals)
    #np.where takes a condition as the first argument and then a value or list of values for the second/third
    #which represents the values if true for each element in the second or false in the third
    vals = np.where(vals-100 < 0, 0, vals-100)
    print(vals)

The way to solve the binomial tree to find the price of an option is something called backwards induction. Essentially, first we need to find the value of all nodes in the final time period, then we go back one level at a time finding the value for each node in the prior time period. To begin with, I am switching to four quarters to make it easier to see and including a function which finds the value only for the last level, and gives a question mark for all other levels.

def backward_induction(node_levels, strike_price):
    opt_values = []
    #The expiration level has values based on the derivative formula
    level = node_levels[-1]
    level = np.array(level)
    level = np.where(level-strike_price < 0, 0, level-strike_price).round(2)
    opt_values.append(level)
    #Go through each level in reverse order after the expiration level
    #For now, each gets ? for the value to represent not knowing
    for level in node_levels[::-1][1:]:
        opt_values.append(["?" for x in level])
    #Reverse the order of the list
    opt_values = opt_values[::-1]
    return opt_values

#Switch to 4 quarters
node_levels = create_nodes(100, 4, .08, 1/4)
print(node_levels)

print()
opt_values = backward_induction(node_levels, 100)
print(opt_values)

We have a few steps to deal with now. First, we can bring in the data.

#Set up the data
node_levels = create_nodes(100, 4, .08, 1/4)
opt_values = backward_induction(node_levels, 100)
node_levels = convert_nodes_to_labels(node_levels)

Next, we need to deal with the way in which our function for labels will fail. If we try to create labels with some of the values as question marks, there becomes issues because we are using something that formats only floating point numbers. To deal with it, let's make a temporary function to use.

#Create a function to convert the option labels while we have the questions marks
def convert_nodes_to_labels_TEMP(node_levels):
    node_labels = []
    for t in range(len(node_levels)):
        #Get rid of the string format for the number
        nodes = ["T={}-${}".format(t, x) for x in node_levels[t]]
        node_labels.append(nodes)
    return node_labels
opt_labels = convert_nodes_to_lab

And the final step... build the position dictionary and draw the binary tree.

pos_dictionary = create_node_positions(node_levels)
draw_binary_tree(node_levels, pos_dictionary,opt_labels)

To find the value of the next level of nodes we do two steps. First, find the expected value which is simply the probability of an up movement multiplied by the option value in an up movement plus the probability of a down movement multiplied by the option value in a down movement. Then with the expected value, we discount by the e raised to the risk-free rate multiplied by the time between node levels (so in this case 1/4 * r).