// OptionPricing.hpp — Black-Scholes, Binomial, Monte Carlo option pricing
// SPDX-License-Identifier: LGPL-3.0-or-later
// Copyright (c) 2026 calx contributors

#ifndef CALX_FINANCE_OPTIONPRICING_HPP
#define CALX_FINANCE_OPTIONPRICING_HPP

#include <cmath>
#include <cstddef>
#include <vector>
#include <algorithm>
#include <random>
#include <stdexcept>

namespace calx::finance {

// ---------------------------------------------------------------
// 正規分布 CDF: Phi(x) = 0.5 * erfc(-x / sqrt(2))
// ---------------------------------------------------------------
namespace detail {
    template<typename T>
    T normalCDF(T x) {
        return T{0.5} * std::erfc(-x / std::sqrt(T{2}));
    }

    template<typename T>
    T normalPDF(T x) {
        static const T INV_SQRT_2PI = T{1} / std::sqrt(T{2} * T{3.14159265358979323846});
        return INV_SQRT_2PI * std::exp(T{-0.5} * x * x);
    }
} // namespace detail

// ---------------------------------------------------------------
// Black-Scholes 結果
// ---------------------------------------------------------------
template<typename T>
struct BSResult {
    T callPrice, putPrice;
    T delta, gamma, theta, vega, rho;
};

// ---------------------------------------------------------------
// Black-Scholes 公式 + Greeks
//   S:       原資産価格
//   K:       行使価格
//   r:       無リスク金利 (年率)
//   sigma:   ボラティリティ (年率)
//   T_years: 満期までの年数
// ---------------------------------------------------------------
template<typename T>
BSResult<T> blackScholes(T S, T K, T r, T sigma, T T_years) {
    if (S <= T{0} || K <= T{0} || sigma <= T{0} || T_years <= T{0})
        throw std::invalid_argument("blackScholes: all parameters must be positive");

    T sqrtT = std::sqrt(T_years);
    T d1 = (std::log(S / K) + (r + T{0.5} * sigma * sigma) * T_years) / (sigma * sqrtT);
    T d2 = d1 - sigma * sqrtT;

    T Nd1 = detail::normalCDF(d1);
    T Nd2 = detail::normalCDF(d2);
    T Nmd1 = detail::normalCDF(-d1);
    T Nmd2 = detail::normalCDF(-d2);
    T discountK = K * std::exp(-r * T_years);

    BSResult<T> res;
    res.callPrice = S * Nd1 - discountK * Nd2;
    res.putPrice  = discountK * Nmd2 - S * Nmd1;

    // Greeks (コールオプション基準)
    T nd1 = detail::normalPDF(d1);
    res.delta = Nd1;
    res.gamma = nd1 / (S * sigma * sqrtT);
    res.theta = -(S * nd1 * sigma) / (T{2} * sqrtT) - r * discountK * Nd2;
    res.vega  = S * nd1 * sqrtT;
    res.rho   = discountK * T_years * Nd2;

    return res;
}

// ---------------------------------------------------------------
// 二項モデル (Cox-Ross-Rubinstein)
// ---------------------------------------------------------------
template<typename T>
T binomialOption(T S, T K, T r, T sigma, T T_years, int steps, bool isCall) {
    if (steps <= 0)
        throw std::invalid_argument("binomialOption: steps must be > 0");

    T dt = T_years / static_cast<T>(steps);
    T u = std::exp(sigma * std::sqrt(dt));
    T d = T{1} / u;
    T disc = std::exp(-r * dt);
    T p = (std::exp(r * dt) - d) / (u - d); // リスク中立確率
    T q = T{1} - p;

    // 満期時点のオプション価値
    std::vector<T> prices(static_cast<size_t>(steps) + 1);
    for (int i = 0; i <= steps; ++i) {
        T ST = S * std::pow(u, static_cast<T>(steps - 2 * i));
        if (isCall)
            prices[static_cast<size_t>(i)] = std::max(ST - K, T{0});
        else
            prices[static_cast<size_t>(i)] = std::max(K - ST, T{0});
    }

    // 逆向き帰納法
    for (int step = steps - 1; step >= 0; --step) {
        for (int i = 0; i <= step; ++i) {
            prices[static_cast<size_t>(i)] =
                disc * (p * prices[static_cast<size_t>(i)] + q * prices[static_cast<size_t>(i) + 1]);
        }
    }

    return prices[0];
}

// ---------------------------------------------------------------
// モンテカルロ法
// 幾何ブラウン運動: S_T = S * exp((r - 0.5*sigma^2)*T + sigma*sqrt(T)*Z)
// ---------------------------------------------------------------
template<typename T>
T monteCarloOption(T S, T K, T r, T sigma, T T_years,
                   size_t nPaths, bool isCall, unsigned seed = 42) {
    if (nPaths == 0)
        throw std::invalid_argument("monteCarloOption: nPaths must be > 0");

    std::mt19937 rng(seed);
    std::normal_distribution<T> dist(T{0}, T{1});

    T drift = (r - T{0.5} * sigma * sigma) * T_years;
    T vol   = sigma * std::sqrt(T_years);
    T disc  = std::exp(-r * T_years);

    T sumPayoff = T{0};
    for (size_t i = 0; i < nPaths; ++i) {
        T z = dist(rng);
        T ST = S * std::exp(drift + vol * z);
        T payoff = isCall ? std::max(ST - K, T{0}) : std::max(K - ST, T{0});
        sumPayoff += payoff;
    }

    return disc * sumPayoff / static_cast<T>(nPaths);
}

} // namespace calx::finance

#endif // CALX_FINANCE_OPTIONPRICING_HPP
