#
# A popular numerical method to value options is the binomial option pricing model
# pioneered by Cox, Ross, and Rubinstein (1979). This method relies on
# representing the possible future evolution of an asset by a (recombining) tree.
#

 

#
# In this model, as in the Black-Scholes-Merton (1973) setup, there is a risky asset,
# an index or stock, and a riskless asset, a bond. The relevant time interval from
# today until the maturity of the option is divided, in general, into equidistant
# subintervals of length dt (Delta t). 
#

​

#
# The code that follows presents a Python implementation that creates a
# recombining tree based on some fixed numerical parameters for the model:
#

​

import math
import numpy as np
import numba
import time

​

S0 = 36. # Initial value of the risky asset.
T = 1.0 # Time horizon for the binomial tree simulation.
r = 0.06 # Constant short rate.

# Eu chamaria de constant risk-free interest rate, sendo em % per annum
sigma = 0.2 # Constant volatility factor. # Também em % per annum !!!

​

def simulate_tree(M):
    dt = T / M # Length of the time intervals.
    u = math.exp(sigma * math.sqrt(dt)) # Factors for the upward and downward movements.
    d = 1 / u # Factors for the upward and downward movements.
    S = np.zeros((M + 1, M + 1))
    S[0, 0] = S0
    z = 1
    for t in range(1, M + 1):
        for i in range(z):
            S[i, t] = S[i, t-1] * u
            S[i+1, t] = S[i, t-1] * d
        z = 1 + z
    return S
 

#
#
# Contrary to what happens in typical tree plots, an upward movement is
# represented in the ndarray object as a sideways movement, which decreases the
# ndarray size considerably:
#
#

​

np.set_printoptions(formatter={'float':lambda x: '%6.2f' % x})

​

simulate_tree(4)

​

#

# Out[99]: array([[ 36.00, 39.79, 43.97, 48.59, 53.71],
#                           [  0.00, 32.57, 36.00, 39.79, 43.97],
#                          [  0.00,   0.00, 29.47, 32.57, 36.00],
#                         [  0.00,   0.00,   0.00, 26.67, 29.47],
#                        [  0.00,   0.00,   0.00,   0.00, 24.13]])

#

​

start_time = time.time()
simulate_tree(500)
print("--- %s seconds ---" % (time.time() - start_time))

# CPU times: user 148 ms, sys: 4.49 ms, total: 152 ms
# Wall time: 154 ms

​

#
# With some trickery, such a binomial tree can be created with NumPy based on
# fully vectorized code:
#
# In the NumPy case, the code is a bit more compact.
#
# However, more importantly, NumPy vectorization achieves a speedup
# of an order of magnitude (Wall time: 154 ms and Wall time: 12.9 ms)

# while not using more memory:
#

​

def simulate_tree_np(M):
    dt = T / M
    up = np.arange(M + 1)
    up = np.resize(up, (M + 1, M + 1))
    down = up.transpose() * 2
    S = S0 * np.exp(sigma * math.sqrt(dt) * (up - down))
    # Tree for four time intervals (upper-right triangle of values)
    # Ou seja, só importa a parte triangular superior direita da Matriz S
    return S

​

start_time = time.time()
simulate_tree_np(500)
print("--- %s seconds ---" % (time.time() - start_time))

​

# %time simulate_tree_np(500)
# Wall time: 3.95 ms

#

 

#

# Numba (dynamic compilation)
# This financial algorithm should be well suited to optimization through Numba
# dynamic compilation. And indeed, another speedup compared to the NumPy
# version of an order of magnitude is observed. This makes the Numba version
# orders of magnitude faster than the Python (or rather hybrid) version:
#

​

simulate_tree_nb = numba.jit(simulate_tree)

​

start_time = time.time()
simulate_tree_nb(500)
print("--- %s seconds ---" % (time.time() - start_time))

​

# %time simulate_tree_nb(500)
# Wall time: 992 µs

​