EuropeanOption.cpp

// EuropeanOption.cpp
//
// Implementation of European option pricing using Black-Scholes model
// with Boost libraries for normal distribution functions.
// 
// 
// // Created by: Andre Boiko, 2025-09-12
//

#include "EuropeanOption.h"
#include <cmath>
#include <boost/math/distributions/normal.hpp>
using namespace boost::math;

//////////// Gaussian functions using Boost /////////////////////////////////
double EuropeanOption::n(double x) const
{
    normal_distribution<> myNormal(0.0, 1.0);
    return pdf(myNormal, x);
}

double EuropeanOption::N(double x) const
{
    normal_distribution<> myNormal(0.0, 1.0);
    return cdf(myNormal, x);
}

/////////////////////////////////////////////////////////////////////////////////////

void EuropeanOption::init()
{
    T = 1.0;
    K = 100.0;
    sig = 0.2;
    r = 0.05;
    b = 0.05;  // Default: stock option (b = r)
    S = 100.0;
}

void EuropeanOption::copy(const EuropeanOption& o2)
{
    T = o2.T;
    K = o2.K;
    sig = o2.sig;
    r = o2.r;
    b = o2.b;
    S = o2.S;
}

EuropeanOption::EuropeanOption() { init(); }
EuropeanOption::EuropeanOption(const EuropeanOption& o2) { copy(o2); }
EuropeanOption::EuropeanOption(double T, double K, double sig, double r, double b, double S)
    : T(T), K(K), sig(sig), r(r), b(b), S(S) {
}
EuropeanOption::~EuropeanOption() {}

EuropeanOption& EuropeanOption::operator=(const EuropeanOption& option2)
{
    if (this == &option2) return *this;
    copy(option2);
    return *this;
}

void EuropeanOption::SetParameters(double T, double K, double sig, double r, double b, double S)
{
    this->T = T; this->K = K; this->sig = sig;
    this->r = r; this->b = b; this->S = S;
}

// CORE PRICING FUNCTIONS 

double EuropeanOption::CallPrice() const
{
    // Cache shared calculations for performance
    const double sqrt_T = sqrt(T);
    const double sig_sqrt_T = sig * sqrt_T;
    const double d1 = (log(S / K) + (b + 0.5 * sig * sig) * T) / sig_sqrt_T;
    const double d2 = d1 - sig_sqrt_T;

    return S * exp((b - r) * T) * N(d1) - K * exp(-r * T) * N(d2);
}

double EuropeanOption::PutPrice() const
{
    const double sqrt_T = sqrt(T);
    const double sig_sqrt_T = sig * sqrt_T;
    const double d1 = (log(S / K) + (b + 0.5 * sig * sig) * T) / sig_sqrt_T;
    const double d2 = d1 - sig_sqrt_T;

    return K * exp(-r * T) * N(-d2) - S * exp((b - r) * T) * N(-d1);
}

double EuropeanOption::CallPrice(double S_new) const
{
    const double sqrt_T = sqrt(T);
    const double sig_sqrt_T = sig * sqrt_T;
    const double d1 = (log(S_new / K) + (b + 0.5 * sig * sig) * T) / sig_sqrt_T;
    const double d2 = d1 - sig_sqrt_T;

    return S_new * exp((b - r) * T) * N(d1) - K * exp(-r * T) * N(d2);
}

double EuropeanOption::PutPrice(double S_new) const
{
    const double sqrt_T = sqrt(T);
    const double sig_sqrt_T = sig * sqrt_T;
    const double d1 = (log(S_new / K) + (b + 0.5 * sig * sig) * T) / sig_sqrt_T;
    const double d2 = d1 - sig_sqrt_T;

    return K * exp(-r * T) * N(-d2) - S_new * exp((b - r) * T) * N(-d1);
}

// VECTOR PRICING FUNCTIONS 

std::vector<double> EuropeanOption::CallPriceVector(const std::vector<double>& S_array) const
{
    std::vector<double> prices;
    prices.reserve(S_array.size());

    for (const auto& S_val : S_array) {
        prices.push_back(CallPrice(S_val));
    }
    return prices;
}

