Leisen Reimer C++ Optimization and Spreadsheet Modeling

Some virtues of Leisen Reimer (1996)

Below, we present three C++ implementations for the Leisen-Reimer binomial tree. The Leisen Reimer (1996) model is noted for its rapid convergence. In the video presented here we implement the VBA code found in the preceding page and the C++ static code found in the first snippet below. It is quite clear Leisen Reimer (1996) converges more quickly to true relative to Cox, Ross and Rubinstein.

C++ Implementation with three level of Optimization

Below, we present three C++ implementations for the Leisen-Reimer binomial tree. The Leisen Reimer (1996) model is noted for its relatively rapid convergence. Below we present three snippets of C++ code: (1) Static, (2) Dynamic and (3) Dynamic- zero region Truncated. The latter being the most optimized. In the video just below, we look at the Leisen Reimer binomial option pricing model. I consider three versions of the model set up in C++ (1) - (3). We go to Table 1 of the original Leisen Reimer (1996) paper. It is clear that the (3) is the fastest for computation and this secured by making use of Dynamic memory and making redundant the zero region. The static model is the slowest and cannot reach 25,000 steps. I also illustrate that the Leisen Reimer tree has a different shape to the Cox Ross and Rubinstein tree. This can be illustrated in excel.

(1) Leisen Reimer Static (C++ Code)

Static Memory storage tends to produce less efficient outcomes. If the entire lattice has to be saved at once this imposes a cost when running any estimation. We include the static code for completeness but would advise making use of the Dynamic approach proposed by Haug

//Static Leisen Reimer from Volopta

//Changed to include static Leisen Reimer

#include <vector>

#include <stdio.h>

#include <math.h>

#include <iostream>

#include <string>

#include <algorithm>

#include<iomanip>

#include <time.h>

using namespace std;

// Returns the sign of a number

double Sgn(double x) {

if (x<0) return -1;

else if (x>0) return 1;

else if (x == 0) return 0;

}

// Function for binomial tree

double Binomial(int n, double Spot, double K, double r, double q, double v, double T, char PutCall, char OpStyle)

{

int i, j; // Counters

if (n % 2 == 0) // Use only odd number of steps

n = n + 1;

//double dt, u, d, p;

int z;

int Method =2;

// Quantities for the tree

double dt = T / n;

double d1 = (log(Spot / K) + (r - q + v*v / 2) * T) / v / sqrt(T); // D Plus

double d2 = (log(Spot / K) + (r - q - v*v / 2) * T) / v / sqrt(T); // D Minus

double rq_n = exp((r - q) * dt);

double Term1 = pow((d1 / (n + 1.0 / 3.0 - (1 - Method) * 0.1 / (n + 1))), 2) * (n + 1.0 / 6.0);

double pp = 0.5 + Sgn(d1) * 0.5 * sqrt(1 - exp(-Term1));

double Term2 = pow((d2 / (n + 1.0 / 3.0 - (1 - Method) * 0.1 / (n + 1))), 2) * (n + 1.0 / 6.0);

double p = 0.5 + Sgn(d2) * 0.5 * sqrt(1 - exp(-Term2));

//changed

double u = rq_n * pp / p;

double d = (rq_n - p * u) / (1 - p);

//Original Static CRR

// int i, j;

vector<vector<double> > S(n + 1, vector<double>(n + 1));

vector<vector<double> > Op(n + 1, vector<double>(n + 1));

// double dt, u, d, p;

// dt = T / n;

// u = exp(v*sqrt(dt));

// d = 1 / u;

// p = (exp((r- q)*dt) - d) / (u - d);

// Build the binomial tree

for (j = 0; j <= n; j++) {

for (i = 0; i <= j; i++) {

S[i][j] = Spot*pow(u, j - i)*pow(d, i);

}

}

// Compute terminal payoffs

for (i = 0; i <= n; i++) {

if (PutCall == 'C')

Op[i][n] = max(S[i][n] - K, 0.0);

else

Op[i][n] = max(K - S[i][n], 0.0);

}

// Backward recursion through the tree

for (j = n - 1; j >= 0; j--) {

for (i = 0; i <= j; i++) {

if (OpStyle == 'E')

Op[i][j] = exp(-r*dt)*(p*(Op[i][j + 1]) + (1 - p)*(Op[i + 1][j + 1]));

else {

if (PutCall == 'C')

Op[i][j] = max(S[i][j] - K, exp(-r*dt)*(p*(Op[i][j + 1]) + (1 - p)*(Op[i + 1][j + 1])));

else

Op[i][j] = max(K - S[i][j], exp(-r*dt)*(p*(Op[i][j + 1]) + (1 - p)*(Op[i + 1][j + 1])));

}

}

}

return Op[0][0];

}

int main()

