14910176578

Committed 08 May 2025 03:28PM UTC coverage: 73.315% (+0.02%) from 73.3%

Build # 14910176578

Build Type

Pull #2195

github

Committed by

web-flow

Commit Message

Merge 3a61f499c into 5d972fb7b

Pull Request Pull Request #2195: Added `Handle<Quote>` for spread in `OISRateHelper`

Run Details

32 of 33 new or added lines in 2 files covered. (96.97%)

277 existing lines in 25 files now uncovered.

56277 of 76761 relevant lines covered (73.31%)

8687029.35 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.35

/ql/pricingengines/vanilla/qdplusamericanengine.cpp

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

/*
 Copyright (C) 2022 Klaus Spanderen

 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
 <http://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.
*/

/*! \file qrplusamericanengine.cpp
*/

#include <algorithm>
#include <ql/exercise.hpp>
#include <ql/utilities/null.hpp>
#include <ql/math/functional.hpp>
#include <ql/math/comparison.hpp>
#include <ql/math/solvers1d/brent.hpp>
#include <ql/math/solvers1d/ridder.hpp>
#include <ql/math/solvers1d/newton.hpp>
#include <ql/math/interpolations/chebyshevinterpolation.hpp>
#include <ql/pricingengines/blackcalculator.hpp>
#include <ql/pricingengines/vanilla/qdplusamericanengine.hpp>
#include <ql/math/integrals/tanhsinhintegral.hpp>
#ifndef QL_BOOST_HAS_TANH_SINH
#include <ql/math/integrals/gausslobattointegral.hpp>
#endif

namespace QuantLib {

    class QdPlusBoundaryEvaluator {
      public:
        QdPlusBoundaryEvaluator(
            Real S, Real strike, Rate rf, Rate dy, Volatility vol, Time t, Time T)
        : tau(t), K(strike), sigma(vol), sigma2(sigma * sigma), v(sigma * std::sqrt(tau)), r(rf),
          q(dy), dr(std::exp(-r * tau)), dq(std::exp(-q * tau)),
          ddr((std::abs(r*tau) > 1e-5)? Real(r/(1-dr)) : Real(1/(tau*(1-0.5*r*tau*(1-r*tau/3))))),
          omega(2 * (r - q) / sigma2),
          lambda(0.5 *
                 (-(omega - 1) - std::sqrt(squared(omega - 1) + 8 * ddr / sigma2))),
          lambdaPrime(2 * ddr*ddr /
                      (sigma2 * std::sqrt(squared(omega - 1) + 8 * ddr / sigma2))),
          alpha(2 * dr / (sigma2 * (2 * lambda + omega - 1))),
          beta(alpha * (ddr + lambdaPrime / (2 * lambda + omega - 1)) - lambda),
          xMax(QdPlusAmericanEngine::xMax(strike, r, q)),
          xMin(QL_EPSILON * 1e4 * std::min(0.5 * (strike + S), xMax)),

          sc(Null<Real>()) {}

        Real operator()(Real S) const {
            ++nrEvaluations;

            if (S != sc)
                preCalculate(S);

            if (close_enough(K-S, npv)) {
                return (1-dq*Phi_dp)*S + alpha*theta/dr;
            }
            else {
                const Real c0 = -beta - lambda + alpha*theta/(dr*(K-S-npv));
                return (1-dq*Phi_dp)*S + (lambda+c0)*(K-S-npv);
            }
        }
        Real derivative(Real S) const {
            if (S != sc)
                preCalculate(S);

            return 1 - dq*Phi_dp + dq/v*phi_dp + beta*(1-dq*Phi_dp)
                    + alpha/dr*charm;
        }
        Real fprime2(Real S) const {
            if (S != sc)
                preCalculate(S);

            const Real gamma = phi_dp*dq/(v*S);
            const Real colour = gamma*(q + (r-q)*dp/v + (1-dp*dm)/(2*tau));

            return dq*(phi_dp/(S*v) - phi_dp*dp/(S*v*v))
                    + beta*gamma + alpha/dr*colour;
        }

        Real xmin() const { return xMin; }
        Real xmax() const { return xMax; }
        Size evaluations() const { return nrEvaluations; }

