Derivative Valuation

American Options: A comparative analysis of the Cox Ross Rubinstein (CRR) and Longstaff & Schwartz (LS) valuation methods.

Paul McAteer, MSc, MBA

pcm353@stern.nyu.edu

1.0 Introduction ¶

This technical note compares the prices obtained from the two most common methods employed for the valuation of American-Style options. We implement the binomial "tree-based" technique of Cox Ross Rubinstein¹ (CRR) and the least-squares Monte Carlo method of Longstaff & Schwartz ² (LS). The design, implementation and output of both are validated by comparing our results to those produced in the academic literature.

¹Cox, J., S. Ross, and M. Rubinstein (1979): “Option Pricing: A Simplified Approach.” Journal of Financial Economics, Volume 7, Issue 3, pp. 229–263.
²Longstaff, Francis & Schwartz, Eduardo. (2001). Valuing American Options by Simulation: A Simple Least-Squares Approach. Review of Financial Studies. 14. 113-47.

2.1 LS Method: Model Theory and Design ¶

The valuation of an American-style option using the Longstaff Shwarz approach requires the Monte Carlo simulation of $N$ realizations of the asset path from now to expiration using time steps of $δt$ . If the random walk is lognormal then we will simulate

S_{j + 1} = S_{j} \exp ((r - \frac{1}{2} σ^{2}) δ t + σ \sqrt{δ} t ϕ)

$S_{j+1}= S_{j}\exp\ (( r-\frac{1}{2}\ \sigma^2)\ \ \delta t+\sigma\sqrt\delta t\ \phi)$

Where $\phi$ is a Normally distributed random variable, $r$ is the risk free rate, $\sigma$ is volatility and and $S_{j}$ is the stock price after $j$ time steps.

For the purposes of exposition, we will follow the example of the Longman Shwartz paper which demonstrates the algorithm over three time steps and 8 simulated price paths, considering a 3-year American put option on a non-dividend-paying stock that can be exercised at the end of year 1, the end of year 2, and the end of year 3. The risk-free rate is 6% per annum (continuously compounded). The current stock price is 1.00 and the strike price is 1.10. In the implementation phase, the timesteps and simulations will obviously be significantly greater. Having obtained the 8 simulations of the asset path, we calculate the intrinsic value for each simulation at each time step. The present value of the average of all of the payoffs at expiration is the value of a European option. See below for the Stock Price Matrix and the Intrinsic Value Matrix.

From the stock price matrix, there are five paths where the option is in the money at time $t+2$ . These are paths 1, 3, 4, 6, and 7. Obviously at the final exercise date, $t+3$ , the optimal exercise strategy for an American option is to exercise the option if it is in the money (ITM). Prior to the final date however, at $t+2$ , the optimal strategy is to compare the immediate exercise value $X-S$ with the expected cash flows from continuing $E[V|S]$ , and then exercise if immediate exercise is more valuable, i.e. where $X-S > E[V|S]$ . Conversely, one would continue to hold where the continuation value is more valuable, i.e. where $E[V|S] > X-S$

The conditional expectation function is obtained at each time step by regressing the present value of subsequent realized cash flows from continuation at t+2 on the simulated stock price at the time step t+2 when the continuation/exercise decision is made. That is we assume an approximate relationship of the form:

E [V | S] = \hat{V} = a + b S + c S^{2}

$E[V|S] =\ \hat{V}= \ a\ +bS\ +\ cS^{2}$

Where $S$ is the stock price at the 2-year point, $\hat{V}$ is the expected value of continuing discounted back to the two year time point, based on the expected subsequent prices generated in the Monte Carlo simulation process. For this example, we visualize below the present value of subsequent realized cash flows and the conditional expectation of the discounted value of subsequent realized cash flows.

Having obtained a vector of continuation values for t+2, we now compare these with a vector of immediate exercise values at t+2 and determine optimal exercise behaviour:

We use this to construct a matrix composed of the cashflows which would result from optimal exercise behaviour, conditional on not exercising prior to t+2. Note that if there is a positive entry in the time t+2 column it means that we should have exercised there (if not earlier) and so the cashflows for later times are set to zero.

Proceeding recursively, we next examine whether the option should be exercised at time t+1. From the stock price matrix, there are again five paths where the option is in the money at time 1. These are paths 1, 4, 6, 7, and 8.

We now construct a cashflow matrix assuming that the options are optimally exercised at each time step. The average exercised value at all time steps is the value of the American option.

The American option valuation problem is an optimal stopping problem, formulated as follows:

V (t, s) = sup_{τ \in [t, T]} E^{Q} [e^{- r (τ - t)} Φ (S_{τ}) | S_{t} = s]

$V(t,s) = \sup_{\tau \in [t,T]} \mathbb{E}^{\mathbb{Q}}\biggl[ e^{-r(\tau-t)} \Phi(S_{\tau}) \bigg| S_t=s \biggr]$

where

Φ (\cdot)

$\Phi(\cdot)$ is the payoff, also called intrinsic value.This formula resembles the formula of a European option, except that in this case the option can be exercised at any time. The problem centres on determining the best exercise time

τ

$\tau$ that maximizes the expectation. Note that the time

τ

$\tau$ is a random variable.

2.2 LS Method: Replication of Canonical Longman Shwartz Example¶

import numpy as np
import scipy.stats as ss
import pandas as pd


import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import display
import sympy; sympy.init_printing()

def display_matrix(m):
    display(sympy.Matrix(m))