std::vector<double> EuropeanOption::PutPriceVector(const std::vector<double>& S_array) const
{
    std::vector<double> prices;
    prices.reserve(S_array.size());

    for (const auto& S_val : S_array) {
        prices.push_back(PutPrice(S_val));
    }
    return prices;
}

std::vector<double> EuropeanOption::CallPriceByT(const std::vector<double>& T_array) const
{
    std::vector<double> prices;
    prices.reserve(T_array.size());
    const double original_T = T;

    for (const auto& T_val : T_array) {
        T = T_val;  // T is mutable for this purpose
        prices.push_back(CallPrice());
    }

    T = original_T;  // Restore original value
    return prices;
}

std::vector<double> EuropeanOption::PutPriceByT(const std::vector<double>& T_array) const
{
    std::vector<double> prices;
    prices.reserve(T_array.size());
    const double original_T = T;

    for (const auto& T_val : T_array) {
        T = T_val;
        prices.push_back(PutPrice());
    }

    T = original_T;
    return prices;
}

std::vector<double> EuropeanOption::CallPriceBySig(const std::vector<double>& sig_array) const
{
    std::vector<double> prices;
    prices.reserve(sig_array.size());
    const double original_sig = sig;

    for (const auto& sig_val : sig_array) {
        sig = sig_val;  // sig is mutable for this purpose
        prices.push_back(CallPrice());
    }

    sig = original_sig;  // Restore original value
    return prices;
}

std::vector<double> EuropeanOption::PutPriceBySig(const std::vector<double>& sig_array) const
{
    std::vector<double> prices;
    prices.reserve(sig_array.size());
    const double original_sig = sig;

    for (const auto& sig_val : sig_array) {
        sig = sig_val;
        prices.push_back(PutPrice());
    }

    sig = original_sig;
    return prices;
}

// PUT-CALL PARITY FUNCTIONS 

double EuropeanOption::PutFromCallParity() const
{
    const double callPrice = CallPrice();
    return callPrice + K * exp(-r * T) - S;
}

bool EuropeanOption::CheckParity(double tolerance) const
{
    const double callPrice = CallPrice();
    const double putPrice = PutPrice();
    const double putFromParity = PutFromCallParity();

    return std::abs(putPrice - putFromParity) < tolerance;
}

// GREEKS - ANALYTICAL FORMULAS 

double EuropeanOption::CallDelta() const
{
    const double sqrt_T = sqrt(T);
    const double d1 = (log(S / K) + (b + 0.5 * sig * sig) * T) / (sig * sqrt_T);
    return exp((b - r) * T) * N(d1);
}

double EuropeanOption::PutDelta() const
{
    const double sqrt_T = sqrt(T);
    const double d1 = (log(S / K) + (b + 0.5 * sig * sig) * T) / (sig * sqrt_T);
    return exp((b - r) * T) * (N(d1) - 1.0);
}

double EuropeanOption::CallDelta(double S_new) const
{
    const double sqrt_T = sqrt(T);
    const double d1 = (log(S_new / K) + (b + 0.5 * sig * sig) * T) / (sig * sqrt_T);
    return exp((b - r) * T) * N(d1);
}

double EuropeanOption::PutDelta(double S_new) const
{
    const double sqrt_T = sqrt(T);
    const double d1 = (log(S_new / K) + (b + 0.5 * sig * sig) * T) / (sig * sqrt_T);
    return exp((b - r) * T) * (N(d1) - 1.0);
}

double EuropeanOption::Gamma() const
{
    const double sqrt_T = sqrt(T);
    const double sig_sqrt_T = sig * sqrt_T;
    const double d1 = (log(S / K) + (b + 0.5 * sig * sig) * T) / sig_sqrt_T;
    return (n(d1) * exp((b - r) * T)) / (S * sig_sqrt_T);
}

