/src/quantlib/ql/pricingengines/vanilla/qdfpamericanengine.cpp

Source (jump to first uncovered line)
/* -*- 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 qrfpamericanengine.cpp
*/

#include <ql/math/distributions/normaldistribution.hpp>
#include <ql/math/functional.hpp>
#include <ql/math/integrals/gaussianquadratures.hpp>
#include <ql/math/integrals/tanhsinhintegral.hpp>
#include <ql/math/interpolations/chebyshevinterpolation.hpp>
#include <ql/pricingengines/blackcalculator.hpp>
#include <ql/pricingengines/vanilla/qdfpamericanengine.hpp>
#include <utility>
#ifndef QL_BOOST_HAS_TANH_SINH
#    include <ql/math/integrals/gausslobattointegral.hpp>
#endif

namespace QuantLib {

    QdFpLegendreScheme::QdFpLegendreScheme(
        Size l, Size m, Size n, Size p):
        m_(m), n_(n),
        fpIntegrator_(ext::make_shared<GaussLegendreIntegrator>(l)),
        exerciseBoundaryIntegrator_(
            ext::make_shared<GaussLegendreIntegrator>(p)) {

        QL_REQUIRE(m_ > 0, "at least one fixed point iteration step is needed");
        QL_REQUIRE(n_ > 0, "at least one interpolation point is needed");
    }

    Size QdFpLegendreScheme::getNumberOfChebyshevInterpolationNodes()
        const {
        return n_;
    }
    Size QdFpLegendreScheme::getNumberOfNaiveFixedPointSteps() const {
        return m_-1;
    }
    Size QdFpLegendreScheme::getNumberOfJacobiNewtonFixedPointSteps()
        const {
        return Size(1);
    }

    ext::shared_ptr<Integrator>
        QdFpLegendreScheme::getFixedPointIntegrator() const {
        return fpIntegrator_;
    }
    ext::shared_ptr<Integrator>
        QdFpLegendreScheme::getExerciseBoundaryToPriceIntegrator()
        const {
        return exerciseBoundaryIntegrator_;
    }

    QdFpTanhSinhIterationScheme::QdFpTanhSinhIterationScheme(
        Size m, Size n, Real eps)
    : m_(m), n_(n),
#ifdef QL_BOOST_HAS_TANH_SINH
      integrator_(ext::make_shared<TanhSinhIntegral>(eps))
#else
      integrator_(ext::make_shared<GaussLobattoIntegral>(
          100000, QL_MAX_REAL, 0.1*eps))
#endif
    {}

    Size QdFpTanhSinhIterationScheme::getNumberOfChebyshevInterpolationNodes()
        const {
        return n_;
    }
    Size QdFpTanhSinhIterationScheme::getNumberOfNaiveFixedPointSteps() const {
        return m_-1;
    }
    Size QdFpTanhSinhIterationScheme::getNumberOfJacobiNewtonFixedPointSteps()
        const {
        return Size(1);
    }

    ext::shared_ptr<Integrator>
    QdFpTanhSinhIterationScheme::getFixedPointIntegrator() const {
        return integrator_;
    }
    ext::shared_ptr<Integrator>
    QdFpTanhSinhIterationScheme::getExerciseBoundaryToPriceIntegrator()
        const {
        return integrator_;
    }

    QdFpLegendreTanhSinhScheme::QdFpLegendreTanhSinhScheme(
        Size l, Size m, Size n, Real eps)
    : QdFpLegendreScheme(l, m, n, 1),
      eps_(eps) {}

    ext::shared_ptr<Integrator>
    QdFpLegendreTanhSinhScheme::getExerciseBoundaryToPriceIntegrator() const {
#ifdef QL_BOOST_HAS_TANH_SINH
            return ext::make_shared<TanhSinhIntegral>(eps_);
#else
            return ext::make_shared<GaussLobattoIntegral>(
                100000, QL_MAX_REAL, 0.1*eps_);
#endif
    }


    class DqFpEquation {
      public:
        DqFpEquation(Rate _r,
                     Rate _q,
                     Volatility _vol,
                     std::function<Real(Real)> B,
                     ext::shared_ptr<Integrator> _integrator)
        : r(_r), q(_q), vol(_vol), B(std::move(B)), integrator(std::move(_integrator)) {
            const auto legendreIntegrator =
                ext::dynamic_pointer_cast<GaussLegendreIntegrator>(integrator);

            if (legendreIntegrator != nullptr) {
                x_i = legendreIntegrator->getIntegration()->x();
                w_i = legendreIntegrator->getIntegration()->weights();
            }
        }

        virtual std::pair<Real, Real> NDd(Real tau, Real b) const = 0;
        virtual std::tuple<Real, Real, Real> f(Real tau, Real b) const = 0;

        virtual ~DqFpEquation() = default;

