36 CumulativeNormalDistribution cumNormalDist;
40 Real lambda = (-rT + gamma * bT + 0.5 * gamma * (gamma - 1.0)
45 return std::exp(lambda) * (cumNormalDist(
d)
46 - std::pow((I /
S),
kappa) *
47 cumNormalDist(
d - 2.0 * std::log(I/
S) / std::sqrt(
variance)));
51 const Real lsh = std::log(
S/H);
52 const Real lis = std::log(I/
S);
53 const Real sv = std::sqrt(
v);
55 return std::exp(bT*gamma - rT + ((-1 +gamma)*gamma*
v)/2.)*((-(std::pow(I/
S,2*(gamma + bT/
v))/(std::exp(
squared(2*bT -
v + 2*gamma*
v + 4*lis + 2*lsh)/(8.*
v))*I))- 1/(std::exp(
squared(2*bT -
v + 2*gamma*
v + 2*lsh)/(8.*
v))*
S))/(
M_SQRT2*
M_SQRTPI*sv) +(std::pow(I/
S,2*(gamma + bT/
v))*(2*bT + (-1 + 2*gamma)*
v)*std::erfc((2*bT-
v + 2*gamma*
v + 4*lis + 2*lsh)/(2.*
M_SQRT2*sv)))/(2.*I*
v));
59 const Real lsh = std::log(
S/H);
60 const Real lis = std::log(I/
S);
61 const Real sv = std::sqrt(
v);
62 const Real ex = std::exp(
squared(2*bT -
v + 2*gamma*
v + 4*lis + 2*lsh)/(8.*
v));
63 const Real ey = std::exp(
squared(2*bT + (-1 + 2*gamma)*
v + 2*lsh)/(8.*
v));
65 return (std::exp(bT*gamma - rT + ((-1 +gamma)*gamma*
v)/2.)*((
M_SQRT2*I*
v*sv)/ey +(2*
M_SQRT2*std::pow(I/
S,2*(gamma + bT/
v))*
S*sv*(2*bT +(-1 + 2*gamma)*
v))/ex -2*std::sqrt(
M_PI)*std::pow(I/
S,2*(gamma + bT/
v))*
S*(bT +gamma*
v)*(2*bT + (-1 + 2*gamma)*
v)*std::erfc((2*bT -
v + 2*gamma*
v +4*lis + 2*lsh)/(2.*
M_SQRT2*sv)) +(
M_SQRT2*I*sv*(bT + (-0.5 + gamma)*
v +lsh))/ey - (std::pow(I/
S,2*(gamma + bT/
v))*
S*sv*(2*bT - 3*
v + 2*gamma*
v + 4*lis +2*lsh))/(
M_SQRT2*ex)))/(2.*I*
M_SQRTPI*
squared(
S*
v));
69 const Real lsh = std::log(
S/H);
70 const Real lis = std::log(I/
S);
71 const Real sv = std::sqrt(
v);
73 return std::exp(bT*gamma - rT + ((-1 + gamma)*gamma*
v)/2)*(((-std::exp(-
squared(2*bT -
v + 2*gamma*
v +2*lsh)/(8*
v)) + std::pow(I/
S,-1 + 2*gamma +(2*bT)/
v)/std::exp(
squared(2*bT -
v + 2*gamma*
v + 4*lis +2*lsh)/(8*
v)))*sv)/(
M_SQRT2*
M_SQRTPI) + ((2*bT+ (-1 + 2*gamma)*
v)*std::erfc((2*bT + (-1 + 2*gamma)*
v +2*lsh)/(2.*
M_SQRT2*sv)))/4. -(std::pow(I/
S,-1 + 2*gamma + (2*bT)/
v)*std::erfc((2*bT -
v + 2*gamma*
v + 4*lis +2*lsh)/(2.*
M_SQRT2*sv))*(2*bT + (-1 + 2*gamma)*
v + 4*lis))/4.);
77 const Real lsh = std::log(
S/H);
79 return (std::exp(bT*gamma - rT + ((-1 + gamma)*gamma*
v)/2.)*(I/std::exp(
squared(2*bT -
v + 2*gamma*
v + 2*lsh)/(8.*
v))- (std::pow(I/
S,2*(gamma + bT/
v))*
S)/std::exp(
squared(2*bT -
v + 2*gamma*
v + 4*std::log(I/
S) + 2*lsh)/(8.*
v))))/(H*I*std::sqrt(2*
M_PI)*std::sqrt(
v));
83 const Real lsh = std::log(
S/H);
84 const Real lis = std::log(I/
S);
85 const Real sv = std::sqrt(
v);
87 return (std::exp(bT*gamma - rT + ((-1 + gamma)*gamma*
v)/2.)*std::pow(I/
S,2*(gamma + bT/
v))*
S*((2*std::sqrt(2/
M_PI))/(std::exp(
squared(2*bT -
v + 2*gamma*
v + 4*lis + 2*lsh)/(8.*
v))*sv) + (1 - 2*gamma - (2*bT)/
v)*std::erfc((2*bT -
v + 2*gamma*
v + 4*lis +2*lsh)/(2.*
M_SQRT2*sv))))/(2.*I*I);
91 const Real lsh = std::log(
S/H);
92 return (std::exp(bT*gamma - rT + ((-1 + gamma)*gamma*
v)/2.)*(-(I*std::erfc((2*bT-
v + 2*gamma*
v + 2*lsh)/(2.*std::sqrt(2*
v)))) +std::pow(I/
S,2*(gamma + bT/
v))*
S*std::erfc((2*bT -
v + 2*gamma*
v + 4*std::log(I/
S) +2*lsh)/(2.*std::sqrt(2*
v)))))/(2.*I);
96 const Real lsh = std::log(
S/H);
97 const Real lis = std::log(I/
S);
98 const Real sv = std::sqrt(
v);
100 return (std::exp(bT*gamma - rT + ((-1 + gamma)*gamma*
v)/2.)*(
M_SQRT2*(-(I/std::exp(
squared(2*bT -
v +2*gamma*
v + 2*lsh)/(8.*
v))) + (std::pow(I/
S,2*(gamma +bT/
v))*
S)/std::exp(
squared(2*bT -
v + 2*gamma*
v + 4*lis +2*lsh)/(8.*
v)))*sv + gamma*I*std::sqrt(
M_PI)*
v*std::erfc((2*bT -
v + 2*gamma*
v +2*lsh)/(2.*
M_SQRT2*sv)) -
M_SQRTPI*std::pow(I/
S,2*(gamma + bT/
v))*
S*std::erfc((2*bT -
v +2*gamma*
v + 4*lis + 2*lsh)/(2.*
M_SQRT2*sv))*(gamma*
v +2*lis)))/(2.*I*std::sqrt(
M_PI)*
v);
104 const Real lsh = std::log(
S/H);
105 const Real lis = std::log(I/
S);
106 const Real sv = std::sqrt(
v);
107 const Real er = std::erfc((2*bT -
v + 2*gamma*
v + 4*lis + 2*lsh)/(2.*
M_SQRT2*sv));
109 return (std::exp(bT*gamma - rT + ((-1 + gamma)*gamma*
v)/2.)*(((-1 +gamma)*gamma*(I*std::erfc((2*bT -
v + 2*gamma*
v + 2*lsh)/(2.*
M_SQRT2*sv)) -std::pow(I/
S,2*(gamma + bT/
v))*
S*er))/(2.*I) +(2*bT*std::pow(I/
S,-1 + 2*gamma + (2*bT)/
v)*er*lis)/(
v*
v)+ (2*bT +
v - 2*gamma*
v + 2*lsh)/(2.*std::exp(std::pow(2*bT + (-1 + 2*gamma)*
v +2*lsh,2)/(8.*
v))*
M_SQRT2*
M_SQRTPI*
v*sv) -(std::pow(I/
S,-1 + 2*gamma + (2*bT)/
v)*(2*bT +
v - 2*gamma*
v +4*lis + 2*lsh))/(2.*std::exp(
squared(2*bT -
v + 2*gamma*
v + 4*lis + 2*lsh)/(8.*
v))*
M_SQRT2*
M_SQRTPI*
v*sv)))/2.;
114 ext::shared_ptr<GeneralizedBlackScholesProcess> process)
115 : process_(
std::move(process)) {
125 const Real forwardPrice =
S * dD/rfD;
138 arguments_.exercise->lastDate());
142 arguments_.exercise->lastDate());
146 arguments_.exercise->lastDate());
153 results.additionalResults[
"exerciseType"] = std::string(
"European");
161 results.value = std::max(0.0,
S - X);
162 results.delta = (
S >= X)? 1.0 : 0.0;
171 results.additionalResults[
"strikeGamma"] =
Real(0.0);
172 results.additionalResults[
"exerciseType"] = std::string(
"Immediate");
186 const Real bT = std::log(dD/rfD);
187 const Real rT = std::log(1.0/rfD);
193 const Real B0 = (bT == rT) ? X : std::max(X, rT / (rT - bT) * X);
194 const Real ht = -(bT + 2.0*std::sqrt(
variance)) * B0 / (BInfinity - B0);
196 const Real I = B0 + (BInfinity - B0) * (1 - std::exp(ht));
198 const Real fwd =
S * dD/rfD;
210 const Real phi_S_beta_I_I_rT_bT_v
212 const Real phi_S_1_I_I_rT_bT_v
214 const Real phi_S_1_X_I_rT_bT_V
217 *(1 - phi_S_beta_I_I_rT_bT_v)
218 +
S * phi_S_1_I_I_rT_bT_v
219 -
S * phi_S_1_X_I_rT_bT_V
220 - X * phi(
S, 0.0, I, I, rT, bT,
variance)
221 + X * phi(
S, 0.0, X, I, rT, bT,
variance);
223 const Real phi_S_S_beta_I_I_rT_bT_v
225 const Real phi_S_S_1_I_I_rT_bT_v
227 const Real phi_S_S_1_X_I_rT_bT_v
230 * (1 - phi_S_beta_I_I_rT_bT_v)
231 - (I - X) * std::pow(
S/I,
beta)
232 * phi_S_S_beta_I_I_rT_bT_v
233 + phi_S_1_I_I_rT_bT_v
234 +
S*phi_S_S_1_I_I_rT_bT_v
235 - phi_S_1_X_I_rT_bT_V
236 -
S*phi_S_S_1_X_I_rT_bT_v
237 - X*phi_S(
S, 0.0, I, I, rT, bT,
variance)
238 + X*phi_S(
S, 0.0, X, I, rT, bT,
variance);
240 const Date refDate =
process_->riskFreeRate()->referenceDate();
241 const Date exerciseDate = arguments_.exercise->lastDate();
249 const Real B0Dq = (dD <= rfD) ?
Real(0.0)
252 const Real htDq = tq * B0 / (BInfinity - B0)
254 *(B0Dq*(BInfinity - B0) - B0*(BInfinityDq - B0Dq))
256 const Real IDq = B0Dq + (BInfinityDq - B0Dq) * (1 - std::exp(ht))
257 - (BInfinity - B0) * std::exp(ht)*htDq;
259 const Real phi_H_S_beta_I_I_rT_bT_v
261 const Real phi_I_S_beta_I_I_rT_bT_v
263 const Real phi_gamma_S_beta_I_I_rT_bT_v
265 const Real phi_bt_S_beta_I_I_rT_bT_v
267 const Real phi_H_S_1_I_I_rT_bT_v
269 const Real phi_I_S_1_I_I_rT_bT_v
271 const Real phi_bt_S_1_I_I_rT_bT_v
272 = phi_bt(
S, 1.0, I, I, rT, bT,
variance);
273 const Real phi_I_S_1_X_I_rT_bT_v
275 const Real phi_bt_S_1_X_I_rT_bT_v
276 = phi_bt(
S, 1.0, X, I, rT, bT,
variance);
277 const Real phi_H_S_0_I_I_rT_bT_v
279 const Real phi_I_S_0_I_I_rT_bT_v
281 const Real phi_bt_S_0_I_I_rT_bT_v
282 = phi_bt(
S, 0.0, I, I, rT, bT,
variance);
283 const Real phi_I_S_0_X_I_rT_bT_v
285 const Real phi_bt_S_0_X_I_rT_bT_v
286 = phi_bt(
S, 0.0, X, I, rT, bT,
variance);
289 (IDq*std::pow(
S/I,
beta)
290 + (I-X)*std::pow(
S/I,
beta)*(betaDq*std::log(
S/I) -
beta*1/I*IDq))
291 * (1 - phi_S_beta_I_I_rT_bT_v)
292 - (I - X) * std::pow(
S/I,
beta)
293 *( phi_H_S_beta_I_I_rT_bT_v*IDq
294 +phi_I_S_beta_I_I_rT_bT_v*IDq
295 +phi_gamma_S_beta_I_I_rT_bT_v*betaDq
296 -phi_bt_S_beta_I_I_rT_bT_v*tq)
297 +
S*( phi_H_S_1_I_I_rT_bT_v*IDq
298 + phi_I_S_1_I_I_rT_bT_v*IDq
299 - phi_bt_S_1_I_I_rT_bT_v*tq)
300 -
S*( phi_I_S_1_X_I_rT_bT_v*IDq
301 - phi_bt_S_1_X_I_rT_bT_v*tq)
302 - X*( phi_H_S_0_I_I_rT_bT_v*IDq
303 + phi_I_S_0_I_I_rT_bT_v*IDq
304 - phi_bt_S_0_I_I_rT_bT_v*tq)
305 + X*( phi_I_S_0_X_I_rT_bT_v*IDq
306 - phi_bt_S_0_X_I_rT_bT_v*tq);
309 const Time tr = rdc.yearFraction(refDate, exerciseDate);
315 const Real B0Dr = (dD <= rfD) ?
Real(0) :
Real(-X*tr/std::log(dD));
316 const Real htDr = -tr * B0 / (BInfinity - B0)
318 *(B0Dr*(BInfinity - B0) - B0*(BInfinityDr - B0Dr))
320 const Real IDr = B0Dr + (BInfinityDr - B0Dr) * (1 - std::exp(ht))
321 - (BInfinity - B0) * std::exp(ht)*htDr;
324 (IDr*std::pow(
S/I,
beta)
325 + (I-X)*std::pow(
S/I,
beta)*(betaDr*std::log(
S/I) -
beta/I*IDr))
326 * (1 - phi_S_beta_I_I_rT_bT_v)
327 - (I - X) * std::pow(
S/I,
beta)
328 *( phi_H_S_beta_I_I_rT_bT_v*IDr
329 + phi_I_S_beta_I_I_rT_bT_v*IDr
330 + phi_gamma_S_beta_I_I_rT_bT_v*betaDr
332 + phi_bt_S_beta_I_I_rT_bT_v*tr)
333 +
S*( phi_H_S_1_I_I_rT_bT_v*IDr
334 + phi_I_S_1_I_I_rT_bT_v*IDr
335 + phi_rt(
S, 1.0, I, I, rT, bT,
variance)*tr
336 + phi_bt_S_1_I_I_rT_bT_v*tr)
337 -
S*( phi_I_S_1_X_I_rT_bT_v*IDr
338 + phi_rt(
S, 1.0, X, I, rT, bT,
variance)*tr
339 + phi_bt_S_1_X_I_rT_bT_v*tr)
340 - X*( phi_H_S_0_I_I_rT_bT_v*IDr
341 + phi_I_S_0_I_I_rT_bT_v*IDr
342 + phi_rt(
S, 0.0, I, I, rT, bT,
variance)*tr
343 + phi_bt_S_0_I_I_rT_bT_v*tr)
344 + X*( phi_I_S_0_X_I_rT_bT_v*IDr
345 + phi_rt(
S, 0.0, X, I, rT, bT,
variance)*tr
346 + phi_bt_S_0_X_I_rT_bT_v*tr);
360 const Real htDv = -1/std::sqrt(
variance)*varianceDv*B0/(BInfinity-B0)
363 const Real IDv = BInfinityDv*(1-std::exp(ht))
364 - (BInfinity-B0)*std::exp(ht)*htDv;
367 (IDv*std::pow(
S/I,
beta)
368 + (I-X)*std::pow(
S/I,
beta)*(betaDv*std::log(
S/I) -
beta/I*IDv))
369 * (1 - phi_S_beta_I_I_rT_bT_v)
370 - (I - X) * std::pow(
S/I,
beta)
371 *( phi_H_S_beta_I_I_rT_bT_v*IDv
372 + phi_I_S_beta_I_I_rT_bT_v*IDv
373 + phi_gamma_S_beta_I_I_rT_bT_v*betaDv
375 +
S*( phi_H_S_1_I_I_rT_bT_v*IDv
376 + phi_I_S_1_I_I_rT_bT_v*IDv
377 + phi_v(
S, 1.0, I, I, rT, bT,
variance)*varianceDv)
378 -
S*( phi_I_S_1_X_I_rT_bT_v*IDv
379 + phi_v(
S, 1.0, X, I, rT, bT,
variance)*varianceDv)
380 - X*( phi_H_S_0_I_I_rT_bT_v*IDv
381 + phi_I_S_0_I_I_rT_bT_v*IDv
382 + phi_v(
S, 0.0, I, I, rT, bT,
variance)*varianceDv)
383 + X*( phi_I_S_0_X_I_rT_bT_v*IDv
384 + phi_v(
S, 0.0, X, I, rT, bT,
variance)*varianceDv);
388 * (1 - phi_S_beta_I_I_rT_bT_v)
389 - 2*(I - X) * std::pow(
S/I,
beta-1)*
beta/I
390 *phi_S_S_beta_I_I_rT_bT_v
391 - (I - X) * std::pow(
S/I,
beta)
394 + 2*phi_S_S_1_I_I_rT_bT_v
397 - 2*phi_S_S_1_X_I_rT_bT_v
400 - X*phi_SS(
S, 0.0, I, I, rT, bT,
variance)
401 + X*phi_SS(
S, 0.0, X, I, rT, bT,
variance);
408 const Time dtr = rdc.yearFraction(refDate, exerciseDate)
409 - rdc.yearFraction(tomorrow, exerciseDate);
414 +
results.rho*rT/(tr*tr)*dtr +
results.dividendRho*(rT-bT)/(tq*tq)*dtq);
420 results.additionalResults[
"exerciseType"] = std::string(
"American");
435 "not an American Option");
437 ext::shared_ptr<AmericanExercise> ex =
438 ext::dynamic_pointer_cast<AmericanExercise>(arguments_.exercise);
439 QL_REQUIRE(ex,
"non-American exercise given");
441 "payoff at expiry not handled");
443 ext::shared_ptr<PlainVanillaPayoff>
payoff =
444 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
448 process_->blackVolatility()->blackVariance(ex->lastDate(),
451 process_->dividendYield()->discount(ex->lastDate());
453 process_->riskFreeRate()->discount(ex->lastDate());
455 QL_REQUIRE(spot > 0.0,
"negative or null underlying given");
460 std::swap(spot, strike);
461 std::swap(riskFreeDiscount, dividendDiscount);
462 payoff = ext::make_shared<PlainVanillaPayoff>(
466 if (dividendDiscount > 1.0 && riskFreeDiscount > dividendDiscount)
467 QL_FAIL(
"double-boundary case r<q<0 for a call given");
470 if (dividendDiscount >= 1.0 && dividendDiscount >= riskFreeDiscount) {
472 spot, strike, riskFreeDiscount, dividendDiscount,
variance);
476 spot, strike, riskFreeDiscount, dividendDiscount,
variance);
484 if (ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff)
491 ext::any_cast<Real>(
results_.additionalResults[
"strikeGamma"]);
492 results_.additionalResults[
"strikeGamma"] = tmp;
496 Time tr =
process_->riskFreeRate()->dayCounter().yearFraction(
497 process_->riskFreeRate()->referenceDate(),
498 arguments_.exercise->lastDate());
499 Time tq =
process_->dividendYield()->dayCounter().yearFraction(
500 process_->dividendYield()->referenceDate(),
501 arguments_.exercise->lastDate());
Maps any to either the boost or std implementation.
Bjerksund and Stensland approximation engine.
Black-formula calculator class.
const Instrument::results * results_
BjerksundStenslandApproximationEngine(ext::shared_ptr< GeneralizedBlackScholesProcess >)
void calculate() const override
OneAssetOption::results americanCallApproximation(Real S, Real X, Real rfD, Real dD, Real variance) const
OneAssetOption::results europeanCallResults(Real S, Real X, Real rfD, Real dD, Real variance) const
OneAssetOption::results immediateExercise(Real S, Real X) const
ext::shared_ptr< GeneralizedBlackScholesProcess > process_
Black 1976 calculator class.
Real dividendRho(Time maturity) const
virtual Real delta(Real spot) const
Real vega(Time maturity) const
Real strikeSensitivity() const
virtual Real gamma(Real spot) const
virtual Real thetaPerDay(Real spot, Time maturity) const
virtual Real theta(Real spot, Time maturity) const
Real rho(Time maturity) const
Time yearFraction(const Date &, const Date &, const Date &refPeriodStart=Date(), const Date &refPeriodEnd=Date()) const
Returns the period between two dates as a fraction of year.
Results from single-asset option calculation
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.
LinearInterpolation variance
Real Time
continuous quantity with 1-year units
Real DiscountFactor
discount factor between dates
Real Volatility
volatility
ext::shared_ptr< QuantLib::Payoff > payoff
functionals and combinators not included in the STL
normal, cumulative and inverse cumulative distributions
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< BlackVolTermStructure > v