double EuropeanOption::Gamma(double S_new) const
{
    const double sqrt_T = sqrt(T);
    const double sig_sqrt_T = sig * sqrt_T;
    const double d1 = (log(S_new / K) + (b + 0.5 * sig * sig) * T) / sig_sqrt_T;
    return (n(d1) * exp((b - r) * T)) / (S_new * sig_sqrt_T);
}

// VECTOR GREEKS 

std::vector<double> EuropeanOption::CallDeltaVector(const std::vector<double>& S_array) const
{
    std::vector<double> deltas;
    deltas.reserve(S_array.size());

    for (const auto& S_val : S_array) {
        deltas.push_back(CallDelta(S_val));
    }
    return deltas;
}

std::vector<double> EuropeanOption::PutDeltaVector(const std::vector<double>& S_array) const
{
    std::vector<double> deltas;
    deltas.reserve(S_array.size());

    for (const auto& S_val : S_array) {
        deltas.push_back(PutDelta(S_val));
    }
    return deltas;
}

std::vector<double> EuropeanOption::GammaVector(const std::vector<double>& S_array) const
{
    std::vector<double> gammas;
    gammas.reserve(S_array.size());

    for (const auto& S_val : S_array) {
        gammas.push_back(Gamma(S_val));
    }
    return gammas;
}

// NUMERICAL GREEKS 

double EuropeanOption::CallDeltaNumerical(double h) const
{
    const double upPrice = CallPrice(S + h);
    const double downPrice = CallPrice(S - h);
    return (upPrice - downPrice) / (2.0 * h);
}

double EuropeanOption::PutDeltaNumerical(double h) const
{
    const double upPrice = PutPrice(S + h);
    const double downPrice = PutPrice(S - h);
    return (upPrice - downPrice) / (2.0 * h);
}

double EuropeanOption::CallDeltaNumerical(double S_new, double h) const
{
    const double upPrice = CallPrice(S_new + h);
    const double downPrice = CallPrice(S_new - h);
    return (upPrice - downPrice) / (2.0 * h);
}

double EuropeanOption::PutDeltaNumerical(double S_new, double h) const
{
    const double upPrice = PutPrice(S_new + h);
    const double downPrice = PutPrice(S_new - h);
    return (upPrice - downPrice) / (2.0 * h);
}

double EuropeanOption::GammaNumerical(double h) const
{
    const double upPrice = CallPrice(S + h);
    const double centerPrice = CallPrice(S);
    const double downPrice = CallPrice(S - h);
    return (upPrice - 2.0 * centerPrice + downPrice) / (h * h);
}

double EuropeanOption::GammaNumerical(double S_new, double h) const
{
    const double upPrice = CallPrice(S_new + h);
    const double centerPrice = CallPrice(S_new);
    const double downPrice = CallPrice(S_new - h);
    return (upPrice - 2.0 * centerPrice + downPrice) / (h * h);
}

// VECTOR NUMERICAL GREEKS 

std::vector<double> EuropeanOption::CallDeltaNumericalVector(const std::vector<double>& S_array, double h) const
{
    std::vector<double> deltas;
    deltas.reserve(S_array.size());

    for (const auto& S_val : S_array) {
        deltas.push_back(CallDeltaNumerical(S_val, h));
    }
    return deltas;
}

std::vector<double> EuropeanOption::GammaNumericalVector(const std::vector<double>& S_array, double h) const
{
    std::vector<double> gammas;
    gammas.reserve(S_array.size());

    for (const auto& S_val : S_array) {
        gammas.push_back(GammaNumerical(S_val, h));
    }
    return gammas;
}

// ADDITIONAL GREEKS 

double EuropeanOption::CallTheta() const
{
    const double sqrt_T = sqrt(T);
    const double sig_sqrt_T = sig * sqrt_T;
    const double d1 = (log(S / K) + (b + 0.5 * sig * sig) * T) / sig_sqrt_T;
    const double d2 = d1 - sig_sqrt_T;

    const double t1 = (S * exp((b - r) * T) * n(d1) * sig) / (2.0 * sqrt_T);
    const double t2 = (b - r) * S * exp((b - r) * T) * N(d1);
    const double t3 = r * K * exp(-r * T) * N(d2);

    return -(t1 + t2 + t3);
}