      private:
        void preCalculate(Real S) const {
            S = std::max(QL_EPSILON, S);
            sc = S;
            dp = std::log(S*dq/(K*dr))/v + 0.5*v;
            dm = dp - v;
            Phi_dp = Phi(-dp);
            Phi_dm = Phi(-dm);
            phi_dp = phi(dp);

            npv = dr*K*Phi_dm - S*dq*Phi_dp;
            theta = r*K*dr*Phi_dm - q*S*dq*Phi_dp - sigma2*S/(2*v)*dq*phi_dp;
            charm = -dq*(phi_dp*((r-q)/v - dm/(2*tau)) +q*Phi_dp);
        }

        const CumulativeNormalDistribution Phi;
        const NormalDistribution phi;
        const Time tau;
        const Real K;
        const Volatility sigma, sigma2, v;
        const Rate r, q;
        const DiscountFactor dr, dq, ddr;
        const Real omega, lambda, lambdaPrime, alpha, beta, xMax, xMin;
        mutable Size nrEvaluations = 0;
        mutable Real sc, dp, dm, Phi_dp, Phi_dm, phi_dp;
        mutable Real npv, theta, charm;
    };


    detail::QdPlusAddOnValue::QdPlusAddOnValue(
        Time T, Real S, Real K, Rate r, Rate q, Volatility vol,
        const Real xmax, ext::shared_ptr<Interpolation> q_z)
    : T_(T), S_(S), K_(K), xmax_(xmax),
      r_(r), q_(q), vol_(vol), q_z_(std::move(q_z)) {}


    Real detail::QdPlusAddOnValue::operator()(Real z) const {
        const Real t = z*z;
        const Real q = (*q_z_)(2*std::sqrt(std::max(0.0, T_-t)/T_) - 1, true);
        const Real b_t = xmax_*std::exp(-std::sqrt(std::max(0.0, q)));

        const Real dr = std::exp(-r_*t);
        const Real dq = std::exp(-q_*t);
        const Real v = vol_*std::sqrt(t);

        Real r;
        if (v >= QL_EPSILON) {
            if (b_t > QL_EPSILON) {
                const Real dp = std::log(S_*dq/(b_t*dr))/v + 0.5*v;
                r = 2*z*(r_*K_*dr*Phi_(-dp+v) - q_*S_*dq*Phi_(-dp));
            }
            else
                r = 0.0;
        }
        else if (close_enough(S_*dq, b_t*dr))
            r = z*(r_*K_*dr - q_*S_*dq);
        else if (b_t*dr > S_*dq)
            r = 2*z*(r_*K_*dr - q_*S_*dq);
        else
            r = 0.0;

        return r;
    }


    detail::QdPutCallParityEngine::QdPutCallParityEngine(
        ext::shared_ptr<GeneralizedBlackScholesProcess> process)
    : process_(std::move(process)) {
        registerWith(process_);
    }

    void detail::QdPutCallParityEngine::calculate() const {
        QL_REQUIRE(arguments_.exercise->type() == Exercise::American,
                   "not an American option");

        const auto payoff =
            ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
        QL_REQUIRE(payoff, "non-striked payoff given");

        const Real spot = process_->x0();
        QL_REQUIRE(spot >= 0.0, "negative underlying given");

        const auto maturity = arguments_.exercise->lastDate();
        const Time T = process_->time(maturity);
        const Real S = process_->x0();
        const Real K = payoff->strike();
        const Rate r = -std::log(process_->riskFreeRate()->discount(maturity))/T;
        const Rate q = -std::log(process_->dividendYield()->discount(maturity))/T;
        const Volatility vol = process_->blackVolatility()->blackVol(T, K);

        QL_REQUIRE(S >= 0, "zero or positive underlying value is required");
        QL_REQUIRE(K >= 0, "zero or positive strike is required");
        QL_REQUIRE(vol >= 0, "zero or positive volatility is required");

        if (payoff->optionType() == Option::Put)
            results_.value = calculatePutWithEdgeCases(S, K, r, q, vol, T);
        else if (payoff->optionType() == Option::Call)
            results_.value = calculatePutWithEdgeCases(K, S, q, r, vol, T);
        else
            QL_FAIL("unknown option type");
    }