N = 4          # number of time steps
r = 0.06       # interest rate
K = 1.1        # strike 
T = 3          # Maturity

dt = T/(N-1)          # time interval
df = np.exp(-r * dt)  # discount factor per time interval

S = np.array([
            [1.00, 1.09, 1.08, 1.34],
            [1.00, 1.16, 1.26, 1.54],
            [1.00, 1.22, 1.07, 1.03],
            [1.00, 0.93, 0.97, 0.92],
            [1.00, 1.11, 1.56, 1.52],
            [1.00, 0.76, 0.77, 0.90],
            [1.00, 0.92, 0.84, 1.01],
            [1.00, 0.88, 1.22, 1.34]])

display_matrix(S)

In the previous cell we defined the stock matrix S.
It has 8 rows that correspond to the number of paths.
The 4 columns correspond to the 4 time steps, i.e. each row is a path with 4 time steps.

H = np.maximum(K - S, 0)           # intrinsic values for put option
V = np.zeros_like(H)               # value matrix
V[:,-1] = H[:,-1]

# Valuation by LS Method
for t in range(N-2, 0, -1):        #Return an object that produces a sequence of integers from start (inclusive) t+1 to 
                                   # to stop (exclusive) t+0 by each step (backwards)

    good_paths = H[:,t] > 0        # Extract paths where the intrinsic value is positive 
                                   # the regression is performed only on these paths 
    
    rg = np.polyfit( S[good_paths, t], V[good_paths, t+1] * df, 2)    # polynomial regression
    C = np.polyval( rg, S[good_paths,t] )                             # evaluation of regression  
    
    exercise = np.zeros( len(good_paths), dtype=bool)    # initialize matrix for exercise rule
    exercise[good_paths] = H[good_paths,t] > C           # paths where it is optimal to exercise
    
    V[exercise,t] = H[exercise,t]                        # set V equal to H where it is optimal to exercise 
    V[exercise,t+1:] = 0                                 # set future cash flows, for that path, equal to zero  
    discount_path = (V[:,t] == 0)                        # paths where we didn't exercise 
    V[discount_path,t] = V[discount_path,t+1] * df       # set V[t] in continuation region
    
V0 = np.mean(V[:,1]) * df  # discounted expectation of V[t=1]
print("Example price= ", V0)

Example price=  0.11443433004505696

The matrix H = np.maximum(K - S, 0), is the matrix of intrinsic values:

display_matrix(H.round(2))

The matrix V contains the cash flows.

display_matrix(V.round(2))

2.3 LS Implementation: A Python Class for American Options pricing using LSMC¶

import IPython
import numpy as np
import warnings
warnings.filterwarnings("ignore")
from sys import version

import numpy as np

