/src/quantlib/ql/pricingengines/blackcalculator.cpp

Source
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2003, 2004, 2005, 2006 Ferdinando Ametrano
 Copyright (C) 2006 StatPro Italia srl

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <https://www.quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include <ql/pricingengines/blackcalculator.hpp>
#include <ql/math/distributions/normaldistribution.hpp>
#include <ql/math/comparison.hpp>

namespace QuantLib {

    class BlackCalculator::Calculator : public AcyclicVisitor,
                                        public Visitor<Payoff>,
                                        public Visitor<PlainVanillaPayoff>,
                                        public Visitor<CashOrNothingPayoff>,
                                        public Visitor<AssetOrNothingPayoff>,
                                        public Visitor<GapPayoff> {
      private:
        BlackCalculator& black_;
      public:
        explicit Calculator(BlackCalculator& black) : black_(black) {}
        void visit(Payoff&) override;
        void visit(PlainVanillaPayoff&) override;
        void visit(CashOrNothingPayoff&) override;
        void visit(AssetOrNothingPayoff&) override;
        void visit(GapPayoff&) override;
    };


    BlackCalculator::BlackCalculator(const ext::shared_ptr<StrikedTypePayoff>& p,
                                     Real forward,
                                     Real stdDev,
                                     Real discount)
    : strike_(p->strike()), forward_(forward), stdDev_(stdDev),
      discount_(discount), variance_(stdDev*stdDev) {
        initialize(p);
    }

    BlackCalculator::BlackCalculator(Option::Type optionType,
                                     Real strike,
                                     Real forward,
                                     Real stdDev,
                                     Real discount)
    : strike_(strike), forward_(forward), stdDev_(stdDev),
      discount_(discount), variance_(stdDev*stdDev) {
        initialize(ext::shared_ptr<StrikedTypePayoff>(new
            PlainVanillaPayoff(optionType, strike)));
    }

    void BlackCalculator::initialize(const ext::shared_ptr<StrikedTypePayoff>& p) {
        QL_REQUIRE(strike_>=0.0,
                   "strike (" << strike_ << ") must be non-negative");
        QL_REQUIRE(forward_>0.0,
                   "forward (" << forward_ << ") must be positive");
        //QL_REQUIRE(displacement_>=0.0,
        //           "displacement (" << displacement_ << ") must be non-negative");
        QL_REQUIRE(stdDev_>=0.0,
                   "stdDev (" << stdDev_ << ") must be non-negative");
        QL_REQUIRE(discount_>0.0,
                   "discount (" << discount_ << ") must be positive");

        if (stdDev_>=QL_EPSILON) {
            if (close(strike_, 0.0)) {
                d1_ = QL_MAX_REAL;
                d2_ = QL_MAX_REAL;
                cum_d1_ = 1.0;
                cum_d2_ = 1.0;
                n_d1_ = 0.0;
                n_d2_ = 0.0;
            } else {
                d1_ = std::log(forward_/strike_)/stdDev_ + 0.5*stdDev_;
                d2_ = d1_-stdDev_;
                CumulativeNormalDistribution f;
                cum_d1_ = f(d1_);
                cum_d2_ = f(d2_);
                n_d1_ = f.derivative(d1_);
                n_d2_ = f.derivative(d2_);
            }
        } else {
            if (close(forward_, strike_)) {
                d1_ = 0;
                d2_ = 0;
                cum_d1_ = 0.5;
                cum_d2_ = 0.5;
                n_d1_ = M_SQRT_2 * M_1_SQRTPI;
                n_d2_ = M_SQRT_2 * M_1_SQRTPI;
            } else if (forward_>strike_) {
                d1_ = QL_MAX_REAL;
                d2_ = QL_MAX_REAL;
                cum_d1_ = 1.0;
                cum_d2_ = 1.0;
                n_d1_ = 0.0;
                n_d2_ = 0.0;
            } else {
                d1_ = QL_MIN_REAL;
                d2_ = QL_MIN_REAL;
                cum_d1_ = 0.0;
                cum_d2_ = 0.0;
                n_d1_ = 0.0;
                n_d2_ = 0.0;
            }
        }

        x_ = strike_;
        DxDstrike_ = 1.0;

        // the following one will probably disappear as soon as
        // super-share will be properly handled
        DxDs_ = 0.0;

        // this part is always executed.
        // in case of plain-vanilla payoffs, it is also the only part
        // which is executed.
        switch (p->optionType()) {
          case Option::Call:
            alpha_     =  cum_d1_;//  N(d1)
            DalphaDd1_ =    n_d1_;//  n(d1)
            beta_      = -cum_d2_;// -N(d2)
            DbetaDd2_  = -  n_d2_;// -n(d2)
            break;
          case Option::Put:
            alpha_     = -1.0+cum_d1_;// -N(-d1)
            DalphaDd1_ =        n_d1_;//  n( d1)
            beta_      =  1.0-cum_d2_;//  N(-d2)
            DbetaDd2_  =     -  n_d2_;// -n( d2)
            break;
          default:
            QL_FAIL("invalid option type");
        }

        // now dispatch on type.

        Calculator calc(*this);
        p->accept(calc);
    }

