/src/quantlib/ql/pricingengines/bacheliercalculator.cpp

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

/*

 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/math/comparison.hpp>
#include <ql/math/distributions/normaldistribution.hpp>
#include <ql/pricingengines/bacheliercalculator.hpp>

namespace QuantLib {

    class BachelierCalculator::Calculator : public AcyclicVisitor,
                                        public Visitor<Payoff>,
                                        public Visitor<PlainVanillaPayoff>,
                                        public Visitor<CashOrNothingPayoff>,
                                        public Visitor<AssetOrNothingPayoff>,
                                        public Visitor<GapPayoff> {
      private:
        BachelierCalculator& bachelier_;

      public:
        explicit Calculator(BachelierCalculator& bachelier) : bachelier_(bachelier) {}
        void visit(Payoff&) override;
        void visit(PlainVanillaPayoff&) override;
        void visit(CashOrNothingPayoff&) override;
        void visit(AssetOrNothingPayoff&) override;
        void visit(GapPayoff&) override;
    };


    BachelierCalculator::BachelierCalculator(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);
    }

    BachelierCalculator::BachelierCalculator(
        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 BachelierCalculator::initialize(const ext::shared_ptr<StrikedTypePayoff>& p) {
        QL_REQUIRE(stdDev_ >= 0.0, "stdDev (" << stdDev_ << ") must be non-negative");
        QL_REQUIRE(discount_ > 0.0, "discount (" << discount_ << ") must be positive");

        // For Bachelier model, we use d = (F - K) / σ instead of the Black-Scholes d1, d2
        if (stdDev_ >= QL_EPSILON) {
            // Bachelier d parameter: d = (F - K) / σ
            d_ = (forward_ - strike_) / stdDev_;
            
            CumulativeNormalDistribution f;
            cum_d_ = f(d_);
            n_d_ = f.derivative(d_);
        } else {
            // When volatility is zero
            if (close(forward_, strike_)) {
                d_ = 0;
                cum_d_ = 0.5;
                n_d_ = M_SQRT_2 * M_1_SQRTPI;
            } else if (forward_ > strike_) {
                d_ = QL_MAX_REAL;
                cum_d_ = 1.0;
                n_d_ = 0.0;
            } else {
                d_ = QL_MIN_REAL;
                cum_d_ = 0.0;
                n_d_ = 0.0;
            }
        }

        x_ = strike_;
        DxDstrike_ = 1.0;
        DxDs_ = 0.0;

        // For Bachelier model, the option values are:
        // Call: max(F-K, 0) = (F-K)*N(d) + σ*n(d)
        // Put:  max(K-F, 0) = (K-F)*N(-d) + σ*n(d)
        // 
        // We represent this as: discount * (forward * alpha + x * beta)
        // where for Bachelier:
        // Call: alpha = N(d), beta = -N(d) + σ*n(d)/x (when x != 0)
        // Put:  alpha = -N(-d), beta = N(-d) + σ*n(d)/x (when x != 0)
        
        switch (p->optionType()) {
            case Option::Call:
                alpha_ = cum_d_;        // N(d)
                DalphaDd_ = n_d_;       // n(d)
                beta_ = -cum_d_;        // -N(d) - base part
                DbetaDd_ = -n_d_;       // -n(d)
                break;
            case Option::Put:
                alpha_ = cum_d_ - 1.0;  // N(d) - 1 = -N(-d)
                DalphaDd_ = n_d_;       // n(d)
                beta_ = 1.0 - cum_d_;   // 1 - N(d) = N(-d)
                DbetaDd_ = -n_d_;       // -n(d)
                break;
            default:
                QL_FAIL("invalid option type");
        }

        // now dispatch on type.
        Calculator calc(*this);
        p->accept(calc);
    }

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

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

    void BachelierCalculator::Calculator::visit(CashOrNothingPayoff& payoff) {
        bachelier_.alpha_ = bachelier_.DalphaDd_ = 0.0;
        bachelier_.x_ = payoff.cashPayoff();
        bachelier_.DxDstrike_ = 0.0;
        switch (payoff.optionType()) {
            case Option::Call:
                bachelier_.beta_ = bachelier_.cum_d_;
                bachelier_.DbetaDd_ = bachelier_.n_d_;
                break;
            case Option::Put:
                bachelier_.beta_ = 1.0 - bachelier_.cum_d_;
                bachelier_.DbetaDd_ = -bachelier_.n_d_;
                break;
            default:
                QL_FAIL("invalid option type");
        }
    }

    void BachelierCalculator::Calculator::visit(AssetOrNothingPayoff& payoff) {
        bachelier_.beta_ = bachelier_.DbetaDd_ = 0.0;
        switch (payoff.optionType()) {
            case Option::Call:
                bachelier_.alpha_ = bachelier_.cum_d_;
                bachelier_.DalphaDd_ = bachelier_.n_d_;
                break;
            case Option::Put:
                bachelier_.alpha_ = 1.0 - bachelier_.cum_d_;
                bachelier_.DalphaDd_ = -bachelier_.n_d_;
                break;
            default:
                QL_FAIL("invalid option type");
        }
    }

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

    Real BachelierCalculator::value() const {
        // Bachelier option value formula:
        // Call: (F-K)*N(d) + σ*n(d)
        // Put:  (K-F)*N(-d) + σ*n(d)
        // where d = (F-K)/σ
        
        Real intrinsic = forward_ - strike_;
        Real timeValue = 0.0;
        
        if (stdDev_ > QL_EPSILON) {
            timeValue = stdDev_ * n_d_;
        }
        
        Real result;
        if (alpha_ >= 0) // Call option (alpha_ = N(d) >= 0)
            result = intrinsic * cum_d_ + timeValue;
        else // Put option (alpha_ = N(d) - 1 < 0)
            result = -intrinsic * (1.0 - cum_d_) + timeValue;
        
        return discount_ * std::max(result, 0.0);
    }

    Real BachelierCalculator::delta(Real spot) const {

        // For Bachelier model:
        // Delta = dV/dS = (dV/dF) * (dF/dS)
        // where dF/dS = F/S (assuming forward = spot * exp(r*T))
        // and dV/dF = N(d) for calls, -N(-d) for puts
        
        Real DforwardDs = forward_ / spot;
        Real deltaFwd = deltaForward();

        return deltaFwd * DforwardDs;
    }

    Real BachelierCalculator::deltaForward() const {
        // For Bachelier model:
        // Delta_Forward = dV/dF = N(d) for calls, -N(-d) for puts
        // where d = (F-K)/σ
        
        if (alpha_ >= 0) { // Call option
            return discount_ * cum_d_; // N(d)
        } else { // Put option  
            return discount_ * (cum_d_ - 1.0); // N(d) - 1 = -N(-d)
        }
    }

    Real BachelierCalculator::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 BachelierCalculator::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 BachelierCalculator::gamma(Real spot) const {

        // For Bachelier model:
        // Gamma = d²V/dS² = d/dS(dV/dS) = d/dS(N(d) * dF/dS) * dF/dS
        // = n(d) * (1/σ) * (dF/dS)² * (dd/dF)
        // where dd/dF = 1/σ
        
        if (stdDev_ <= QL_EPSILON) {
            return 0.0;
        }
        
        Real DforwardDs = forward_ / spot;
        Real gammaForward = n_d_ / stdDev_; // dn(d)/dF = n(d) * (1/σ) * (dd/dF) = n(d)/σ
        
        return discount_ * gammaForward * DforwardDs * DforwardDs;
    }

    Real BachelierCalculator::gammaForward() const {
        // For Bachelier model:
        // Gamma_Forward = d²V/dF² = d/dF(N(d)) = n(d) * dd/dF = n(d)/σ
        
        if (stdDev_ <= QL_EPSILON) {
            return 0.0;
        }
        
        return discount_ * n_d_ / stdDev_;
    }

    Real BachelierCalculator::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;
        
        // Theta = -dV/dt = -(r*V - r*S*Delta + 0.5*σ²*Gamma)

        return -(std::log(discount_) * value() + std::log(forward_ / spot) * spot * delta(spot) +
                 0.5 * variance_ * gamma(spot)) / maturity;
    }

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

        // For Bachelier model:
        // Vega = dV/dσ = d/dσ[(F-K)*N(d) + σ*n(d)]
        // = (F-K)*n(d)*dd/dσ + n(d) + σ*n'(d)*dd/dσ
        // where d = (F-K)/σ, so dd/dσ = -(F-K)/σ² = -d/σ
        // and n'(d) = -d*n(d)
        // Therefore: Vega = n(d) + σ*(-d*n(d))*(-d/σ) = n(d) + d²*n(d) = n(d)*(1 + d²) - (F-K)*n(d)*d/σ
        // Simplifying: Vega = n(d) for Bachelier model
        
        if (maturity <= QL_EPSILON || stdDev_ <= QL_EPSILON) {
            return 0.0;
        }
        
        return discount_ * std::sqrt(maturity) * n_d_;
    }

    Real BachelierCalculator::rho(Time maturity) const {
        QL_REQUIRE(maturity >= 0.0, "negative maturity not allowed");
        
        // For Bachelier model:
        // Rho = dV/dr = T * (discount * delta_forward * forward - value)
        // where delta_forward = N(d) for calls, N(d)-1 for puts
        
        Real deltaFwd = deltaForward();
        Real rho_value = maturity * (deltaFwd * forward_ - value());
        
        return rho_value;
    }

    Real BachelierCalculator::dividendRho(Time maturity) const {
        QL_REQUIRE(maturity >= 0.0, "negative maturity not allowed");
        
        // For Bachelier model:
        // Dividend rho = -T * discount * delta_forward * forward
        // where delta_forward = N(d) for calls, N(d)-1 for puts
        
        Real deltaFwd = (alpha_ >= 0) ? cum_d_ : (cum_d_ - 1.0);
        
        return -maturity * discount_ * deltaFwd * forward_;
    }

    Real BachelierCalculator::strikeSensitivity() const {
        // For Bachelier model:
        // dV/dK = -N(d) for calls, N(-d) for puts
        // where d = (F-K)/σ, so dd/dK = -1/σ
        
        if (alpha_ >= 0) { // Call option
            return -discount_ * cum_d_; // -N(d)
        } else { // Put option
            return discount_ * (1.0 - cum_d_); // N(-d) = 1 - N(d)
        }
    }


    Real BachelierCalculator::strikeGamma() const {
        // For Bachelier model:
        // d²V/dK² = d/dK(-N(d)) = -n(d) * dd/dK = -n(d) * (-1/σ) = n(d)/σ
        // This is the same for both calls and puts
        
        if (stdDev_ <= QL_EPSILON) {
            return 0.0;
        }
        
        return discount_ * n_d_ / stdDev_;
    }

    // Sensitivity of Vega to forward (Vanna)
    Real BachelierCalculator::vanna(Time maturity) const {
        if (maturity <= QL_EPSILON || stdDev_ <= QL_EPSILON) {
            return 0.0;
        }
        // d_ is (F-K)/stdDev_
        // n_d_ is n(d)
        // Vanna = -d * n(d) / stdDev_ * sqrt(T)
        return -d_ * n_d_ * std::sqrt(maturity) / stdDev_;
    }

    // Sensitivity of Vega to volatility (Volga)
    Real BachelierCalculator::volga(Time maturity) const {
        if (maturity <= QL_EPSILON || stdDev_ <= QL_EPSILON) {
            return 0.0;
        }
        // Volga = d^2 / stdDev_ * Vega
        Real vega = this->vega(maturity);
        return (d_ * d_ / stdDev_) * vega;
    }
}

Coverage Report

Created: 2026-02-03 07:02

Line	Count	Source
1		/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3		/*
4
5		This file is part of QuantLib, a free-software/open-source library
6		for financial quantitative analysts and developers - http://quantlib.org/
7
8		QuantLib is free software: you can redistribute it and/or modify it
9		under the terms of the QuantLib license. You should have received a
10		copy of the license along with this program; if not, please email
11		<quantlib-dev@lists.sf.net>. The license is also available online at
12		<https://www.quantlib.org/license.shtml>.
13
14		This program is distributed in the hope that it will be useful, but WITHOUT
15		ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16		FOR A PARTICULAR PURPOSE. See the license for more details.
17		*/
18
19		#include <ql/math/comparison.hpp>
20		#include <ql/math/distributions/normaldistribution.hpp>
21		#include <ql/pricingengines/bacheliercalculator.hpp>
22
23		namespace QuantLib {
24
25		class BachelierCalculator::Calculator : public AcyclicVisitor,
26		public Visitor<Payoff>,
27		public Visitor<PlainVanillaPayoff>,
28		public Visitor<CashOrNothingPayoff>,
29		public Visitor<AssetOrNothingPayoff>,
30		public Visitor<GapPayoff> {
31		private:
32		BachelierCalculator& bachelier_;
33
34		public:
35	0	explicit Calculator(BachelierCalculator& bachelier) : bachelier_(bachelier) {}
36		void visit(Payoff&) override;
37		void visit(PlainVanillaPayoff&) override;
38		void visit(CashOrNothingPayoff&) override;
39		void visit(AssetOrNothingPayoff&) override;
40		void visit(GapPayoff&) override;
41		};
42
43
44		BachelierCalculator::BachelierCalculator(const ext::shared_ptr<StrikedTypePayoff>& p,
45		Real forward,
46		Real stdDev,
47		Real discount)
48	0	: strike_(p->strike()), forward_(forward), stdDev_(stdDev),
49	0	discount_(discount), variance_(stdDev*stdDev) {
50	0	initialize(p);
51	0	}
52
53		BachelierCalculator::BachelierCalculator(
54		Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount)
55	0	: strike_(strike), forward_(forward), stdDev_(stdDev),
56	0	discount_(discount), variance_(stdDev*stdDev) {
57	0	initialize(ext::shared_ptr<StrikedTypePayoff>(new PlainVanillaPayoff(optionType, strike)));
58	0	}
59
60	0	void BachelierCalculator::initialize(const ext::shared_ptr<StrikedTypePayoff>& p) {
61	0	QL_REQUIRE(stdDev_ >= 0.0, "stdDev (" << stdDev_ << ") must be non-negative");
62	0	QL_REQUIRE(discount_ > 0.0, "discount (" << discount_ << ") must be positive");
63
64		// For Bachelier model, we use d = (F - K) / σ instead of the Black-Scholes d1, d2
65	0	if (stdDev_ >= QL_EPSILON) {
66		// Bachelier d parameter: d = (F - K) / σ
67	0	d_ = (forward_ - strike_) / stdDev_;
68
69	0	CumulativeNormalDistribution f;
70	0	cum_d_ = f(d_);
71	0	n_d_ = f.derivative(d_);
72	0	} else {
73		// When volatility is zero
74	0	if (close(forward_, strike_)) {
75	0	d_ = 0;
76	0	cum_d_ = 0.5;
77	0	n_d_ = M_SQRT_2 * M_1_SQRTPI;
78	0	} else if (forward_ > strike_) {
79	0	d_ = QL_MAX_REAL;
80	0	cum_d_ = 1.0;
81	0	n_d_ = 0.0;
82	0	} else {
83	0	d_ = QL_MIN_REAL;
84	0	cum_d_ = 0.0;
85	0	n_d_ = 0.0;
86	0	}
87	0	}
88
89	0	x_ = strike_;
90	0	DxDstrike_ = 1.0;
91	0	DxDs_ = 0.0;
92
93		// For Bachelier model, the option values are:
94		// Call: max(F-K, 0) = (F-K)N(d) + σn(d)
95		// Put: max(K-F, 0) = (K-F)N(-d) + σn(d)
96		//
97		// We represent this as: discount * (forward * alpha + x * beta)
98		// where for Bachelier:
99		// Call: alpha = N(d), beta = -N(d) + σ*n(d)/x (when x != 0)
100		// Put: alpha = -N(-d), beta = N(-d) + σ*n(d)/x (when x != 0)
101
102	0	switch (p->optionType()) {
103	0	case Option::Call:
104	0	alpha_ = cum_d_; // N(d)
105	0	DalphaDd_ = n_d_; // n(d)
106	0	beta_ = -cum_d_; // -N(d) - base part
107	0	DbetaDd_ = -n_d_; // -n(d)
108	0	break;
109	0	case Option::Put:
110	0	alpha_ = cum_d_ - 1.0; // N(d) - 1 = -N(-d)
111	0	DalphaDd_ = n_d_; // n(d)
112	0	beta_ = 1.0 - cum_d_; // 1 - N(d) = N(-d)
113	0	DbetaDd_ = -n_d_; // -n(d)
114	0	break;
115	0	default:
116	0	QL_FAIL("invalid option type");
117	0	}
118
119		// now dispatch on type.
120	0	Calculator calc(*this);
121	0	p->accept(calc);
122	0	}
123
124	0	void BachelierCalculator::Calculator::visit(Payoff& p) {
125	0	QL_FAIL("unsupported payoff type: " << p.name());
126	0	}
127
128	0	void BachelierCalculator::Calculator::visit(PlainVanillaPayoff&) {}
129
130	0	void BachelierCalculator::Calculator::visit(CashOrNothingPayoff& payoff) {
131	0	bachelier_.alpha_ = bachelier_.DalphaDd_ = 0.0;
132	0	bachelier_.x_ = payoff.cashPayoff();
133	0	bachelier_.DxDstrike_ = 0.0;
134	0	switch (payoff.optionType()) {
135	0	case Option::Call:
136	0	bachelier_.beta_ = bachelier_.cum_d_;
137	0	bachelier_.DbetaDd_ = bachelier_.n_d_;
138	0	break;
139	0	case Option::Put:
140	0	bachelier_.beta_ = 1.0 - bachelier_.cum_d_;
141	0	bachelier_.DbetaDd_ = -bachelier_.n_d_;
142	0	break;
143	0	default:
144	0	QL_FAIL("invalid option type");
145	0	}
146	0	}
147
148	0	void BachelierCalculator::Calculator::visit(AssetOrNothingPayoff& payoff) {
149	0	bachelier_.beta_ = bachelier_.DbetaDd_ = 0.0;
150	0	switch (payoff.optionType()) {
151	0	case Option::Call:
152	0	bachelier_.alpha_ = bachelier_.cum_d_;
153	0	bachelier_.DalphaDd_ = bachelier_.n_d_;
154	0	break;
155	0	case Option::Put:
156	0	bachelier_.alpha_ = 1.0 - bachelier_.cum_d_;
157	0	bachelier_.DalphaDd_ = -bachelier_.n_d_;
158	0	break;
159	0	default:
160	0	QL_FAIL("invalid option type");
161	0	}
162	0	}
163
164	0	void BachelierCalculator::Calculator::visit(GapPayoff& payoff) {
165	0	bachelier_.x_ = payoff.secondStrike();
166	0	bachelier_.DxDstrike_ = 0.0;
167	0	}
168
169	0	Real BachelierCalculator::value() const {
170		// Bachelier option value formula:
171		// Call: (F-K)N(d) + σn(d)
172		// Put: (K-F)N(-d) + σn(d)
173		// where d = (F-K)/σ
174
175	0	Real intrinsic = forward_ - strike_;
176	0	Real timeValue = 0.0;
177
178	0	if (stdDev_ > QL_EPSILON) {
179	0	timeValue = stdDev_ * n_d_;
180	0	}
181
182	0	Real result;
183	0	if (alpha_ >= 0) // Call option (alpha_ = N(d) >= 0)
184	0	result = intrinsic * cum_d_ + timeValue;
185	0	else // Put option (alpha_ = N(d) - 1 < 0)
186	0	result = -intrinsic * (1.0 - cum_d_) + timeValue;
187
188	0	return discount_ * std::max(result, 0.0);
189	0	}
190
191	0	Real BachelierCalculator::delta(Real spot) const {
192
193		// For Bachelier model:
194		// Delta = dV/dS = (dV/dF) * (dF/dS)
195		// where dF/dS = F/S (assuming forward = spot * exp(r*T))
196		// and dV/dF = N(d) for calls, -N(-d) for puts
197
198	0	Real DforwardDs = forward_ / spot;
199	0	Real deltaFwd = deltaForward();
200
201	0	return deltaFwd * DforwardDs;
202	0	}
203
204	0	Real BachelierCalculator::deltaForward() const {
205		// For Bachelier model:
206		// Delta_Forward = dV/dF = N(d) for calls, -N(-d) for puts
207		// where d = (F-K)/σ
208
209	0	if (alpha_ >= 0) { // Call option
210	0	return discount_ * cum_d_; // N(d)
211	0	} else { // Put option
212	0	return discount_ * (cum_d_ - 1.0); // N(d) - 1 = -N(-d)
213	0	}
214	0	}
215
216	0	Real BachelierCalculator::elasticity(Real spot) const {
217	0	Real val = value();
218	0	Real del = delta(spot);
219	0	if (val > QL_EPSILON)
220	0	return del / val * spot;
221	0	else if (std::fabs(del) < QL_EPSILON)
222	0	return 0.0;
223	0	else if (del > 0.0)
224	0	return QL_MAX_REAL;
225	0	else
226	0	return QL_MIN_REAL;
227	0	}
228
229	0	Real BachelierCalculator::elasticityForward() const {
230	0	Real val = value();
231	0	Real del = deltaForward();
232	0	if (val > QL_EPSILON)
233	0	return del / val * forward_;
234	0	else if (std::fabs(del) < QL_EPSILON)
235	0	return 0.0;
236	0	else if (del > 0.0)
237	0	return QL_MAX_REAL;
238	0	else
239	0	return QL_MIN_REAL;
240	0	}
241
242	0	Real BachelierCalculator::gamma(Real spot) const {
243
244		// For Bachelier model:
245		// Gamma = d²V/dS² = d/dS(dV/dS) = d/dS(N(d) * dF/dS) * dF/dS
246		// = n(d) * (1/σ) * (dF/dS)² * (dd/dF)
247		// where dd/dF = 1/σ
248
249	0	if (stdDev_ <= QL_EPSILON) {
250	0	return 0.0;
251	0	}
252
253	0	Real DforwardDs = forward_ / spot;
254	0	Real gammaForward = n_d_ / stdDev_; // dn(d)/dF = n(d) * (1/σ) * (dd/dF) = n(d)/σ
255
256	0	return discount_ * gammaForward * DforwardDs * DforwardDs;
257	0	}
258
259	0	Real BachelierCalculator::gammaForward() const {
260		// For Bachelier model:
261		// Gamma_Forward = d²V/dF² = d/dF(N(d)) = n(d) * dd/dF = n(d)/σ
262
263	0	if (stdDev_ <= QL_EPSILON) {
264	0	return 0.0;
265	0	}
266
267	0	return discount_ * n_d_ / stdDev_;
268	0	}
269
270	0	Real BachelierCalculator::theta(Real spot, Time maturity) const {
271
272	0	QL_REQUIRE(maturity >= 0.0, "maturity (" << maturity << ") must be non-negative");
273	0	if (close(maturity, 0.0)) return 0.0;
274
275		// Theta = -dV/dt = -(rV - rSDelta + 0.5σ²*Gamma)
276
277	0	return -(std::log(discount_) * value() + std::log(forward_ / spot) * spot * delta(spot) +
278	0	0.5 * variance_ * gamma(spot)) / maturity;
279	0	}
280
281	0	Real BachelierCalculator::vega(Time maturity) const {
282	0	QL_REQUIRE(maturity >= 0.0, "negative maturity not allowed");
283
284		// For Bachelier model:
285		// Vega = dV/dσ = d/dσ[(F-K)N(d) + σn(d)]
286		// = (F-K)n(d)dd/dσ + n(d) + σn'(d)dd/dσ
287		// where d = (F-K)/σ, so dd/dσ = -(F-K)/σ² = -d/σ
288		// and n'(d) = -d*n(d)
289		// Therefore: Vega = n(d) + σ(-dn(d))(-d/σ) = n(d) + d²n(d) = n(d)(1 + d²) - (F-K)n(d)*d/σ
290		// Simplifying: Vega = n(d) for Bachelier model
291
292	0	if (maturity <= QL_EPSILON \|\| stdDev_ <= QL_EPSILON) {
293	0	return 0.0;
294	0	}
295
296	0	return discount_ * std::sqrt(maturity) * n_d_;
297	0	}
298
299	0	Real BachelierCalculator::rho(Time maturity) const {
300	0	QL_REQUIRE(maturity >= 0.0, "negative maturity not allowed");
301
302		// For Bachelier model:
303		// Rho = dV/dr = T * (discount * delta_forward * forward - value)
304		// where delta_forward = N(d) for calls, N(d)-1 for puts
305
306	0	Real deltaFwd = deltaForward();
307	0	Real rho_value = maturity * (deltaFwd * forward_ - value());
308
309	0	return rho_value;
310	0	}
311
312	0	Real BachelierCalculator::dividendRho(Time maturity) const {
313	0	QL_REQUIRE(maturity >= 0.0, "negative maturity not allowed");
314
315		// For Bachelier model:
316		// Dividend rho = -T * discount * delta_forward * forward
317		// where delta_forward = N(d) for calls, N(d)-1 for puts
318
319	0	Real deltaFwd = (alpha_ >= 0) ? cum_d_ : (cum_d_ - 1.0);
320
321	0	return -maturity * discount_ * deltaFwd * forward_;
322	0	}
323
324	0	Real BachelierCalculator::strikeSensitivity() const {
325		// For Bachelier model:
326		// dV/dK = -N(d) for calls, N(-d) for puts
327		// where d = (F-K)/σ, so dd/dK = -1/σ
328
329	0	if (alpha_ >= 0) { // Call option
330	0	return -discount_ * cum_d_; // -N(d)
331	0	} else { // Put option
332	0	return discount_ * (1.0 - cum_d_); // N(-d) = 1 - N(d)
333	0	}
334	0	}
335
336
337	0	Real BachelierCalculator::strikeGamma() const {
338		// For Bachelier model:
339		// d²V/dK² = d/dK(-N(d)) = -n(d) * dd/dK = -n(d) * (-1/σ) = n(d)/σ
340		// This is the same for both calls and puts
341
342	0	if (stdDev_ <= QL_EPSILON) {
343	0	return 0.0;
344	0	}
345
346	0	return discount_ * n_d_ / stdDev_;
347	0	}
348
349		// Sensitivity of Vega to forward (Vanna)
350	0	Real BachelierCalculator::vanna(Time maturity) const {
351	0	if (maturity <= QL_EPSILON \|\| stdDev_ <= QL_EPSILON) {
352	0	return 0.0;
353	0	}
354		// d_ is (F-K)/stdDev_
355		// n_d_ is n(d)
356		// Vanna = -d * n(d) / stdDev_ * sqrt(T)
357	0	return -d_ * n_d_ * std::sqrt(maturity) / stdDev_;
358	0	}
359
360		// Sensitivity of Vega to volatility (Volga)
361	0	Real BachelierCalculator::volga(Time maturity) const {
362	0	if (maturity <= QL_EPSILON \|\| stdDev_ <= QL_EPSILON) {
363	0	return 0.0;
364	0	}
365		// Volga = d^2 / stdDev_ * Vega
366	0	Real vega = this->vega(maturity);
367	0	return (d_ * d_ / stdDev_) * vega;
368	0	}
369		}