American Call Pricing with Multiple Discrete Dividends
Article Quant Q&A · Author: Sarah Van Distel
Summary
The document describes a Clenshaw–Curtis quadrature implementation for pricing an American call on a stock that pays multiple discrete dividends. The author compares the computed price with a published result and with a Vellekoop–Nieuwenhuis implementation, which agrees with the paper. Reproducing the same discrepancy in Octave suggests a logic error rather than a language-specific numerical issue.
The accepted answer locates the error in the backward recursion: the original loop applies exercise and continuation operations at the wrong points relative to dividend dates. The corrected sequence checks early exercise at a dividend date, then integrates the continuation value over the next interval, and finally applies the exercise test at the first dividend date before integrating back to time zero. The example illustrates how event timing matters in numerical option pricing. It provides a correction for the presented implementation, but does not develop the quadrature method or establish how the correction behaves across other parameter choices.
Key ideas
- The backward recursion must align continuation calculations with the intervals between dividend dates.
- Early exercise should be checked at the appropriate dividend dates before propagating values backward.
- A matching European option result does not guarantee that American exercise logic is correctly implemented.
- Agreement with another method can help diagnose errors in an option-pricing implementation.
Tags
Full text
# 79380
# 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.
## Answer by David Walthour (score 1, accepted)
https://quant.stackexchange.com/a/79573
I fixed the last part of the code as follows and now it gives the correct answers:
```
for (int k = tD.size() - 1; k >= 1; --k) {
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.array()).cwiseMax(V_Ac_exercise);
}
xiLUplus = xiLU - D[k-1];
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);
}
V_Ac = exp(-r * tau[k]) * 0.5 * (U - L) * fnxkLU; // This is the continuation value
}
if (D[0] > Khat[0]) {
Eigen::ArrayXd V_Ac_exercise = (xiLUplus + D[0] - K).cwiseMax(0.0); // This is the exercise value
V_Ac = (V_Ac.array()).cwiseMax(V_Ac_exercise);
}
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;
`
```Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.