    void BlackCalculator::Calculator::visit(Payoff& p) {
        QL_FAIL("unsupported payoff type: " << p.name());
    }

    void BlackCalculator::Calculator::visit(PlainVanillaPayoff&) {}

    void BlackCalculator::Calculator::visit(CashOrNothingPayoff& payoff) {
        black_.alpha_ = black_.DalphaDd1_ = 0.0;
        black_.x_ = payoff.cashPayoff();
        black_.DxDstrike_ = 0.0;
        switch (payoff.optionType()) {
          case Option::Call:
            black_.beta_     = black_.cum_d2_;
            black_.DbetaDd2_ = black_.n_d2_;
            break;
          case Option::Put:
            black_.beta_     = 1.0-black_.cum_d2_;
            black_.DbetaDd2_ =    -black_.n_d2_;
            break;
          default:
            QL_FAIL("invalid option type");
        }
    }

    void BlackCalculator::Calculator::visit(AssetOrNothingPayoff& payoff) {
        black_.beta_ = black_.DbetaDd2_ = 0.0;
        switch (payoff.optionType()) {
          case Option::Call:
            black_.alpha_     = black_.cum_d1_;
            black_.DalphaDd1_ = black_.n_d1_;
            break;
          case Option::Put:
            black_.alpha_     = 1.0-black_.cum_d1_;
            black_.DalphaDd1_ = -black_.n_d1_;
            break;
          default:
            QL_FAIL("invalid option type");
        }
    }

    void BlackCalculator::Calculator::visit(GapPayoff& payoff) {
        black_.x_ = payoff.secondStrike();
        black_.DxDstrike_ = 0.0;
    }

    Real BlackCalculator::value() const {
        Real result = discount_ * (forward_ * alpha_ + x_ * beta_);
        return result;
    }

    Real BlackCalculator::delta(Real spot) const {

        QL_REQUIRE(spot > 0.0, "positive spot value required: " <<
                   spot << " not allowed");

        // Handle zero volatility case
        if (stdDev_ <= QL_EPSILON) {
            // For zero volatility, delta is:
            // ITM Call: 1.0, OTM Call: 0.0, ATM Call: 0.5
            // ITM Put: -1.0, OTM Put: 0.0, ATM Put: -0.5
            Real DforwardDs = forward_ / spot;
            if (close(forward_, strike_)) {
                // ATM case
                if (alpha_ >= 0) { // Call
                    return discount_ * 0.5 * DforwardDs;
                } else { // Put
                    return discount_ * (-0.5) * DforwardDs;
                }
            } else if (forward_ > strike_) {
                // ITM Call, OTM Put
                if (alpha_ >= 0) { // Call
                    return discount_ * 1.0 * DforwardDs;
                } else { // Put  
                    return 0.0;
                }
            } else {
                // OTM Call, ITM Put
                if (alpha_ >= 0) { // Call
                    return 0.0;
                } else { // Put
                    return discount_ * (-1.0) * DforwardDs;
                }
            }
        }

        Real DforwardDs = forward_ / spot;

        Real temp = stdDev_*spot;
        Real DalphaDs = DalphaDd1_/temp;
        Real DbetaDs  = DbetaDd2_/temp;
        Real temp2 = DalphaDs * forward_ + alpha_ * DforwardDs
                      +DbetaDs  * x_       + beta_  * DxDs_;

        return discount_ * temp2;
    }

    Real BlackCalculator::deltaForward() const {

        // Handle zero volatility case
        if (stdDev_ <= QL_EPSILON) {
            // For zero volatility, forward delta is:
            // ITM Call: 1.0, OTM Call: 0.0, ATM Call: 0.5
            // ITM Put: -1.0, OTM Put: 0.0, ATM Put: -0.5
            if (close(forward_, strike_)) {
                // ATM case
                if (alpha_ >= 0) { // Call
                    return discount_ * 0.5;
                } else { // Put
                    return discount_ * (-0.5);
                }
            } else if (forward_ > strike_) {
                // ITM Call, OTM Put
                if (alpha_ >= 0) { // Call
                    return discount_ * 1.0;
                } else { // Put
                    return 0.0;
                }
            } else {
                // OTM Call, ITM Put
                if (alpha_ >= 0) { // Call
                    return 0.0;
                } else { // Put
                    return discount_ * (-1.0);
                }
            }
        }

        Real temp = stdDev_*forward_;
        Real DalphaDforward = DalphaDd1_/temp;
        Real DbetaDforward  = DbetaDd2_/temp;
        Real temp2 = DalphaDforward * forward_ + alpha_
                      +DbetaDforward  * x_; // DXDforward = 0.0

        return discount_ * temp2;
    }

