31#ifndef QL_BOOST_HAS_TANH_SINH
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);
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);
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);
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;
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;
294 -(0.5*vol2 +
r-
q) / (
b*vol*vol2*sqTau) + 1 / (
b*vol*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;
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);
Black-formula calculator class.
chebyshev interpolation between discrete Chebyshev nodes
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_
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
#define QL_FAIL(message)
throw an error (possibly with file and line information)
ext::function< Real(Real)> b
Integral of a 1-dimensional function using the Gauss quadratures.
integral of a one-dimensional function using the adaptive Gauss-Lobatto integral
Real Time
continuous quantity with 1-year units
Real Volatility
volatility
QL_INTEGER Integer
integer number
std::size_t Size
size of a container
functionals and combinators not included in the STL
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)
normal, cumulative and inverse cumulative distributions
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< YieldTermStructure > r
ext::shared_ptr< BlackVolTermStructure > v