Leisen Reimer

The Leisen Reimer (1996)

The Leisen Reimer (1996) model is noted for its rapid convergence. (Reference here). Valuation of European options with the Black-Scholes formula is a long-solved problem . Efficient and accurate analogues valuation for American options using a Black-Scholes style closed form solution remains a challenging problem. Zhu (2006) developed an analytical formula for the value of such an option. The proposed method however incorporates two infinite sums of infinite double integrals which slows down computational time. The difficulties associated with obtaining analytical results have rendered numerical methods the natural choice to pricing American options. Monte Carlo simulation is of less use for pricing American options and the rate of convergence for Monte Carlo simulation can be relatively slow. Amin and Khanna (1994) demonstrated that the Cox, Ross and Rubinstein (1979) tree converges to the correct price. The rate of convergence was later studied by Leisen (1998) who showed that the CRR tree converges with order 1. Staunton (2005) benchmarked several of these approximations to price American options, including: two analytic approximations due to Ju (1998) and Ju and Zhong (1999); the explicit finite difference method; the implicit finite difference method; and the Leisen and Reimer (1996) binomial tree with and without truncation. Where applicable these techniques were tested with and without Richardson extrapolation. Considering time versus error, Staunton concluded that the Leisen-Reimer tree with Richardson extrapolation and truncation was the most effective method, in preference to other competing methods, including those based on numerical integration. Joshi (2009) provides an overview of varying formulae, worth taking a look at before proceeding here. We initially present VBA code for a relatively faster Dynamic- zero Truncated Leisen Reimer binomial estimation based on Espen Haug. Python code and R code is also based on Espen haug tree Design which incorporates dynamic memory usage important for speed.

VBA Code for Leisen Reimer with Truncated zero region and dynamic memory

You might consider this this video to confirm code below is consistent with results from original Leisen Reimer (1996) paper results. Also, I provide some additional analysis of the zero region in the Leisen Lattice lattice here and here.

Public Function LeisenReimerTrunc(AmeEurFlag As String, CallPutFlag As String, S As Double, X As Double, _

T As Double, r As Double, b As Double, v As Double, n As Integer) As Variant

Dim OptionValue() As Double

'Dim ReturnValue(3) As Double

Dim d1 As Double, d2 As Double

Dim hd1 As Double, hd2 As Double

Dim u As Double, d As Double, p As Double

Dim dt As Double, Df As Double

Dim i As Integer, j As Integer, z As Integer

n = Application.Odd(n)

ReDim OptionValue(0 To n)

If CallPutFlag = "c" Then

z = 1

ElseIf CallPutFlag = "p" Then

z = -1

End If

d1 = (Log(S / X) + (b + v ^ 2 / 2) * T) / (v * Sqr(T))

d2 = d1 - v * Sqr(T)

'// Using Preizer-Pratt inversion method 2

hd1 = 0.5 + Sgn(d1) * (0.25 - 0.25 * Exp(-(d1 / (n + 1 / 3 + 0.1 / (n + 1))) ^ 2 * (n + 1 / 6))) ^ 0.5

hd2 = 0.5 + Sgn(d2) * (0.25 - 0.25 * Exp(-(d2 / (n + 1 / 3 + 0.1 / (n + 1))) ^ 2 * (n + 1 / 6))) ^ 0.5

dt = T / n

p = hd2

u = Exp(b * dt) * hd1 / hd2

d = (Exp(b * dt) - p * u) / (1 - p)

Df = Exp(-r * dt)

For i = 0 To n

OptionValue(i) = Application.Max(0, z * (S * u ^ i * d ^ (n - i) - X))

Next

If z = 1 Then

For j = n - 1 To ((n + 1) / 2) Step -1

For i = (j - ((n - 1) / 2)) To j

If AmeEurFlag = "e" Then

OptionValue(i) = (p * OptionValue(i + 1) + (1 - p) * OptionValue(i)) * Df

ElseIf AmeEurFlag = "a" Then

OptionValue(i) = Application.Max((z * (S * u ^ i * d ^ (j - i) - X)), _

(p * OptionValue(i + 1) + (1 - p) * OptionValue(i)) * Df)

End If

Next

Next

ElseIf z = -1 Then

For j = n - 1 To ((n + 1) / 2) Step -1

For i = 0 To ((n - 1) / 2)

OptionValue((n + 1) / 2) = 0

If AmeEurFlag = "e" Then

OptionValue(i) = (p * OptionValue(i + 1) + (1 - p) * OptionValue(i)) * Df

ElseIf AmeEurFlag = "a" Then

