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>
31#ifndef QL_BOOST_HAS_TANH_SINH
32# include <ql/math/integrals/gausslobattointegral.hpp>
41 exerciseBoundaryIntegrator_(
44 QL_REQUIRE(
m_ > 0,
"at least one fixed point iteration step is needed");
45 QL_REQUIRE(
n_ > 0,
"at least one interpolation point is needed");
60 ext::shared_ptr<Integrator>
64 ext::shared_ptr<Integrator>
73#ifdef QL_BOOST_HAS_TANH_SINH
93 ext::shared_ptr<Integrator>
97 ext::shared_ptr<Integrator>
108 ext::shared_ptr<Integrator>
110#ifdef QL_BOOST_HAS_TANH_SINH
111 return ext::make_shared<TanhSinhIntegral>(
eps_);
113 return ext::make_shared<GaussLobattoIntegral>(
121 DqFpEquation(
Rate _r,
125 ext::shared_ptr<Integrator> _integrator)
126 : r(_r), q(_q), vol(_vol), B(
std::move(B)), integrator(
std::move(_integrator)) {
127 const auto legendreIntegrator =
128 ext::dynamic_pointer_cast<GaussLegendreIntegrator>(integrator);
130 if (legendreIntegrator !=
nullptr) {
131 x_i = legendreIntegrator->getIntegration()->x();
132 w_i = legendreIntegrator->getIntegration()->weights();
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;
139 virtual ~DqFpEquation() =
default;
142 std::pair<Real, Real> d(
Time t,
Real z)
const {
143 const Real v = vol * std::sqrt(t);
144 const Real m = (std::log(z) + (r-q)*t)/v + 0.5*v;
146 return std::make_pair(m, m-v);
154 const ext::shared_ptr<Integrator> integrator;
156 const NormalDistribution phi;
157 const CumulativeNormalDistribution Phi;
160 class DqFpEquation_B:
public DqFpEquation {
162 DqFpEquation_B(
Real K,
167 ext::shared_ptr<Integrator> _integrator);
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;
176 class DqFpEquation_A:
public DqFpEquation {
178 DqFpEquation_A(
Real K,
183 ext::shared_ptr<Integrator> _integrator);
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;
192 DqFpEquation_A::DqFpEquation_A(
Real K,
197 ext::shared_ptr<Integrator> _integrator)
198 : DqFpEquation(_r, _q, _vol,
std::move(B),
std::move(_integrator)), K(K) {}
200 std::tuple<Real, Real, Real> DqFpEquation_A::f(
Real tau,
Real b)
const {
201 const Real v = vol * std::sqrt(tau);
206 N = 1/(M_SQRT2*M_SQRTPI * v);
211 D = (b > K)? 1.0: 0.0;
215 const Real stv = std::sqrt(tau)/vol;
221 for (
Integer i = x_i.size()-1; i >= 0; --i) {
222 const Real y = x_i[i];
224 const std::pair<Real, Real> dpm = d(m, b/B(tau-m));
226 K12 += w_i[i] * std::exp(q*tau - q*m)
227 *(0.5*tau*(y+1)*Phi(dpm.first) + stv*phi(dpm.first));
228 K3 += w_i[i] * stv*std::exp(r*tau-r*m)*phi(dpm.second);
231 K12 = (*integrator)([&,
this](
Real y) ->
Real {
233 const Real df = std::exp(q*tau - q*m);
237 return df*stv/(M_SQRT2*M_SQRTPI);
242 const Real dp = d(m, b/B(tau-m)).first;
243 return df*(0.5*tau*(y+1)*Phi(dp) + stv*phi(dp));
247 K3 = (*integrator)([&,
this](
Real y) ->
Real {
249 const Real df = std::exp(r*tau-r*m);
253 return df*stv/(M_SQRT2*M_SQRTPI);
258 return df*stv*phi(d(m, b/B(tau-m)).second);
261 const std::pair<Real, Real> dpm = d(tau, b/K);
262 N = phi(dpm.second)/v + r*K3;
263 D = phi(dpm.first)/v + Phi(dpm.first) + q*K12;
266 const Real alpha = K*std::exp(-(r-q)*tau);
275 fv = alpha*r*((q < 0)? -1.0 : 1.0)/
QL_EPSILON;
283 return std::make_tuple(N, D, fv);
286 std::pair<Real, Real> DqFpEquation_A::NDd(
Real tau,
Real b)
const {
291 const Real sqTau = std::sqrt(tau);
292 const Real vol2 = vol*vol;
293 Dd = M_1_SQRTPI*M_SQRT_2*(
294 -(0.5*vol2 + r-q) / (b*vol*vol2*sqTau) + 1 / (b*vol*sqTau));
295 Nd = M_1_SQRTPI*M_SQRT_2 * (-0.5*vol2 + r-q) / (b*vol*vol2*sqTau);
301 const std::pair<Real, Real> dpm = d(tau, b/K);
303 Dd = -phi(dpm.first) * dpm.first / (b*vol*vol*tau) +
304 phi(dpm.first) / (b*vol * std::sqrt(tau));
305 Nd = -phi(dpm.second) * dpm.second / (b*vol*vol*tau);
308 return std::make_pair(Nd, Dd);
312 DqFpEquation_B::DqFpEquation_B(
Real K,
317 ext::shared_ptr<Integrator> _integrator)
318 : DqFpEquation(_r, _q, _vol,
std::move(B),
std::move(_integrator)), K(K) {}
321 std::tuple<Real, Real, Real> DqFpEquation_B::f(
Real tau,
Real b)
const {
334 const Real c = 0.5*tau;
337 for (
Integer i = x_i.size()-1; i >= 0; --i) {
338 const Real u = c*x_i[i] + c;
339 const std::pair<Real, Real> dpm = d(tau - u, b/B(u));
340 ni += w_i[i] * std::exp(r*u)*Phi(dpm.second);
341 di += w_i[i] * std::exp(q*u)*Phi(dpm.first);
346 ni = (*integrator)([&,
this](
Real u) ->
Real {
347 const Real df = std::exp(r*u);
352 return df*((b < B(u)? 0.0: 1.0));
355 return df*Phi(d(tau - u, b/B(u)).second);
357 di = (*integrator)([&,
this](
Real u) ->
Real {
358 const Real df = std::exp(q*u);
363 return df*((b < B(u)? 0.0: 1.0));
366 return df*Phi(d(tau - u, b/B(u)).first);
370 const std::pair<Real, Real> dpm = d(tau, b/K);
372 N = Phi(dpm.second) + r*ni;
373 D = Phi(dpm.first) + q*di;
377 const Real alpha = K*std::exp(-(r-q)*tau);
383 fv = alpha*r*((q < 0)? -1.0 : 1.0)/
QL_EPSILON;
391 return std::make_tuple(N, D, fv);
394 std::pair<Real, Real> DqFpEquation_B::NDd(
Real tau,
Real b)
const {
395 const std::pair<Real, Real> dpm = d(tau, b/K);
396 return std::make_pair(
397 phi(dpm.second) / (b*vol*std::sqrt(tau)),
398 phi(dpm.first) / (b*vol*std::sqrt(tau))
403 ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess,
404 ext::shared_ptr<QdFpIterationScheme> iterationScheme,
406 : detail::QdPutCallParityEngine(
std::move(bsProcess)),
407 iterationScheme_(
std::move(iterationScheme)),
408 fpEquation_(fpEquation) {
411 ext::shared_ptr<QdFpIterationScheme>
413 static auto scheme = ext::make_shared<QdFpLegendreScheme>(7, 2, 7, 27);
417 ext::shared_ptr<QdFpIterationScheme>
419 static auto scheme = ext::make_shared<QdFpLegendreTanhSinhScheme>(25, 5, 13, 1e-8);
423 ext::shared_ptr<QdFpIterationScheme>
425 static auto scheme = ext::make_shared<QdFpTanhSinhIterationScheme>(10, 30, 1e-10);
432 if (r < 0.0 && q < r)
433 QL_FAIL(
"double-boundary case q<r<0 for a put option is given");
438 const ext::shared_ptr<ChebyshevInterpolation> interp =
443 const Array z = interp->nodes();
444 const Array x = 0.5*std::sqrt(T)*(1.0+z);
446 const auto B = [xmax, T, &interp](
Real tau) ->
Real {
447 const Real z = 2*std::sqrt(std::abs(tau)/T)-1;
448 return xmax*std::exp(-std::sqrt(std::max(
Real(0), (*interp)(z,
true))));
451 const auto h = [=](
Real fv) ->
Real {
452 return squared(std::log(fv/xmax));
455 const ext::shared_ptr<DqFpEquation> eqn
458 ext::shared_ptr<DqFpEquation>(
new DqFpEquation_A(
461 : ext::shared_ptr<DqFpEquation>(
new DqFpEquation_B(
470 for (
Size k=0; k < n_newton; ++k) {
473 const Real b = B(tau);
475 const std::tuple<Real, Real, Real>
results = eqn->f(tau, b);
483 const std::pair<Real, Real> ndd = eqn->NDd(tau, b);
484 const Real Nd = std::get<0>(ndd);
485 const Real Dd = std::get<1>(ndd);
487 const Real fd = K*std::exp(-(r-q)*tau) * (Nd/D - Dd*N/(D*D));
489 y[i] = h(b - (fv - b)/ (fd-1));
496 for (
Size k=0; k < n_fp; ++k) {
499 const Real fv = std::get<2>(eqn->f(tau, B(tau)));
509 aov, 0.0, std::sqrt(T));
513 vol*std::sqrt(T), std::exp(-r*T)).
value();
515 return std::max(europeanValue, 0.0) + std::max(0.0, addOn);
1-D array used in linear algebra.
Size size() const
dimension of the array
Black 1976 calculator class.
Integral of a one-dimensional function.
static ext::shared_ptr< QdFpIterationScheme > highPrecisionScheme()
static ext::shared_ptr< QdFpIterationScheme > fastScheme()
QdFpAmericanEngine(ext::shared_ptr< GeneralizedBlackScholesProcess > bsProcess, ext::shared_ptr< QdFpIterationScheme > iterationScheme=accurateScheme(), FixedPointEquation fpEquation=Auto)
const FixedPointEquation fpEquation_
const ext::shared_ptr< QdFpIterationScheme > iterationScheme_
static ext::shared_ptr< QdFpIterationScheme > accurateScheme()
Real calculatePut(Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const override
Gauss-Legendre (l,m,n)-p Scheme.
Size getNumberOfChebyshevInterpolationNodes() const override
QdFpLegendreScheme(Size l, Size m, Size n, Size p)
ext::shared_ptr< Integrator > getFixedPointIntegrator() const override
Size getNumberOfNaiveFixedPointSteps() const override
const ext::shared_ptr< Integrator > exerciseBoundaryIntegrator_
ext::shared_ptr< Integrator > getExerciseBoundaryToPriceIntegrator() const override
Size getNumberOfJacobiNewtonFixedPointSteps() const override
const ext::shared_ptr< Integrator > fpIntegrator_
QdFpLegendreTanhSinhScheme(Size l, Size m, Size n, Real eps)
ext::shared_ptr< Integrator > getExerciseBoundaryToPriceIntegrator() const override
Size getNumberOfChebyshevInterpolationNodes() const override
QdFpTanhSinhIterationScheme(Size m, Size n, Real eps)
ext::shared_ptr< Integrator > getFixedPointIntegrator() const override
Size getNumberOfNaiveFixedPointSteps() const override
ext::shared_ptr< Integrator > getExerciseBoundaryToPriceIntegrator() const override
Size getNumberOfJacobiNewtonFixedPointSteps() const override
const ext::shared_ptr< Integrator > integrator_
American engine based on the QD+ approximation to the exercise boundary.
ext::shared_ptr< ChebyshevInterpolation > getPutExerciseBoundary(Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const
static Real xMax(Real K, Rate r, Rate q)
const ext::shared_ptr< GeneralizedBlackScholesProcess > process_
Real Time
continuous quantity with 1-year units
Real Volatility
volatility
QL_INTEGER Integer
integer number
std::size_t Size
size of a container
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)