34#ifndef QL_BOOST_HAS_TANH_SINH
40 class QdPlusBoundaryEvaluator {
42 QdPlusBoundaryEvaluator(
45 q(dy), dr(
std::exp(-
r * tau)), dq(
std::exp(-
q * tau)),
46 ddr((
std::abs(
r*tau) > 1e-5)?
Real(
r/(1-dr)) :
Real(1/(tau*(1-0.5*
r*tau*(1-
r*tau/3))))),
47 omega(2 * (
r -
q) / sigma2),
49 (-(omega - 1) -
std::sqrt(
squared(omega - 1) + 8 * ddr / sigma2))),
50 lambdaPrime(2 * ddr*ddr /
51 (sigma2 *
std::sqrt(
squared(omega - 1) + 8 * ddr / sigma2))),
52 alpha(2 * dr / (sigma2 * (2 * lambda + omega - 1))),
53 beta(
alpha * (ddr + lambdaPrime / (2 * lambda + omega - 1)) - lambda),
54 xMax(QdPlusAmericanEngine::xMax(strike,
r,
q)),
70 return (1-dq*Phi_dp)*
S + (lambda+c0)*(K-
S-npv);
77 return 1 - dq*Phi_dp + dq/
v*phi_dp +
beta*(1-dq*Phi_dp)
84 const Real gamma = phi_dp*dq/(
v*
S);
85 const Real colour = gamma*(
q + (
r-
q)*dp/v + (1-dp*dm)/(2*tau));
87 return dq*(phi_dp/(
S*
v) - phi_dp*dp/(
S*v*v))
88 + beta*gamma + alpha/dr*colour;
91 Real xmin()
const {
return xMin; }
92 Real xmax()
const {
return xMax; }
93 Size evaluations()
const {
return nrEvaluations; }
96 void preCalculate(
Real S)
const {
99 dp = std::log(
S*dq/(K*dr))/
v + 0.5*
v;
105 npv = dr*K*Phi_dm -
S*dq*Phi_dp;
106 theta =
r*K*dr*Phi_dm -
q*
S*dq*Phi_dp - sigma2*
S/(2*
v)*dq*phi_dp;
107 charm = -dq*(phi_dp*((
r-
q)/v - dm/(2*tau)) +q*Phi_dp);
110 const CumulativeNormalDistribution Phi;
111 const NormalDistribution phi;
117 const Real omega, lambda, lambdaPrime,
alpha,
beta, xMax, xMin;
118 mutable Size nrEvaluations = 0;
119 mutable Real sc, dp, dm, Phi_dp, Phi_dm, phi_dp;
126 const Real xmax, ext::shared_ptr<Interpolation> q_z)
127 :
T_(
T), S_(
S), K_(K), xmax_(xmax),
128 r_(
r), q_(
q),
vol_(vol), q_z_(
std::move(q_z)) {}
133 const Real q = (*q_z_)(2*std::sqrt(std::max(0.0,
T_-
t)/
T_) - 1,
true);
134 const Real b_t = xmax_*std::exp(-std::sqrt(std::max(0.0,
q)));
136 const Real dr = std::exp(-r_*
t);
137 const Real dq = std::exp(-q_*
t);
143 const Real dp = std::log(S_*dq/(b_t*dr))/
v + 0.5*
v;
144 r = 2*z*(r_*K_*dr*Phi_(-dp+
v) - q_*S_*dq*Phi_(-dp));
150 r = z*(r_*K_*dr - q_*S_*dq);
151 else if (b_t*dr > S_*dq)
152 r = 2*z*(r_*K_*dr - q_*S_*dq);
161 ext::shared_ptr<GeneralizedBlackScholesProcess> process)
162 : process_(
std::move(process)) {
168 "not an American option");
171 ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
174 const Real spot = process_->x0();
175 QL_REQUIRE(spot >= 0.0,
"negative underlying given");
177 const auto maturity = arguments_.exercise->lastDate();
178 const Time T = process_->time(maturity);
179 const Real S = process_->x0();
181 const Rate r = -std::log(process_->riskFreeRate()->discount(maturity))/
T;
182 const Rate q = -std::log(process_->dividendYield()->discount(maturity))/
T;
183 const Volatility vol = process_->blackVolatility()->blackVol(
T, K);
185 QL_REQUIRE(
S >= 0,
"zero or positive underlying value is required");
186 QL_REQUIRE(K >= 0,
"zero or positive strike is required");
187 QL_REQUIRE(vol >= 0,
"zero or positive volatility is required");
190 results_.value = calculatePutWithEdgeCases(
S, K,
r,
q, vol,
T);
192 results_.value = calculatePutWithEdgeCases(K,
S,
q,
r, vol,
T);
194 QL_FAIL(
"unknown option type");
204 return std::max(K, K*std::exp(-
r*
T));
206 if (
r <= 0.0 &&
r <=
q)
209 vol*std::sqrt(
T), std::exp(-
r*
T)).value());
211 if (
close(vol, 0.0)) {
212 const auto intrinsic = [&](
Real t) ->
Real {
213 return std::max(0.0, K*std::exp(-
r*
t)-
S*std::exp(-
q*
t));
215 const Real npv0 = intrinsic(0.0);
216 const Real npvT = intrinsic(
T);
220 if (extremT > 0.0 && extremT <
T)
221 return std::max(npv0, std::max(npvT, intrinsic(extremT)));
223 return std::max(npv0, npvT);
226 return calculatePut(
S, K,
r,
q, vol,
T);
234 if (
r > 0.0 &&
q > 0.0)
235 return K*std::min(1.0,
r/
q);
236 else if (
r > 0.0 &&
q <= 0.0)
238 else if (
r == 0.0 &&
q < 0.0)
240 else if (
r == 0.0 &&
q >= 0.0)
242 else if (r < 0.0 && q >= 0.0)
244 else if (
r < 0.0 &&
q <
r)
246 else if (
r < 0.0 &&
r <=
q &&
q < 0.0)
253 ext::shared_ptr<GeneralizedBlackScholesProcess> process,
254 Size interpolationPoints,
257 : detail::QdPutCallParityEngine(
std::move(process)),
258 interpolationPoints_(interpolationPoints),
259 solverType_(solverType),
261 maxIter_((maxIter ==
Null<
Size>()) ? (
263 || solverType ==
Brent || solverType==
Ridder)? 100 : 10)
266 template <
class Solver>
268 const QdPlusBoundaryEvaluator& eval,
271 solver.setMaxEvaluations(maxIter);
272 solver.setLowerBound(eval.xmin());
274 const Real fxmin = eval(eval.xmin());
275 Real xmax = std::max(0.5*(eval.xmax() +
S), eval.xmax());
276 while (eval(xmax)*fxmin > 0.0 && eval.evaluations() <
maxIter_)
280 guess = 0.5*(xmax +
S);
283 guess = std::nextafter(xmax,
Real(-1));
284 else if (guess <= eval.xmin())
287 return solver.solve(eval,
eps_, guess, eval.xmin(), xmax);
295 return std::pair<Size, Real>(
298 const QdPlusBoundaryEvaluator eval(
S, K,
r,
q, vol, tau,
T);
314 bool resultCloseEnough;
317 const Real xmin = eval.xmin();
322 const Real fPrime = eval.derivative(x);
323 const Real lf = fx*eval.fprime2(x)/(fPrime*fPrime);
325 ?
Real(1/(1 - 0.5*lf)*fx/fPrime)
326 :
Real((1 + 0.5*lf/(1-lf))*fx/fPrime);
328 x = std::max(xmin, x - step);
329 resultCloseEnough = std::fabs(x-xOld) < 0.5*
eps_;
331 while (!resultCloseEnough && eval.evaluations() <
maxIter_);
333 if (!resultCloseEnough && !
close(std::fabs(fx), 0.0)) {
339 QL_FAIL(
"unknown solver type");
342 return std::pair<Size, Real>(eval.evaluations(), x);
345 ext::shared_ptr<ChebyshevInterpolation>
351 return ext::make_shared<ChebyshevInterpolation>(
366 if (
r < 0.0 &&
q <
r)
367 QL_FAIL(
"double-boundary case q<r<0 for a put option is given");
369 const ext::shared_ptr<Interpolation> q_z
375#ifdef QL_BOOST_HAS_TANH_SINH
379 (aov, 0.0, std::sqrt(
T));
383 "negative early exercise value " << addOn);
385 const Real europeanValue = std::max(
390 vol*std::sqrt(
T), std::exp(-
r*
T)).value()
393 return europeanValue + std::max(0.0, addOn);
Black-formula calculator class.
ext::shared_ptr< SimpleQuote > vol_
const Instrument::results * results_
chebyshev interpolation between discrete Chebyshev nodes
Black 1976 calculator class.
Integral of a one-dimensional function.
template class providing a null value for a given type.
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
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 SolverType solverType_
const Size interpolationPoints_
std::pair< Size, Real > putExerciseBoundaryAtTau(Real S, Real K, Rate r, Rate q, Volatility vol, Time T, Time tau) const
Real calculatePut(Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const override
QdPlusAmericanEngine(ext::shared_ptr< GeneralizedBlackScholesProcess >, Size interpolationPoints=8, SolverType solverType=Halley, Real eps=1e-6, Size maxIter=Null< Size >())
Real buildInSolver(const QdPlusBoundaryEvaluator &eval, Solver solver, Real S, Real strike, Size maxIter, Real guess=Null< Real >()) const
Real operator()(Real z) const
QdPlusAddOnValue(Time T, Real S, Real K, Rate r, Rate q, Volatility vol, Real xmax, ext::shared_ptr< Interpolation > q_z)
void calculate() const override
QdPutCallParityEngine(ext::shared_ptr< GeneralizedBlackScholesProcess > process)
const ext::shared_ptr< GeneralizedBlackScholesProcess > process_
Real calculatePutWithEdgeCases(Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const
floating-point comparisons
#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)
Option exercise classes and payoff function.
integral of a one-dimensional function using the adaptive Gauss-Lobatto integral
Real Time
continuous quantity with 1-year units
Real DiscountFactor
discount factor between dates
Real Volatility
volatility
std::size_t Size
size of a container
ext::shared_ptr< QuantLib::Payoff > payoff
functionals and combinators not included in the STL
bool close(const Quantity &m1, const Quantity &m2, Size n)
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< YieldTermStructure > r
ext::shared_ptr< BlackVolTermStructure > v