class AmericanOptionsLSMC(object):
    """ We define a new class to price American options pricing using Longstaff-Schwartz (2001):
    
    The AmericanOptionsLSMC class describes a data object with the following attributes:
    
    option_type: str: 'call' or 'put'
    S0 : float : initial stock/index level
    strike : float : strike price
    T : float : time to maturity (in year fractions)
    M : int : grid or granularity for time (in number of total points). Time to Expiry or Maturity
    r : float : constant risk-free short rate
    div :    float : dividend yield
    sigma :  float : volatility factor in diffusion term 
    
    Unitest(doctest): 
    >>> AmericanPUT = AmericanOptionsLSMC('put', 36., 40., 1., 50, 0.06, 0.06, 0.2, 10000  )
    >>> AmericanPUT.price
    4.4731177017712209
    """
    
    
    def __init__(self, option_type, S0, strike, T, M, r, div, sigma, simulations):
        try:
            self.option_type = option_type
            assert isinstance(option_type, str)
            self.S0 = float(S0)
            self.strike = float(strike)
            assert T > 0
            self.T = float(T)
            assert M > 0
            self.M = int(M)
            assert r >= 0
            self.r = float(r)
            assert div >= 0
            self.div = float(div)
            assert sigma > 0
            self.sigma = float(sigma)
            assert simulations > 0
            self.simulations = int(simulations)
        except ValueError:
            print('Error passing Options parameters')


        if option_type != 'call' and option_type != 'put':
            raise ValueError("Error: option type not valid. Enter 'call' or 'put'")
        if S0 < 0 or strike < 0 or T <= 0 or r < 0 or div < 0 or sigma < 0:
            raise ValueError('Error: Negative inputs not allowed')

        self.time_unit = self.T / float(self.M)
        self.discount = np.exp(-self.r * self.time_unit)

    # Stock Price Matrix generated by Geometric Brownian Motion. Input is instance of AMLS with defined attributes + seed
    # for reproducibility
    @property 
    def MCprice_matrix(self, seed = 123):
        """ Returns MC price matrix rows: time columns: price-path simulation """
        np.random.seed(seed) 
        # Initialize Matrix with number of time steps (+ 1 to have t_0 values incorporated) and no. of sims
        MCprice_matrix = np.zeros((self.M + 1, self.simulations), dtype=np.float64)
        # Fill First row all columns with initial stock price
        MCprice_matrix[0,:] = self.S0
        # For each time step from first (inclusive) to last (last +1 exclusive)
        for t in range(1, self.M + 1):
            # Brownian motion with antithetic variance
            brownian = np.random.standard_normal(self.simulations // 2)
            brownian = np.concatenate((brownian, -brownian))
            MCprice_matrix[t, :] = (MCprice_matrix[t - 1, :]
                                  * np.exp((self.r - self.div - self.sigma ** 2 / 2.) * self.time_unit
                                  + self.sigma * brownian * np.sqrt(self.time_unit)))
        return MCprice_matrix
    
    #Intrinsic Value Matrix
    @property
    def MCpayoff(self):
        """Returns the inner-value of American Option"""
        if self.option_type == 'call':
            payoff = np.maximum(self.MCprice_matrix - self.strike,
                           np.zeros((self.M + 1, self.simulations),dtype=np.float64))
        else:
            payoff = np.maximum(self.strike - self.MCprice_matrix,
                            np.zeros((self.M + 1, self.simulations),
                            dtype=np.float64))
        return payoff

    @property
    def value_vector(self):
        value_matrix = np.zeros_like(self.MCpayoff)
        value_matrix[-1, :] = self.MCpayoff[-1, :]
        for t in range(self.M - 1, 0 , -1):
            regression = np.polyfit(self.MCprice_matrix[t, :], value_matrix[t + 1, :] * self.discount, 5)
            continuation_value = np.polyval(regression, self.MCprice_matrix[t, :])
            value_matrix[t, :] = np.where(self.MCpayoff[t, :] > continuation_value,
                                          self.MCpayoff[t, :],
                                          value_matrix[t + 1, :] * self.discount)

        return value_matrix[1,:] * self.discount


    @property
    def price(self): return np.sum(self.value_vector) / float(self.simulations)
    
    @property
    def delta(self):
        diff = self.S0 * 0.01
        myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0 + diff, 
                                       self.strike, self.T, self.M, 
                                       self.r, self.div, self.sigma, self.simulations)
        myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0 - diff, 
                                       self.strike, self.T, self.M, 
                                       self.r, self.div, self.sigma, self.simulations)
        return (myCall_1.price - myCall_2.price) / float(2. * diff)
    
    @property
    def gamma(self):
        diff = self.S0 * 0.01
        myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0 + diff, 
                                       self.strike, self.T, self.M, 
                                       self.r, self.div, self.sigma, self.simulations)
        myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0 - diff, 
                                       self.strike, self.T, self.M, 
                                       self.r, self.div, self.sigma, self.simulations)
        return (myCall_1.delta - myCall_2.delta) / float(2. * diff)
    
    @property
    def vega(self):
        diff = self.sigma * 0.01
        myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0, 
                                       self.strike, self.T, self.M, 
                                       self.r, self.div, self.sigma + diff, 
                                       self.simulations)
        myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0,
                                       self.strike, self.T, self.M, 
                                       self.r, self.div, self.sigma - diff, 
                                       self.simulations)
        return (myCall_1.price - myCall_2.price) / float(2. * diff)    
    
    @property
    def rho(self):        
        diff = self.r * 0.01
        if (self.r - diff) < 0:        
            myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0, 
                                       self.strike, self.T, self.M, 
                                       self.r + diff, self.div, self.sigma, 
                                       self.simulations)
            myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0, 
                                       self.strike, self.T, self.M, 
                                       self.r, self.div, self.sigma, 
                                       self.simulations)
            return (myCall_1.price - myCall_2.price) / float(diff)
        else:
            myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0, 
                                       self.strike, self.T, self.M, 
                                       self.r + diff, self.div, self.sigma, 
                                       self.simulations)
            myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0, 
                                       self.strike, self.T, self.M, 
                                       self.r - diff, self.div, self.sigma, 
                                       self.simulations)
            return (myCall_1.price - myCall_2.price) / float(2. * diff)
    
    @property
    def theta(self): 
        diff = 1 / 252.
        myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0, 
                                       self.strike, self.T + diff, self.M, 
                                       self.r, self.div, self.sigma, 
                                       self.simulations)
        myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0, 
                                       self.strike, self.T - diff, self.M, 
                                       self.r, self.div, self.sigma, 
                                       self.simulations)
        return (myCall_2.price - myCall_1.price) / float(2. * diff)

OPTION_TYPE = 'put'
S_0 = 44
strike = 40
T = 2
M = 50
r = 0.06
div =  0.0
sigma = 0.4
simulations = 10000

AmericanPUT = AmericanOptionsLSMC(OPTION_TYPE, S_0, strike, T, M, r, div, sigma, simulations)

print ('Price: ', AmericanPUT.price)
print ('Delta: ', AmericanPUT.delta)
print ('Gamma: ', AmericanPUT.gamma)
print ('Vega:  ', AmericanPUT.vega)
print ('Rho:   ', AmericanPUT.rho)
print ('Theta: ', AmericanPUT.theta)

Price:  5.454866604491297
Delta:  -0.29914743601896093
Gamma:  -0.0025164728062225276
Vega:   20.273079995041577
Rho:    -21.10592022499747
Theta:  -1.0181598018009304

# FUNCTION TO CREATE COMPLETE SURFACE 

# THREE VARIABLES: 
# 1. CHANGE IN PRICE
# 2. CHANGE IN SIGMA (VOLATILITY)
# 3. CHANGE IN TIME TO EXPIRY


def LS_Surface(px,sig,tte):
    
    
    list_px = []
    list_vol = []
    list_op = []
    list_tte = []

    
    for S_0 in (px): 
        for sigma in (sig):
            for T in (tte):
                list_px.append(S_0)
                arr_px = np.array(list_px)
                
                list_vol.append(sigma)
                arr_vol = np.array(list_vol)
                
                list_tte.append(T)
                arr_tte = np.array(list_tte)
                
                AmericanOPT = AmericanOptionsLSMC(OPTION_TYPE, S_0, strike, T, M, r, div, sigma, simulations)
                list_op.append(AmericanOPT.price)
                arr_op = np.array(list_op)
    
    return arr_px, arr_vol, arr_tte, arr_op