{

clock_t start_time, end_time;

start_time = clock();

//for (double ii = 0; ii <= 20; ii++) {

int n = 15;

//double Spot = ii * 10;

double Spot = 100;

double K = 100;

double r = 0.07;

double q = 0.0;

double v = 0.3;

double T = 0.5;

char PutCall = 'C';

char OpStyle = 'A';

cout << setprecision(10);

cout << Spot << " " << "The binomial price is " << Binomial(n, Spot, K, r, q, v, T, PutCall, OpStyle) << endl;

end_time = clock();

cout << " " << start_time << " " << end_time << " " << (end_time - start_time) << endl;

cout << " " << start_time << " " << end_time << " " << (end_time - start_time) / (double)CLOCKS_PER_SEC << " seconds" << endl;

system("PAUSE");

}

(2) Leisen Reimer Dynamic (C++ Code)

//This is code for Leisen Reimer but based on dynamic programming

//The code does not truncate the zero region

//Make sure both methods are set to Method 1 to get same results

#include <stdio.h>

#include <math.h>

#include <iostream>

#include <algorithm>

#include <iostream>

#include <vector>

#include<iomanip>

#include <string>

//using namespace System;

using namespace std;

// Returns the sign of a number

double Sgn(double x) {

if (x<0) return -1;

else if (x>0) return 1;

else if (x == 0) return 0;

}

// Function for binomial tree

double Binomial(int n, double S, double K, double r, double q, double v, double T, char PutCall, char OpStyle) {

int i, j; // Counters

if (n % 2 == 0) // Use only odd number of steps

n = n + 1;

//double dt, u, d, p;

int z;

int Method =2;

// Quantities for the tree

double dt = T / n;

double d1 = (log(S / K) + (r - q + v*v / 2) * T) / v / sqrt(T); // D Plus

double d2 = (log(S / K) + (r - q - v*v / 2) * T) / v / sqrt(T); // D Minus

double rq_n = exp((r - q) * dt);

double Term1 = pow((d1 / (n + 1.0 / 3.0 - (1 - Method) * 0.1 / (n + 1))), 2) * (n + 1.0 / 6.0);

double pp = 0.5 + Sgn(d1) * 0.5 * sqrt(1 - exp(-Term1));

double Term2 = pow((d2 / (n + 1.0 / 3.0 - (1 - Method) * 0.1 / (n + 1))), 2) * (n + 1.0 / 6.0);

double p = 0.5 + Sgn(d2) * 0.5 * sqrt(1 - exp(-Term2));

//changed

double u = rq_n * pp / p;

double d = (rq_n - p * u) / (1 - p);

//u = exp(v*sqrt(dt));

//d = 1.0 / u;

//p = (exp((r - q)*dt) - d) / (u - d);

if (PutCall == 'C')

{

z = 1;

}

else if (PutCall == 'P')

{

z = -1;

}

vector<double> OptionValue;

//resize the row of the sub-vector: (1,2,3..)

OptionValue.resize(n + 1);

// Calculating Intrinsic Value

for (i = 0; i <= n; i++) {

OptionValue[i] = max(z*(S*pow(u, i)*pow(d, n - i) - K), 0.0);

}

// Backward recursion through the tree

for (j = n - 1; j >= 0; j--)

for (i = 0; i <= j; i++) {

if (OpStyle == 'E')

OptionValue[i] = exp(-r*dt)*(p*(OptionValue[i + 1]) + (1.0 - p)*(OptionValue[i]));

else {

OptionValue[i] = max(z*(S*pow(u, i)*pow(d, j - i) - K), exp(-r*dt)*(p*(OptionValue[i + 1]) + (1.0 - p)*(OptionValue[i])));

}

}

// Return the option price

return OptionValue[0];

}

int main() {

clock_t start_time, end_time;

start_time = clock();

//double S, K, T, v, r;

//char PutCall, OpStyle;

//int n;

int n = 25000; // Number of steps

double S = 100; // Spot Price

double K = 110.0; // Strike Price

double T = 0.5; // Years to maturity

double r = 0.07; // Risk Free Rate

double q = 0.0; // q added

double v = 0.30;

char PutCall = 'C';

char OpStyle = 'A';

cout << setprecision(10);

cout << "The binomial price is " << Binomial(n, S, K, r, q, v, T, PutCall, OpStyle) << endl;

end_time = clock();

cout << " " << start_time << " " << end_time << " " << (end_time - start_time) << endl;

cout << " " << start_time << " " << end_time << " " << (end_time - start_time) / (double)CLOCKS_PER_SEC << " seconds" << endl;

//system("PAUSE");

}

(3) Leisen Reimer Dynamic and Truncated (C++ Code)

// Leisen Reimer dynamic and truncated zeros

// Seems to tally with Leisen Reimer paper

#include <stdio.h>

#include <math.h>

#include <iostream>

#include <algorithm>

#include <iostream>

#include <vector>

#include<iomanip>

#include <string>

//using namespace System;

using namespace std;

// Returns the sign of a number

double Sgn(double x) {

if (x<0) return -1;

else if (x>0) return 1;

else if (x == 0) return 0;

}

// Function for binomial tree

double Binomial(int n, double S, double K, double r, double q, double v, double T, char PutCall, char OpStyle) {

int i, j; // Counters

if (n % 2 == 0) // Use only odd number of steps

n = n + 1;

//double dt, u, d, p;

int z;

int Method =2;

// Quantities for the tree

double dt = T / n;

double d1 = (log(S / K) + (r - q + v*v / 2) * T) / v / sqrt(T); // D Plus

double d2 = (log(S / K) + (r - q - v*v / 2) * T) / v / sqrt(T); // D Minus

double rq_n = exp((r - q) * dt);

double Term1 = pow((d1 / (n + 1.0 / 3.0 - (1 - Method) * 0.1 / (n + 1))), 2) * (n + 1.0 / 6.0);

double pp = 0.5 + Sgn(d1) * 0.5 * sqrt(1 - exp(-Term1));

double Term2 = pow((d2 / (n + 1.0 / 3.0 - (1 - Method) * 0.1 / (n + 1))), 2) * (n + 1.0 / 6.0);

double p = 0.5 + Sgn(d2) * 0.5 * sqrt(1 - exp(-Term2));

//changed

double u = rq_n * pp / p;

double d = (rq_n - p * u) / (1 - p);

//CRR

//u = exp(v*sqrt(dt));

//d = 1.0 / u;

//p = (exp((r - q)*dt) - d) / (u - d);

if (PutCall == 'C')

{

z = 1;

}

else if (PutCall == 'P')

{

z = -1;

}

vector<double> OptionValue;

//resize the row of the sub-vector: (1,2,3..)

OptionValue.resize(n + 1);

// Calculating Intrinsic Value

for (i = 0; i <= n; i++) {

OptionValue[i] = max(z*(S*pow(u, i)*pow(d, n - i) - K), 0.0);

}

if (z ==1)

// Backward recursion through the tree

// The tree is segmented for truncating zeros

for (j = n - 1; j >= (n+1)/2; j--)

for (i = (j - ((n+1)/2)); i <= j; i++) {

if (OpStyle == 'E')

OptionValue[i] = exp(-r*dt)*(p*(OptionValue[i + 1]) + (1.0 - p)*(OptionValue[i]));

else {

OptionValue[i] = max(z*(S*pow(u, i)*pow(d, j - i) - K), exp(-r*dt)*(p*(OptionValue[i + 1]) + (1.0 - p)*(OptionValue[i])));

}

}

else if (z == -1)

for (j = n - 1; j >= (n+1)/2; j--)

for (i = 0; i <= (n+1)/2; i++) {

OptionValue[(n+1)/2]=0;

if (OpStyle == 'E')

OptionValue[i] = exp(-r*dt)*(p*(OptionValue[i + 1]) + (1.0 - p)*(OptionValue[i]));

else {

OptionValue[i] = max(z*(S*pow(u, i)*pow(d, j - i) - K), exp(-r*dt)*(p*(OptionValue[i + 1]) + (1.0 - p)*(OptionValue[i])));

}

}

for (j = ((n+1)/2)-1; j >= 0; j--)

for (i =0; i <= j + 1; i++) {

if (OpStyle == 'E')

OptionValue[i] = exp(-r*dt)*(p*(OptionValue[i + 1]) + (1.0 - p)*(OptionValue[i]));

else {

OptionValue[i] = max(z*(S*pow(u, i)*pow(d, j - i) - K), exp(-r*dt)*(p*(OptionValue[i + 1]) + (1.0 - p)*(OptionValue[i])));

}

}

// Return the option price

return OptionValue[0];

}

int main() {

clock_t start_time, end_time;

start_time = clock();

//double S, K, T, v, r;

//char PutCall, OpStyle;

//int n;

int n = 10000; // Number of steps

double S = 100; // Spot Price

double K = 110.0; // Strike Price

double T = 0.5; // Years to maturity

double r = 0.07; // Risk Free Rate

double q = 0.0; // q added

double v = 0.30;

char PutCall = 'C';

char OpStyle = 'E';

cout << setprecision(10);

cout << "The binomial price is " << Binomial(n, S, K, r, q, v, T, PutCall, OpStyle) << endl;

end_time = clock();

cout << " " << start_time << " " << end_time << " " << (end_time - start_time) << endl;

cout << " " << start_time << " " << end_time << " " << (end_time - start_time) / (double)CLOCKS_PER_SEC << " seconds" << endl;

//system("PAUSE");

}

A Deeper Dive in Leisen Reimer (1996) Modeling

In the video below, we develop a little bit more the main features that distinguish Leisen Reimer from Cox, Ross and Rubinstein. The zero region in the Leisen Reimer tree is fixed. It does not adjust when we adjust the parameter values. This is different to Cox, Ross and Rubinstein. I use the manual trees set out in excel to capture the main differences. Making the zero region redundant through truncation will all else being equal speed up lattice estimation. In the exposition here it is possible to unfurl some properties linked to the early exercise boundaries for both Leisen reimer and Cox Ross and Rubinstein. This is relevant for understanding Intelligent lattice Search.