3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ This code is for American options or European? What is the difference between the two codes? $\endgroup$
    – nbbo2
    May 14 at 18:50
  • 2
    $\begingroup$ @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. $\endgroup$ May 15 at 7:33

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.