    Real BlackCalculator::elasticity(Real spot) const {
        Real val = value();
        Real del = delta(spot);
        if (val>QL_EPSILON)
            return del/val*spot;
        else if (std::fabs(del)<QL_EPSILON)
            return 0.0;
        else if (del>0.0)
            return QL_MAX_REAL;
        else
            return QL_MIN_REAL;
    }

    Real BlackCalculator::elasticityForward() const {
        Real val = value();
        Real del = deltaForward();
        if (val>QL_EPSILON)
            return del/val*forward_;
        else if (std::fabs(del)<QL_EPSILON)
            return 0.0;
        else if (del>0.0)
            return QL_MAX_REAL;
        else
            return QL_MIN_REAL;
    }

    Real BlackCalculator::gamma(Real spot) const {

        QL_REQUIRE(spot > 0.0, "positive spot value required: " <<
                   spot << " not allowed");

        // Handle zero volatility case
        if (stdDev_ <= QL_EPSILON) {
            // For zero volatility, gamma is always 0 (no convexity)
            return 0.0;
        }

        Real DforwardDs = forward_ / spot;

        Real temp = stdDev_*spot;
        Real DalphaDs = DalphaDd1_/temp;
        Real DbetaDs  = DbetaDd2_/temp;

        Real D2alphaDs2 = - DalphaDs/spot*(1+d1_/stdDev_);
        Real D2betaDs2  = - DbetaDs /spot*(1+d2_/stdDev_);

        Real temp2 = D2alphaDs2 * forward_ + 2.0 * DalphaDs * DforwardDs
                      +D2betaDs2  * x_       + 2.0 * DbetaDs  * DxDs_;

        return  discount_ * temp2;
    }

    Real BlackCalculator::gammaForward() const {

        // Handle zero volatility case
        if (stdDev_ <= QL_EPSILON) {
            // For zero volatility, gamma is always 0 (no convexity)
            return 0.0;
        }

        Real temp = stdDev_*forward_;
        Real DalphaDforward = DalphaDd1_/temp;
        Real DbetaDforward  = DbetaDd2_/temp;

        Real D2alphaDforward2 = - DalphaDforward/forward_*(1+d1_/stdDev_);
        Real D2betaDforward2  = - DbetaDforward /forward_*(1+d2_/stdDev_);

        Real temp2 = D2alphaDforward2 * forward_ + 2.0 * DalphaDforward
                      +D2betaDforward2  * x_; // DXDforward = 0.0

        return discount_ * temp2;
    }

    Real BlackCalculator::theta(Real spot,
                                Time maturity) const {

        QL_REQUIRE(maturity>=0.0,
                   "maturity (" << maturity << ") must be non-negative");
        if (close(maturity, 0.0)) return 0.0;
        return -( std::log(discount_)            * value()
                 +std::log(forward_/spot) * spot * delta(spot)
                 +0.5*variance_ * spot  * spot * gamma(spot))/maturity;
    }

    Real BlackCalculator::vega(Time maturity) const {
        QL_REQUIRE(maturity>=0.0,
                   "negative maturity not allowed");

        // Handle zero volatility case
        if (stdDev_ <= QL_EPSILON) {
            return 0.0;
        }

        Real temp = std::log(strike_/forward_)/variance_;
        // actually DalphaDsigma / SQRT(T)
        Real DalphaDsigma = DalphaDd1_*(temp+0.5);
        Real DbetaDsigma  = DbetaDd2_ *(temp-0.5);

        Real temp2 = DalphaDsigma * forward_ + DbetaDsigma * x_;

        return discount_ * std::sqrt(maturity) * temp2;
    }

    Real BlackCalculator::rho(Time maturity) const {
        QL_REQUIRE(maturity>=0.0,
                   "negative maturity not allowed");

        // Handle zero volatility case
        if (stdDev_ <= QL_EPSILON) {
            // For zero volatility, rho = T * (delta_forward * forward - value/discount)
            Real deltaFwd = deltaForward();
            return maturity * (deltaFwd * forward_ - value());
        }

        // actually DalphaDr / T
        Real DalphaDr = DalphaDd1_/stdDev_;
        Real DbetaDr  = DbetaDd2_/stdDev_;
        Real temp = DalphaDr * forward_ + alpha_ * forward_ + DbetaDr * x_;

        return maturity * (discount_ * temp - value());
    }

    Real BlackCalculator::dividendRho(Time maturity) const {
        QL_REQUIRE(maturity>=0.0,
                   "negative maturity not allowed");

        // Handle zero volatility case
        if (stdDev_ <= QL_EPSILON) {
            // For zero volatility, dividend rho = -T * discount * delta_forward * forward
            Real deltaFwd = deltaForward() / discount_;  // Remove discount to get pure delta
            return -maturity * discount_ * deltaFwd * forward_;
        }

        // actually DalphaDq / T
        Real DalphaDq = -DalphaDd1_/stdDev_;
        Real DbetaDq  = -DbetaDd2_/stdDev_;

        Real temp = DalphaDq * forward_ - alpha_ * forward_ + DbetaDq * x_;

        return maturity * discount_ * temp;
    }