      protected:
        std::pair<Real, Real> d(Time t, Real z) const {
            const Real v = vol * std::sqrt(t);
            const Real m = (std::log(z) + (r-q)*t)/v + 0.5*v;

            return std::make_pair(m, m-v);
        }

        Array x_i, w_i;
        const Rate r, q;
        const Volatility vol;

        const std::function<Real(Real)> B;
        const ext::shared_ptr<Integrator> integrator;

        const NormalDistribution phi;
        const CumulativeNormalDistribution Phi;
    };

    class DqFpEquation_B: public DqFpEquation {
      public:
        DqFpEquation_B(Real K,
                       Rate _r,
                       Rate _q,
                       Volatility _vol,
                       std::function<Real(Real)> B,
                       ext::shared_ptr<Integrator> _integrator);

        std::pair<Real, Real> NDd(Real tau, Real b) const override;
        std::tuple<Real, Real, Real> f(Real tau, Real b) const override;

      private:
          const Real K;
    };

    class DqFpEquation_A: public DqFpEquation {
      public:
        DqFpEquation_A(Real K,
                       Rate _r,
                       Rate _q,
                       Volatility _vol,
                       std::function<Real(Real)> B,
                       ext::shared_ptr<Integrator> _integrator);

        std::pair<Real, Real> NDd(Real tau, Real b) const override;
        std::tuple<Real, Real, Real> f(Real tau, Real b) const override;

      private:
          const Real K;
    };

    DqFpEquation_A::DqFpEquation_A(Real K,
                                   Rate _r,
                                   Rate _q,
                                   Volatility _vol,
                                   std::function<Real(Real)> B,
                                   ext::shared_ptr<Integrator> _integrator)
    : DqFpEquation(_r, _q, _vol, std::move(B), std::move(_integrator)), K(K) {}

    std::tuple<Real, Real, Real> DqFpEquation_A::f(Real tau, Real b) const {
        const Real v = vol * std::sqrt(tau);

        Real N, D;
        if (tau < squared(QL_EPSILON)) {
            if (close_enough(b, K)) {
                N = 1/(M_SQRT2*M_SQRTPI * v);
                D = N + 0.5;
            }
            else {
                N = 0.0;
                D = (b > K)? 1.0: 0.0;
            }
        }
        else {
            const Real stv = std::sqrt(tau)/vol;

            Real K12, K3;
            if (!x_i.empty()) {
                K12 = K3 = 0.0;

                for (Integer i = x_i.size()-1; i >= 0; --i) {
                    const Real y = x_i[i];
                    const Real m = 0.25*tau*squared(1+y);
                    const std::pair<Real, Real> dpm = d(m, b/B(tau-m));

                    K12 += w_i[i] * std::exp(q*tau - q*m)
                        *(0.5*tau*(y+1)*Phi(dpm.first) + stv*phi(dpm.first));
                    K3 += w_i[i] * stv*std::exp(r*tau-r*m)*phi(dpm.second);
                }
            } else {
                K12 = (*integrator)([&, this](Real y) -> Real {
                    const Real m = 0.25*tau*squared(1+y);
                    const Real df = std::exp(q*tau - q*m);

                    if (y <= 5*QL_EPSILON - 1) {
                        if (close_enough(b, B(tau-m)))
                            return df*stv/(M_SQRT2*M_SQRTPI);
                        else
                            return 0.0;
                    }
                    else {
                        const Real dp = d(m, b/B(tau-m)).first;
                        return df*(0.5*tau*(y+1)*Phi(dp) + stv*phi(dp));
                    }
                }, -1, 1);

                K3 = (*integrator)([&, this](Real y) -> Real {
                    const Real m = 0.25*tau*squared(1+y);
                    const Real df = std::exp(r*tau-r*m);

                    if (y <= 5*QL_EPSILON - 1) {
                        if (close_enough(b, B(tau-m)))
                            return df*stv/(M_SQRT2*M_SQRTPI);
                        else
                            return 0.0;
                    }
                    else
                        return df*stv*phi(d(m, b/B(tau-m)).second);
                }, -1, 1);
            }
            const std::pair<Real, Real> dpm = d(tau, b/K);
            N = phi(dpm.second)/v + r*K3;
            D = phi(dpm.first)/v + Phi(dpm.first) + q*K12;
        }

        const Real alpha = K*std::exp(-(r-q)*tau);
        Real fv;
        if (tau < squared(QL_EPSILON)) {
            if (close_enough(b, K))
                fv = alpha;
            else if (b > K)
                fv = 0.0;
            else {
                if (close_enough(q, Real(0)))
                    fv = alpha*r*((q < 0)? -1.0 : 1.0)/QL_EPSILON;
                else
                    fv = alpha*r/q;
            }
        }
        else
            fv = alpha*N/D;

        return std::make_tuple(N, D, fv);
    }

