# Unable to correctly implement the pricing of an American call with multiple discrete dividends using the Clenshaw-Curtis quadrature

I'm not a quant, just an enthusiast. I am trying to implement in C++ the methodology published in the paper "Fast Quadrature Methods for Options with Discrete Dividends", by Thakoor and Bhuruth.

I have succeeded in implementing the pricing of an European call with multiple dividends and my results are identical to the author's. However, in trying to implement the pricing of an American call with multiple dividends, I am getting results which are slightly off. Using the Vellekoop & Nieuwenhuis method, I am able to get exactly the results published by Thakoor & Bhuruth (they use this method as a comparison).

This is my code:

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

double normalCDF(double value)
{
return 0.5 * erfc(-value * M_SQRT1_2);
}

// Utility function to compute differences between consecutive elements in a vector
VectorXd diff(const VectorXd& vec, double T) {
ArrayXd extendedVec(vec.size() + 2);  // Extend the size to accommodate 0.0 and T
extendedVec << 0.0, vec, T;          // Prepend 0.0 and append T
return extendedVec.tail(extendedVec.size() - 1) - extendedVec.head(extendedVec.size() - 1);
}

// Computes Clenshaw-Curtis (aka Fejér) nodes and weights
void clenshaw_curtis_compute(int order, VectorXd& x, VectorXd& w) {
if (order < 1) {
cerr << "clenshaw_curtis_compute - Fatal error! Illegal value of ORDER = " << order << "\n";
exit(1);
}
x.resize(order);
w.resize(order);

if (order == 1) {
x[0] = 0.0;
w[0] = 2.0;
} else {
#pragma omp parallel for
for (int i = 0; i < order; i++) {
x[i] = cos((double)(order - 1 - i) * M_PIl / (double)(order - 1));
}

x[0] = -1.0;
if (order % 2 == 1) {
x[(order - 1) / 2] = 0.0;
}
x[order - 1] = 1.0;

#pragma omp parallel for
for (int i = 0; i < order; i++) {
double theta = (double)i * M_PIl / (double)(order - 1);
w[i] = 1.0;

#pragma omp parallel for
for (int j = 1; j <= (order - 1) / 2; j++) {
double b = (2 * j == (order - 1)) ? 1.0 : 2.0;
w[i] -= b * cos(2.0 * j * theta) / (4.0 * j * j - 1);
}
}

w[0] /= (double)(order - 1);
for (int i = 1; i < order - 1; i++) {
w[i] *= 2.0 / (double)(order - 1);
}
w[order - 1] /= (double)(order - 1);
}
}