double EuropeanOption::PutTheta() const
{
    const double sqrt_T = sqrt(T);
    const double sig_sqrt_T = sig * sqrt_T;
    const double d1 = (log(S / K) + (b + 0.5 * sig * sig) * T) / sig_sqrt_T;
    const double d2 = d1 - sig_sqrt_T;

    const double t1 = (S * exp((b - r) * T) * n(d1) * sig) / (2.0 * sqrt_T);
    const double t2 = (b - r) * S * exp((b - r) * T) * N(-d1);
    const double t3 = r * K * exp(-r * T) * N(-d2);

    return t2 + t3 - t1;
}

double EuropeanOption::CallVega() const
{
    const double sqrt_T = sqrt(T);
    const double d1 = (log(S / K) + (b + 0.5 * sig * sig) * T) / (sig * sqrt_T);
    return S * exp((b - r) * T) * n(d1) * sqrt_T;
}

double EuropeanOption::PutVega() const
{
    return CallVega();  // Same for calls and puts
}

double EuropeanOption::CallRho() const
{
    const double sqrt_T = sqrt(T);
    const double d1 = (log(S / K) + (b + 0.5 * sig * sig) * T) / (sig * sqrt_T);
    const double d2 = d1 - sig * sqrt_T;

    if (b != 0.0)
        return T * K * exp(-r * T) * N(d2);
    else
        return -T * CallPrice();
}

double EuropeanOption::PutRho() const
{
    const double sqrt_T = sqrt(T);
    const double d1 = (log(S / K) + (b + 0.5 * sig * sig) * T) / (sig * sqrt_T);
    const double d2 = d1 - sig * sqrt_T;

    if (b != 0.0)
        return -T * K * exp(-r * T) * N(-d2);
    else
        return -T * PutPrice();
}

// MATRIX FUNCTIONALITY (CLASS STATICS - WRAPPERS) 

std::vector<std::vector<double>> EuropeanOption::CallPriceMatrix(
    const std::vector<std::vector<double>>& S_matrix,
    const std::vector<std::vector<EuropeanOption>>& option_matrix)
{
    std::vector<std::vector<double>> result;
    result.reserve(S_matrix.size());

    for (size_t i = 0; i < S_matrix.size(); ++i) {
        std::vector<double> row;
        row.reserve(S_matrix[i].size());

        for (size_t j = 0; j < S_matrix[i].size(); ++j) {
            row.push_back(option_matrix[i][j].CallPrice(S_matrix[i][j]));
        }
        result.push_back(std::move(row));
    }
    return result;
}

std::vector<std::vector<double>> EuropeanOption::DeltaMatrix(
    const std::vector<std::vector<double>>& S_matrix,
    const std::vector<std::vector<EuropeanOption>>& option_matrix)
{
    std::vector<std::vector<double>> result;
    result.reserve(S_matrix.size());

    for (size_t i = 0; i < S_matrix.size(); ++i) {
        std::vector<double> row;
        row.reserve(S_matrix[i].size());

        for (size_t j = 0; j < S_matrix[i].size(); ++j) {
            row.push_back(option_matrix[i][j].CallDelta(S_matrix[i][j]));
        }
        result.push_back(std::move(row));
    }
    return result;
}

std::vector<std::vector<double>> EuropeanOption::GammaMatrix(
    const std::vector<std::vector<double>>& S_matrix,
    const std::vector<std::vector<EuropeanOption>>& option_matrix)
{
    std::vector<std::vector<double>> result;
    result.reserve(S_matrix.size());

    for (size_t i = 0; i < S_matrix.size(); ++i) {
        std::vector<double> row;
        row.reserve(S_matrix[i].size());

        for (size_t j = 0; j < S_matrix[i].size(); ++j) {
            row.push_back(option_matrix[i][j].Gamma(S_matrix[i][j]));
        }
        result.push_back(std::move(row));
    }
    return result;
}