OptionValue(i) = Application.Max((z * (S * u ^ i * d ^ (j - i) - X)), _

(p * OptionValue(i + 1) + (1 - p) * OptionValue(i)) * Df)

End If

Next

Next

End If

For j = ((n + 1) / 2) - 1 To 0 Step -1

For i = 0 To j + 1

If AmeEurFlag = "e" Then

OptionValue(i) = (p * OptionValue(i + 1) + (1 - p) * OptionValue(i)) * Df

ElseIf AmeEurFlag = "a" Then

OptionValue(i) = Application.Max((z * (S * u ^ i * d ^ (j - i) - X)), _

(p * OptionValue(i + 1) + (1 - p) * OptionValue(i)) * Df)

End If

Next

Next

LeisenReimerTrunc = OptionValue(0)

End Function

Python Code for Leisen Reimer

# https://github.com/tienusss/Option_Calculations

# Import packages

import numpy as np

# Defined functions

def LeisenReimerBinomial(OutputFlag, AmeEurFlag, CallPutFlag, S, X, T, r, c, v, n):

# This code is based on "The complete guide to Option Pricing Formulas" by Espen Gaarder Haug (2007)

# Translated from a VBA code

# OutputFlag:

# "P" Returns the options price

# "d" Returns the options delta

# "a" Returns an array containing the option value, delta and gamma

# AmeEurFlag:

# "a" Returns the American option value

# "e" Returns the European option value

# CallPutFlag:

# "C" Returns the call value

# "P" Returns the put value

# S is the share price at time t

# X is the strike price

# T is the time to maturity in years (days/365)

# r is the risk-free interest rate

# c is the cost of carry rate

# v is the volatility

# n determines the stepsize

# Start of the code

# rounds n up tot the nearest odd integer (the function is displayed below the LeisenReimerBinomial function in line x)

n = round_up_to_odd(n)

# Creates a list with values from 0 up to n (which will be used to determine to exercise or not)

n_list = np.arange(0, (n + 1), 1)

# Checks if the input option is a put or a call, if not it returns an error

if CallPutFlag == 'C':

z = 1

elif CallPutFlag == 'P':

z = -1

else:

return 'Call or put not defined'

# d1 and d2 formulas of the Black-Scholes formula for European options

d1 = (np.log(S / X) + (c + v ** 2 / 2) * T) / (v * np.sqrt(T))

d2 = d1 - v * np.sqrt(T)

# The Preizer-Pratt inversion method 1

hd1 = 0.5 + np.sign(d1) * (0.25 - 0.25 * np.exp(-(d1 / (n + 1 / 3 + 0.1 / (n + 1))) ** 2 * (n + 1 / 6))) ** 0.5

# The Preizer-Prat inversion method 2

hd2 = 0.5 + np.sign(d2) * (0.25 - 0.25 * np.exp(-(d2 / (n + 1 / 3 + 0.1 / (n + 1))) ** 2 * (n + 1 / 6))) ** 0.5

# Calculates the stepsize in years

dt = T / n

# The probability of going up

p = hd2

# The up and down factors

u = np.exp(c * dt) * hd1 / hd2

d = (np.exp(c * dt) - p * u) / (1 - p)

df = np.exp(-r * dt)

# Creates the most right column of the tree

max_pay_off_list = []

for i in n_list:

i = i.astype('int')

max_pay_off = np.maximum(0, z * (S * u ** i * d ** (n - i) - X))

max_pay_off_list.append(max_pay_off)

# The binominal tree

for j in np.arange(n - 1, 0 - 1, -1):

for i in np.arange(0, j + 1, 1):

i = i.astype(int) # Need to be converted to a integer

if AmeEurFlag == 'e':

max_pay_off_list[i] = (p * max_pay_off_list[i + 1] + (1 - p) * max_pay_off_list[i]) * df

elif AmeEurFlag == 'a':

max_pay_off_list[i] = np.maximum((z * (S * u ** i * d ** (j - i) - X)),

(p * max_pay_off_list[i + 1] + (1 - p) * max_pay_off_list[i]) * df)

if j == 2:

gamma = ((max_pay_off_list[2] - max_pay_off_list[1]) / (S * u ** 2 - S * u * d) - (

max_pay_off_list[1] - max_pay_off_list[0]) / (S * u * d - S * d ** 2)) / (

0.5 * (S * u ** 2 - S * d ** 2))

if j == 1:

delta = ((max_pay_off_list[1] - max_pay_off_list[0])) / (S * u - S * d)

price = max_pay_off_list[0]

# Put all the variables in the list

variable_list = [delta, gamma, price]

# Return values

if OutputFlag == 'P':

return price

elif OutputFlag == 'd':