int main() {

// Parameters for option pricing
double r = 0.05;     // Risk-free rate
double sigma = 0.2;  // Volatility
double S0 = 100.0;   // Initial stock price
double K = 100.0;    // Strike price
double T = 1.0;      // Time to maturity
const double Xi = 6.5;     // Number of standard deviations

// Dividend details
ArrayXd D(4); D << 5.0, 5.0, 5.0, 5.0;  // Dividend values
ArrayXd tD(4); tD << 0.2, 0.4, 0.6, 0.8;  // Dividend payment times

const double N = 1000;  // Number of steps or nodes in numerical method

// Validation of input parameters
if (sigma <= 0 || S0 <= 0 || K <= 0 || T <= 0) {
cerr << "Error: Invalid parameters!\n";
return -1;
}

// Check for consistent input sizes
if (tD.size() != D.size()) {
cerr << "Error: Mismatch between dividend times and amounts.\n";
return -1;
}

// Precompute constants used in the pricing formula
const double sigma_sqr = sigma * sigma;
const double twoPiVariance = 2.0 * M_PIl * sigma_sqr;
const double mu = r - 0.5 * sigma_sqr;  // Drift term adjusted for volatility

// Compute time intervals for option pricing steps
const u_int64_t m = tD.size();
const ArrayXd tau = diff(tD, T);

//ArrayXd time(tD.size()+2); time << 0.0, tD.array(), T;
ArrayXd fnxkLU = VectorXd::Zero(N);

// Compute Clenshaw-Curtis nodes and weights
VectorXd xi, wi;
clenshaw_curtis_compute(N, xi, wi);

// Compute maximum and minimum stock prices for range of integration or approximation
double U = S0 * exp(mu * T + Xi * sigma * sqrt(T));
double L = max(S0 * exp(mu * T - Xi * sigma * sqrt(T)), *min_element(D.begin(), D.end()));

// Calculate xiLU for cum-dividend stock prices
ArrayXd xiLU = 0.5 * (U - L) * xi.array() + 0.5 * (U + L);

double c1 = (r + sigma_sqr / 2) * tau(tau.size() - 1);
double c2 = sigma * sqrt(tau(tau.size() - 1));
double c3 = K * exp(-r * tau(tau.size() - 1));
ArrayXd xiLUplus = xiLU - D(D.size() - 1);

ArrayXd d1 = (Eigen::log(xiLUplus / K) + c1) / c2;
ArrayXd d2 = d1 - c2;

ArrayXd normal_cdf_d1 = d1.unaryExpr(&normalCDF);
ArrayXd normal_cdf_d2 = d2.unaryExpr(&normalCDF);
VectorXd V_Ac = xiLUplus * normal_cdf_d1 - c3 * normal_cdf_d2;

// Calculate the critical values Khat
ArrayXd Khat(tD.size());
for (int k = 0; k < tD.size() - 1; ++k) {
Khat(k) = K * (1 - exp(-r * (tD(k+1) - tD(k))));
}
Khat(tD.size() - 1) = K * (1 - exp(-r * (T - tD(tD.size() - 1))));

for (int k = tD.size() - 1; k >= 0; --k) {

xiLUplus = xiLU - D[k];

ArrayXd fnxkLU(xiLUplus.size());

// Computing the discounted lognormal probability density function of the stock price
#pragma omp parallel for
for (int i = 0; i < xiLU.size(); ++i) {
Eigen::ArrayXd term = (log(xiLU / xiLUplus[i]) - mu * tau[k]).square() / (2 * sigma_sqr * tau[k]);
Eigen::VectorXd weights = V_Ac.array() / sqrt(twoPiVariance * tau[k]) / xiLU * exp(-term);
fnxkLU[i] = weights.dot(wi);
}

Eigen::ArrayXd V_Ac_continue = exp(-r * tau[k]) * 0.5 * (U - L) * fnxkLU;   // This is the continuation value

if (D[k] > Khat[k]) {
Eigen::ArrayXd V_Ac_exercise = (xiLUplus + D[k] - K).cwiseMax(0.0); // This is the exercise value
V_Ac = (V_Ac_continue).cwiseMax(V_Ac_exercise);
}
else {
V_Ac = V_Ac_continue;
}
}

V_Ac = V_Ac.array() * 1 / (xiLU * sqrt(twoPiVariance * tau[0])) * Eigen::exp(-Eigen::pow(Eigen::log(xiLU / S0) - mu * tau[0], 2) / (2.0 * sigma_sqr * tau[0]));
double V_0 = std::exp(-r * tau[0]) * 0.5 * (U - L) * V_Ac.dot(wi);

cout << "Price of the option at t=0: " << V_0 << endl;

return 0;
}

// Compilation: g++ -std=c++20 -O3 -march=native -o fejer_AC fejer_AC.cpp -I/usr/local/include/eigen3


With the parameters above ($$r=0.05$$, $$\sigma=0.2$$, $$S_0=100.0$$, $$K=100.0$$, $$T=1.0$$, $$\xi=6.5$$, $$N=1000$$, $$dividends = \{5.0, 5.0, 5.0, 5.0\}$$ and $$dividend\:dates = \{0.2, 0.4, 0.6, 0.8\}$$ ), I am getting a price of the option at $$t=0$$ of 4.8881, while the original paper obtains 4.8629.

My implementation of the Vellekoop & Nieuwenhuis method leads to 4.8629, which is in agreement with the author's implementation of that method.

I'm stuck and if anyone can point me in the right direction I'd be very grateful!

Cheers!

Adenda 1: I tried the exact same algorithm in Octave and I got the exact same small error. This leads me to believe that it's not a numerical problem, but rather a misunderstanding in the logic.

Adenda 2: I read and re-read Hull's chapter on pricing of American options with dividends, as well as the seminal paper by Andricopoulos et al (Universal option pricing using quadrature methods), but i still can's figure what's wrong with the code.

• This code is for American options or European? What is the difference between the two codes? May 14 at 18:50
• @nbbo2 This code is an attempt at pricing an American call. The European call version doesn't need to compute Khat (i.e. the boundary condition) nor the if (D[k] > Khat[k]) condition. May 15 at 7:33