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