# FUNCTION INPUTS AND OBJECT INSTANTIATION 

OPTION_TYPE = 'put'
S_0 = 36
strike = 40
T = 1 
M = 50
r = 0.06
div =  0.0
sigma = 0.2
simulations = 10000



#px = np.linspace(36, 44, 10)
px = (36., 38., 40., 42., 44.)
#sig = np.linspace(0.2, 0.4, 10)
sig =(0.2, 0.4)
tte = (1.0, 2.0)

result_LS = LS_Surface(px,sig,tte)

# MODEL OUTPUT
cols=['Stock Price', 'Volatility', 'Expiration','Option Price']
df_surf = pd.DataFrame(result_LS)
df_surf = df_surf.T
df_surf.columns = cols
df_surf

2.4 LS Model Validation: Comparing our model results for LSM American Put valuations in Longstaff & Schwartz (2001) Paper¶

The model results are aligned with "Simulated American" values produced by Longman & Schwartz, indicating that the model is correctly constructed and specified.

2.5 LS Model Validation: Comparing our model results for LSM American Call valuations in Zhao (2018) Paper¶

The model results are aligned with "LSM" values produced by Zhao³, indicating that the model is correctly constructed and specified.

³Zhao, J.n (2018): “American Option Valuation Methods” International Journal of Economics and Finance; Vol. 10, No. 5; 2018

# FUNCTION INPUTS AND OBJECT INSTANTIATION 

OPTION_TYPE = 'call'
S_0 = 100
strike = 100
T = 1 
M = 4
r = 0.05
div =  0.1
sigma = 0.2
simulations = 50000


px = (70,80,90,100,110,120,130)
sig =(0.2, 0.2)
tte = (1.0, 1.0)

result_LS_call = LS_Surface(px,sig,tte)

Our results:

# MODEL OUTPUT
cols=['Stock Price', 'Volatility', 'Expiration','Option Price']
df_surf = pd.DataFrame(result_LS_call)
df_surf = df_surf.T
df_surf.columns = cols
df_surf
df_surf = df_surf.iloc[[0,4,8,12,16,22,26]]
df_surf

Zhao's results:

2.6 LS Model Results: Visualizing the Value Surface¶

# FUNCTION TO CREATE OPTION VALUE SURFACE 

# TWO VARIABLES: 
# 1. CHANGE IN PRICE
# 2. CHANGE IN SIGMA (VOLATILITY)

def LS_Surface_val(px,sig):
    
    
    list_px = []
    list_vol = []
    list_op = []
    
    
    for S_0 in (px): 
        for sigma in (sig):
                list_px.append(S_0)
                arr_px = np.array(list_px)
                
                list_vol.append(sigma)
                arr_vol = np.array(list_vol)
                
                
                AmericanOPT = AmericanOptionsLSMC(OPTION_TYPE, S_0, strike, T, M, r, div, sigma, simulations)
                list_op.append(AmericanOPT.price)
                arr_op = np.array(list_op)
    
    return arr_px, arr_vol, arr_op

OPTION_TYPE = 'put'
S_0 = 36
strike = 40
T = 1 
M = 50
r = 0.06
div =  0.06
sigma = 0.2
simulations = 10000

px = np.linspace(36, 44, 10)
sig = np.linspace(0.2, 0.4, 10)

result_LS_val = LS_Surface_val(px,sig)

pd.set_option('display.max_rows', None)
cols_val=['Stock Price', 'Volatility','Option Price']
df_surf_val = pd.DataFrame(result_LS_val)
df_surf_val = df_surf_val.T
df_surf_val.columns = cols_val
df_surf_val

# Preparation of output for 3D surface chart, requiring 3 2-D matrices of equal dimensions (here 10x10)

arr_surf_px = np.array(df_surf_val['Stock Price'])
arr_surf_px = np.reshape(arr_surf_px, (10, 10))
arr_surf_vol = np.array(df_surf_val['Volatility'])
arr_surf_vol = np.reshape(arr_surf_vol, (10, 10))
arr_surf_opt = np.array(df_surf_val['Option Price'])
arr_surf_opt = np.reshape(arr_surf_opt, (10, 10))

from pylab import plt
plt.style.use('seaborn')
import matplotlib as mpl
mpl.rcParams['font.family'] = 'serif'
mpl.rcParams['axes.labelpad'] = 15
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline



from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(14, 10))
ax = fig.gca(projection='3d')

surf = ax.plot_surface(arr_surf_px, arr_surf_vol, arr_surf_opt, rstride=2, cstride=2,
                       cmap=plt.cm.coolwarm, linewidth=0.5,
                       antialiased=True)

ax.set_xlabel('Price')
ax.set_ylabel('Volatility')
ax.set_zlabel('OPtion Price')
ax.set(facecolor='white')
ax.dist = 12

fig.colorbar(surf, shrink=0.5, aspect=5)

<matplotlib.colorbar.Colorbar at 0x1cf421b1550>

3.1 CRR Method: Model Theory and Design ¶

The probability of reaching a particular node in the binomial tree depends on the numbers of distinct paths to that node and the probabilities of the up and down moves.
The binomial grid of 4 time steps below shows all asset price paths to each node. "u" represents an up-move and "d" a down-move.

The probability of reaching a particular node in the binomial tree depends on the numbers of distinct paths to that node and the probabilities of the up and down moves. The tree below displays the probability of reaching each node. "p" denotes the probability of an up-move and 1-p the probability of a down-move.