    std::pair<Real, Real> DqFpEquation_A::NDd(Real tau, Real b) const {
        Real Dd, Nd;

        if (tau < squared(QL_EPSILON)) {
            if (close_enough(b, K)) {
                const Real sqTau = std::sqrt(tau);
                const Real vol2 = vol*vol;
                Dd = M_1_SQRTPI*M_SQRT_2*(
                    -(0.5*vol2 + r-q) / (b*vol*vol2*sqTau) + 1 / (b*vol*sqTau));
                Nd = M_1_SQRTPI*M_SQRT_2 * (-0.5*vol2 + r-q)  / (b*vol*vol2*sqTau);
            }
            else
                Dd = Nd = 0.0;
        }
        else {
            const std::pair<Real, Real> dpm = d(tau, b/K);

            Dd = -phi(dpm.first) * dpm.first / (b*vol*vol*tau) +
                    phi(dpm.first) / (b*vol * std::sqrt(tau));
            Nd = -phi(dpm.second) * dpm.second / (b*vol*vol*tau);
        }

        return std::make_pair(Nd, Dd);
    }


    DqFpEquation_B::DqFpEquation_B(Real K,
                                   Rate _r,
                                   Rate _q,
                                   Volatility _vol,
                                   std::function<Real(Real)> B,
                                   ext::shared_ptr<Integrator> _integrator)
    : DqFpEquation(_r, _q, _vol, std::move(B), std::move(_integrator)), K(K) {}


    std::tuple<Real, Real, Real> DqFpEquation_B::f(Real tau, Real b) const {
        Real N, D;
        if (tau < squared(QL_EPSILON)) {
            if (close_enough(b, K))
                N = D = 0.5;
            else if (b < K)
                N = D = 0.0;
            else
                N = D = 1.0;
        }
        else {
            Real ni, di;
            if (!x_i.empty()) {
                const Real c = 0.5*tau;

                ni = di = 0.0;
                for (Integer i = x_i.size()-1; i >= 0; --i) {
                    const Real u = c*x_i[i] + c;
                    const std::pair<Real, Real> dpm = d(tau - u, b/B(u));
                    ni += w_i[i] * std::exp(r*u)*Phi(dpm.second);
                    di += w_i[i] * std::exp(q*u)*Phi(dpm.first);
                }
                ni *= c;
                di *= c;
            } else {
                ni = (*integrator)([&, this](Real u) -> Real {
                  const Real df = std::exp(r*u);
                    if (u >= tau*(1 - 5*QL_EPSILON)) {
                        if (close_enough(b, B(u)))
                            return 0.5*df;
                        else
                            return df*((b < B(u)? 0.0: 1.0));
                    }
                    else
                        return df*Phi(d(tau - u, b/B(u)).second);
                }, 0, tau);
                di = (*integrator)([&, this](Real u) -> Real {
                  const Real df = std::exp(q*u);
                    if (u >= tau*(1 - 5*QL_EPSILON)) {
                        if (close_enough(b, B(u)))
                            return 0.5*df;
                        else
                            return df*((b < B(u)? 0.0: 1.0));
                    }
                    else
                        return df*Phi(d(tau - u, b/B(u)).first);
                    }, 0, tau);
            }

            const std::pair<Real, Real> dpm = d(tau, b/K);

            N = Phi(dpm.second) + r*ni;
            D = Phi(dpm.first) + q*di;
        }

        Real fv;
        const Real alpha = K*std::exp(-(r-q)*tau);
        if (tau < squared(QL_EPSILON)) {
            if (close_enough(b, K) || b > K)
                fv = alpha;
            else {
                if (close_enough(q, Real(0)))
                    fv = alpha*r*((q < 0)? -1.0 : 1.0)/QL_EPSILON;
                else
                    fv = alpha*r/q;
            }
        }
        else
            fv = alpha*N/D;

        return std::make_tuple(N, D, fv);
    }

    std::pair<Real, Real> DqFpEquation_B::NDd(Real tau, Real b) const {
        const std::pair<Real, Real> dpm = d(tau, b/K);
        return std::make_pair(
            phi(dpm.second) / (b*vol*std::sqrt(tau)),
            phi(dpm.first)  / (b*vol*std::sqrt(tau))
        );
    }

    QdFpAmericanEngine::QdFpAmericanEngine(
        ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess,
        ext::shared_ptr<QdFpIterationScheme> iterationScheme,
        FixedPointEquation fpEquation)
    : detail::QdPutCallParityEngine(std::move(bsProcess)),
      iterationScheme_(std::move(iterationScheme)),
      fpEquation_(fpEquation) {
    }

    ext::shared_ptr<QdFpIterationScheme>
    QdFpAmericanEngine::fastScheme() {
        static auto scheme = ext::make_shared<QdFpLegendreScheme>(7, 2, 7, 27);
        return scheme;
    }

    ext::shared_ptr<QdFpIterationScheme>
    QdFpAmericanEngine::accurateScheme() {
        static auto scheme = ext::make_shared<QdFpLegendreTanhSinhScheme>(25, 5, 13, 1e-8);
        return scheme;
    }

    ext::shared_ptr<QdFpIterationScheme>
    QdFpAmericanEngine::highPrecisionScheme() {
        static auto scheme = ext::make_shared<QdFpTanhSinhIterationScheme>(10, 30, 1e-10);
        return scheme;
    }

    Real QdFpAmericanEngine::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 Real xmax =  QdPlusAmericanEngine::xMax(K, r, q);
        const Size n = iterationScheme_->getNumberOfChebyshevInterpolationNodes();

        const ext::shared_ptr<ChebyshevInterpolation> interp =
            QdPlusAmericanEngine(
                    process_, n+1, QdPlusAmericanEngine::Halley, 1e-8)
                .getPutExerciseBoundary(S, K, r, q, vol, T);

        const Array z = interp->nodes();
        const Array x = 0.5*std::sqrt(T)*(1.0+z);

        const auto B = [xmax, T, &interp](Real tau) -> Real {
            const Real z = 2*std::sqrt(std::abs(tau)/T)-1;
            return xmax*std::exp(-std::sqrt(std::max(Real(0), (*interp)(z, true))));
        };

        const auto h = [=](Real fv) -> Real {
            return squared(std::log(fv/xmax));
        };

        const ext::shared_ptr<DqFpEquation> eqn
            = (fpEquation_ == FP_A
               || (fpEquation_ == Auto && std::abs(r-q) < 0.001))?
              ext::shared_ptr<DqFpEquation>(new DqFpEquation_A(
                  K, r, q, vol, B,
                  iterationScheme_->getFixedPointIntegrator()))
            : ext::shared_ptr<DqFpEquation>(new DqFpEquation_B(
                    K, r, q, vol, B,
                    iterationScheme_->getFixedPointIntegrator()));

        Array y(x.size());
        y[0] = 0.0;

        const Size n_newton
            = iterationScheme_->getNumberOfJacobiNewtonFixedPointSteps();
        for (Size k=0; k < n_newton; ++k) {
            for (Size i=1; i < x.size(); ++i) {
                const Real tau = squared(x[i]);
                const Real b = B(tau);

                const std::tuple<Real, Real, Real> results = eqn->f(tau, b);
                const Real N = std::get<0>(results);
                const Real D = std::get<1>(results);
                const Real fv = std::get<2>(results);

                if (tau < QL_EPSILON)
                    y[i] = h(fv);
                else {
                    const std::pair<Real, Real> ndd = eqn->NDd(tau, b);
                    const Real Nd = std::get<0>(ndd);
                    const Real Dd = std::get<1>(ndd);

                    const Real fd = K*std::exp(-(r-q)*tau) * (Nd/D - Dd*N/(D*D));

                    y[i] = h(b - (fv - b)/ (fd-1));
                }
            }
            interp->updateY(y);
        }

        const Size n_fp = iterationScheme_->getNumberOfNaiveFixedPointSteps();
        for (Size k=0; k < n_fp; ++k) {
            for (Size i=1; i < x.size(); ++i) {
                const Real tau = squared(x[i]);
                const Real fv = std::get<2>(eqn->f(tau, B(tau)));

                y[i] = h(fv);
            }
            interp->updateY(y);
        }

        const detail::QdPlusAddOnValue aov(T, S, K, r, q, vol, xmax, interp);
        const Real addOn =
           (*iterationScheme_->getExerciseBoundaryToPriceIntegrator())(
               aov, 0.0, std::sqrt(T));

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

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

}





Coverage Report

Created: 2025-08-11 06:28

