QuantLib: a free/open-source library for quantitative finance
Fully annotated sources - version 1.32
Loading...
Searching...
No Matches
blackformula.cpp
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2001, 2002, 2003 Sadruddin Rejeb
5 Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2012 Ferdinando Ametrano
6 Copyright (C) 2006 Mark Joshi
7 Copyright (C) 2006 StatPro Italia srl
8 Copyright (C) 2007 Cristina Duminuco
9 Copyright (C) 2007 Chiara Fornarola
10 Copyright (C) 2013 Gary Kennedy
11 Copyright (C) 2015 Peter Caspers
12 Copyright (C) 2017 Klaus Spanderen
13 Copyright (C) 2019 Wojciech Ślusarski
14 Copyright (C) 2020 Marcin Rybacki
15
16 This file is part of QuantLib, a free-software/open-source library
17 for financial quantitative analysts and developers - http://quantlib.org/
18
19 QuantLib is free software: you can redistribute it and/or modify it
20 under the terms of the QuantLib license. You should have received a
21 copy of the license along with this program; if not, please email
22 <quantlib-dev@lists.sf.net>. The license is also available online at
23 <http://quantlib.org/license.shtml>.
24
25 This program is distributed in the hope that it will be useful, but WITHOUT
26 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
27 FOR A PARTICULAR PURPOSE. See the license for more details.
28*/
29
30#include <ql/pricingengines/blackformula.hpp>
31#include <ql/math/functional.hpp>
32#include <ql/math/solvers1d/newtonsafe.hpp>
33#include <ql/math/distributions/normaldistribution.hpp>
34#include <boost/math/special_functions/fpclassify.hpp>
35#include <boost/math/special_functions/atanh.hpp>
36#include <boost/math/special_functions/sign.hpp>
37
38namespace {
39 void checkParameters(QuantLib::Real strike,
40 QuantLib::Real forward,
41 QuantLib::Real displacement)
42 {
43 QL_REQUIRE(displacement >= 0.0, "displacement ("
44 << displacement
45 << ") must be non-negative");
46 QL_REQUIRE(strike + displacement >= 0.0,
47 "strike + displacement (" << strike << " + " << displacement
48 << ") must be non-negative");
49 QL_REQUIRE(forward + displacement > 0.0, "forward + displacement ("
50 << forward << " + "
51 << displacement
52 << ") must be positive");
53 }
54}
55
56namespace QuantLib {
57
58 Real blackFormula(Option::Type optionType,
59 Real strike,
60 Real forward,
61 Real stdDev,
62 Real discount,
63 Real displacement)
64 {
65 checkParameters(strike, forward, displacement);
66 QL_REQUIRE(stdDev>=0.0,
67 "stdDev (" << stdDev << ") must be non-negative");
68 QL_REQUIRE(discount>0.0,
69 "discount (" << discount << ") must be positive");
70
71 auto sign = Integer(optionType);
72
73 if (stdDev == 0.0)
74 return std::max((forward-strike) * sign, Real(0.0)) * discount;
75
76 forward = forward + displacement;
77 strike = strike + displacement;
78
79 // since displacement is non-negative strike==0 iff displacement==0
80 // so returning forward*discount is OK
81 if (strike==0.0)
82 return (optionType==Option::Call ? Real(forward*discount) : 0.0);
83
84 Real d1 = std::log(forward/strike)/stdDev + 0.5*stdDev;
85 Real d2 = d1 - stdDev;
86 CumulativeNormalDistribution phi;
87 Real nd1 = phi(sign * d1);
88 Real nd2 = phi(sign * d2);
89 Real result = discount * sign * (forward*nd1 - strike*nd2);
90 QL_ENSURE(result>=0.0,
91 "negative value (" << result << ") for " <<
92 stdDev << " stdDev, " <<
93 optionType << " option, " <<
94 strike << " strike , " <<
95 forward << " forward");
96 return result;
97 }
98
99 Real blackFormula(const ext::shared_ptr<PlainVanillaPayoff>& payoff,
100 Real forward,
101 Real stdDev,
102 Real discount,
103 Real displacement) {
104 return blackFormula(payoff->optionType(),
105 payoff->strike(), forward, stdDev, discount, displacement);
106 }
107
108 Real blackFormulaForwardDerivative(Option::Type optionType,
109 Real strike,
110 Real forward,
111 Real stdDev,
112 Real discount,
113 Real displacement)
114 {
115 checkParameters(strike, forward, displacement);
116 QL_REQUIRE(stdDev>=0.0,
117 "stdDev (" << stdDev << ") must be non-negative");
118 QL_REQUIRE(discount>0.0,
119 "discount (" << discount << ") must be positive");
120
121 auto sign = Integer(optionType);
122
123 if (stdDev == 0.0)
124 return sign * std::max(1.0 * boost::math::sign((forward - strike) * sign), 0.0) * discount;
125
126 forward = forward + displacement;
127 strike = strike + displacement;
128
129 if (strike == 0.0)
130 return (optionType == Option::Call ? discount : 0.0);
131
132 Real d1 = std::log(forward/strike)/stdDev + 0.5*stdDev;
133 CumulativeNormalDistribution phi;
134 return sign * phi(sign * d1) * discount;
135 }
136
137 Real blackFormulaForwardDerivative(const ext::shared_ptr<PlainVanillaPayoff>& payoff,
138 Real forward,
139 Real stdDev,
140 Real discount,
141 Real displacement)
142 {
143 return blackFormulaForwardDerivative(payoff->optionType(),
144 payoff->strike(), forward, stdDev, discount, displacement);
145 }
146
147 Real blackFormulaImpliedStdDevApproximation(Option::Type optionType,
148 Real strike,
149 Real forward,
150 Real blackPrice,
151 Real discount,
152 Real displacement)
153 {
154 checkParameters(strike, forward, displacement);
155 QL_REQUIRE(blackPrice>=0.0,
156 "blackPrice (" << blackPrice << ") must be non-negative");
157 QL_REQUIRE(discount>0.0,
158 "discount (" << discount << ") must be positive");
159
160 Real stdDev;
161 forward = forward + displacement;
162 strike = strike + displacement;
163 if (strike==forward)
164 // Brenner-Subrahmanyan (1988) and Feinstein (1988) ATM approx.
165 stdDev = blackPrice/discount*std::sqrt(2.0 * M_PI)/forward;
166 else {
167 // Corrado and Miller extended moneyness approximation
168 Real moneynessDelta = Integer(optionType) * (forward-strike);
169 Real moneynessDelta_2 = moneynessDelta/2.0;
170 Real temp = blackPrice/discount - moneynessDelta_2;
171 Real moneynessDelta_PI = moneynessDelta*moneynessDelta/M_PI;
172 Real temp2 = temp*temp-moneynessDelta_PI;
173 if (temp2<0.0) // approximation breaks down, 2 alternatives:
174 // 1. zero it
175 temp2=0.0;
176 // 2. Manaster-Koehler (1982) efficient Newton-Raphson seed
177 //return std::fabs(std::log(forward/strike))*std::sqrt(2.0);
178 temp2 = std::sqrt(temp2);
179 temp += temp2;
180 temp *= std::sqrt(2.0 * M_PI);
181 stdDev = temp/(forward+strike);
182 }
183 QL_ENSURE(stdDev>=0.0,
184 "stdDev (" << stdDev << ") must be non-negative");
185 return stdDev;
186 }
187
188 Real blackFormulaImpliedStdDevApproximation(
189 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
190 Real forward,
191 Real blackPrice,
192 Real discount,
193 Real displacement) {
194 return blackFormulaImpliedStdDevApproximation(payoff->optionType(),
195 payoff->strike(), forward, blackPrice, discount, displacement);
196 }
197
198 Real blackFormulaImpliedStdDevChambers(Option::Type optionType,
199 Real strike,
200 Real forward,
201 Real blackPrice,
202 Real blackAtmPrice,
203 Real discount,
204 Real displacement) {
205 checkParameters(strike, forward, displacement);
206 QL_REQUIRE(blackPrice >= 0.0,
207 "blackPrice (" << blackPrice << ") must be non-negative");
208 QL_REQUIRE(blackAtmPrice >= 0.0, "blackAtmPrice ("
209 << blackAtmPrice
210 << ") must be non-negative");
211 QL_REQUIRE(discount > 0.0, "discount (" << discount
212 << ") must be positive");
213
214 Real stdDev;
215
216 forward = forward + displacement;
217 strike = strike + displacement;
218 blackPrice /= discount;
219 blackAtmPrice /= discount;
220
221 Real s0 = M_SQRT2 * M_SQRTPI * blackAtmPrice /
222 forward; // Brenner-Subrahmanyam formula
223 Real priceAtmVol =
224 blackFormula(optionType, strike, forward, s0, 1.0, 0.0);
225 Real dc = blackPrice - priceAtmVol;
226
227 if (close(dc, 0.0)) {
228 stdDev = s0;
229 } else {
230 Real d1 =
231 blackFormulaStdDevDerivative(strike, forward, s0, 1.0, 0.0);
232 Real d2 = blackFormulaStdDevSecondDerivative(strike, forward, s0,
233 1.0, 0.0);
234 Real ds = 0.0;
235 Real tmp = d1 * d1 + 2.0 * d2 * dc;
236 if (std::fabs(d2) > 1E-10 && tmp >= 0.0)
237 ds = (-d1 + std::sqrt(tmp)) / d2; // second order approximation
238 else
239 if(std::fabs(d1) > 1E-10)
240 ds = dc / d1; // first order approximation
241 stdDev = s0 + ds;
242 }
243
244 QL_ENSURE(stdDev >= 0.0, "stdDev (" << stdDev
245 << ") must be non-negative");
246 return stdDev;
247 }
248
249 Real blackFormulaImpliedStdDevChambers(
250 const ext::shared_ptr<PlainVanillaPayoff> &payoff,
251 Real forward,
252 Real blackPrice,
253 Real blackAtmPrice,
254 Real discount,
255 Real displacement) {
256 return blackFormulaImpliedStdDevChambers(
257 payoff->optionType(), payoff->strike(), forward, blackPrice,
258 blackAtmPrice, discount, displacement);
259 }
260
261 namespace {
262 Real Af(Real x) {
263 return 0.5*(1.0+boost::math::sign(x)
264 *std::sqrt(1.0-std::exp(-M_2_PI*x*x)));
265 }
266 }
267
268 Real blackFormulaImpliedStdDevApproximationRS(
269 Option::Type type, Real K, Real F,
270 Real marketValue, Real df, Real displacement) {
271
272 checkParameters(K, F, displacement);
273 QL_REQUIRE(marketValue >= 0.0,
274 "blackPrice (" << marketValue << ") must be non-negative");
275 QL_REQUIRE(df > 0.0, "discount (" << df << ") must be positive");
276
277 F = F + displacement;
278 K = K + displacement;
279
280 const Real ey = F/K;
281 const Real ey2 = ey*ey;
282 const Real y = std::log(ey);
283 const Real alpha = marketValue/(K*df);
284 const Real R = 2 * alpha + ((type == Option::Call) ? Real(-ey + 1.0) : ey - 1.0);
285 const Real R2 = R*R;
286
287 const Real a = std::exp((1.0-M_2_PI)*y);
288 const Real A = squared(a - 1.0/a);
289 const Real b = std::exp(M_2_PI*y);
290 const Real B = 4.0*(b + 1/b)
291 - 2*K/F*(a + 1.0/a)*(ey2 + 1 - R2);
292 const Real C = (R2-squared(ey-1))*(squared(ey+1)-R2)/ey2;
293
294 const Real beta = 2*C/(B+std::sqrt(B*B+4*A*C));
295 const Real gamma = -M_PI_2*std::log(beta);
296
297 if (y >= 0.0) {
298 const Real M0 = K*df*(
299 (type == Option::Call) ? Real(ey*Af(std::sqrt(2*y)) - 0.5)
300 : 0.5-ey*Af(-std::sqrt(2*y)));
301
302 if (marketValue <= M0)
303 return std::sqrt(gamma+y)-std::sqrt(gamma-y);
304 else
305 return std::sqrt(gamma+y)+std::sqrt(gamma-y);
306 }
307 else {
308 const Real M0 = K*df*(
309 (type == Option::Call) ? Real(0.5*ey - Af(-std::sqrt(-2*y)))
310 : Af(std::sqrt(-2*y)) - 0.5*ey);
311
312 if (marketValue <= M0)
313 return std::sqrt(gamma-y)-std::sqrt(gamma+y);
314 else
315 return std::sqrt(gamma+y)+std::sqrt(gamma-y);
316 }
317 }
318
319 Real blackFormulaImpliedStdDevApproximationRS(
320 const ext::shared_ptr<PlainVanillaPayoff> &payoff,
321 Real F, Real marketValue,
322 Real df, Real displacement) {
323
324 return blackFormulaImpliedStdDevApproximationRS(
325 payoff->optionType(), payoff->strike(),
326 F, marketValue, df, displacement);
327 }
328
329 class BlackImpliedStdDevHelper {
330 public:
331 BlackImpliedStdDevHelper(Option::Type optionType,
332 Real strike,
333 Real forward,
334 Real undiscountedBlackPrice,
335 Real displacement = 0.0)
336 : halfOptionType_(0.5 * Integer(optionType)),
337 signedStrike_(Integer(optionType) * (strike+displacement)),
338 signedForward_(Integer(optionType) * (forward+displacement)),
339 undiscountedBlackPrice_(undiscountedBlackPrice)
340 {
341 checkParameters(strike, forward, displacement);
342 QL_REQUIRE(undiscountedBlackPrice>=0.0,
343 "undiscounted Black price (" <<
344 undiscountedBlackPrice << ") must be non-negative");
345 signedMoneyness_ = Integer(optionType) * std::log((forward+displacement)/(strike+displacement));
346 }
347
348 Real operator()(Real stdDev) const {
349 #if defined(QL_EXTRA_SAFETY_CHECKS)
350 QL_REQUIRE(stdDev>=0.0,
351 "stdDev (" << stdDev << ") must be non-negative");
352 #endif
353 if (stdDev==0.0)
354 return std::max(signedForward_-signedStrike_, Real(0.0))
355 - undiscountedBlackPrice_;
356 Real temp = halfOptionType_*stdDev;
357 Real d = signedMoneyness_/stdDev;
358 Real signedD1 = d + temp;
359 Real signedD2 = d - temp;
360 Real result = signedForward_ * N_(signedD1)
361 - signedStrike_ * N_(signedD2);
362 // numerical inaccuracies can yield a negative answer
363 return std::max(Real(0.0), result) - undiscountedBlackPrice_;
364 }
365 Real derivative(Real stdDev) const {
366 #if defined(QL_EXTRA_SAFETY_CHECKS)
367 QL_REQUIRE(stdDev>=0.0,
368 "stdDev (" << stdDev << ") must be non-negative");
369 #endif
370 Real signedD1 = signedMoneyness_/stdDev + halfOptionType_*stdDev;
371 return signedForward_*N_.derivative(signedD1);
372 }
373 private:
374 Real halfOptionType_;
375 Real signedStrike_, signedForward_;
376 Real undiscountedBlackPrice_, signedMoneyness_;
377 CumulativeNormalDistribution N_;
378 };
379
380
381 Real blackFormulaImpliedStdDev(Option::Type optionType,
382 Real strike,
383 Real forward,
384 Real blackPrice,
385 Real discount,
386 Real displacement,
387 Real guess,
388 Real accuracy,
389 Natural maxIterations)
390 {
391 checkParameters(strike, forward, displacement);
392
393 QL_REQUIRE(discount>0.0,
394 "discount (" << discount << ") must be positive");
395
396 QL_REQUIRE(blackPrice>=0.0,
397 "option price (" << blackPrice << ") must be non-negative");
398 // check the price of the "other" option implied by put-call paity
399 Real otherOptionPrice = blackPrice - Integer(optionType) * (forward-strike)*discount;
400 QL_REQUIRE(otherOptionPrice>=0.0,
401 "negative " << Option::Type(-1*optionType) <<
402 " price (" << otherOptionPrice <<
403 ") implied by put-call parity. No solution exists for " <<
404 optionType << " strike " << strike <<
405 ", forward " << forward <<
406 ", price " << blackPrice <<
407 ", deflator " << discount);
408
409 // solve for the out-of-the-money option which has
410 // greater vega/price ratio, i.e.
411 // it is numerically more robust for implied vol calculations
412 if (optionType==Option::Put && strike>forward) {
413 optionType = Option::Call;
414 blackPrice = otherOptionPrice;
415 }
416 if (optionType==Option::Call && strike<forward) {
417 optionType = Option::Put;
418 blackPrice = otherOptionPrice;
419 }
420
421 strike = strike + displacement;
422 forward = forward + displacement;
423
424 if (guess==Null<Real>())
425 guess = blackFormulaImpliedStdDevApproximation(
426 optionType, strike, forward, blackPrice, discount, displacement);
427 else
428 QL_REQUIRE(guess>=0.0,
429 "stdDev guess (" << guess << ") must be non-negative");
430 BlackImpliedStdDevHelper f(optionType, strike, forward,
431 blackPrice/discount);
432 NewtonSafe solver;
433 solver.setMaxEvaluations(maxIterations);
434 Real minSdtDev = 0.0, maxStdDev = 24.0; // 24 = 300% * sqrt(60)
435 Real stdDev = solver.solve(f, accuracy, guess, minSdtDev, maxStdDev);
436 QL_ENSURE(stdDev>=0.0,
437 "stdDev (" << stdDev << ") must be non-negative");
438 return stdDev;
439 }
440
441 Real blackFormulaImpliedStdDev(
442 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
443 Real forward,
444 Real blackPrice,
445 Real discount,
446 Real displacement,
447 Real guess,
448 Real accuracy,
449 Natural maxIterations) {
450 return blackFormulaImpliedStdDev(payoff->optionType(), payoff->strike(),
451 forward, blackPrice, discount, displacement, guess, accuracy, maxIterations);
452 }
453
454
455 namespace {
456 Real Np(Real x, Real v) {
457 return CumulativeNormalDistribution()(x/v + 0.5*v);
458 }
459 Real Nm(Real x, Real v) {
460 return std::exp(-x)*CumulativeNormalDistribution()(x/v - 0.5*v);
461 }
462 Real phi(Real x, Real v) {
463 const Real ax = 2*std::fabs(x);
464 const Real v2 = v*v;
465 return (v2-ax)/(v2+ax);
466 }
467 Real F(Real v, Real x, Real cs, Real w) {
468 return cs+Nm(x,v)+w*Np(x,v);
469 }
470 Real G(Real v, Real x, Real cs, Real w) {
471 const Real q = F(v,x,cs,w)/(1+w);
472
473 // Acklam's inverse w/o Halley's refinement step
474 // does not provide enough accuracy. But both together are
475 // slower than the boost replacement.
476 const Real k = MaddockInverseCumulativeNormal()(q);
477
478 return k + std::sqrt(k*k + 2*std::fabs(x));
479 }
480 }
481
482 Real blackFormulaImpliedStdDevLiRS(
483 Option::Type optionType,
484 Real strike,
485 Real forward,
486 Real blackPrice,
487 Real discount,
488 Real displacement,
489 Real guess,
490 Real w,
491 Real accuracy,
492 Natural maxIterations) {
493
494 QL_REQUIRE(discount>0.0,
495 "discount (" << discount << ") must be positive");
496
497 QL_REQUIRE(blackPrice>=0.0,
498 "option price (" << blackPrice << ") must be non-negative");
499
500 strike = strike + displacement;
501 forward = forward + displacement;
502
503 if (guess == Null<Real>()) {
504 guess = blackFormulaImpliedStdDevApproximationRS(
505 optionType, strike, forward,
506 blackPrice, discount, displacement);
507 }
508 else {
509 QL_REQUIRE(guess>=0.0,
510 "stdDev guess (" << guess << ") must be non-negative");
511 }
512
513 Real x = std::log(forward/strike);
514 Real cs = (optionType == Option::Call)
515 ? Real(blackPrice / (forward*discount))
516 : (blackPrice/ (forward*discount) + 1.0 - strike/forward);
517
518 QL_REQUIRE(cs >= 0.0, "normalized call price (" << cs
519 << ") must be positive");
520
521 if (x > 0) {
522 // use in-out duality
523 cs = forward/strike*cs + 1.0 - forward/strike;
524 QL_REQUIRE(cs >= 0.0, "negative option price from in-out duality");
525 x = -x;
526 }
527
528 Size nIter = 0;
529 Real dv, vk, vkp1 = guess;
530
531 do {
532 vk = vkp1;
533 const Real alphaK = (1+w)/(1+phi(x,vk));
534 vkp1 = alphaK*G(vk,x,cs,w) + (1-alphaK)*vk;
535 dv = std::fabs(vkp1 - vk);
536 } while (dv > accuracy && ++nIter < maxIterations);
537
538 QL_REQUIRE(dv <= accuracy, "max iterations exceeded");
539 QL_REQUIRE(vk >= 0.0, "stdDev (" << vk << ") must be non-negative");
540
541 return vk;
542 }
543
544 Real blackFormulaImpliedStdDevLiRS(
545 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
546 Real forward,
547 Real blackPrice,
548 Real discount,
549 Real displacement,
550 Real guess,
551 Real omega,
552 Real accuracy,
553 Natural maxIterations) {
554
555 return blackFormulaImpliedStdDevLiRS(
556 payoff->optionType(), payoff->strike(),
557 forward, blackPrice, discount, displacement,
558 guess, omega, accuracy, maxIterations);
559 }
560
561
562 Real blackFormulaCashItmProbability(Option::Type optionType,
563 Real strike,
564 Real forward,
565 Real stdDev,
566 Real displacement) {
567 checkParameters(strike, forward, displacement);
568
569 auto sign = Integer(optionType);
570
571 if (stdDev==0.0)
572 return (forward * sign > strike * sign ? 1.0 : 0.0);
573
574 forward = forward + displacement;
575 strike = strike + displacement;
576 if (strike==0.0)
577 return (optionType == Option::Call ? 1.0 : 0.0);
578 Real d2 = std::log(forward/strike)/stdDev - 0.5*stdDev;
579 CumulativeNormalDistribution phi;
580 return phi(sign * d2);
581 }
582
583 Real blackFormulaCashItmProbability(
584 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
585 Real forward,
586 Real stdDev,
587 Real displacement) {
588 return blackFormulaCashItmProbability(payoff->optionType(),
589 payoff->strike(), forward, stdDev , displacement);
590 }
591
592 Real blackFormulaAssetItmProbability(
593 Option::Type optionType,
594 Real strike,
595 Real forward,
596 Real stdDev,
597 Real displacement) {
598 checkParameters(strike, forward, displacement);
599
600 auto sign = Integer(optionType);
601
602 if (stdDev==0.0)
603 return (forward * sign < strike * sign ? 1.0 : 0.0);
604
605 forward = forward + displacement;
606 strike = strike + displacement;
607 if (strike == 0.0)
608 return (optionType == Option::Call ? 1.0 : 0.0);
609 Real d1 = std::log(forward/strike)/stdDev + 0.5*stdDev;
610 CumulativeNormalDistribution phi;
611 return phi(sign * d1);
612 }
613
614 Real blackFormulaAssetItmProbability(
615 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
616 Real forward,
617 Real stdDev,
618 Real displacement) {
619 return blackFormulaAssetItmProbability(payoff->optionType(),
620 payoff->strike(), forward, stdDev , displacement);
621 }
622
623 Real blackFormulaVolDerivative(Rate strike,
624 Rate forward,
625 Real stdDev,
626 Real expiry,
627 Real discount,
628 Real displacement)
629 {
630 return blackFormulaStdDevDerivative(strike,
631 forward,
632 stdDev,
633 discount,
634 displacement)*std::sqrt(expiry);
635 }
636
637 Real blackFormulaStdDevDerivative(Rate strike,
638 Rate forward,
639 Real stdDev,
640 Real discount,
641 Real displacement)
642 {
643 checkParameters(strike, forward, displacement);
644 QL_REQUIRE(stdDev>=0.0,
645 "stdDev (" << stdDev << ") must be non-negative");
646 QL_REQUIRE(discount>0.0,
647 "discount (" << discount << ") must be positive");
648
649 forward = forward + displacement;
650 strike = strike + displacement;
651
652 if (stdDev==0.0 || strike==0.0)
653 return 0.0;
654
655 Real d1 = std::log(forward/strike)/stdDev + .5*stdDev;
656 return discount * forward *
657 CumulativeNormalDistribution().derivative(d1);
658 }
659
660 Real blackFormulaStdDevDerivative(
661 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
662 Real forward,
663 Real stdDev,
664 Real discount,
665 Real displacement) {
666 return blackFormulaStdDevDerivative(payoff->strike(), forward,
667 stdDev, discount, displacement);
668 }
669
670 Real blackFormulaStdDevSecondDerivative(Rate strike,
671 Rate forward,
672 Real stdDev,
673 Real discount,
674 Real displacement)
675 {
676 checkParameters(strike, forward, displacement);
677 QL_REQUIRE(stdDev>=0.0,
678 "stdDev (" << stdDev << ") must be non-negative");
679 QL_REQUIRE(discount>0.0,
680 "discount (" << discount << ") must be positive");
681
682 forward = forward + displacement;
683 strike = strike + displacement;
684
685 if (stdDev==0.0 || strike==0.0)
686 return 0.0;
687
688 Real d1 = std::log(forward/strike)/stdDev + .5*stdDev;
689 Real d1p = -std::log(forward/strike)/(stdDev*stdDev) + .5;
690 return discount * forward *
691 NormalDistribution().derivative(d1) * d1p;
692 }
693
694 Real blackFormulaStdDevSecondDerivative(
695 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
696 Real forward,
697 Real stdDev,
698 Real discount,
699 Real displacement) {
700 return blackFormulaStdDevSecondDerivative(payoff->strike(), forward,
701 stdDev, discount, displacement);
702 }
703
704 Real bachelierBlackFormula(Option::Type optionType,
705 Real strike,
706 Real forward,
707 Real stdDev,
708 Real discount)
709 {
710 QL_REQUIRE(stdDev>=0.0,
711 "stdDev (" << stdDev << ") must be non-negative");
712 QL_REQUIRE(discount>0.0,
713 "discount (" << discount << ") must be positive");
714 Real d = (forward-strike) * Integer(optionType), h = d / stdDev;
715 if (stdDev==0.0)
716 return discount*std::max(d, 0.0);
717 CumulativeNormalDistribution phi;
718 Real result = discount*(stdDev*phi.derivative(h) + d*phi(h));
719 QL_ENSURE(result>=0.0,
720 "negative value (" << result << ") for " <<
721 stdDev << " stdDev, " <<
722 optionType << " option, " <<
723 strike << " strike , " <<
724 forward << " forward");
725 return result;
726 }
727
728 Real bachelierBlackFormula(
729 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
730 Real forward,
731 Real stdDev,
732 Real discount) {
733 return bachelierBlackFormula(payoff->optionType(),
734 payoff->strike(), forward, stdDev, discount);
735 }
736
737 Real bachelierBlackFormulaForwardDerivative(
738 Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount)
739 {
740 QL_REQUIRE(stdDev>=0.0,
741 "stdDev (" << stdDev << ") must be non-negative");
742 QL_REQUIRE(discount>0.0,
743 "discount (" << discount << ") must be positive");
744 auto sign = Integer(optionType);
745 if (stdDev == 0.0)
746 return sign * std::max(1.0 * boost::math::sign((forward - strike) * sign), 0.0) * discount;
747 Real d = (forward - strike) * sign, h = d / stdDev;
748 CumulativeNormalDistribution phi;
749 return sign * phi(h) * discount;
750 }
751
752 Real bachelierBlackFormulaForwardDerivative(
753 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
754 Real forward,
755 Real stdDev,
756 Real discount)
757 {
758 return bachelierBlackFormulaForwardDerivative(payoff->optionType(),
759 payoff->strike(), forward, stdDev, discount);
760 }
761
762 static Real h(Real eta) {
763
764 const static Real A0 = 3.994961687345134e-1;
765 const static Real A1 = 2.100960795068497e+1;
766 const static Real A2 = 4.980340217855084e+1;
767 const static Real A3 = 5.988761102690991e+2;
768 const static Real A4 = 1.848489695437094e+3;
769 const static Real A5 = 6.106322407867059e+3;
770 const static Real A6 = 2.493415285349361e+4;
771 const static Real A7 = 1.266458051348246e+4;
772
773 const static Real B0 = 1.000000000000000e+0;
774 const static Real B1 = 4.990534153589422e+1;
775 const static Real B2 = 3.093573936743112e+1;
776 const static Real B3 = 1.495105008310999e+3;
777 const static Real B4 = 1.323614537899738e+3;
778 const static Real B5 = 1.598919697679745e+4;
779 const static Real B6 = 2.392008891720782e+4;
780 const static Real B7 = 3.608817108375034e+3;
781 const static Real B8 = -2.067719486400926e+2;
782 const static Real B9 = 1.174240599306013e+1;
783
784 QL_REQUIRE(eta>=0.0,
785 "eta (" << eta << ") must be non-negative");
786
787 const Real num = A0 + eta * (A1 + eta * (A2 + eta * (A3 + eta * (A4 + eta
788 * (A5 + eta * (A6 + eta * A7))))));
789
790 const Real den = B0 + eta * (B1 + eta * (B2 + eta * (B3 + eta * (B4 + eta
791 * (B5 + eta * (B6 + eta * (B7 + eta * (B8 + eta * B9))))))));
792
793 return std::sqrt(eta) * (num / den);
794
795 }
796
797 Real bachelierBlackFormulaImpliedVol(Option::Type optionType,
798 Real strike,
799 Real forward,
800 Real tte,
801 Real bachelierPrice,
802 Real discount) {
803
804 const static Real SQRT_QL_EPSILON = std::sqrt(QL_EPSILON);
805
806 QL_REQUIRE(tte>0.0,
807 "tte (" << tte << ") must be positive");
808
809 Real forwardPremium = bachelierPrice/discount;
810
811 Real straddlePremium;
812 if (optionType==Option::Call){
813 straddlePremium = 2.0 * forwardPremium - (forward - strike);
814 } else {
815 straddlePremium = 2.0 * forwardPremium + (forward - strike);
816 }
817
818 Real nu = (forward - strike) / straddlePremium;
819 QL_REQUIRE(nu<1.0 || close_enough(nu,1.0),
820 "nu (" << nu << ") must be <= 1.0");
821 QL_REQUIRE(nu>-1.0 || close_enough(nu,-1.0),
822 "nu (" << nu << ") must be >= -1.0");
823
824 nu = std::max(-1.0 + QL_EPSILON, std::min(nu,1.0 - QL_EPSILON));
825
826 // nu / arctanh(nu) -> 1 as nu -> 0
827 Real eta = (std::fabs(nu) < SQRT_QL_EPSILON) ? 1.0 : Real(nu / boost::math::atanh(nu));
828
829 Real heta = h(eta);
830
831 Real impliedBpvol = std::sqrt(M_PI / (2 * tte)) * straddlePremium * heta;
832
833 return impliedBpvol;
834 }
835
836
837 Real bachelierBlackFormulaStdDevDerivative(Rate strike,
838 Rate forward,
839 Real stdDev,
840 Real discount)
841 {
842 QL_REQUIRE(stdDev>=0.0,
843 "stdDev (" << stdDev << ") must be non-negative");
844 QL_REQUIRE(discount>0.0,
845 "discount (" << discount << ") must be positive");
846
847 if (stdDev==0.0)
848 return 0.0;
849
850 Real d1 = (forward - strike)/stdDev;
851 return discount *
852 CumulativeNormalDistribution().derivative(d1);
853 }
854
855 Real bachelierBlackFormulaStdDevDerivative(
856 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
857 Real forward,
858 Real stdDev,
859 Real discount) {
860 return bachelierBlackFormulaStdDevDerivative(payoff->strike(), forward,
861 stdDev, discount);
862 }
863
864 Real bachelierBlackFormulaAssetItmProbability(
865 Option::Type optionType,
866 Real strike,
867 Real forward,
868 Real stdDev) {
869 QL_REQUIRE(stdDev>=0.0,
870 "stdDev (" << stdDev << ") must be non-negative");
871 Real d = (forward - strike) * Integer(optionType), h = d / stdDev;
872 if (stdDev==0.0)
873 return std::max(d, 0.0);
874 CumulativeNormalDistribution phi;
875 Real result = phi(h);
876 return result;
877 }
878
879 Real bachelierBlackFormulaAssetItmProbability(
880 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
881 Real forward,
882 Real stdDev) {
883 return bachelierBlackFormulaAssetItmProbability(payoff->optionType(),
884 payoff->strike(), forward, stdDev);
885 }
886}
QuantLib::CumulativeNormalDistribution
Cumulative normal distribution function.
Definition: normaldistribution.hpp:80
QuantLib::NewtonSafe
safe Newton 1-D solver
Definition: newtonsafe.hpp:40
QuantLib::NormalDistribution
Normal distribution function.
Definition: normaldistribution.hpp:44
QuantLib::Null
template class providing a null value for a given type.
Definition: null.hpp:76
QuantLib::Solver1D::setMaxEvaluations
void setMaxEvaluations(Size evaluations)
Definition: solver1d.hpp:238
QuantLib::Solver1D::solve
Real solve(const F &f, Real accuracy, Real guess, Real step) const
Definition: solver1d.hpp:84
QL_EPSILON
#define QL_EPSILON
Definition: qldefines.hpp:178
QuantLib::Real
QL_REAL Real
real number
Definition: types.hpp:50
QuantLib::Natural
unsigned QL_INTEGER Natural
positive integer
Definition: types.hpp:43
QuantLib::Integer
QL_INTEGER Integer
integer number
Definition: types.hpp:35
QuantLib::Rate
Real Rate
interest rates
Definition: types.hpp:70
QuantLib::Size
std::size_t Size
size of a container
Definition: types.hpp:58
QuantLib
Definition: any.hpp:35
QuantLib::blackFormulaAssetItmProbability
Real blackFormulaAssetItmProbability(Option::Type optionType, Real strike, Real forward, Real stdDev, Real displacement)
Definition: blackformula.cpp:592
QuantLib::blackFormulaCashItmProbability
Real blackFormulaCashItmProbability(Option::Type optionType, Real strike, Real forward, Real stdDev, Real displacement)
Definition: blackformula.cpp:562
QuantLib::blackFormulaImpliedStdDevLiRS
Real blackFormulaImpliedStdDevLiRS(Option::Type optionType, Real strike, Real forward, Real blackPrice, Real discount, Real displacement, Real guess, Real w, Real accuracy, Natural maxIterations)
Definition: blackformula.cpp:482
QuantLib::blackFormulaImpliedStdDevChambers
Real blackFormulaImpliedStdDevChambers(Option::Type optionType, Real strike, Real forward, Real blackPrice, Real blackAtmPrice, Real discount, Real displacement)
Definition: blackformula.cpp:198
QuantLib::bachelierBlackFormulaAssetItmProbability
Real bachelierBlackFormulaAssetItmProbability(Option::Type optionType, Real strike, Real forward, Real stdDev)
Definition: blackformula.cpp:864
QuantLib::bachelierBlackFormulaStdDevDerivative
Real bachelierBlackFormulaStdDevDerivative(Rate strike, Rate forward, Real stdDev, Real discount)
Definition: blackformula.cpp:837
QuantLib::bachelierBlackFormula
Real bachelierBlackFormula(Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount)
Definition: blackformula.cpp:704
QuantLib::blackFormulaImpliedStdDev
Real blackFormulaImpliedStdDev(Option::Type optionType, Real strike, Real forward, Real blackPrice, Real discount, Real displacement, Real guess, Real accuracy, Natural maxIterations)
Definition: blackformula.cpp:381
QuantLib::squared
T squared(T x)
Definition: functional.hpp:37
QuantLib::close
bool close(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:163
QuantLib::blackFormulaForwardDerivative
Real blackFormulaForwardDerivative(Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount, Real displacement)
Definition: blackformula.cpp:108
QuantLib::bachelierBlackFormulaForwardDerivative
Real bachelierBlackFormulaForwardDerivative(Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount)
Definition: blackformula.cpp:737
QuantLib::blackFormulaVolDerivative
Real blackFormulaVolDerivative(Rate strike, Rate forward, Real stdDev, Real expiry, Real discount, Real displacement)
Definition: blackformula.cpp:623
QuantLib::blackFormula
Real blackFormula(Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount, Real displacement)
Definition: blackformula.cpp:58
QuantLib::close_enough
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:182
QuantLib::bachelierBlackFormulaImpliedVol
Real bachelierBlackFormulaImpliedVol(Option::Type optionType, Real strike, Real forward, Real tte, Real bachelierPrice, Real discount)
Definition: blackformula.cpp:797
QuantLib::blackFormulaStdDevSecondDerivative
Real blackFormulaStdDevSecondDerivative(Rate strike, Rate forward, Real stdDev, Real discount, Real displacement)
Definition: blackformula.cpp:670
QuantLib::blackFormulaImpliedStdDevApproximation
Real blackFormulaImpliedStdDevApproximation(Option::Type optionType, Real strike, Real forward, Real blackPrice, Real discount, Real displacement)
Definition: blackformula.cpp:147
QuantLib::blackFormulaImpliedStdDevApproximationRS
Real blackFormulaImpliedStdDevApproximationRS(Option::Type type, Real K, Real F, Real marketValue, Real df, Real displacement)
Definition: blackformula.cpp:268
QuantLib::blackFormulaStdDevDerivative
Real blackFormulaStdDevDerivative(Rate strike, Rate forward, Real stdDev, Real discount, Real displacement)
Definition: blackformula.cpp:637
Generated by Doxygen 1.9.5