Derivation of Parameters¶

Recall that the binomial tree model is formed by building forward the “tree” of interim and terminal asset prices then performing backward induction, which values the option at each node as the probability-adjusted, discounted value of its child nodes. The option price is the sum of these discounted expected values.

The underlying asset price path is determined by the number of time steps and the magnitude of “up-moves” and “down-moves” at each node.

The magnitude of these moves is an exponential function of volatility scaled by square root of time:

UP-MOVE:

u = e^{σ \sqrt{Δ t}}

$u=e^{\sigma \sqrt{\Delta t}}$

DOWN-MOVE:

d = \frac{1}{u} = e^{- σ \sqrt{Δ t}}

$d=\frac{1}{u}= e^{-\sigma \sqrt{\Delta t}}$

Risk neutral valuation allows us to infer the probability of an “up-move”. Risk neutral probabilities are the implied probabilities, which equate the current market price of the option to the discounted expected payoff. They are the outputs of a function which takes the market value, future expected value and the discount rate as the inputs. The fundamental assumption is that the call option seller can form a riskless portfolio by going long a fraction of the share - determined by delta - and writing a call option ( $f$ ):

S_{0} Δ - f

$S_0\Delta -f$

The payoff at each up and down node - the sum of (delta*stock) PLUS the intrinsic value of option - will be the same.

S_{0} u Δ - f_{u} = S_{0} d Δ - f_{d}

$S_0u\Delta -f_u=\ S_0d\Delta -f_d$

Delta is the fraction which equates the payoff at up and down nodes.

Δ = \frac{f_{u} - f_{d}}{S_{0} u - S_{0} d}

$\Delta \ =\ \frac{f_u-f_d}{S_0u-\ S_0d}$

The seller of the option has delta hedged his exposure. The portfolio is riskless and, for there to be no arbitrage opportunities, it must earn the risk-free interest rate. The payoff at each node - the future portfolio value, composed of option and fraction of share - is discounted to the present value at the risk-free rate:

S_{0} Δ - f = (S_{0} u Δ - f_{u}) e^{- r Δ t}

$S_0\Delta -f=\ \left(S_0u\Delta -f_u\right)e^{-r\Delta t}$

Rearranging, we find the option value at each node to be the discounted value of the fractional share position and the discounted intrinsic value of the option:

f = S_{0} Δ (1 - u e^{- r Δ t}) + f_{u} e^{- r Δ t}

$f=\ \left.S_0\Delta \ \left(1-ue^{-r\Delta t}\right)+f_u\right.e^{-r\Delta t}$

Or, equivalently:

f = S_{0} \frac{f_{u} - f_{d}}{S_{0} u - S_{0} d} (1 - u e^{- r Δ t}) + f_{u} e^{- r Δ t}

$f=\ \left.S_0\frac{f_u-f_d}{S_0u-\ S_0d}\ \left(1-ue^{-r\Delta t}\right)+f_u\right.e^{-r\Delta t}$

Which, upon removing the $S_0$ terms, simplifies to:

f = \frac{f_{u} (1 - d e^{- r Δ t}) + f_{d} (u e^{- r Δ t} - 1)}{u - d}

$f=\ \frac{f_u\left.\left(1-de^{-r\Delta t}\right)+f_d\right.\left(ue^{-r\Delta t}-1\right)}{u-d}$

Which can be expressed in terms of a risk neutral probability:

f = e^{- r Δ t} [f_{u} p + f_{d} (1 - p)]

$f\mathrm{=\ }e^{-r\Delta t}\left[f_up+\ f_d(1-p)\ \right]\mathrm{\ }$

Where:

RISK-NEUTRAL PROBABILITY OF UP-MOVE, p:

p = \frac{e^{- r Δ t} - d}{u - d}

$p=\ \frac{e^{-r\Delta t}-d}{u-d}$

3.2. CRR Method: Implementation ¶

import numpy as np


def Binomial(n, S, K, r,q, v, t, PutCall):  

    # Solved Parameters

    At = t/n                                #Timestep

    u = np.exp(v*np.sqrt(At))               # Upmove

    d = 1./u                                # Downmove

    p = (np.exp((r-q)*At)-d) / (u-d)         # Risk-Neutral Probability



    #Binomial price tree

    stockvalue = np.zeros((n+1,n+1))                #Create symmetrical zero-filled array with dimensions of Timesteps + 1

    stockvalue[0,0] = S                             #Location of current stock price in matrix (1st element in 1st column)
                                                    # Note the rows express time progression and columns tree level
    
    for i in range(1,n+1):                          #Outer loop: For each timestep...

        stockvalue[i,0] = stockvalue[i-1,0]*u       # Compute future stockprice as product of prev. price and upmoves

        for j in range(1,i+1):                      # Inner Loop: For each time step at each tree level...
            
            stockvalue[i,j] = stockvalue[i-1,j-1]*d # Compute associated downmove stock value for each upmove
                                                    # (Nested for loop."Inner loop" will be executed one time 
                                                    # for each iteration of the "outer loop").
                                                     
             

    #Option value at terminal nodes   

    optionvalue = np.zeros((n+1,n+1))               #Create symmetrical zero-filled matrix with dimensions of Timesteps + 1

    for j in range(n+1):                            # For each node at the terminal time step...
                                                    # Compute instrinsic value based on put or call position
        if PutCall=="C": # Call

            optionvalue[n,j] = max(0, stockvalue[n,j]-K)      

        elif PutCall=="P": #Put

            optionvalue[n,j] = max(0, K-stockvalue[n,j])

    

    #Backward induction process for option values at interim nodes    

    for i in range(n-1,-1,-1):                    # For each time step, working backwards, from penultimate to first...

        for j in range(i+1):                      # For each node at each of above time step...
                                                  # Compute the maximum of the expected value and the current intrinsic val.

                if PutCall=="P":                 

                    optionvalue[i,j] = max(0, 
                                           K-stockvalue[i,j],np.exp(-r*At)*(p*optionvalue[i+1,j]+(1-p)*optionvalue[i+1,j+1]))

                elif PutCall=="C":

                    optionvalue[i,j] = max(0, 
                                           stockvalue[i,j]-K,np.exp(-r*At)*(p*optionvalue[i+1,j]+(1-p)*optionvalue[i+1,j+1]))
    


    #return optionvalue[0,0]
    return optionvalue[0,0], stockvalue, optionvalue, n, S, K, r, q, v,t

from tabulate import tabulate

# Inputs

n = 500       #number of steps

S = 80       #initial underlying asset price

r = 0.08     #risk-free interest rate

q = 0.12     #risk-free interest rate

K = 100       #strike price

v = 0.2      #volatility

t = 0.5        #time to expiry

PutCall="C"  #option type



result = Binomial(n, S, K, r, q, v, t, PutCall)


header3 = ['Opt. Val.(Basic Model)','Spot', 'Strike', 'Rfr','Div. Yld.', 'Vol.',
           'Time to Exp.', 'Time Steps']
table3 = [[result[0], result[4],result[5],result[6],result[7],result[8],result[9],result[3]]]

print(tabulate(table3,header3))

  Opt. Val.(Basic Model)    Spot    Strike    Rfr    Div. Yld.    Vol.    Time to Exp.    Time Steps
------------------------  ------  --------  -----  -----------  ------  --------------  ------------
                0.214788      80       100   0.08         0.12     0.2             0.5           500

3.3 CRR Model Validation: Comparing our model results to Binomial American Call and Put valuations in Zhao (2018) Paper¶

The model results are aligned with "Binomial" values produced by Zhao, indicating that the model is correctly constructed and specified.

# FUNCTION TO CREATE COMPLETE SURFACE 

# THREE VARIABLES: 
# 1. CHANGE IN PRICE
# 2. CHANGE IN SIGMA (VOLATILITY)
# 3. CHANGE IN TIME TO EXPIRY


def CRR_Surface(px,sig,tte):
    
    
    list_px = []
    list_vol = []
    list_op = []
    list_tte = []

    
    for S in (px): 
        for v in (sig):
            for t in (tte):
                list_px.append(S)
                arr_px = np.array(list_px)
                
                list_vol.append(v)
                arr_vol = np.array(list_vol)
                
                list_tte.append(t)
                arr_tte = np.array(list_tte)
                
                AmericanOPT_bin = Binomial(n, S, K, r, q, v, t, PutCall)
                bin_prices = AmericanOPT_bin[0]
                
                list_op.append(bin_prices)
                arr_op = np.array(list_op)
    
    return arr_px, arr_vol, arr_tte, arr_op

Our results:

# FUNCTION INPUTS AND OBJECT INSTANTIATION 


n = 500       #number of steps

S = 80       #initial underlying asset price

r = 0.08     #risk-free interest rate

q = 0.12     #risk-free interest rate

K = 100       #strike price

v = 0.2      #volatility

t = 0.5        #time to expiry

PutCall="C"  #option type


px = (80, 90, 100, 110, 120)
sig =(0.2, 0.2)
tte = (0.5, 0.5)

result_CRR = CRR_Surface(px,sig,tte)

# MODEL OUTPUT
cols=['Stock Price', 'Volatility', 'Expiration','Option Price']
df_surf_CRR = pd.DataFrame(result_CRR)
df_surf_CRR = df_surf_CRR.T
df_surf_CRR.columns = cols
df_surf_CRR = df_surf_CRR.iloc[[0,4,8,12,16]]
df_surf_CRR

# FUNCTION INPUTS AND OBJECT INSTANTIATION 


n = 500       #number of steps

S = 80       #initial underlying asset price

r = 0.08     #risk-free interest rate

q = 0.12     #risk-free interest rate

K = 100       #strike price

v = 0.2      #volatility

t = 0.5        #time to expiry

PutCall="P"  #option type


px = (80, 90, 100, 110, 120)
sig =(0.2, 0.2)
tte = (0.5, 0.5)

result_CRR = CRR_Surface(px,sig,tte)

# MODEL OUTPUT
cols=['Stock Price', 'Volatility', 'Expiration','Option Price']
df_surf_CRR = pd.DataFrame(result_CRR)
df_surf_CRR = df_surf_CRR.T
df_surf_CRR.columns = cols
df_surf_CRR = df_surf_CRR.iloc[[0,4,8,12,16]]
df_surf_CRR

Zhao's results:

4.1. Model Review: Comparing valuations produced by the validated CRR and LS models¶

# Inputs - CRR

n = 1000       #number of steps
S = 36       #initial underlying asset price
r = 0.06     #risk-free interest rate
q = 0.03     #constant dividend yield
K = 40       #strike price
v = 0.2      #volatility
t = 1        #time to expiry
PutCall="P"  #option type

px = (36., 38., 40., 42., 44.)
sig =(0.2, 0.4)
tte = (1.0, 2.0)

result_CRR = CRR_Surface(px,sig,tte)

# MODEL OUTPUT - CRR
cols=['Stock Price', 'Volatility', 'Expiration','Option Price_CRR']
df_surf_CRR = pd.DataFrame(result_CRR)
df_surf_CRR = df_surf_CRR.T
df_surf_CRR.columns = cols
df_surf_CRR

# Inputs - LS

OPTION_TYPE = 'put'
S_0 = 36
strike = 40
T = 1 
M = 50
r = 0.06
div =  0.03
sigma = 0.2
simulations = 10000




px = (36., 38., 40., 42., 44.)
sig =(0.2, 0.4)
tte = (1.0, 2.0)

result_LS = LS_Surface(px,sig,tte)

# MODEL OUTPUT - LS

cols=['Stock Price', 'Volatility', 'Expiration','Option Price_LS']
df_surf = pd.DataFrame(result_LS)
df_surf = df_surf.T
df_surf.columns = cols
df_surf
df_surf['Option Price_CRR'] =  df_surf_CRR['Option Price_CRR']
df_surf['Difference ($)'] = df_surf['Option Price_LS'] - df_surf['Option Price_CRR']
df_surf['Difference (%)'] = df_surf['Option Price_LS'] / df_surf['Option Price_CRR'] -1
df_surf

mean_abs_diff = df_surf['Difference (%)'].abs().mean()
mean_abs_diff = round(mean_abs_diff,3)
print('The mean absolute difference in option values is', mean_abs_diff)

The mean absolute difference in option values is 0.029

df_scatter = df_surf['Option Price_LS'] 
df_scatter = pd.DataFrame(df_scatter)
df_scatter['Option Price_CRR'] = df_surf['Option Price_CRR']

def render_scatter_plot(data, x_stock_name, y_stock_name, xlim=None, ylim=None):
    fig = plt.figure(figsize=(12,8))
    ax = fig.add_subplot(111)
    ax.scatter(data[x_stock_name], data[y_stock_name])
    if xlim is not None: ax.set_xlim(xlim)
    ax.autoscale(False)
    ax.vlines(0, -10, 10)
    ax.hlines(0, -10, 10)
    ax.plot((-10, 10), (-10, 10))
    ax.set_xlabel(x_stock_name)
    ax.set_ylabel(y_stock_name)

render_scatter_plot(df_scatter,'Option Price_LS', 'Option Price_CRR')

5.0 Conclusion¶

The mean absolute difference in the valuations produced by the two modelling techniques is 0.29% As we can see from the scatter plot, in which the blue diagonal line represents equality, the difference is relatively constant and not significantly impacted by time, volatility or moneyness. The CRR model offers some advantage in terms of computational expense given that model accuracy is not significantly impacted if the number of time steps is set at approximately 500. The L&S method requires some 10,000 simulations to achieve some degree of convergence and stability.

	Stock Price	Volatility	Option Price
0	36.000000	0.200000	5.206054
1	36.000000	0.222222	5.471741
2	36.000000	0.244444	5.752070
3	36.000000	0.266667	6.056894
4	36.000000	0.288889	6.322004
5	36.000000	0.311111	6.599891
6	36.000000	0.333333	6.896866
7	36.000000	0.355556	7.204511
8	36.000000	0.377778	7.492173
9	36.000000	0.400000	7.782835
10	36.888889	0.200000	4.627043
11	36.888889	0.222222	4.920564
12	36.888889	0.244444	5.202559
13	36.888889	0.266667	5.506364
14	36.888889	0.288889	5.806709
15	36.888889	0.311111	6.113080
16	36.888889	0.333333	6.412738
17	36.888889	0.355556	6.710785
18	36.888889	0.377778	7.008733
19	36.888889	0.400000	7.308940
20	37.777778	0.200000	4.087386
21	37.777778	0.222222	4.396794
22	37.777778	0.244444	4.717167
23	37.777778	0.266667	5.021173
24	37.777778	0.288889	5.342402
25	37.777778	0.311111	5.655916
26	37.777778	0.333333	5.965341
27	37.777778	0.355556	6.279945
28	37.777778	0.377778	6.580187
29	37.777778	0.400000	6.885706
30	38.666667	0.200000	3.621215
31	38.666667	0.222222	3.945228
32	38.666667	0.244444	4.269696
33	38.666667	0.266667	4.581571
34	38.666667	0.288889	4.911046
35	38.666667	0.311111	5.230309
36	38.666667	0.333333	5.551210
37	38.666667	0.355556	5.846468
38	38.666667	0.377778	6.163644
39	38.666667	0.400000	6.486886
40	39.555556	0.200000	3.179831
41	39.555556	0.222222	3.507200
42	39.555556	0.244444	3.835114
43	39.555556	0.266667	4.159987
44	39.555556	0.288889	4.481262
45	39.555556	0.311111	4.810664
46	39.555556	0.333333	5.131173
47	39.555556	0.355556	5.454154
48	39.555556	0.377778	5.787542
49	39.555556	0.400000	6.109407
50	40.444444	0.200000	2.772369
51	40.444444	0.222222	3.101532
52	40.444444	0.244444	3.427957
53	40.444444	0.266667	3.761857
54	40.444444	0.288889	4.073630
55	40.444444	0.311111	4.403973
56	40.444444	0.333333	4.731690
57	40.444444	0.355556	5.052291
58	40.444444	0.377778	5.397449
59	40.444444	0.400000	5.719805
60	41.333333	0.200000	2.416657
61	41.333333	0.222222	2.750773
62	41.333333	0.244444	3.070639
63	41.333333	0.266667	3.396295
64	41.333333	0.288889	3.717889
65	41.333333	0.311111	4.040805
66	41.333333	0.333333	4.360816
67	41.333333	0.355556	4.696933
68	41.333333	0.377778	5.028388
69	41.333333	0.400000	5.358597
70	42.222222	0.200000	2.086653
71	42.222222	0.222222	2.411249
72	42.222222	0.244444	2.738304
73	42.222222	0.266667	3.065173
74	42.222222	0.288889	3.393410
75	42.222222	0.311111	3.716457
76	42.222222	0.333333	4.035015
77	42.222222	0.355556	4.360892
78	42.222222	0.377778	4.694451
79	42.222222	0.400000	5.021857
80	43.111111	0.200000	1.784809
81	43.111111	0.222222	2.107221
82	43.111111	0.244444	2.431824
83	43.111111	0.266667	2.756291
84	43.111111	0.288889	3.084635
85	43.111111	0.311111	3.405475
86	43.111111	0.333333	3.728438
87	43.111111	0.355556	4.053395
88	43.111111	0.377778	4.368599
89	43.111111	0.400000	4.695098
90	44.000000	0.200000	1.520289
91	44.000000	0.222222	1.830534
92	44.000000	0.244444	2.141651
93	44.000000	0.266667	2.470506
94	44.000000	0.288889	2.803661
95	44.000000	0.311111	3.120689
96	44.000000	0.333333	3.437835
97	44.000000	0.355556	3.769759
98	44.000000	0.377778	4.082568
99	44.000000	0.400000	4.409799

	Stock Price	Volatility	Expiration	Option Price_LS	Option Price_CRR	Difference ($)	Difference (%)
