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 Rubinstein1 (CRR) and the least-squares Monte Carlo method of Longstaff & Schwartz 2 (LS). The design, implementation and output of both are validated by comparing our results to those produced in the academic literature.
1 Cox, J., S. Ross, and M. Rubinstein (1979): “Option Pricing: A Simplified Approach.” Journal
of Financial Economics, Volume 7, Issue 3, pp. 229–263.
2 Longstaff, Francis & Schwartz, Eduardo. (2001). Valuing American Options by Simulation: A Simple Least-Squares Approach. Review of Financial Studies. 14. 113-47.
The valuation of an American-style option using the Longstaff Shwarz approach requires the Monte Carlo simulation of realizations of the asset path from now to expiration using time steps of . If the random walk is lognormal then we will simulate
Where is a Normally distributed random variable, is the risk free rate, is volatility and and is the stock price after 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 . These are paths 1, 3, 4, 6, and 7. Obviously at the final exercise date, , 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 , the optimal strategy is to compare the immediate exercise value with the expected cash flows from continuing , and then exercise if immediate exercise is more valuable, i.e. where . Conversely, one would continue to hold where the continuation value is more valuable, i.e. where
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:
Where is the stock price at the 2-year point, 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:
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)
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))
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)
# 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
The model results are aligned with "Simulated American" values produced by Longman & Schwartz, indicating that the model is correctly constructed and specified.
The model results are aligned with "LSM" values produced by Zhao3 , indicating that the model is correctly constructed and specified.
3 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:
# 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)
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.
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:
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 ():
The payoff at each up and down node - the sum of (delta*stock) PLUS the intrinsic value of option - will be the same.
Delta is the fraction which equates the payoff at up and down nodes.
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:
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:
Or, equivalently:
Which, upon removing the terms, simplifies to:
Which can be expressed in terms of a risk neutral probability:
Where:
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))
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:
# 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)
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')
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.