/*
 Niels Elken Sønderby 2003
 NesQuant - http://www.nielses.dk/quantlib/nesquant
 NesQuant is an extension to QuantLib - http://quantlib.org/

 Modification of EuropeanOption example
*/

#include <nq/nesquant.hpp>

using namespace QuantLib;
using namespace NesQuant;

using QuantLib::Pricers::EuropeanOption;
using QuantLib::Math::SegmentIntegral;

using NesQuant::Math::KronrodIntegral;

using QuantLib::PlainVanillaPayoff;
using QuantLib::Option;

// This will be included in the library after a bit of redesign
class Payoff2 {
    public:
        Payoff2(Time maturity,
               double strike,
               double s0,
               double sigma,
               Rate r)
        : maturity_(maturity),
        strike_(strike),
        s0_(s0),
        sigma_(sigma),r_(r){}

        double operator()(double x) const {
           PlainVanillaPayoff payoff(Option::Call, strike_);

           double nuT = (r_-0.5*sigma_*sigma_)*maturity_;
           return QL_EXP(-r_*maturity_)
               *payoff(s0_*QL_EXP(x))
               *QL_EXP(-(x - nuT)*(x -nuT)/(2*sigma_*sigma_*maturity_))
               /QL_SQRT(2.0*3.141592*sigma_*sigma_*maturity_);
        }
private:
    Time maturity_;
    double strike_;
    double s0_;
    double sigma_;
    Rate r_;
};

int main(int argc, char* argv[])
{
    try {
        // our option
        double underlying = 102;
        double strike = 100;      // at the money
        Spread dividendYield = 0.0; // no dividends
        Rate riskFreeRate = 0.05; // 5%
        Time maturity = 0.25;      // 3 months
        double volatility = 0.20; // 20%
        std::cout << "Time to maturity = "        << maturity
                  << std::endl;
        std::cout << "Underlying price = "        << underlying
                  << std::endl;
        std::cout << "Strike = "                  << strike
                  << std::endl;
        std::cout << "Risk-free interest rate = " << riskFreeRate
                  << std::endl;
        std::cout << "Volatility = "              << volatility
                  << std::endl;
        std::cout << std::endl;

        // write column headings
        std::cout << "Method\t\tValue\tEstimatedError\tDiscrepancy"
            "\tRel. Discr." << std::endl;

        // Black Scholes analytic solution
        std::string method ="Black Scholes";
        double value = EuropeanOption(Option::Call, underlying, strike,
            dividendYield, riskFreeRate, maturity, volatility).value();
        double estimatedError = 0.0;
        double discrepancy = 0.0;
        double relativeDiscrepancy = 0.0;
        std::cout << method << "\t"
             << DoubleFormatter::toString(value, 4) << "\t"
             << DoubleFormatter::toString(estimatedError, 4) << "\t\t"
             << DoubleFormatter::toString(discrepancy, 6) << "\t"
             << DoubleFormatter::toString(relativeDiscrepancy, 6)
             << std::endl;

        // store the Black Scholes value as the correct one
        double rightValue = value;


        // Segment Integral
        method ="SegmentIntegral";
        Payoff2 po(maturity, strike, underlying, volatility, riskFreeRate);

        double nuT = (riskFreeRate - 0.5*volatility*volatility)*maturity;
        double infinity = 10.0*volatility*QL_SQRT(maturity);

		int intervals = 1000;

        SegmentIntegral sintegrator(intervals);

        value = sintegrator(po, nuT-infinity, nuT+infinity);
        discrepancy = QL_FABS(value-rightValue);
        relativeDiscrepancy = discrepancy/rightValue;
        std::cout << method << "\t"
             << DoubleFormatter::toString(value, 4) << "\t"
             << "N/A\t\t"
             << DoubleFormatter::toString(discrepancy, 6) << "\t"
             << DoubleFormatter::toString(relativeDiscrepancy, 6)
             << std::endl;
		std::cout << "Function evaluations : " << intervals + 1 << std::endl;

        // Gauss-Kronrod Integral
        method ="KronrodIntegral";

        KronrodIntegral kintegrator(0.0001);

        value = kintegrator(po, nuT-infinity, nuT+infinity);
        discrepancy = QL_FABS(value-rightValue);
        relativeDiscrepancy = discrepancy/rightValue;
        std::cout << method << "\t"
             << DoubleFormatter::toString(value, 4) << "\t"
             << "N/A\t\t"
             << DoubleFormatter::toString(discrepancy, 6) << "\t"
             << DoubleFormatter::toString(relativeDiscrepancy, 6)
             << std::endl;
		std::cout << "Function evaluations : " 
			      << kintegrator.functionEvaluations() << std::endl;

        return 0;
    } catch (std::exception& e) {
        std::cout << e.what() << std::endl;
        return 1;
    } catch (...) {
        std::cout << "unknown error" << std::endl;
        return 1;
    }
}