    Real detail::QdPutCallParityEngine::calculatePutWithEdgeCases(
        Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {

        if (close(K, 0.0))
            return 0.0;

        if (close(S, 0.0))
            return std::max(K, K*std::exp(-r*T));

        if (r <= 0.0 && r <= q)
            return std::max(0.0,
                BlackCalculator(Option::Put, K, S*std::exp((r-q)*T),
                                vol*std::sqrt(T), std::exp(-r*T)).value());

        if (close(vol, 0.0)) {
            const auto intrinsic = [&](Real t)  -> Real {
                return std::max(0.0, K*std::exp(-r*t)-S*std::exp(-q*t));
            };
            const Real npv0 = intrinsic(0.0);
            const Real npvT = intrinsic(T);
            const Real extremT
                = close_enough(r, q)? QL_MAX_REAL : Real(std::log(r*K/(q*S))/(r-q));

            if (extremT > 0.0 && extremT < T)
                return std::max({npv0, npvT, intrinsic(extremT)});
            else
                return std::max(npv0, npvT);
        }

        return calculatePut(S, K, r, q, vol, T);
    }


    Real QdPlusAmericanEngine::xMax(Real K, Rate r, Rate q) {
        //Table 2 from Leif Andersen, Mark Lake (2021)
        //"Fast American Option Pricing: The Double-Boundary Case"

        if (r > 0.0 && q > 0.0)
            return K*std::min(1.0, r/q);
        else if (r > 0.0 && q <= 0.0)
            return K;
        else if (r == 0.0 && q < 0.0)
            return K;
        else if (r == 0.0 && q >= 0.0)
            return 0.0; // Eurpoean case
        else if (r < 0.0 && q >= 0.0)
            return 0.0; // European case
        else if (r < 0.0 && q < r)
            return K; // double boundary case
        else if (r < 0.0 && r <= q && q < 0.0)
            return 0; // European case
        else
            QL_FAIL("internal error");
    }

    QdPlusAmericanEngine::QdPlusAmericanEngine(
        ext::shared_ptr<GeneralizedBlackScholesProcess> process,
        Size interpolationPoints,
        QdPlusAmericanEngine::SolverType solverType,
        Real eps, Size maxIter)
    : detail::QdPutCallParityEngine(std::move(process)),
      interpolationPoints_(interpolationPoints),
      solverType_(solverType),
      eps_(eps),
      maxIter_((maxIter == Null<Size>()) ? (
          (solverType == Newton
           || solverType == Brent || solverType== Ridder)? 100 : 10)
          : maxIter ) { }

    template <class Solver>
    Real QdPlusAmericanEngine::buildInSolver(
        const QdPlusBoundaryEvaluator& eval,
        Solver solver, Real S, Real strike, Size maxIter, Real guess) const {

        solver.setMaxEvaluations(maxIter);
        solver.setLowerBound(eval.xmin());

        const Real fxmin = eval(eval.xmin());
        Real xmax = std::max(0.5*(eval.xmax() + S), eval.xmax());
        while (eval(xmax)*fxmin > 0.0 && eval.evaluations() < maxIter_)
            xmax*=2;

        if (guess == Null<Real>())
            guess = 0.5*(xmax + S);

        if (guess >= xmax)
            guess = std::nextafter(xmax, Real(-1));
        else if (guess <= eval.xmin())
            guess = std::nextafter(eval.xmin(), QL_MAX_REAL);

        return solver.solve(eval, eps_, guess, eval.xmin(), xmax);
    }

    std::pair<Size, Real> QdPlusAmericanEngine::putExerciseBoundaryAtTau(
            Real S, Real K, Rate r, Rate q,
            Volatility vol, Time T, Time tau) const {

        if (tau < QL_EPSILON)
            return std::pair<Size, Real>(
                Size(0), xMax(K, r, q));

        const QdPlusBoundaryEvaluator eval(S, K, r, q, vol, tau, T);

        Real x;
        switch (solverType_) {
          case Brent:
            x = buildInSolver(eval, QuantLib::Brent(), S, K, maxIter_);
            break;
          case Newton:
            x = buildInSolver(eval, QuantLib::Newton(), S, K, maxIter_);
            break;
          case Ridder:
            x = buildInSolver(eval, QuantLib::Ridder(), S, K, maxIter_);
            break;
          case Halley:
          case SuperHalley:
            {
                bool resultCloseEnough;
                x = eval.xmax();
                Real xOld, fx;
                const Real xmin = eval.xmin();

                do {
                    xOld = x;
                    fx = eval(x);
                    const Real fPrime = eval.derivative(x);
                    const Real lf = fx*eval.fprime2(x)/(fPrime*fPrime);
                    const Real step = (solverType_ == Halley)
                        ? Real(1/(1 - 0.5*lf)*fx/fPrime)
                        : Real((1 + 0.5*lf/(1-lf))*fx/fPrime);

                    x = std::max(xmin, x - step);
                    resultCloseEnough = std::fabs(x-xOld) < 0.5*eps_;
                }
                while (!resultCloseEnough && eval.evaluations() < maxIter_);

                if (!resultCloseEnough && !close(std::fabs(fx), 0.0)) {
                    x = buildInSolver(eval, QuantLib::Brent(), S, K, 10*maxIter_, x);
                }
            }
            break;
          default:
            QL_FAIL("unknown solver type");
        }

        return std::pair<Size, Real>(eval.evaluations(), x);
    }

    ext::shared_ptr<ChebyshevInterpolation>
        QdPlusAmericanEngine::getPutExerciseBoundary(
        Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {

        const Real xmax = xMax(K, r, q);

        return ext::make_shared<ChebyshevInterpolation>(
            interpolationPoints_,
            [&, this](Real z) {
                const Real x_sq = 0.25*T*squared(1+z);
                return squared(std::log(
                    this->putExerciseBoundaryAtTau(S, K, r, q, vol, T, x_sq)
                        .second/xmax));
            },
            ChebyshevInterpolation::SecondKind
        );
    }

    Real QdPlusAmericanEngine::calculatePut(
        Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {

        if (r < 0.0 && q < r)
            QL_FAIL("double-boundary case q<r<0 for a put option is given");

        const ext::shared_ptr<Interpolation> q_z
            = getPutExerciseBoundary(S, K, r, q, vol, T);

        const Real xmax = xMax(K, r, q);
        const detail::QdPlusAddOnValue aov(T, S, K, r, q, vol, xmax, q_z);

#ifdef QL_BOOST_HAS_TANH_SINH
        const Real addOn = TanhSinhIntegral(eps_)(aov, 0.0, std::sqrt(T));
#else
        const Real addOn = GaussLobattoIntegral(100000, QL_MAX_REAL, 0.1*eps_)
                (aov, 0.0, std::sqrt(T));
#endif

        QL_REQUIRE(addOn > -10*eps_,
            "negative early exercise value " << addOn);

        const Real europeanValue = std::max(
            0.0,
            BlackCalculator(
                Option::Put, K,
                S*std::exp((r-q)*T),
                vol*std::sqrt(T), std::exp(-r*T)).value()
        );

        return europeanValue + std::max(0.0, addOn);
    }
}

1	/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3	/*
4	Copyright (C) 2022 Klaus Spanderen
5
6	This file is part of QuantLib, a free-software/open-source library
7	for financial quantitative analysts and developers - http://quantlib.org/
8
9	QuantLib is free software: you can redistribute it and/or modify it
10	under the terms of the QuantLib license. You should have received a
11	copy of the license along with this program; if not, please email
12	<quantlib-dev@lists.sf.net>. The license is also available online at
13	<http://quantlib.org/license.shtml>.
14
15	This program is distributed in the hope that it will be useful, but WITHOUT
16	ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17	FOR A PARTICULAR PURPOSE. See the license for more details.
18	*/
19
20	/*! \file qrplusamericanengine.cpp
21	*/
22
23	#include <algorithm>
24	#include <ql/exercise.hpp>
25	#include <ql/utilities/null.hpp>
26	#include <ql/math/functional.hpp>
27	#include <ql/math/comparison.hpp>
28	#include <ql/math/solvers1d/brent.hpp>
29	#include <ql/math/solvers1d/ridder.hpp>
30	#include <ql/math/solvers1d/newton.hpp>
31	#include <ql/math/interpolations/chebyshevinterpolation.hpp>
32	#include <ql/pricingengines/blackcalculator.hpp>
33	#include <ql/pricingengines/vanilla/qdplusamericanengine.hpp>
34	#include <ql/math/integrals/tanhsinhintegral.hpp>
35	#ifndef QL_BOOST_HAS_TANH_SINH
36	#include <ql/math/integrals/gausslobattointegral.hpp>
37	#endif
38
39	namespace QuantLib {
40
41	class QdPlusBoundaryEvaluator {
42	public:
43	QdPlusBoundaryEvaluator(	27,664✔
44	Real S, Real strike, Rate rf, Rate dy, Volatility vol, Time t, Time T)
45	: tau(t), K(strike), sigma(vol), sigma2(sigma * sigma), v(sigma * std::sqrt(tau)), r(rf),	55,328✔
46	q(dy), dr(std::exp(-r * tau)), dq(std::exp(-q * tau)),	27,664✔
47	ddr((std::abs(rtau) > 1e-5)? Real(r/(1-dr)) : Real(1/(tau(1-0.5rtau(1-rtau/3))))),	27,664✔
48	omega(2 * (r - q) / sigma2),	27,664✔
49	lambda(0.5 *	27,664✔
50	(-(omega - 1) - std::sqrt(squared(omega - 1) + 8 * ddr / sigma2))),	27,664✔
51	lambdaPrime(2 * ddr*ddr /	27,664✔
52	(sigma2 * std::sqrt(squared(omega - 1) + 8 * ddr / sigma2))),	27,664✔
53	alpha(2 * dr / (sigma2 * (2 * lambda + omega - 1))),	27,664✔
54	beta(alpha * (ddr + lambdaPrime / (2 * lambda + omega - 1)) - lambda),	27,664✔
55	xMax(QdPlusAmericanEngine::xMax(strike, r, q)),	27,664✔
56	xMin(QL_EPSILON * 1e4 * std::min(0.5 * (strike + S), xMax)),	27,664✔
57
58	sc(Null<Real>()) {}	27,664✔
59
60	Real operator()(Real S) const {	130,227✔
61	++nrEvaluations;	130,227✔
62
63	if (S != sc)	130,227✔
64	preCalculate(S);	129,957✔
65
66	if (close_enough(K-S, npv)) {	130,227✔
67	return (1-dqPhi_dp)S + alpha*theta/dr;	22✔
68	}
69	else {
70	const Real c0 = -beta - lambda + alphatheta/(dr(K-S-npv));	130,205✔
71	return (1-dqPhi_dp)S + (lambda+c0)*(K-S-npv);	130,205✔
72	}
73	}
74	Real derivative(Real S) const {	106,784✔
75	if (S != sc)	106,784✔
UNCOV 76	preCalculate(S);	×
77
78	return 1 - dqPhi_dp + dq/vphi_dp + beta(1-dqPhi_dp)	106,784✔
79	+ alpha/dr*charm;	106,784✔
80	}
81	Real fprime2(Real S) const {	103,392✔
82	if (S != sc)	103,392✔
UNCOV 83	preCalculate(S);	×
84
85	const Real gamma = phi_dpdq/(vS);	103,392✔
86	const Real colour = gamma(q + (r-q)dp/v + (1-dpdm)/(2tau));	103,392✔
87
88	return dq(phi_dp/(Sv) - phi_dpdp/(Sv*v))	103,392✔
89	+ betagamma + alpha/drcolour;	103,392✔
90	}
91
92	Real xmin() const { return xMin; }	2,143✔
93	Real xmax() const { return xMax; }	27,707✔
94	Size evaluations() const { return nrEvaluations; }	105,535✔
95
96	private:
97	void preCalculate(Real S) const {	129,957✔
98	S = std::max(QL_EPSILON, S);	129,957✔
99	sc = S;	129,957✔
100	dp = std::log(Sdq/(Kdr))/v + 0.5*v;	129,957✔
101	dm = dp - v;	129,957✔
102	Phi_dp = Phi(-dp);	129,957✔
103	Phi_dm = Phi(-dm);	129,957✔
104	phi_dp = phi(dp);	129,957✔
105
106	npv = drKPhi_dm - SdqPhi_dp;	129,957✔
107	theta = rKdrPhi_dm - qSdqPhi_dp - sigma2S/(2v)dqphi_dp;	129,957✔
108	charm = -dq(phi_dp((r-q)/v - dm/(2tau)) +qPhi_dp);	129,957✔
109	}	129,957✔
110
111	const CumulativeNormalDistribution Phi;
112	const NormalDistribution phi;
113	const Time tau;
114	const Real K;
115	const Volatility sigma, sigma2, v;
116	const Rate r, q;
117	const DiscountFactor dr, dq, ddr;
118	const Real omega, lambda, lambdaPrime, alpha, beta, xMax, xMin;
119	mutable Size nrEvaluations = 0;
120	mutable Real sc, dp, dm, Phi_dp, Phi_dm, phi_dp;
121	mutable Real npv, theta, charm;
122	};
123
124
125	detail::QdPlusAddOnValue::QdPlusAddOnValue(	2,732✔
126	Time T, Real S, Real K, Rate r, Rate q, Volatility vol,
127	const Real xmax, ext::shared_ptr<Interpolation> q_z)	2,732✔
128	: T_(T), S_(S), K_(K), xmax_(xmax),	2,732✔
129	r_(r), q_(q), vol_(vol), q_z_(std::move(q_z)) {}	2,732✔
130
131
132	Real detail::QdPlusAddOnValue::operator()(Real z) const {	1,591,269✔
133	const Real t = z*z;	1,591,269✔
134	const Real q = (q_z_)(2std::sqrt(std::max(0.0, T_-t)/T_) - 1, true);	3,181,078✔
135	const Real b_t = xmax_*std::exp(-std::sqrt(std::max(0.0, q)));	1,591,269✔
136
137	const Real dr = std::exp(-r_*t);	1,591,269✔
138	const Real dq = std::exp(-q_*t);	1,591,269✔
139	const Real v = vol_*std::sqrt(t);	1,591,269✔
140
141	Real r;
142	if (v >= QL_EPSILON) {	1,591,269✔
143	if (b_t > QL_EPSILON) {	1,219,489✔
144	const Real dp = std::log(S_dq/(b_tdr))/v + 0.5*v;	1,219,489✔
145	r = 2z(r_K_drPhi_(-dp+v) - q_S_dqPhi_(-dp));	1,219,489✔
146	}
147	else
148	r = 0.0;
149	}
150	else if (close_enough(S_dq, b_tdr))	371,780✔
UNCOV 151	r = z(r_K_dr - q_S_*dq);	×
152	else if (b_tdr > S_dq)	371,780✔
153	r = 2z(r_K_dr - q_S_dq);	15,218✔
154	else
155	r = 0.0;
156
157	return r;	1,591,269✔
158	}
159
160
161	detail::QdPutCallParityEngine::QdPutCallParityEngine(	2,287✔
162	ext::shared_ptr<GeneralizedBlackScholesProcess> process)	2,287✔
163	: process_(std::move(process)) {	2,287✔
164	registerWith(process_);	4,574✔
165	}	2,287✔
166
167	void detail::QdPutCallParityEngine::calculate() const {	2,750✔
168	QL_REQUIRE(arguments_.exercise->type() == Exercise::American,	2,750✔
169	"not an American option");
170
171	const auto payoff =
172	ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);	2,750✔
173	QL_REQUIRE(payoff, "non-striked payoff given");	2,750✔
174
175	const Real spot = process_->x0();	2,750✔
176	QL_REQUIRE(spot >= 0.0, "negative underlying given");	2,750✔
177
178	const auto maturity = arguments_.exercise->lastDate();	2,750✔
179	const Time T = process_->time(maturity);	2,750✔
180	const Real S = process_->x0();	2,750✔
181	const Real K = payoff->strike();	2,750✔
182	const Rate r = -std::log(process_->riskFreeRate()->discount(maturity))/T;	2,750✔
183	const Rate q = -std::log(process_->dividendYield()->discount(maturity))/T;	2,750✔
184	const Volatility vol = process_->blackVolatility()->blackVol(T, K);	2,750✔
185
186	QL_REQUIRE(S >= 0, "zero or positive underlying value is required");	2,750✔
187	QL_REQUIRE(K >= 0, "zero or positive strike is required");	2,750✔
188	QL_REQUIRE(vol >= 0, "zero or positive volatility is required");	2,750✔
189
190	if (payoff->optionType() == Option::Put)	2,750✔
191	results_.value = calculatePutWithEdgeCases(S, K, r, q, vol, T);	2,482✔
192	else if (payoff->optionType() == Option::Call)	268✔
193	results_.value = calculatePutWithEdgeCases(K, S, q, r, vol, T);	268✔
194	else
UNCOV 195	QL_FAIL("unknown option type");	×
196	}	2,750✔
197
198	Real detail::QdPutCallParityEngine::calculatePutWithEdgeCases(	2,750✔
199	Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {
200
201	if (close(K, 0.0))	2,750✔
202	return 0.0;
203
204	if (close(S, 0.0))	2,746✔
205	return std::max(K, Kstd::exp(-rT));	5✔
206
207	if (r <= 0.0 && r <= q)	2,742✔
208	return std::max(0.0,	12✔
209	BlackCalculator(Option::Put, K, Sstd::exp((r-q)T),	6✔
210	volstd::sqrt(T), std::exp(-rT)).value());	12✔
211
212	if (close(vol, 0.0)) {	2,736✔
213	const auto intrinsic = [&](Real t) -> Real {	9✔
214	return std::max(0.0, Kstd::exp(-rt)-Sstd::exp(-qt));	9✔
215	};	4✔
216	const Real npv0 = intrinsic(0.0);	4✔
217	const Real npvT = intrinsic(T);	4✔
218	const Real extremT
219	= close_enough(r, q)? QL_MAX_REAL : Real(std::log(rK/(qS))/(r-q));	4✔
220
221	if (extremT > 0.0 && extremT < T)	4✔
222	return std::max({npv0, npvT, intrinsic(extremT)});	1✔
223	else
224	return std::max(npv0, npvT);	3✔
225	}
226
227	return calculatePut(S, K, r, q, vol, T);	2,732✔
228	}
229
230
231	Real QdPlusAmericanEngine::xMax(Real K, Rate r, Rate q) {	35,895✔
232	//Table 2 from Leif Andersen, Mark Lake (2021)
233	//"Fast American Option Pricing: The Double-Boundary Case"
234
235	if (r > 0.0 && q > 0.0)	35,895✔
236	return K*std::min(1.0, r/q);	43,841✔
237	else if (r > 0.0 && q <= 0.0)	8,464✔
238	return K;
239	else if (r == 0.0 && q < 0.0)	78✔
240	return K;
UNCOV 241	else if (r == 0.0 && q >= 0.0)	×
242	return 0.0; // Eurpoean case
UNCOV 243	else if (r < 0.0 && q >= 0.0)	×
244	return 0.0; // European case
UNCOV 245	else if (r < 0.0 && q < r)	×
246	return K; // double boundary case
UNCOV 247	else if (r < 0.0 && r <= q && q < 0.0)	×
248	return 0; // European case
249	else
UNCOV 250	QL_FAIL("internal error");	×
251	}
252
253	QdPlusAmericanEngine::QdPlusAmericanEngine(	2,256✔
254	ext::shared_ptr<GeneralizedBlackScholesProcess> process,
255	Size interpolationPoints,
256	QdPlusAmericanEngine::SolverType solverType,
257	Real eps, Size maxIter)	2,256✔
258	: detail::QdPutCallParityEngine(std::move(process)),
259	interpolationPoints_(interpolationPoints),	2,256✔
260	solverType_(solverType),	2,256✔
261	eps_(eps),	2,256✔
262	maxIter_((maxIter == Null<Size>()) ? (	4,512✔
263	(solverType == Newton
264	\|\| solverType == Brent \|\| solverType== Ridder)? 100 : 10)	2,256✔
265	: maxIter ) { }	6,768✔
266
267	template <class Solver>
268	Real QdPlusAmericanEngine::buildInSolver(	2,143✔
269	const QdPlusBoundaryEvaluator& eval,
270	Solver solver, Real S, Real strike, Size maxIter, Real guess) const {
271
272	solver.setMaxEvaluations(maxIter);
273	solver.setLowerBound(eval.xmin());
274
275	const Real fxmin = eval(eval.xmin());	2,143✔
276	Real xmax = std::max(0.5*(eval.xmax() + S), eval.xmax());	2,143✔
277	while (eval(xmax)*fxmin > 0.0 && eval.evaluations() < maxIter_)	2,143✔
UNCOV 278	xmax*=2;	×
279
280	if (guess == Null<Real>())	2,143✔
281	guess = 0.5*(xmax + S);	2,100✔
282
283	if (guess >= xmax)	2,143✔
284	guess = std::nextafter(xmax, Real(-1));	1,200✔
285	else if (guess <= eval.xmin())	943✔
UNCOV 286	guess = std::nextafter(eval.xmin(), QL_MAX_REAL);	×
287
288	return solver.solve(eval, eps_, guess, eval.xmin(), xmax);	2,143✔
289	}
290
291	std::pair<Size, Real> QdPlusAmericanEngine::putExerciseBoundaryAtTau(	30,431✔
292	Real S, Real K, Rate r, Rate q,
293	Volatility vol, Time T, Time tau) const {
294
295	if (tau < QL_EPSILON)	30,431✔
296	return std::pair<Size, Real>(
297	Size(0), xMax(K, r, q));	2,767✔
298
299	const QdPlusBoundaryEvaluator eval(S, K, r, q, vol, tau, T);	27,664✔
300
301	Real x;
302	switch (solverType_) {	27,664✔
303	case Brent:
304	x = buildInSolver(eval, QuantLib::Brent(), S, K, maxIter_);	700✔
305	break;	700✔
306	case Newton:
307	x = buildInSolver(eval, QuantLib::Newton(), S, K, maxIter_);	700✔
308	break;	700✔
309	case Ridder:
310	x = buildInSolver(eval, QuantLib::Ridder(), S, K, maxIter_);	700✔
311	break;	700✔
312	case Halley:
313	case SuperHalley:
314	{
315	bool resultCloseEnough;
316	x = eval.xmax();
317	Real xOld, fx;
318	const Real xmin = eval.xmin();	25,564✔
319
320	do {
321	xOld = x;
322	fx = eval(x);	103,392✔
323	const Real fPrime = eval.derivative(x);	103,392✔
324	const Real lf = fxeval.fprime2(x)/(fPrimefPrime);	103,392✔
325	const Real step = (solverType_ == Halley)	103,392✔
326	? Real(1/(1 - 0.5lf)fx/fPrime)	103,392✔
327	: Real((1 + 0.5lf/(1-lf))fx/fPrime);	2,371✔
328
329	x = std::max(xmin, x - step);	103,392✔
330	resultCloseEnough = std::fabs(x-xOld) < 0.5*eps_;	103,392✔
331	}
332	while (!resultCloseEnough && eval.evaluations() < maxIter_);	103,392✔
333
334	if (!resultCloseEnough && !close(std::fabs(fx), 0.0)) {	25,564✔
335	x = buildInSolver(eval, QuantLib::Brent(), S, K, 10*maxIter_, x);	43✔
336	}
337	}
338	break;	25,564✔
339	default:	×
UNCOV 340	QL_FAIL("unknown solver type");	×
341	}
342
343	return std::pair<Size, Real>(eval.evaluations(), x);	27,664✔
344	}
345
346	ext::shared_ptr<ChebyshevInterpolation>
347	QdPlusAmericanEngine::getPutExerciseBoundary(	2,732✔
348	Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {
349
350	const Real xmax = xMax(K, r, q);	2,732✔
351
352	return ext::make_shared<ChebyshevInterpolation>(
353	interpolationPoints_,	2,732✔
354	[&, this](Real z) {	26,891✔
355	const Real x_sq = 0.25Tsquared(1+z);	26,891✔
356	return squared(std::log(	26,891✔
357	this->putExerciseBoundaryAtTau(S, K, r, q, vol, T, x_sq)	26,891✔
358	.second/xmax));	26,891✔
359	},
360	ChebyshevInterpolation::SecondKind	5,464✔
361	);	2,732✔
362	}
363
364	Real QdPlusAmericanEngine::calculatePut(	514✔
365	Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {
366
367	if (r < 0.0 && q < r)	514✔
UNCOV 368	QL_FAIL("double-boundary case q<r<0 for a put option is given");	×
369
370	const ext::shared_ptr<Interpolation> q_z
371	= getPutExerciseBoundary(S, K, r, q, vol, T);	514✔
372
373	const Real xmax = xMax(K, r, q);	514✔
374	const detail::QdPlusAddOnValue aov(T, S, K, r, q, vol, xmax, q_z);	514✔
375
376	#ifdef QL_BOOST_HAS_TANH_SINH
377	const Real addOn = TanhSinhIntegral(eps_)(aov, 0.0, std::sqrt(T));	514✔
378	#else
379	const Real addOn = GaussLobattoIntegral(100000, QL_MAX_REAL, 0.1*eps_)
380	(aov, 0.0, std::sqrt(T));
381	#endif
382
383	QL_REQUIRE(addOn > -10*eps_,	514✔
384	"negative early exercise value " << addOn);
385
386	const Real europeanValue = std::max(
387	0.0,	1,028✔
UNCOV 388	BlackCalculator(	×
389	Option::Put, K,
390	Sstd::exp((r-q)T),	514✔
391	volstd::sqrt(T), std::exp(-rT)).value()	514✔
392	);	514✔
393
394	return europeanValue + std::max(0.0, addOn);	1,519✔
395	}
396	}

lballabio / QuantLib / 14910176578

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous