Optimizing the Cox Ross and Rubinstein Algorithm

Optimizing Cox, Ross and Rubinstein in C++

Cox, Ross and Rubinstein (1979) (CRR) developed a binomial approach. It is widely accepted that greater accuracy can be introduced by making the CRR lattice mesh finer or by, more prosaically, incorporating a larger number of steps. Broadie and Detemple (1996) used a 15,000-step binomial model to obtain “True values”. The binomial model is therefore acknowledged as a reliable workhorse and serves to benchmark other techniques. This lattice method, however, is distinctly viewed as having substantial computational workloads, runtimes and gluttonous memory requirements. Making use of Dynamic Memory Storage can produce an important reduction in computational cost. Using the three segments of code below, we try to reveal how computer memory can be used more efficiently. A conventional two-dimensional static n-step binomial model requires (n + 1)(n + 2) / 2 nodes to be memorized. Broadie and Detemple (1996) and Haug (1997) propose using a one-dimensional dynamic binomial tree. This approach takes the option values at the last column and stores them in a dynamic vector Opt(j) for j = 0, 1, … , n. After moving one step back, the values in the re-dimensioned Opt(j) for j = 0, 1, …, n - 1 will be replaced by the option values of the corresponding nodes at the penultimate column. Similarly, the values of Opt(j) for j = 0, 1, …, k - 1 at k^th column will always be substituted by the option values at (k - 1)^th column for 1 ≤ k ≤ n. Therefore, a dynamic binomial tree only requires n + 1 contemporaneous storage spaces. To make the best use of computer memory we will compare the performance of three segments of code that incorporate the following

(1) CRR Lattice with Static Memory Storage design

(2) CRR Lattice with Dynamic Memory design

(3) CRR Lattice with Dynamic Memory design and a truncated zero region

In the following 3 snippets of C++ code we will examine the accuracy and performance in terms of execution speed. It should become clear that optimizing C++ code is important for making better use of computing resources. For a deep dive down the rabbit hole explaining (1), (2) and (3) CRR algorithms with a fairly dense navigation of array mapping, backward recursion and truncation - please follow link. Otherwise brief outline below may suffice.

(1) CRR Lattice with Static Memory design

This code was adapted from Yinqiu Zhang and is based on Fabrice Rouah static storage Lattice design.

//This code must have come from Yinqiu Zhang and is based on Fabrice Rouah static Lattice design

#include <vector>

#include <stdio.h>

#include <math.h>

#include <iostream>

#include <string>

#include <algorithm>

#include<iomanip>

#include <time.h>

using namespace std;

// 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;

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 = 5000;

//double Spot = ii * 10;

double Spot = 100;

double K = 100;

double r = 0.03;

double q = 0.07;

double v = 0.2;

double T = 3;

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");

}

OnlineGBD and CPU time

The OnlineGBD C++ estimation time can vary. The values I will present below for CPU time are representative. To get some sense of speed you need to run C++ snippets on your own machine. Table 2 below presents a CRR valuation for an American Call at 15,000 step size.

Table 2 from from Broadie and Detemple (1996)

(2) CRR Lattice with Dynamic Memory design

//C++ code based on Espen Haug and Broadie and Detemple (1996) Dynamic Memory design.

//Tested for American Call Option and giving the same results for Broadie Detemple (1996) 9.066 p. 1224 Table 2

//The binomial price is 9.06594 where

//int n = 15000; // Number of steps

//double S = 100.0; // Spot Price

//double K = 100.0; // Strike Price

//double T = 3; // Years to maturity

//double r = 0.03; // Risk Free Rate

//double q = 0.07; //double q = 0.03;

//double v = 0.20;

//char PutCall = 'C';

//char OpStyle = 'A';

#include <stdio.h>

#include <math.h>

#include <iostream>

#include <algorithm>

#include <iostream>

#include <vector>

#include<iomanip>

#include <string>

#include <time.h>

//using namespace System;

using namespace std;

// 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;

double dt, u, d, p;

int z;

// Quantities for the tree

dt = T / n;

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 column

OptionValue.resize(n + 1);

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

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

//OptionValue[i] = max(z*(S*pow(u, (2 * i - n)) - 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])));

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

}

}

// Return the option price

return OptionValue[0];

}

int main() {

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

//char PutCall, OpStyle;

//int n;

clock_t start_time, end_time;

start_time = clock();

int n = 5000; // Number of steps

double S = 100.0; // Spot Price

double K = 100.0; // Strike Price

double T = 3; // Years to maturity

double r = 0.03; // Risk Free Rate

double q = 0.07; //double q = 0.03; // q added

double v = 0.20;

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) CRR Lattice with Dynamic Memory design and a truncated zero region

//Dynamic Binomial with Truncated zeroes

//Faster with power term

//Got rid of z

//Working.

//Tested for American Call Option and giving the same results for Broadie Detemple (1996) 9.066 p. 1224 Table 2

//The binomial price is 9.06594 where

//int n = 15000; // Number of steps

//double S = 100.0; // Spot Price

//double K = 100.0; // Strike Price

//double T = 3; // Years to maturity

//double r = 0.03; // Risk Free Rate

//double q = 0.07; //double q = 0.03;

//double v = 0.20;

//char PutCall = 'C';

//char OpStyle = 'A';

#include <stdio.h>

#include <math.h>

#include <iostream>

#include <algorithm>

#include <iostream>

#include <vector>

#include<iomanip>

#include <string>

#include <time.h>

//using namespace System;

using namespace std;

// 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;

double dt, u, d, p, Transfer;

int Num, Start;

// Quantities for the tree

dt = T / n;

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

d = 1.0 / u;

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

Num = int(log(K / (S * pow(d, n))) / log(u / d));

//Qianru Shang adapted this section for preliminary tests

//MacDonald Schroeder

if (PutCall == 'P')

{

Transfer = S; //Exchange the value of S and K4

S = K;

K = Transfer;

Transfer = r; //Exchange the value of Interest and Dividends

r = q;

q = Transfer;

PutCall = 'C';

}

if (PutCall == 'C'){

vector<double> OptionValue;

//resize the column

OptionValue.resize(n + 1);

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

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

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

}

// Backward recursion through the tree

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

//Qianru Shang adapted this section for preliminary tests

Start = Num - (n - 1 - j); //Truncation: omit the zero are

if (j < (n - Num))

{

Start = 0;

}

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

if (OpStyle == 'E')

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

else {

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

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

}

}

}

// Return the option price

return OptionValue[0];

}

}

int main() {

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

//char PutCall, OpStyle;

//int n;

clock_t start_time, end_time;

start_time = clock();

int n = 5000; // Number of steps

double S = 100.0; // Spot Price

double K = 100.0; // Strike Price

double T = 3; // Years to maturity

double r = 0.03; // Risk Free Rate

double q = 0.07; //double q = 0.03; // q added

double v = 0.20;

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");

}