Line	Count	Source (jump to first uncovered line)
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 qrfpamericanengine.cpp
21		*/
22
23		#include <ql/math/distributions/normaldistribution.hpp>
24		#include <ql/math/functional.hpp>
25		#include <ql/math/integrals/gaussianquadratures.hpp>
26		#include <ql/math/integrals/tanhsinhintegral.hpp>
27		#include <ql/math/interpolations/chebyshevinterpolation.hpp>
28		#include <ql/pricingengines/blackcalculator.hpp>
29		#include <ql/pricingengines/vanilla/qdfpamericanengine.hpp>
30		#include <utility>
31		#ifndef QL_BOOST_HAS_TANH_SINH
32		# include <ql/math/integrals/gausslobattointegral.hpp>
33		#endif
34
35		namespace QuantLib {
36
37		QdFpLegendreScheme::QdFpLegendreScheme(
38		Size l, Size m, Size n, Size p):
39	0	m_(m), n_(n),
40	0	fpIntegrator_(ext::make_shared<GaussLegendreIntegrator>(l)),
41	0	exerciseBoundaryIntegrator_(
42	0	ext::make_shared<GaussLegendreIntegrator>(p)) {
43
44	0	QL_REQUIRE(m_ > 0, "at least one fixed point iteration step is needed");
45	0	QL_REQUIRE(n_ > 0, "at least one interpolation point is needed");
46	0	}
47
48		Size QdFpLegendreScheme::getNumberOfChebyshevInterpolationNodes()
49	0	const {
50	0	return n_;
51	0	}
52	0	Size QdFpLegendreScheme::getNumberOfNaiveFixedPointSteps() const {
53	0	return m_-1;
54	0	}
55		Size QdFpLegendreScheme::getNumberOfJacobiNewtonFixedPointSteps()
56	0	const {
57	0	return Size(1);
58	0	}
59
60		ext::shared_ptr<Integrator>
61	0	QdFpLegendreScheme::getFixedPointIntegrator() const {
62	0	return fpIntegrator_;
63	0	}
64		ext::shared_ptr<Integrator>
65		QdFpLegendreScheme::getExerciseBoundaryToPriceIntegrator()
66	0	const {
67	0	return exerciseBoundaryIntegrator_;
68	0	}
69
70		QdFpTanhSinhIterationScheme::QdFpTanhSinhIterationScheme(
71		Size m, Size n, Real eps)
72	0	: m_(m), n_(n),
73		#ifdef QL_BOOST_HAS_TANH_SINH
74	0	integrator_(ext::make_shared<TanhSinhIntegral>(eps))
75		#else
76		integrator_(ext::make_shared<GaussLobattoIntegral>(
77		100000, QL_MAX_REAL, 0.1*eps))
78		#endif
79	0	{}
80
81		Size QdFpTanhSinhIterationScheme::getNumberOfChebyshevInterpolationNodes()
82	0	const {
83	0	return n_;
84	0	}
85	0	Size QdFpTanhSinhIterationScheme::getNumberOfNaiveFixedPointSteps() const {
86	0	return m_-1;
87	0	}
88		Size QdFpTanhSinhIterationScheme::getNumberOfJacobiNewtonFixedPointSteps()
89	0	const {
90	0	return Size(1);
91	0	}
92
93		ext::shared_ptr<Integrator>
94	0	QdFpTanhSinhIterationScheme::getFixedPointIntegrator() const {
95	0	return integrator_;
96	0	}
97		ext::shared_ptr<Integrator>
98		QdFpTanhSinhIterationScheme::getExerciseBoundaryToPriceIntegrator()
99	0	const {
100	0	return integrator_;
101	0	}
102
103		QdFpLegendreTanhSinhScheme::QdFpLegendreTanhSinhScheme(
104		Size l, Size m, Size n, Real eps)
105	0	: QdFpLegendreScheme(l, m, n, 1),
106	0	eps_(eps) {}
107
108		ext::shared_ptr<Integrator>
109	0	QdFpLegendreTanhSinhScheme::getExerciseBoundaryToPriceIntegrator() const {
110	0	#ifdef QL_BOOST_HAS_TANH_SINH
111	0	return ext::make_shared<TanhSinhIntegral>(eps_);
112		#else
113		return ext::make_shared<GaussLobattoIntegral>(
114		100000, QL_MAX_REAL, 0.1*eps_);
115		#endif
116	0	}
117
118
119		class DqFpEquation {
120		public:
121		DqFpEquation(Rate _r,
122		Rate _q,
123		Volatility _vol,
124		std::function<Real(Real)> B,
125		ext::shared_ptr<Integrator> _integrator)
126	0	: r(_r), q(_q), vol(_vol), B(std::move(B)), integrator(std::move(_integrator)) {
127	0	const auto legendreIntegrator =
128	0	ext::dynamic_pointer_cast<GaussLegendreIntegrator>(integrator);
129
130	0	if (legendreIntegrator != nullptr) {
131	0	x_i = legendreIntegrator->getIntegration()->x();
132	0	w_i = legendreIntegrator->getIntegration()->weights();
133	0	}
134	0	}
135
136		virtual std::pair<Real, Real> NDd(Real tau, Real b) const = 0;
137		virtual std::tuple<Real, Real, Real> f(Real tau, Real b) const = 0;
138
139	0	virtual ~DqFpEquation() = default;
140
141		protected:
142	0	std::pair<Real, Real> d(Time t, Real z) const {
143	0	const Real v = vol * std::sqrt(t);
144	0	const Real m = (std::log(z) + (r-q)t)/v + 0.5v;
145
146	0	return std::make_pair(m, m-v);
147	0	}
148
149		Array x_i, w_i;
150		const Rate r, q;
151		const Volatility vol;
152
153		const std::function<Real(Real)> B;
154		const ext::shared_ptr<Integrator> integrator;
155
156		const NormalDistribution phi;
157		const CumulativeNormalDistribution Phi;
158		};
159
160		class DqFpEquation_B: public DqFpEquation {
161		public:
162		DqFpEquation_B(Real K,
163		Rate _r,
164		Rate _q,
165		Volatility _vol,
166		std::function<Real(Real)> B,
167		ext::shared_ptr<Integrator> _integrator);
168
169		std::pair<Real, Real> NDd(Real tau, Real b) const override;
170		std::tuple<Real, Real, Real> f(Real tau, Real b) const override;
171
172		private:
173		const Real K;
174		};
175
176		class DqFpEquation_A: public DqFpEquation {
177		public:
178		DqFpEquation_A(Real K,
179		Rate _r,
180		Rate _q,
181		Volatility _vol,
182		std::function<Real(Real)> B,
183		ext::shared_ptr<Integrator> _integrator);
184
185		std::pair<Real, Real> NDd(Real tau, Real b) const override;
186		std::tuple<Real, Real, Real> f(Real tau, Real b) const override;
187
188		private:
189		const Real K;
190		};
191
192		DqFpEquation_A::DqFpEquation_A(Real K,
193		Rate _r,
194		Rate _q,
195		Volatility _vol,
196		std::function<Real(Real)> B,
197		ext::shared_ptr<Integrator> _integrator)
198	0	: DqFpEquation(_r, _q, _vol, std::move(B), std::move(_integrator)), K(K) {}
199
200	0	std::tuple<Real, Real, Real> DqFpEquation_A::f(Real tau, Real b) const {
201	0	const Real v = vol * std::sqrt(tau);
202
203	0	Real N, D;
204	0	if (tau < squared(QL_EPSILON)) {
205	0	if (close_enough(b, K)) {
206	0	N = 1/(M_SQRT2M_SQRTPI v);
207	0	D = N + 0.5;
208	0	}
209	0	else {
210	0	N = 0.0;
211	0	D = (b > K)? 1.0: 0.0;
212	0	}
213	0	}
214	0	else {
215	0	const Real stv = std::sqrt(tau)/vol;
216
217	0	Real K12, K3;
218	0	if (!x_i.empty()) {
219	0	K12 = K3 = 0.0;
220
221	0	for (Integer i = x_i.size()-1; i >= 0; --i) {
222	0	const Real y = x_i[i];
223	0	const Real m = 0.25tausquared(1+y);
224	0	const std::pair<Real, Real> dpm = d(m, b/B(tau-m));
225
226	0	K12 += w_i[i] * std::exp(qtau - qm)
227	0	(0.5tau(y+1)Phi(dpm.first) + stv*phi(dpm.first));
228	0	K3 += w_i[i] * stvstd::exp(rtau-rm)phi(dpm.second);
229	0	}
230	0	} else {
231	0	K12 = (*integrator)([&, this](Real y) -> Real {
232	0	const Real m = 0.25tausquared(1+y);
233	0	const Real df = std::exp(qtau - qm);
234
235	0	if (y <= 5*QL_EPSILON - 1) {
236	0	if (close_enough(b, B(tau-m)))
237	0	return dfstv/(M_SQRT2M_SQRTPI);
238	0	else
239	0	return 0.0;
240	0	}
241	0	else {
242	0	const Real dp = d(m, b/B(tau-m)).first;
243	0	return df(0.5tau(y+1)Phi(dp) + stv*phi(dp));
244	0	}
245	0	}, -1, 1);
246
247	0	K3 = (*integrator)([&, this](Real y) -> Real {
248	0	const Real m = 0.25tausquared(1+y);
249	0	const Real df = std::exp(rtau-rm);
250
251	0	if (y <= 5*QL_EPSILON - 1) {
252	0	if (close_enough(b, B(tau-m)))
253	0	return dfstv/(M_SQRT2M_SQRTPI);
254	0	else
255	0	return 0.0;
256	0	}
257	0	else
258	0	return dfstvphi(d(m, b/B(tau-m)).second);
259	0	}, -1, 1);
260	0	}
261	0	const std::pair<Real, Real> dpm = d(tau, b/K);
262	0	N = phi(dpm.second)/v + r*K3;
263	0	D = phi(dpm.first)/v + Phi(dpm.first) + q*K12;
264	0	}
265
266	0	const Real alpha = Kstd::exp(-(r-q)tau);
267	0	Real fv;
268	0	if (tau < squared(QL_EPSILON)) {
269	0	if (close_enough(b, K))
270	0	fv = alpha;
271	0	else if (b > K)
272	0	fv = 0.0;
273	0	else {
274	0	if (close_enough(q, Real(0)))
275	0	fv = alphar((q < 0)? -1.0 : 1.0)/QL_EPSILON;
276	0	else
277	0	fv = alpha*r/q;
278	0	}
279	0	}
280	0	else
281	0	fv = alpha*N/D;
282
283	0	return std::make_tuple(N, D, fv);
284	0	}
285
286	0	std::pair<Real, Real> DqFpEquation_A::NDd(Real tau, Real b) const {
287	0	Real Dd, Nd;
288
289	0	if (tau < squared(QL_EPSILON)) {
290	0	if (close_enough(b, K)) {
291	0	const Real sqTau = std::sqrt(tau);
292	0	const Real vol2 = vol*vol;
293	0	Dd = M_1_SQRTPIM_SQRT_2(
294	0	-(0.5vol2 + r-q) / (bvolvol2sqTau) + 1 / (bvolsqTau));
295	0	Nd = M_1_SQRTPIM_SQRT_2 (-0.5vol2 + r-q) / (bvolvol2sqTau);
296	0	}
297	0	else
298	0	Dd = Nd = 0.0;
299	0	}
300	0	else {
301	0	const std::pair<Real, Real> dpm = d(tau, b/K);
302
303	0	Dd = -phi(dpm.first) * dpm.first / (bvolvol*tau) +
304	0	phi(dpm.first) / (bvol std::sqrt(tau));
305	0	Nd = -phi(dpm.second) * dpm.second / (bvolvol*tau);
306	0	}
307
308	0	return std::make_pair(Nd, Dd);
309	0	}
310
311
312		DqFpEquation_B::DqFpEquation_B(Real K,
313		Rate _r,
314		Rate _q,
315		Volatility _vol,
316		std::function<Real(Real)> B,
317		ext::shared_ptr<Integrator> _integrator)
318	0	: DqFpEquation(_r, _q, _vol, std::move(B), std::move(_integrator)), K(K) {}
319
320
321	0	std::tuple<Real, Real, Real> DqFpEquation_B::f(Real tau, Real b) const {
322	0	Real N, D;
323	0	if (tau < squared(QL_EPSILON)) {
324	0	if (close_enough(b, K))
325	0	N = D = 0.5;
326	0	else if (b < K)
327	0	N = D = 0.0;
328	0	else
329	0	N = D = 1.0;
330	0	}
331	0	else {
332	0	Real ni, di;
333	0	if (!x_i.empty()) {
334	0	const Real c = 0.5*tau;
335
336	0	ni = di = 0.0;
337	0	for (Integer i = x_i.size()-1; i >= 0; --i) {
338	0	const Real u = c*x_i[i] + c;
339	0	const std::pair<Real, Real> dpm = d(tau - u, b/B(u));
340	0	ni += w_i[i] * std::exp(ru)Phi(dpm.second);
341	0	di += w_i[i] * std::exp(qu)Phi(dpm.first);
342	0	}
343	0	ni *= c;
344	0	di *= c;
345	0	} else {
346	0	ni = (*integrator)([&, this](Real u) -> Real {
347	0	const Real df = std::exp(r*u);
348	0	if (u >= tau(1 - 5QL_EPSILON)) {
349	0	if (close_enough(b, B(u)))
350	0	return 0.5*df;
351	0	else
352	0	return df*((b < B(u)? 0.0: 1.0));
353	0	}
354	0	else
355	0	return df*Phi(d(tau - u, b/B(u)).second);
356	0	}, 0, tau);
357	0	di = (*integrator)([&, this](Real u) -> Real {
358	0	const Real df = std::exp(q*u);
359	0	if (u >= tau(1 - 5QL_EPSILON)) {
360	0	if (close_enough(b, B(u)))
361	0	return 0.5*df;
362	0	else
363	0	return df*((b < B(u)? 0.0: 1.0));
364	0	}
365	0	else
366	0	return df*Phi(d(tau - u, b/B(u)).first);
367	0	}, 0, tau);
368	0	}
369
370	0	const std::pair<Real, Real> dpm = d(tau, b/K);
371
372	0	N = Phi(dpm.second) + r*ni;
373	0	D = Phi(dpm.first) + q*di;
374	0	}
375
376	0	Real fv;
377	0	const Real alpha = Kstd::exp(-(r-q)tau);
378	0	if (tau < squared(QL_EPSILON)) {
379	0	if (close_enough(b, K) \|\| b > K)
380	0	fv = alpha;
381	0	else {
382	0	if (close_enough(q, Real(0)))
383	0	fv = alphar((q < 0)? -1.0 : 1.0)/QL_EPSILON;
384	0	else
385	0	fv = alpha*r/q;
386	0	}
387	0	}
388	0	else
389	0	fv = alpha*N/D;
390
391	0	return std::make_tuple(N, D, fv);
392	0	}
393
394	0	std::pair<Real, Real> DqFpEquation_B::NDd(Real tau, Real b) const {
395	0	const std::pair<Real, Real> dpm = d(tau, b/K);
396	0	return std::make_pair(
397	0	phi(dpm.second) / (bvolstd::sqrt(tau)),
398	0	phi(dpm.first) / (bvolstd::sqrt(tau))
399	0	);
400	0	}
401
402		QdFpAmericanEngine::QdFpAmericanEngine(
403		ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess,
404		ext::shared_ptr<QdFpIterationScheme> iterationScheme,
405		FixedPointEquation fpEquation)
406	0	: detail::QdPutCallParityEngine(std::move(bsProcess)),
407	0	iterationScheme_(std::move(iterationScheme)),
408	0	fpEquation_(fpEquation) {
409	0	}
410
411		ext::shared_ptr<QdFpIterationScheme>
412	0	QdFpAmericanEngine::fastScheme() {
413	0	static auto scheme = ext::make_shared<QdFpLegendreScheme>(7, 2, 7, 27);
414	0	return scheme;
415	0	}
416
417		ext::shared_ptr<QdFpIterationScheme>
418	0	QdFpAmericanEngine::accurateScheme() {
419	0	static auto scheme = ext::make_shared<QdFpLegendreTanhSinhScheme>(25, 5, 13, 1e-8);
420	0	return scheme;
421	0	}
422
423		ext::shared_ptr<QdFpIterationScheme>
424	0	QdFpAmericanEngine::highPrecisionScheme() {
425	0	static auto scheme = ext::make_shared<QdFpTanhSinhIterationScheme>(10, 30, 1e-10);
426	0	return scheme;
427	0	}
428
429		Real QdFpAmericanEngine::calculatePut(
430	0	Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {
431
432	0	if (r < 0.0 && q < r)
433	0	QL_FAIL("double-boundary case q<r<0 for a put option is given");
434
435	0	const Real xmax = QdPlusAmericanEngine::xMax(K, r, q);
436	0	const Size n = iterationScheme_->getNumberOfChebyshevInterpolationNodes();
437
438	0	const ext::shared_ptr<ChebyshevInterpolation> interp =
439	0	QdPlusAmericanEngine(
440	0	process_, n+1, QdPlusAmericanEngine::Halley, 1e-8)
441	0	.getPutExerciseBoundary(S, K, r, q, vol, T);
442
443	0	const Array z = interp->nodes();
444	0	const Array x = 0.5std::sqrt(T)(1.0+z);
445
446	0	const auto B = [xmax, T, &interp](Real tau) -> Real {
447	0	const Real z = 2*std::sqrt(std::abs(tau)/T)-1;
448	0	return xmaxstd::exp(-std::sqrt(std::max(Real(0), (interp)(z, true))));
449	0	};
450
451	0	const auto h = [=](Real fv) -> Real {
452	0	return squared(std::log(fv/xmax));
453	0	};
454
455	0	const ext::shared_ptr<DqFpEquation> eqn
456	0	= (fpEquation_ == FP_A
457	0	\|\| (fpEquation_ == Auto && std::abs(r-q) < 0.001))?
458	0	ext::shared_ptr<DqFpEquation>(new DqFpEquation_A(
459	0	K, r, q, vol, B,
460	0	iterationScheme_->getFixedPointIntegrator()))
461	0	: ext::shared_ptr<DqFpEquation>(new DqFpEquation_B(
462	0	K, r, q, vol, B,
463	0	iterationScheme_->getFixedPointIntegrator()));
464
465	0	Array y(x.size());
466	0	y[0] = 0.0;
467
468	0	const Size n_newton
469	0	= iterationScheme_->getNumberOfJacobiNewtonFixedPointSteps();
470	0	for (Size k=0; k < n_newton; ++k) {
471	0	for (Size i=1; i < x.size(); ++i) {
472	0	const Real tau = squared(x[i]);
473	0	const Real b = B(tau);
474
475	0	const std::tuple<Real, Real, Real> results = eqn->f(tau, b);
476	0	const Real N = std::get<0>(results);
477	0	const Real D = std::get<1>(results);
478	0	const Real fv = std::get<2>(results);
479
480	0	if (tau < QL_EPSILON)
481	0	y[i] = h(fv);
482	0	else {
483	0	const std::pair<Real, Real> ndd = eqn->NDd(tau, b);
484	0	const Real Nd = std::get<0>(ndd);
485	0	const Real Dd = std::get<1>(ndd);
486
487	0	const Real fd = Kstd::exp(-(r-q)tau) * (Nd/D - DdN/(DD));
488
489	0	y[i] = h(b - (fv - b)/ (fd-1));
490	0	}
491	0	}
492	0	interp->updateY(y);
493	0	}
494
495	0	const Size n_fp = iterationScheme_->getNumberOfNaiveFixedPointSteps();
496	0	for (Size k=0; k < n_fp; ++k) {
497	0	for (Size i=1; i < x.size(); ++i) {
498	0	const Real tau = squared(x[i]);
499	0	const Real fv = std::get<2>(eqn->f(tau, B(tau)));
500
501	0	y[i] = h(fv);
502	0	}
503	0	interp->updateY(y);
504	0	}
505
506	0	const detail::QdPlusAddOnValue aov(T, S, K, r, q, vol, xmax, interp);
507	0	const Real addOn =
508	0	(*iterationScheme_->getExerciseBoundaryToPriceIntegrator())(
509	0	aov, 0.0, std::sqrt(T));
510
511	0	const Real europeanValue = BlackCalculator(
512	0	Option::Put, K, Sstd::exp((r-q)T),
513	0	volstd::sqrt(T), std::exp(-rT)).value();
514
515	0	return std::max(europeanValue, 0.0) + std::max(0.0, addOn);
516	0	}
517
518		}
519
520
521
522