0	36.0	0.2	1.0	4.810223	4.804027	0.006196	0.001290
1	36.0	0.2	2.0	5.343238	5.375088	-0.031850	-0.005925
2	36.0	0.4	1.0	7.344409	7.478321	-0.133912	-0.017907
3	36.0	0.4	2.0	8.938887	9.112739	-0.173853	-0.019078
4	38.0	0.2	1.0	3.558451	3.603704	-0.045253	-0.012557
5	38.0	0.2	2.0	4.234528	4.312842	-0.078314	-0.018158
6	38.0	0.4	1.0	6.378546	6.518742	-0.140196	-0.021507
7	38.0	0.4	2.0	8.081265	8.276763	-0.195499	-0.023620
8	40.0	0.2	1.0	2.582319	2.647844	-0.065525	-0.024747
9	40.0	0.2	2.0	3.322020	3.437734	-0.115714	-0.033660
10	40.0	0.4	1.0	5.539348	5.670761	-0.131413	-0.023174
11	40.0	0.4	2.0	7.333656	7.516210	-0.182554	-0.024288
12	42.0	0.2	1.0	1.828810	1.907586	-0.078776	-0.041296
13	42.0	0.2	2.0	2.613911	2.724794	-0.110882	-0.040694
14	42.0	0.4	1.0	4.741304	4.924880	-0.183575	-0.037275
15	42.0	0.4	2.0	6.670912	6.833726	-0.162814	-0.023825
16	44.0	0.2	1.0	1.251240	1.348671	-0.097431	-0.072243
17	44.0	0.2	2.0	2.028939	2.146882	-0.117943	-0.054937
18	44.0	0.4	1.0	4.075738	4.270971	-0.195232	-0.045711
19	44.0	0.4	2.0	6.016795	6.214218	-0.197423	-0.031770

	Stock Price	Volatility	Expiration	Option Price
0	36.0	0.2	1.0	4.473118
1	36.0	0.2	2.0	4.804318
2	36.0	0.4	1.0	6.981569
3	36.0	0.4	2.0	8.380087
4	38.0	0.2	1.0	3.223645
5	38.0	0.2	2.0	3.656794
6	38.0	0.4	1.0	6.007105
7	38.0	0.4	2.0	7.515862
8	40.0	0.2	1.0	2.227884
9	40.0	0.2	2.0	2.763844
10	40.0	0.4	1.0	5.157563
11	40.0	0.4	2.0	6.750051
12	42.0	0.2	1.0	1.532825
13	42.0	0.2	2.0	2.108770
14	42.0	0.4	1.0	4.417692
15	42.0	0.4	2.0	6.077077
16	44.0	0.2	1.0	1.029997
17	44.0	0.2	2.0	1.595208
18	44.0	0.4	1.0	3.775644
19	44.0	0.4	2.0	5.454867

	Stock Price	Volatility	Expiration	Option Price
0	70.0	0.2	1.0	0.118699
4	80.0	0.2	1.0	0.650706
8	90.0	0.2	1.0	2.284714
12	100.0	0.2	1.0	5.732011
16	110.0	0.2	1.0	11.494249
22	120.0	0.2	1.0	19.230460
26	130.0	0.2	1.0	28.225184

	Stock Price	Volatility	Expiration	Option Price
0	80.0	0.2	0.5	0.214788
4	90.0	0.2	0.5	1.360938
8	100.0	0.2	0.5	4.708581
12	110.0	0.2	0.5	10.999584
16	120.0	0.2	0.5	20.000000

	Stock Price	Volatility	Expiration	Option Price
0	80.0	0.2	0.5	20.956660
4	90.0	0.2	0.5	12.633988
8	100.0	0.2	0.5	6.364592
12	110.0	0.2	0.5	2.650154
16	120.0	0.2	0.5	0.918096