    Real BlackCalculator::strikeSensitivity() const {

        // Handle zero volatility case  
        if (stdDev_ <= QL_EPSILON) {
            // For zero volatility, strike sensitivity is:
            // Call: -N(d2) where d2 -> 1 (ITM), 0 (OTM), 0.5 (ATM)
            // Put: N(-d2) = 1 - N(d2)
            if (close(forward_, strike_)) {
                // ATM case
                if (alpha_ >= 0) { // Call
                    return -discount_ * 0.5;
                } else { // Put
                    return discount_ * 0.5;
                }
            } else if (forward_ > strike_) {
                // ITM Call, OTM Put
                if (alpha_ >= 0) { // Call
                    return -discount_ * 1.0;
                } else { // Put
                    return discount_ * 0.0;
                }
            } else {
                // OTM Call, ITM Put
                if (alpha_ >= 0) { // Call
                    return -discount_ * 0.0;
                } else { // Put
                    return discount_ * 1.0;
                }
            }
        }

        Real temp = stdDev_*strike_;
        Real DalphaDstrike = -DalphaDd1_/temp;
        Real DbetaDstrike  = -DbetaDd2_/temp;

        Real temp2 =
            DalphaDstrike * forward_ + DbetaDstrike * x_ + beta_ * DxDstrike_;

        return discount_ * temp2;
    }


    Real BlackCalculator::strikeGamma() const {

        // Handle zero volatility case
        if (stdDev_ <= QL_EPSILON) {
            // For zero volatility, strike gamma is 0 (no convexity)
            return 0.0;
        }

        Real temp = stdDev_*strike_;
        Real DalphaDstrike = -DalphaDd1_/temp;
        Real DbetaDstrike  = -DbetaDd2_/temp;

        Real D2alphaD2strike = -DalphaDstrike/strike_*(1-d1_/stdDev_);
        Real D2betaD2strike  = -DbetaDstrike /strike_*(1-d2_/stdDev_);

        Real temp2 =
                D2alphaD2strike * forward_ + D2betaD2strike * x_
                + 2.0*DbetaDstrike *DxDstrike_;

        return discount_ * temp2;
    }

    //! Sensitivity of vega to spot (Vanna)
    Real BlackCalculator::vanna(Real spot, Time maturity) const {
        QL_REQUIRE(spot > 0.0, "positive spot value required: " << spot << " not allowed");
        QL_REQUIRE(maturity >= 0.0, "negative maturity not allowed");

        if (stdDev_ <= QL_EPSILON) {
            return 0.0;
        }

        // Calculate Vega
        Real vega = this->vega(maturity);

        // d2 = d1 - stdDev_
        // Vanna = -d2 / (spot * stdDev_) * Vega
        return -d2_ / (spot * stdDev_) * vega;
    }

    //! Sensitivity of vega to volatility (Volga/Vomma)
    Real BlackCalculator::volga(Time maturity) const {
        QL_REQUIRE(maturity >= 0.0, "negative maturity not allowed");

        if (stdDev_ <= QL_EPSILON) {
            return 0.0;
        }

        // Calculate Vega
        Real vega = this->vega(maturity);

        // Volga = Vega * d1 * d2 / stdDev_
        return vega * d1_ * d2_ / stdDev_;
    }
}

Coverage Report

Created: 2026-03-11 06:44