return delta

elif OutputFlag == 'g':

return gamma

elif OutputFlag == 'a':

return variable_list

else:

return 'Indicate if you want to return P, d, g or a'

def round_up_to_odd(n):

# This function returns a number rounded up to the nearest odd integer

# For example when n = 100, the function returns 101

return np.ceil(n) // 2 * 2 + 1

#Body of the script

# LeisenReimerBinomial(OutputFlag, AmeEurFlag, CallPutFlag, S, X, T, r, c, v, n)

#Eur_call_result = LeisenReimerBinomial('P', 'e', 'C', 100, 100, 0.5, 0.07, 0.07, 0.3, 25)

American_call_result = LeisenReimerBinomial('P', 'a', 'C', 100, 100, 0.5, 0.07, 0.07, 0.3, 25)

#Eur_put_result = LeisenReimerBinomial('P', 'e', 'P', 100, 100, 0.5, 0.07, 0.07, 0.3, 25)

#American_put_result = LeisenReimerBinomial('P', 'a', 'P', 100, 100, 0.5, 0.07, 0.07, 0.3, 25)

#Print the output of the results

#print('The price of the European call option is equal to ' +str(Eur_call_result))

print('The price of the American call option is equal to ' +str(American_call_result))

#print('The price of the European put option is equal to ' +str(Eur_put_result))

#print('The price of the American put option is equal to ' +str(American_put_result))

Leisen Reimer R Code

# This code was originally taken from foptions rpackage

# fOptions however did not have a function for LR

# https://rdrr.io/rforge/fOptions/src/R/BinomialTreeOptions.R

# results check out with Haug

LRBinomialTreeOption <- function(TypeFlag, S, X, Time, r, b, sigma, n)

{ # A function implemented by Diethelm Wuertz - similar to haug style

# Description:

# Tian's Modification to the Binomial Tree Option

# FUNCTION:

# Check Flags:

TypeFlag = TypeFlag[1]

if (TypeFlag == "ce" || TypeFlag == "ca") z = +1

if (TypeFlag == "pe" || TypeFlag == "pa") z = -1

# Parameters:

# M = exp ( b*dt )

# V = exp ( sigma^2 * dt )

# u = (M*V/2) * ( V + 1 + sqrt(V*V + 2*V - 3) )

# d = (M*V/2) * ( V + 1 - sqrt(V*V + 2*V - 3) )

# p = (M-d)/(u-d)

if (n %% 2 == 0){

n <- (n + 1)

}

dt = Time/n

Method <- 2

d1 <- (logb(S/X) + (b + sigma*sigma / 2) * Time) / sigma / sqrt(Time)

d2 <- (d1 - sigma * sqrt(Time))

rq_n <- exp((b) * dt);

Term1 <- ((d1 / (n + 1.0 / 3.0 - (1 - Method) * 0.1 / (n + 1))) ^ 2) * (n + 1.0 / 6.0)

pp <- (0.5 + sign(d1) * 0.5 * sqrt(1 - exp(-Term1)))

Term2 <- ((d2 / (n + 1.0 / 3.0 - (1 - Method) * 0.1 / (n + 1))) ^ 2) * (n + 1.0 / 6.0)

p <- (0.5 + sign(d2) * 0.5 * sqrt(1 - exp(-Term2)))

u <- (rq_n * pp / p)

d <- ((rq_n - p * u) / (1 - p))

Df = exp(-r*dt)

# Iteration:

OptionValue = z*(S*u^(0:n)*d^(n:0) - X)

OptionValue = (abs(OptionValue) + OptionValue) / 2

# European Option:

if (TypeFlag == "ce" || TypeFlag == "pe") {

for ( j in seq(from = n-1, to = 0, by = -1) )

for ( i in 0:j )

OptionValue[i+1] =

(p*OptionValue[i+2] + (1-p)*OptionValue[i+1]) * Df }

# American Option:

if (TypeFlag == "ca" || TypeFlag == "pa") {

for ( j in seq(from = n-1, to = 0, by = -1) )

for ( i in 0:j )

OptionValue[i+1] = max((z * (S*u^i*d^(abs(i-j)) - X)),

(p*OptionValue[i+2] + (1-p)*OptionValue[i+1]) * Df) }

# Return Value:

OptionValue[1]

}

LRBinomialTreeOption("pa", 100, 100, 0.5, 0.07, 0.07, 0.3, 25)

Testing for True

Broadie and Detemple (1996) used a 15,000-step binomial model to obtain “True values”. Amin and Khanna (1994) demonstrated that the Cox, Ross and Rubinstein (1979) tree converges to the correct price.