Line	Count	Source
1		/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3		/*
4		Copyright (C) 2003, 2004, 2005, 2006 Ferdinando Ametrano
5		Copyright (C) 2006 StatPro Italia srl
6
7		This file is part of QuantLib, a free-software/open-source library
8		for financial quantitative analysts and developers - http://quantlib.org/
9
10		QuantLib is free software: you can redistribute it and/or modify it
11		under the terms of the QuantLib license. You should have received a
12		copy of the license along with this program; if not, please email
13		<quantlib-dev@lists.sf.net>. The license is also available online at
14		<https://www.quantlib.org/license.shtml>.
15
16		This program is distributed in the hope that it will be useful, but WITHOUT
17		ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18		FOR A PARTICULAR PURPOSE. See the license for more details.
19		*/
20
21		#include <ql/pricingengines/blackcalculator.hpp>
22		#include <ql/math/distributions/normaldistribution.hpp>
23		#include <ql/math/comparison.hpp>
24
25		namespace QuantLib {
26
27		class BlackCalculator::Calculator : public AcyclicVisitor,
28		public Visitor<Payoff>,
29		public Visitor<PlainVanillaPayoff>,
30		public Visitor<CashOrNothingPayoff>,
31		public Visitor<AssetOrNothingPayoff>,
32		public Visitor<GapPayoff> {
33		private:
34		BlackCalculator& black_;
35		public:
36	0	explicit Calculator(BlackCalculator& black) : black_(black) {}
37		void visit(Payoff&) override;
38		void visit(PlainVanillaPayoff&) override;
39		void visit(CashOrNothingPayoff&) override;
40		void visit(AssetOrNothingPayoff&) override;
41		void visit(GapPayoff&) override;
42		};
43
44
45		BlackCalculator::BlackCalculator(const ext::shared_ptr<StrikedTypePayoff>& p,
46		Real forward,
47		Real stdDev,
48		Real discount)
49	0	: strike_(p->strike()), forward_(forward), stdDev_(stdDev),
50	0	discount_(discount), variance_(stdDev*stdDev) {
51	0	initialize(p);
52	0	}
53
54		BlackCalculator::BlackCalculator(Option::Type optionType,
55		Real strike,
56		Real forward,
57		Real stdDev,
58		Real discount)
59	0	: strike_(strike), forward_(forward), stdDev_(stdDev),
60	0	discount_(discount), variance_(stdDev*stdDev) {
61	0	initialize(ext::shared_ptr<StrikedTypePayoff>(new
62	0	PlainVanillaPayoff(optionType, strike)));
63	0	}
64
65	0	void BlackCalculator::initialize(const ext::shared_ptr<StrikedTypePayoff>& p) {
66	0	QL_REQUIRE(strike_>=0.0,
67	0	"strike (" << strike_ << ") must be non-negative");
68	0	QL_REQUIRE(forward_>0.0,
69	0	"forward (" << forward_ << ") must be positive");
70		//QL_REQUIRE(displacement_>=0.0,
71		// "displacement (" << displacement_ << ") must be non-negative");
72	0	QL_REQUIRE(stdDev_>=0.0,
73	0	"stdDev (" << stdDev_ << ") must be non-negative");
74	0	QL_REQUIRE(discount_>0.0,
75	0	"discount (" << discount_ << ") must be positive");
76
77	0	if (stdDev_>=QL_EPSILON) {
78	0	if (close(strike_, 0.0)) {
79	0	d1_ = QL_MAX_REAL;
80	0	d2_ = QL_MAX_REAL;
81	0	cum_d1_ = 1.0;
82	0	cum_d2_ = 1.0;
83	0	n_d1_ = 0.0;
84	0	n_d2_ = 0.0;
85	0	} else {
86	0	d1_ = std::log(forward_/strike_)/stdDev_ + 0.5*stdDev_;
87	0	d2_ = d1_-stdDev_;
88	0	CumulativeNormalDistribution f;
89	0	cum_d1_ = f(d1_);
90	0	cum_d2_ = f(d2_);
91	0	n_d1_ = f.derivative(d1_);
92	0	n_d2_ = f.derivative(d2_);
93	0	}
94	0	} else {
95	0	if (close(forward_, strike_)) {
96	0	d1_ = 0;
97	0	d2_ = 0;
98	0	cum_d1_ = 0.5;
99	0	cum_d2_ = 0.5;
100	0	n_d1_ = M_SQRT_2 * M_1_SQRTPI;
101	0	n_d2_ = M_SQRT_2 * M_1_SQRTPI;
102	0	} else if (forward_>strike_) {
103	0	d1_ = QL_MAX_REAL;
104	0	d2_ = QL_MAX_REAL;
105	0	cum_d1_ = 1.0;
106	0	cum_d2_ = 1.0;
107	0	n_d1_ = 0.0;
108	0	n_d2_ = 0.0;
109	0	} else {
110	0	d1_ = QL_MIN_REAL;
111	0	d2_ = QL_MIN_REAL;
112	0	cum_d1_ = 0.0;
113	0	cum_d2_ = 0.0;
114	0	n_d1_ = 0.0;
115	0	n_d2_ = 0.0;
116	0	}
117	0	}
118
119	0	x_ = strike_;
120	0	DxDstrike_ = 1.0;
121
122		// the following one will probably disappear as soon as
123		// super-share will be properly handled
124	0	DxDs_ = 0.0;
125
126		// this part is always executed.
127		// in case of plain-vanilla payoffs, it is also the only part
128		// which is executed.
129	0	switch (p->optionType()) {
130	0	case Option::Call:
131	0	alpha_ = cum_d1_;// N(d1)
132	0	DalphaDd1_ = n_d1_;// n(d1)
133	0	beta_ = -cum_d2_;// -N(d2)
134	0	DbetaDd2_ = - n_d2_;// -n(d2)
135	0	break;
136	0	case Option::Put:
137	0	alpha_ = -1.0+cum_d1_;// -N(-d1)
138	0	DalphaDd1_ = n_d1_;// n( d1)
139	0	beta_ = 1.0-cum_d2_;// N(-d2)
140	0	DbetaDd2_ = - n_d2_;// -n( d2)
141	0	break;
142	0	default:
143	0	QL_FAIL("invalid option type");
144	0	}
145
146		// now dispatch on type.
147
148	0	Calculator calc(*this);
149	0	p->accept(calc);
150	0	}
151
152	0	void BlackCalculator::Calculator::visit(Payoff& p) {
153	0	QL_FAIL("unsupported payoff type: " << p.name());
154	0	}
155
156	0	void BlackCalculator::Calculator::visit(PlainVanillaPayoff&) {}
157
158	0	void BlackCalculator::Calculator::visit(CashOrNothingPayoff& payoff) {
159	0	black_.alpha_ = black_.DalphaDd1_ = 0.0;
160	0	black_.x_ = payoff.cashPayoff();
161	0	black_.DxDstrike_ = 0.0;
162	0	switch (payoff.optionType()) {
163	0	case Option::Call:
164	0	black_.beta_ = black_.cum_d2_;
165	0	black_.DbetaDd2_ = black_.n_d2_;
166	0	break;
167	0	case Option::Put:
168	0	black_.beta_ = 1.0-black_.cum_d2_;
169	0	black_.DbetaDd2_ = -black_.n_d2_;
170	0	break;
171	0	default:
172	0	QL_FAIL("invalid option type");
173	0	}
174	0	}
175
176	0	void BlackCalculator::Calculator::visit(AssetOrNothingPayoff& payoff) {
177	0	black_.beta_ = black_.DbetaDd2_ = 0.0;
178	0	switch (payoff.optionType()) {
179	0	case Option::Call:
180	0	black_.alpha_ = black_.cum_d1_;
181	0	black_.DalphaDd1_ = black_.n_d1_;
182	0	break;
183	0	case Option::Put:
184	0	black_.alpha_ = 1.0-black_.cum_d1_;
185	0	black_.DalphaDd1_ = -black_.n_d1_;
186	0	break;
187	0	default:
188	0	QL_FAIL("invalid option type");
189	0	}
190	0	}
191
192	0	void BlackCalculator::Calculator::visit(GapPayoff& payoff) {
193	0	black_.x_ = payoff.secondStrike();
194	0	black_.DxDstrike_ = 0.0;
195	0	}
196
197	0	Real BlackCalculator::value() const {
198	0	Real result = discount_ * (forward_ * alpha_ + x_ * beta_);
199	0	return result;
200	0	}
201
202	0	Real BlackCalculator::delta(Real spot) const {
203
204	0	QL_REQUIRE(spot > 0.0, "positive spot value required: " <<
205	0	spot << " not allowed");
206
207		// Handle zero volatility case
208	0	if (stdDev_ <= QL_EPSILON) {
209		// For zero volatility, delta is:
210		// ITM Call: 1.0, OTM Call: 0.0, ATM Call: 0.5
211		// ITM Put: -1.0, OTM Put: 0.0, ATM Put: -0.5
212	0	Real DforwardDs = forward_ / spot;
213	0	if (close(forward_, strike_)) {
214		// ATM case
215	0	if (alpha_ >= 0) { // Call
216	0	return discount_ * 0.5 * DforwardDs;
217	0	} else { // Put
218	0	return discount_ * (-0.5) * DforwardDs;
219	0	}
220	0	} else if (forward_ > strike_) {
221		// ITM Call, OTM Put
222	0	if (alpha_ >= 0) { // Call
223	0	return discount_ * 1.0 * DforwardDs;
224	0	} else { // Put
225	0	return 0.0;
226	0	}
227	0	} else {
228		// OTM Call, ITM Put
229	0	if (alpha_ >= 0) { // Call
230	0	return 0.0;
231	0	} else { // Put
232	0	return discount_ * (-1.0) * DforwardDs;
233	0	}
234	0	}
235	0	}
236
237	0	Real DforwardDs = forward_ / spot;
238
239	0	Real temp = stdDev_*spot;
240	0	Real DalphaDs = DalphaDd1_/temp;
241	0	Real DbetaDs = DbetaDd2_/temp;
242	0	Real temp2 = DalphaDs * forward_ + alpha_ * DforwardDs
243	0	+DbetaDs * x_ + beta_ * DxDs_;
244
245	0	return discount_ * temp2;
246	0	}
247
248	0	Real BlackCalculator::deltaForward() const {
249
250		// Handle zero volatility case
251	0	if (stdDev_ <= QL_EPSILON) {
252		// For zero volatility, forward delta is:
253		// ITM Call: 1.0, OTM Call: 0.0, ATM Call: 0.5
254		// ITM Put: -1.0, OTM Put: 0.0, ATM Put: -0.5
255	0	if (close(forward_, strike_)) {
256		// ATM case
257	0	if (alpha_ >= 0) { // Call
258	0	return discount_ * 0.5;
259	0	} else { // Put
260	0	return discount_ * (-0.5);
261	0	}
262	0	} else if (forward_ > strike_) {
263		// ITM Call, OTM Put
264	0	if (alpha_ >= 0) { // Call
265	0	return discount_ * 1.0;
266	0	} else { // Put
267	0	return 0.0;
268	0	}
269	0	} else {
270		// OTM Call, ITM Put
271	0	if (alpha_ >= 0) { // Call
272	0	return 0.0;
273	0	} else { // Put
274	0	return discount_ * (-1.0);
275	0	}
276	0	}
277	0	}
278
279	0	Real temp = stdDev_*forward_;
280	0	Real DalphaDforward = DalphaDd1_/temp;
281	0	Real DbetaDforward = DbetaDd2_/temp;
282	0	Real temp2 = DalphaDforward * forward_ + alpha_
283	0	+DbetaDforward * x_; // DXDforward = 0.0
284
285	0	return discount_ * temp2;
286	0	}
287
288	0	Real BlackCalculator::elasticity(Real spot) const {
289	0	Real val = value();
290	0	Real del = delta(spot);
291	0	if (val>QL_EPSILON)
292	0	return del/val*spot;
293	0	else if (std::fabs(del)<QL_EPSILON)
294	0	return 0.0;
295	0	else if (del>0.0)
296	0	return QL_MAX_REAL;
297	0	else
298	0	return QL_MIN_REAL;
299	0	}
300
301	0	Real BlackCalculator::elasticityForward() const {
302	0	Real val = value();
303	0	Real del = deltaForward();
304	0	if (val>QL_EPSILON)
305	0	return del/val*forward_;
306	0	else if (std::fabs(del)<QL_EPSILON)
307	0	return 0.0;
308	0	else if (del>0.0)
309	0	return QL_MAX_REAL;
310	0	else
311	0	return QL_MIN_REAL;
312	0	}
313
314	0	Real BlackCalculator::gamma(Real spot) const {
315
316	0	QL_REQUIRE(spot > 0.0, "positive spot value required: " <<
317	0	spot << " not allowed");
318
319		// Handle zero volatility case
320	0	if (stdDev_ <= QL_EPSILON) {
321		// For zero volatility, gamma is always 0 (no convexity)
322	0	return 0.0;
323	0	}
324
325	0	Real DforwardDs = forward_ / spot;
326
327	0	Real temp = stdDev_*spot;
328	0	Real DalphaDs = DalphaDd1_/temp;
329	0	Real DbetaDs = DbetaDd2_/temp;
330
331	0	Real D2alphaDs2 = - DalphaDs/spot*(1+d1_/stdDev_);
332	0	Real D2betaDs2 = - DbetaDs /spot*(1+d2_/stdDev_);
333
334	0	Real temp2 = D2alphaDs2 * forward_ + 2.0 * DalphaDs * DforwardDs
335	0	+D2betaDs2 * x_ + 2.0 * DbetaDs * DxDs_;
336
337	0	return discount_ * temp2;
338	0	}
339
340	0	Real BlackCalculator::gammaForward() const {
341
342		// Handle zero volatility case
343	0	if (stdDev_ <= QL_EPSILON) {
344		// For zero volatility, gamma is always 0 (no convexity)
345	0	return 0.0;
346	0	}
347
348	0	Real temp = stdDev_*forward_;
349	0	Real DalphaDforward = DalphaDd1_/temp;
350	0	Real DbetaDforward = DbetaDd2_/temp;
351
352	0	Real D2alphaDforward2 = - DalphaDforward/forward_*(1+d1_/stdDev_);
353	0	Real D2betaDforward2 = - DbetaDforward /forward_*(1+d2_/stdDev_);
354
355	0	Real temp2 = D2alphaDforward2 * forward_ + 2.0 * DalphaDforward
356	0	+D2betaDforward2 * x_; // DXDforward = 0.0
357
358	0	return discount_ * temp2;
359	0	}
360
361		Real BlackCalculator::theta(Real spot,
362	0	Time maturity) const {
363
364	0	QL_REQUIRE(maturity>=0.0,
365	0	"maturity (" << maturity << ") must be non-negative");
366	0	if (close(maturity, 0.0)) return 0.0;
367	0	return -( std::log(discount_) * value()
368	0	+std::log(forward_/spot) * spot * delta(spot)
369	0	+0.5variance_ spot * spot * gamma(spot))/maturity;
370	0	}
371
372	0	Real BlackCalculator::vega(Time maturity) const {
373	0	QL_REQUIRE(maturity>=0.0,
374	0	"negative maturity not allowed");
375
376		// Handle zero volatility case
377	0	if (stdDev_ <= QL_EPSILON) {
378	0	return 0.0;
379	0	}
380
381	0	Real temp = std::log(strike_/forward_)/variance_;
382		// actually DalphaDsigma / SQRT(T)
383	0	Real DalphaDsigma = DalphaDd1_*(temp+0.5);
384	0	Real DbetaDsigma = DbetaDd2_ *(temp-0.5);
385
386	0	Real temp2 = DalphaDsigma * forward_ + DbetaDsigma * x_;
387
388	0	return discount_ * std::sqrt(maturity) * temp2;
389	0	}
390
391	0	Real BlackCalculator::rho(Time maturity) const {
392	0	QL_REQUIRE(maturity>=0.0,
393	0	"negative maturity not allowed");
394
395		// Handle zero volatility case
396	0	if (stdDev_ <= QL_EPSILON) {
397		// For zero volatility, rho = T * (delta_forward * forward - value/discount)
398	0	Real deltaFwd = deltaForward();
399	0	return maturity * (deltaFwd * forward_ - value());
400	0	}
401
402		// actually DalphaDr / T
403	0	Real DalphaDr = DalphaDd1_/stdDev_;
404	0	Real DbetaDr = DbetaDd2_/stdDev_;
405	0	Real temp = DalphaDr * forward_ + alpha_ * forward_ + DbetaDr * x_;
406
407	0	return maturity * (discount_ * temp - value());
408	0	}
409
410	0	Real BlackCalculator::dividendRho(Time maturity) const {
411	0	QL_REQUIRE(maturity>=0.0,
412	0	"negative maturity not allowed");
413
414		// Handle zero volatility case
415	0	if (stdDev_ <= QL_EPSILON) {
416		// For zero volatility, dividend rho = -T * discount * delta_forward * forward
417	0	Real deltaFwd = deltaForward() / discount_; // Remove discount to get pure delta
418	0	return -maturity * discount_ * deltaFwd * forward_;
419	0	}
420
421		// actually DalphaDq / T
422	0	Real DalphaDq = -DalphaDd1_/stdDev_;
423	0	Real DbetaDq = -DbetaDd2_/stdDev_;
424
425	0	Real temp = DalphaDq * forward_ - alpha_ * forward_ + DbetaDq * x_;
426
427	0	return maturity * discount_ * temp;
428	0	}
429
430	0	Real BlackCalculator::strikeSensitivity() const {
431
432		// Handle zero volatility case
433	0	if (stdDev_ <= QL_EPSILON) {
434		// For zero volatility, strike sensitivity is:
435		// Call: -N(d2) where d2 -> 1 (ITM), 0 (OTM), 0.5 (ATM)
436		// Put: N(-d2) = 1 - N(d2)
437	0	if (close(forward_, strike_)) {
438		// ATM case
439	0	if (alpha_ >= 0) { // Call
440	0	return -discount_ * 0.5;
441	0	} else { // Put
442	0	return discount_ * 0.5;
443	0	}
444	0	} else if (forward_ > strike_) {
445		// ITM Call, OTM Put
446	0	if (alpha_ >= 0) { // Call
447	0	return -discount_ * 1.0;
448	0	} else { // Put
449	0	return discount_ * 0.0;
450	0	}
451	0	} else {
452		// OTM Call, ITM Put
453	0	if (alpha_ >= 0) { // Call
454	0	return -discount_ * 0.0;
455	0	} else { // Put
456	0	return discount_ * 1.0;
457	0	}
458	0	}
459	0	}
460
461	0	Real temp = stdDev_*strike_;
462	0	Real DalphaDstrike = -DalphaDd1_/temp;
463	0	Real DbetaDstrike = -DbetaDd2_/temp;
464
465	0	Real temp2 =
466	0	DalphaDstrike * forward_ + DbetaDstrike * x_ + beta_ * DxDstrike_;
467
468	0	return discount_ * temp2;
469	0	}
470
471
472	0	Real BlackCalculator::strikeGamma() const {
473
474		// Handle zero volatility case
475	0	if (stdDev_ <= QL_EPSILON) {
476		// For zero volatility, strike gamma is 0 (no convexity)
477	0	return 0.0;
478	0	}
479
480	0	Real temp = stdDev_*strike_;
481	0	Real DalphaDstrike = -DalphaDd1_/temp;
482	0	Real DbetaDstrike = -DbetaDd2_/temp;
483
484	0	Real D2alphaD2strike = -DalphaDstrike/strike_*(1-d1_/stdDev_);
485	0	Real D2betaD2strike = -DbetaDstrike /strike_*(1-d2_/stdDev_);
486
487	0	Real temp2 =
488	0	D2alphaD2strike * forward_ + D2betaD2strike * x_
489	0	+ 2.0DbetaDstrike DxDstrike_;
490
491	0	return discount_ * temp2;
492	0	}
493
494		//! Sensitivity of vega to spot (Vanna)
495	0	Real BlackCalculator::vanna(Real spot, Time maturity) const {
496	0	QL_REQUIRE(spot > 0.0, "positive spot value required: " << spot << " not allowed");
497	0	QL_REQUIRE(maturity >= 0.0, "negative maturity not allowed");
498
499	0	if (stdDev_ <= QL_EPSILON) {
500	0	return 0.0;
501	0	}
502
503		// Calculate Vega
504	0	Real vega = this->vega(maturity);
505
506		// d2 = d1 - stdDev_
507		// Vanna = -d2 / (spot * stdDev_) * Vega
508	0	return -d2_ / (spot * stdDev_) * vega;
509	0	}
510
511		//! Sensitivity of vega to volatility (Volga/Vomma)
512	0	Real BlackCalculator::volga(Time maturity) const {
513	0	QL_REQUIRE(maturity >= 0.0, "negative maturity not allowed");
514
515	0	if (stdDev_ <= QL_EPSILON) {
516	0	return 0.0;
517	0	}
518
519		// Calculate Vega
520	0	Real vega = this->vega(maturity);
521
522		// Volga = Vega * d1 * d2 / stdDev_
523	0	return vega * d1_ * d2_ / stdDev_;
524	0	}
525		}