QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.38
Loading...
Searching...
No Matches
blackformula.cpp
Go to the documentation of this file.
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, 2024 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/distributions/normal.hpp>
35#include <boost/math/special_functions/fpclassify.hpp>
36#include <boost/math/special_functions/atanh.hpp>
37#include <boost/math/special_functions/sign.hpp>
38
39namespace {
40 void checkParameters(QuantLib::Real strike,
41 QuantLib::Real forward,
42 QuantLib::Real displacement)
43 {
44 QL_REQUIRE(displacement >= 0.0, "displacement ("
45 << displacement
46 << ") must be non-negative");
47 QL_REQUIRE(strike + displacement >= 0.0,
48 "strike + displacement (" << strike << " + " << displacement
49 << ") must be non-negative");
50 QL_REQUIRE(forward + displacement > 0.0, "forward + displacement ("
51 << forward << " + "
52 << displacement
53 << ") must be positive");
54 }
55}
56
57namespace QuantLib {
58
59 Real blackFormula(Option::Type optionType,
60 Real strike,
61 Real forward,
62 Real stdDev,
63 Real discount,
64 Real displacement)
65 {
66 checkParameters(strike, forward, displacement);
67 QL_REQUIRE(stdDev>=0.0,
68 "stdDev (" << stdDev << ") must be non-negative");
69 QL_REQUIRE(discount>0.0,
70 "discount (" << discount << ") must be positive");
71
72 auto sign = Integer(optionType);
73
74 if (stdDev == 0.0)
75 return std::max((forward-strike) * sign, Real(0.0)) * discount;
76
77 forward = forward + displacement;
78 strike = strike + displacement;
79
80 // since displacement is non-negative strike==0 iff displacement==0
81 // so returning forward*discount is OK
82 if (strike==0.0)
83 return (optionType==Option::Call ? Real(forward*discount) : 0.0);
84
85 Real d1 = std::log(forward/strike)/stdDev + 0.5*stdDev;
86 Real d2 = d1 - stdDev;
87 CumulativeNormalDistribution phi;
88 Real nd1 = phi(sign * d1);
89 Real nd2 = phi(sign * d2);
90 Real result = discount * sign * (forward*nd1 - strike*nd2);
91 QL_ENSURE(result>=0.0,
92 "negative value (" << result << ") for " <<
93 stdDev << " stdDev, " <<
94 optionType << " option, " <<
95 strike << " strike , " <<
96 forward << " forward");
97 return result;
98 }
99
100 Real blackFormula(const ext::shared_ptr<PlainVanillaPayoff>& payoff,
101 Real forward,
102 Real stdDev,
103 Real discount,
104 Real displacement) {
105 return blackFormula(payoff->optionType(),
106 payoff->strike(), forward, stdDev, discount, displacement);
107 }
108
109 Real blackFormulaForwardDerivative(Option::Type optionType,
110 Real strike,
111 Real forward,
112 Real stdDev,
113 Real discount,
114 Real displacement)
115 {
116 checkParameters(strike, forward, displacement);
117 QL_REQUIRE(stdDev>=0.0,
118 "stdDev (" << stdDev << ") must be non-negative");
119 QL_REQUIRE(discount>0.0,
120 "discount (" << discount << ") must be positive");
121
122 auto sign = Integer(optionType);
123
124 if (stdDev == 0.0)
125 return sign * std::max(1.0 * boost::math::sign((forward - strike) * sign), 0.0) * discount;
126
127 forward = forward + displacement;
128 strike = strike + displacement;
129
130 if (strike == 0.0)
131 return (optionType == Option::Call ? discount : 0.0);
132
133 Real d1 = std::log(forward/strike)/stdDev + 0.5*stdDev;
134 CumulativeNormalDistribution phi;
135 return sign * phi(sign * d1) * discount;
136 }
137
138 Real blackFormulaForwardDerivative(const ext::shared_ptr<PlainVanillaPayoff>& payoff,
139 Real forward,
140 Real stdDev,
141 Real discount,
142 Real displacement)
143 {
144 return blackFormulaForwardDerivative(payoff->optionType(),
145 payoff->strike(), forward, stdDev, discount, displacement);
146 }
147
148 Real blackFormulaImpliedStdDevApproximation(Option::Type optionType,
149 Real strike,
150 Real forward,
151 Real blackPrice,
152 Real discount,
153 Real displacement)
154 {
155 checkParameters(strike, forward, displacement);
156 QL_REQUIRE(blackPrice>=0.0,
157 "blackPrice (" << blackPrice << ") must be non-negative");
158 QL_REQUIRE(discount>0.0,
159 "discount (" << discount << ") must be positive");
160
161 Real stdDev;
162 forward = forward + displacement;
163 strike = strike + displacement;
164 if (strike==forward)
165 // Brenner-Subrahmanyan (1988) and Feinstein (1988) ATM approx.
166 stdDev = blackPrice/discount*std::sqrt(2.0 * M_PI)/forward;
167 else {
168 // Corrado and Miller extended moneyness approximation
169 Real moneynessDelta = Integer(optionType) * (forward-strike);
170 Real moneynessDelta_2 = moneynessDelta/2.0;
171 Real temp = blackPrice/discount - moneynessDelta_2;
172 Real moneynessDelta_PI = moneynessDelta*moneynessDelta/M_PI;
173 Real temp2 = temp*temp-moneynessDelta_PI;
174 if (temp2<0.0) // approximation breaks down, 2 alternatives:
175 // 1. zero it
176 temp2=0.0;
177 // 2. Manaster-Koehler (1982) efficient Newton-Raphson seed
178 //return std::fabs(std::log(forward/strike))*std::sqrt(2.0);
179 temp2 = std::sqrt(temp2);
180 temp += temp2;
181 temp *= std::sqrt(2.0 * M_PI);
182 stdDev = temp/(forward+strike);
183 }
184 QL_ENSURE(stdDev>=0.0,
185 "stdDev (" << stdDev << ") must be non-negative");
186 return stdDev;
187 }
188
189 Real blackFormulaImpliedStdDevApproximation(
190 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
191 Real forward,
192 Real blackPrice,
193 Real discount,
194 Real displacement) {
195 return blackFormulaImpliedStdDevApproximation(payoff->optionType(),
196 payoff->strike(), forward, blackPrice, discount, displacement);
197 }
198
199 Real blackFormulaImpliedStdDevChambers(Option::Type optionType,
200 Real strike,
201 Real forward,
202 Real blackPrice,
203 Real blackAtmPrice,
204 Real discount,
205 Real displacement) {
206 checkParameters(strike, forward, displacement);
207 QL_REQUIRE(blackPrice >= 0.0,
208 "blackPrice (" << blackPrice << ") must be non-negative");
209 QL_REQUIRE(blackAtmPrice >= 0.0, "blackAtmPrice ("
210 << blackAtmPrice
211 << ") must be non-negative");
212 QL_REQUIRE(discount > 0.0, "discount (" << discount
213 << ") must be positive");
214
215 Real stdDev;
216
217 forward = forward + displacement;
218 strike = strike + displacement;
219 blackPrice /= discount;
220 blackAtmPrice /= discount;
221
222 Real s0 = M_SQRT2 * M_SQRTPI * blackAtmPrice /
223 forward; // Brenner-Subrahmanyam formula
224 Real priceAtmVol =
225 blackFormula(optionType, strike, forward, s0, 1.0, 0.0);
226 Real dc = blackPrice - priceAtmVol;
227
228 if (close(dc, 0.0)) {
229 stdDev = s0;
230 } else {
231 Real d1 =
232 blackFormulaStdDevDerivative(strike, forward, s0, 1.0, 0.0);
233 Real d2 = blackFormulaStdDevSecondDerivative(strike, forward, s0,
234 1.0, 0.0);
235 Real ds = 0.0;
236 Real tmp = d1 * d1 + 2.0 * d2 * dc;
237 if (std::fabs(d2) > 1E-10 && tmp >= 0.0)
238 ds = (-d1 + std::sqrt(tmp)) / d2; // second order approximation
239 else
240 if(std::fabs(d1) > 1E-10)
241 ds = dc / d1; // first order approximation
242 stdDev = s0 + ds;
243 }
244
245 QL_ENSURE(stdDev >= 0.0, "stdDev (" << stdDev
246 << ") must be non-negative");
247 return stdDev;
248 }
249
250 Real blackFormulaImpliedStdDevChambers(
251 const ext::shared_ptr<PlainVanillaPayoff> &payoff,
252 Real forward,
253 Real blackPrice,
254 Real blackAtmPrice,
255 Real discount,
256 Real displacement) {
257 return blackFormulaImpliedStdDevChambers(
258 payoff->optionType(), payoff->strike(), forward, blackPrice,
259 blackAtmPrice, discount, displacement);
260 }
261
262 namespace {
263 Real Af(Real x) {
264 return 0.5*(1.0+boost::math::sign(x)
265 *std::sqrt(1.0-std::exp(-M_2_PI*x*x)));
266 }
267 }
268
269 Real blackFormulaImpliedStdDevApproximationRS(
270 Option::Type type, Real K, Real F,
271 Real marketValue, Real df, Real displacement) {
272
273 checkParameters(K, F, displacement);
274 QL_REQUIRE(marketValue >= 0.0,
275 "blackPrice (" << marketValue << ") must be non-negative");
276 QL_REQUIRE(df > 0.0, "discount (" << df << ") must be positive");
277
278 F = F + displacement;
279 K = K + displacement;
280
281 const Real ey = F/K;
282 const Real ey2 = ey*ey;
283 const Real y = std::log(ey);
284 const Real alpha = marketValue/(K*df);
285 const Real R = 2 * alpha + ((type == Option::Call) ? Real(-ey + 1.0) : ey - 1.0);
286 const Real R2 = R*R;
287
288 const Real a = std::exp((1.0-M_2_PI)*y);
289 const Real A = squared(a - 1.0/a);
290 const Real b = std::exp(M_2_PI*y);
291 const Real B = 4.0*(b + 1/b)
292 - 2*K/F*(a + 1.0/a)*(ey2 + 1 - R2);
293 const Real C = (R2-squared(ey-1))*(squared(ey+1)-R2)/ey2;
294
295 const Real beta = 2*C/(B+std::sqrt(B*B+4*A*C));
296 const Real gamma = -M_PI_2*std::log(beta);
297
298 if (y >= 0.0) {
299 const Real M0 = K*df*(
300 (type == Option::Call) ? Real(ey*Af(std::sqrt(2*y)) - 0.5)
301 : 0.5-ey*Af(-std::sqrt(2*y)));
302
303 if (marketValue <= M0)
304 return std::sqrt(gamma+y)-std::sqrt(gamma-y);
305 else
306 return std::sqrt(gamma+y)+std::sqrt(gamma-y);
307 }
308 else {
309 const Real M0 = K*df*(
310 (type == Option::Call) ? Real(0.5*ey - Af(-std::sqrt(-2*y)))
311 : Af(std::sqrt(-2*y)) - 0.5*ey);
312
313 if (marketValue <= M0)
314 return std::sqrt(gamma-y)-std::sqrt(gamma+y);
315 else
316 return std::sqrt(gamma+y)+std::sqrt(gamma-y);
317 }
318 }
319
320 Real blackFormulaImpliedStdDevApproximationRS(
321 const ext::shared_ptr<PlainVanillaPayoff> &payoff,
322 Real F, Real marketValue,
323 Real df, Real displacement) {
324
325 return blackFormulaImpliedStdDevApproximationRS(
326 payoff->optionType(), payoff->strike(),
327 F, marketValue, df, displacement);
328 }
329
330 class BlackImpliedStdDevHelper {
331 public:
332 BlackImpliedStdDevHelper(Option::Type optionType,
333 Real strike,
334 Real forward,
335 Real undiscountedBlackPrice,
336 Real displacement = 0.0)
337 : halfOptionType_(0.5 * Integer(optionType)),
338 signedStrike_(Integer(optionType) * (strike+displacement)),
339 signedForward_(Integer(optionType) * (forward+displacement)),
340 undiscountedBlackPrice_(undiscountedBlackPrice)
341 {
342 checkParameters(strike, forward, displacement);
343 QL_REQUIRE(undiscountedBlackPrice>=0.0,
344 "undiscounted Black price (" <<
345 undiscountedBlackPrice << ") must be non-negative");
346 signedMoneyness_ = Integer(optionType) * std::log((forward+displacement)/(strike+displacement));
347 }
348
349 Real operator()(Real stdDev) const {
350 #if defined(QL_EXTRA_SAFETY_CHECKS)
351 QL_REQUIRE(stdDev>=0.0,
352 "stdDev (" << stdDev << ") must be non-negative");
353 #endif
354 if (stdDev==0.0)
355 return std::max(signedForward_-signedStrike_, Real(0.0))
356 - undiscountedBlackPrice_;
357 Real temp = halfOptionType_*stdDev;
358 Real d = signedMoneyness_/stdDev;
359 Real signedD1 = d + temp;
360 Real signedD2 = d - temp;
361 Real result = signedForward_ * N_(signedD1)
362 - signedStrike_ * N_(signedD2);
363 // numerical inaccuracies can yield a negative answer
364 return std::max(Real(0.0), result) - undiscountedBlackPrice_;
365 }
366 Real derivative(Real stdDev) const {
367 #if defined(QL_EXTRA_SAFETY_CHECKS)
368 QL_REQUIRE(stdDev>=0.0,
369 "stdDev (" << stdDev << ") must be non-negative");
370 #endif
371 Real signedD1 = signedMoneyness_/stdDev + halfOptionType_*stdDev;
372 return signedForward_*N_.derivative(signedD1);
373 }
374 private:
375 Real halfOptionType_;
376 Real signedStrike_, signedForward_;
377 Real undiscountedBlackPrice_, signedMoneyness_;
378 CumulativeNormalDistribution N_;
379 };
380
381
382 Real blackFormulaImpliedStdDev(Option::Type optionType,
383 Real strike,
384 Real forward,
385 Real blackPrice,
386 Real discount,
387 Real displacement,
388 Real guess,
389 Real accuracy,
390 Natural maxIterations)
391 {
392 checkParameters(strike, forward, displacement);
393
394 QL_REQUIRE(discount>0.0,
395 "discount (" << discount << ") must be positive");
396
397 QL_REQUIRE(blackPrice>=0.0,
398 "option price (" << blackPrice << ") must be non-negative");
399 // check the price of the "other" option implied by put-call paity
400 Real otherOptionPrice = blackPrice - Integer(optionType) * (forward-strike)*discount;
401 QL_REQUIRE(otherOptionPrice>=0.0,
402 "negative " << Option::Type(-1*optionType) <<
403 " price (" << otherOptionPrice <<
404 ") implied by put-call parity. No solution exists for " <<
405 optionType << " strike " << strike <<
406 ", forward " << forward <<
407 ", price " << blackPrice <<
408 ", deflator " << discount);
409
410 // solve for the out-of-the-money option which has
411 // greater vega/price ratio, i.e.
412 // it is numerically more robust for implied vol calculations
413 if (optionType==Option::Put && strike>forward) {
414 optionType = Option::Call;
415 blackPrice = otherOptionPrice;
416 }
417 if (optionType==Option::Call && strike<forward) {
418 optionType = Option::Put;
419 blackPrice = otherOptionPrice;
420 }
421
422 strike = strike + displacement;
423 forward = forward + displacement;
424
425 if (guess==Null<Real>())
426 guess = blackFormulaImpliedStdDevApproximation(
427 optionType, strike, forward, blackPrice, discount, displacement);
428 else
429 QL_REQUIRE(guess>=0.0,
430 "stdDev guess (" << guess << ") must be non-negative");
431 BlackImpliedStdDevHelper f(optionType, strike, forward,
432 blackPrice/discount);
433 NewtonSafe solver;
434 solver.setMaxEvaluations(maxIterations);
435 Real minSdtDev = 0.0, maxStdDev = 24.0; // 24 = 300% * sqrt(60)
436 Real stdDev = solver.solve(f, accuracy, guess, minSdtDev, maxStdDev);
437 QL_ENSURE(stdDev>=0.0,
438 "stdDev (" << stdDev << ") must be non-negative");
439 return stdDev;
440 }
441
442 Real blackFormulaImpliedStdDev(
443 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
444 Real forward,
445 Real blackPrice,
446 Real discount,
447 Real displacement,
448 Real guess,
449 Real accuracy,
450 Natural maxIterations) {
451 return blackFormulaImpliedStdDev(payoff->optionType(), payoff->strike(),
452 forward, blackPrice, discount, displacement, guess, accuracy, maxIterations);
453 }
454
455
456 namespace {
457 Real Np(Real x, Real v) {
458 return CumulativeNormalDistribution()(x/v + 0.5*v);
459 }
460 Real Nm(Real x, Real v) {
461 return std::exp(-x)*CumulativeNormalDistribution()(x/v - 0.5*v);
462 }
463 Real phi(Real x, Real v) {
464 const Real ax = 2*std::fabs(x);
465 const Real v2 = v*v;
466 return (v2-ax)/(v2+ax);
467 }
468 Real F(Real v, Real x, Real cs, Real w) {
469 return cs+Nm(x,v)+w*Np(x,v);
470 }
471 Real G(Real v, Real x, Real cs, Real w) {
472 const Real q = F(v,x,cs,w)/(1+w);
473
474 // Acklam's inverse w/o Halley's refinement step
475 // does not provide enough accuracy. But both together are
476 // slower than the boost replacement.
477 const Real k = MaddockInverseCumulativeNormal()(q);
478
479 return k + std::sqrt(k*k + 2*std::fabs(x));
480 }
481 }
482
483 Real blackFormulaImpliedStdDevLiRS(
484 Option::Type optionType,
485 Real strike,
486 Real forward,
487 Real blackPrice,
488 Real discount,
489 Real displacement,
490 Real guess,
491 Real w,
492 Real accuracy,
493 Natural maxIterations) {
494
495 QL_REQUIRE(discount>0.0,
496 "discount (" << discount << ") must be positive");
497
498 QL_REQUIRE(blackPrice>=0.0,
499 "option price (" << blackPrice << ") must be non-negative");
500
501 strike = strike + displacement;
502 forward = forward + displacement;
503
504 if (guess == Null<Real>()) {
505 guess = blackFormulaImpliedStdDevApproximationRS(
506 optionType, strike, forward,
507 blackPrice, discount, displacement);
508 }
509 else {
510 QL_REQUIRE(guess>=0.0,
511 "stdDev guess (" << guess << ") must be non-negative");
512 }
513
514 Real x = std::log(forward/strike);
515 Real cs = (optionType == Option::Call)
516 ? Real(blackPrice / (forward*discount))
517 : (blackPrice/ (forward*discount) + 1.0 - strike/forward);
518
519 QL_REQUIRE(cs >= 0.0, "normalized call price (" << cs
520 << ") must be positive");
521
522 if (x > 0) {
523 // use in-out duality
524 cs = forward/strike*cs + 1.0 - forward/strike;
525 QL_REQUIRE(cs >= 0.0, "negative option price from in-out duality");
526 x = -x;
527 }
528
529 Size nIter = 0;
530 Real dv, vk, vkp1 = guess;
531
532 do {
533 vk = vkp1;
534 const Real alphaK = (1+w)/(1+phi(x,vk));
535 vkp1 = alphaK*G(vk,x,cs,w) + (1-alphaK)*vk;
536 dv = std::fabs(vkp1 - vk);
537 } while (dv > accuracy && ++nIter < maxIterations);
538
539 QL_REQUIRE(dv <= accuracy, "max iterations exceeded");
540 QL_REQUIRE(vk >= 0.0, "stdDev (" << vk << ") must be non-negative");
541
542 return vk;
543 }
544
545 Real blackFormulaImpliedStdDevLiRS(
546 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
547 Real forward,
548 Real blackPrice,
549 Real discount,
550 Real displacement,
551 Real guess,
552 Real omega,
553 Real accuracy,
554 Natural maxIterations) {
555
556 return blackFormulaImpliedStdDevLiRS(
557 payoff->optionType(), payoff->strike(),
558 forward, blackPrice, discount, displacement,
559 guess, omega, accuracy, maxIterations);
560 }
561
562
563 Real blackFormulaCashItmProbability(Option::Type optionType,
564 Real strike,
565 Real forward,
566 Real stdDev,
567 Real displacement) {
568 checkParameters(strike, forward, displacement);
569
570 auto sign = Integer(optionType);
571
572 if (stdDev==0.0)
573 return (forward * sign > strike * sign ? 1.0 : 0.0);
574
575 forward = forward + displacement;
576 strike = strike + displacement;
577 if (strike==0.0)
578 return (optionType == Option::Call ? 1.0 : 0.0);
579 Real d2 = std::log(forward/strike)/stdDev - 0.5*stdDev;
580 CumulativeNormalDistribution phi;
581 return phi(sign * d2);
582 }
583
584 Real blackFormulaCashItmProbability(
585 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
586 Real forward,
587 Real stdDev,
588 Real displacement) {
589 return blackFormulaCashItmProbability(payoff->optionType(),
590 payoff->strike(), forward, stdDev , displacement);
591 }
592
593 Real blackFormulaAssetItmProbability(
594 Option::Type optionType,
595 Real strike,
596 Real forward,
597 Real stdDev,
598 Real displacement) {
599 checkParameters(strike, forward, displacement);
600
601 auto sign = Integer(optionType);
602
603 if (stdDev==0.0)
604 return (forward * sign < strike * sign ? 1.0 : 0.0);
605
606 forward = forward + displacement;
607 strike = strike + displacement;
608 if (strike == 0.0)
609 return (optionType == Option::Call ? 1.0 : 0.0);
610 Real d1 = std::log(forward/strike)/stdDev + 0.5*stdDev;
611 CumulativeNormalDistribution phi;
612 return phi(sign * d1);
613 }
614
615 Real blackFormulaAssetItmProbability(
616 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
617 Real forward,
618 Real stdDev,
619 Real displacement) {
620 return blackFormulaAssetItmProbability(payoff->optionType(),
621 payoff->strike(), forward, stdDev , displacement);
622 }
623
624 Real blackFormulaVolDerivative(Rate strike,
625 Rate forward,
626 Real stdDev,
627 Real expiry,
628 Real discount,
629 Real displacement)
630 {
631 return blackFormulaStdDevDerivative(strike,
632 forward,
633 stdDev,
634 discount,
635 displacement)*std::sqrt(expiry);
636 }
637
638 Real blackFormulaStdDevDerivative(Rate strike,
639 Rate forward,
640 Real stdDev,
641 Real discount,
642 Real displacement)
643 {
644 checkParameters(strike, forward, displacement);
645 QL_REQUIRE(stdDev>=0.0,
646 "stdDev (" << stdDev << ") must be non-negative");
647 QL_REQUIRE(discount>0.0,
648 "discount (" << discount << ") must be positive");
649
650 forward = forward + displacement;
651 strike = strike + displacement;
652
653 if (stdDev==0.0 || strike==0.0)
654 return 0.0;
655
656 Real d1 = std::log(forward/strike)/stdDev + .5*stdDev;
657 return discount * forward *
658 CumulativeNormalDistribution().derivative(d1);
659 }
660
661 Real blackFormulaStdDevDerivative(
662 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
663 Real forward,
664 Real stdDev,
665 Real discount,
666 Real displacement) {
667 return blackFormulaStdDevDerivative(payoff->strike(), forward,
668 stdDev, discount, displacement);
669 }
670
671 Real blackFormulaStdDevSecondDerivative(Rate strike,
672 Rate forward,
673 Real stdDev,
674 Real discount,
675 Real displacement)
676 {
677 checkParameters(strike, forward, displacement);
678 QL_REQUIRE(stdDev>=0.0,
679 "stdDev (" << stdDev << ") must be non-negative");
680 QL_REQUIRE(discount>0.0,
681 "discount (" << discount << ") must be positive");
682
683 forward = forward + displacement;
684 strike = strike + displacement;
685
686 if (stdDev==0.0 || strike==0.0)
687 return 0.0;
688
689 Real d1 = std::log(forward/strike)/stdDev + .5*stdDev;
690 Real d1p = -std::log(forward/strike)/(stdDev*stdDev) + .5;
691 return discount * forward *
692 NormalDistribution().derivative(d1) * d1p;
693 }
694
695 Real blackFormulaStdDevSecondDerivative(
696 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
697 Real forward,
698 Real stdDev,
699 Real discount,
700 Real displacement) {
701 return blackFormulaStdDevSecondDerivative(payoff->strike(), forward,
702 stdDev, discount, displacement);
703 }
704
705 Real bachelierBlackFormula(Option::Type optionType,
706 Real strike,
707 Real forward,
708 Real stdDev,
709 Real discount)
710 {
711 QL_REQUIRE(stdDev>=0.0,
712 "stdDev (" << stdDev << ") must be non-negative");
713 QL_REQUIRE(discount>0.0,
714 "discount (" << discount << ") must be positive");
715 Real d = (forward-strike) * Integer(optionType), h = d / stdDev;
716 if (stdDev==0.0)
717 return discount*std::max(d, 0.0);
718 CumulativeNormalDistribution phi;
719 Real result = discount*(stdDev*phi.derivative(h) + d*phi(h));
720 QL_ENSURE(result>=0.0,
721 "negative value (" << result << ") for " <<
722 stdDev << " stdDev, " <<
723 optionType << " option, " <<
724 strike << " strike , " <<
725 forward << " forward");
726 return result;
727 }
728
729 Real bachelierBlackFormula(
730 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
731 Real forward,
732 Real stdDev,
733 Real discount) {
734 return bachelierBlackFormula(payoff->optionType(),
735 payoff->strike(), forward, stdDev, discount);
736 }
737
738 Real bachelierBlackFormulaForwardDerivative(
739 Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount)
740 {
741 QL_REQUIRE(stdDev>=0.0,
742 "stdDev (" << stdDev << ") must be non-negative");
743 QL_REQUIRE(discount>0.0,
744 "discount (" << discount << ") must be positive");
745 auto sign = Integer(optionType);
746 if (stdDev == 0.0)
747 return sign * std::max(1.0 * boost::math::sign((forward - strike) * sign), 0.0) * discount;
748 Real d = (forward - strike) * sign, h = d / stdDev;
749 CumulativeNormalDistribution phi;
750 return sign * phi(h) * discount;
751 }
752
753 Real bachelierBlackFormulaForwardDerivative(
754 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
755 Real forward,
756 Real stdDev,
757 Real discount)
758 {
759 return bachelierBlackFormulaForwardDerivative(payoff->optionType(),
760 payoff->strike(), forward, stdDev, discount);
761 }
762
763 static Real h(Real eta) {
764
765 const static Real A0 = 3.994961687345134e-1;
766 const static Real A1 = 2.100960795068497e+1;
767 const static Real A2 = 4.980340217855084e+1;
768 const static Real A3 = 5.988761102690991e+2;
769 const static Real A4 = 1.848489695437094e+3;
770 const static Real A5 = 6.106322407867059e+3;
771 const static Real A6 = 2.493415285349361e+4;
772 const static Real A7 = 1.266458051348246e+4;
773
774 const static Real B0 = 1.000000000000000e+0;
775 const static Real B1 = 4.990534153589422e+1;
776 const static Real B2 = 3.093573936743112e+1;
777 const static Real B3 = 1.495105008310999e+3;
778 const static Real B4 = 1.323614537899738e+3;
779 const static Real B5 = 1.598919697679745e+4;
780 const static Real B6 = 2.392008891720782e+4;
781 const static Real B7 = 3.608817108375034e+3;
782 const static Real B8 = -2.067719486400926e+2;
783 const static Real B9 = 1.174240599306013e+1;
784
785 QL_REQUIRE(eta>=0.0,
786 "eta (" << eta << ") must be non-negative");
787
788 const Real num = A0 + eta * (A1 + eta * (A2 + eta * (A3 + eta * (A4 + eta
789 * (A5 + eta * (A6 + eta * A7))))));
790
791 const Real den = B0 + eta * (B1 + eta * (B2 + eta * (B3 + eta * (B4 + eta
792 * (B5 + eta * (B6 + eta * (B7 + eta * (B8 + eta * B9))))))));
793
794 return std::sqrt(eta) * (num / den);
795
796 }
797
798 Real bachelierBlackFormulaImpliedVolChoi(Option::Type optionType,
799 Real strike,
800 Real forward,
801 Real tte,
802 Real bachelierPrice,
803 Real discount) {
804
805 const static Real SQRT_QL_EPSILON = std::sqrt(QL_EPSILON);
806
807 QL_REQUIRE(tte>0.0,
808 "tte (" << tte << ") must be positive");
809
810 Real forwardPremium = bachelierPrice/discount;
811
812 Real straddlePremium;
813 if (optionType==Option::Call){
814 straddlePremium = 2.0 * forwardPremium - (forward - strike);
815 } else {
816 straddlePremium = 2.0 * forwardPremium + (forward - strike);
817 }
818
819 Real nu = (forward - strike) / straddlePremium;
820 QL_REQUIRE(nu<1.0 || close_enough(nu,1.0),
821 "nu (" << nu << ") must be <= 1.0");
822 QL_REQUIRE(nu>-1.0 || close_enough(nu,-1.0),
823 "nu (" << nu << ") must be >= -1.0");
824
825 nu = std::max(-1.0 + QL_EPSILON, std::min(nu,1.0 - QL_EPSILON));
826
827 // nu / arctanh(nu) -> 1 as nu -> 0
828 Real eta = (std::fabs(nu) < SQRT_QL_EPSILON) ? 1.0 : Real(nu / boost::math::atanh(nu));
829
830 Real heta = h(eta);
831
832 Real impliedBpvol = std::sqrt(M_PI / (2 * tte)) * straddlePremium * heta;
833
834 return impliedBpvol;
835 }
836
837 namespace {
838
839 boost::math::normal_distribution<Real> normal_dist;
840
841 Real phi(const Real x) {
842 return boost::math::pdf(normal_dist, x);
843 }
844
845 Real Phi(const Real x) {
846 return boost::math::cdf(normal_dist, x);
847 }
848
849 Real PhiTilde(const Real x) {
850 return Phi(x) + phi(x) / x;
851 }
852
853 Real inversePhiTilde(const Real PhiTildeStar) {
854 QL_REQUIRE(PhiTildeStar < 0.0,
855 "inversePhiTilde(" << PhiTildeStar << "): negative argument required");
856 Real xbar;
857 if (PhiTildeStar < -0.001882039271) {
858 Real g = 1.0 / (PhiTildeStar - 0.5);
859 Real xibar =
860 (0.032114372355 -
861 g * g *
862 (0.016969777977 - g * g * (2.6207332461E-3 - 9.6066952861E-5 * g * g))) /
863 (1.0 -
864 g * g * (0.6635646938 - g * g * (0.14528712196 - 0.010472855461 * g * g)));
865 xbar = g * (0.3989422804014326 + xibar * g * g);
866 } else {
867 Real h = std::sqrt(-std::log(-PhiTildeStar));
868 xbar =
869 (9.4883409779 - h * (9.6320903635 - h * (0.58556997323 + 2.1464093351 * h))) /
870 (1.0 - h * (0.65174820867 + h * (1.5120247828 + 6.6437847132E-5 * h)));
871 }
872 Real q = (PhiTilde(xbar) - PhiTildeStar) / phi(xbar);
873 Real xstar =
874 xbar +
875 3.0 * q * xbar * xbar * (2.0 - q * xbar * (2.0 + xbar * xbar)) /
876 (6.0 + q * xbar *
877 (-12.0 +
878 xbar * (6.0 * q + xbar * (-6.0 + q * xbar * (3.0 + xbar * xbar)))));
879 return xstar;
880 }
881
882 } // namespace
883
884 Real bachelierBlackFormulaImpliedVol(Option::Type optionType,
885 Real strike,
886 Real forward,
887 Real tte,
888 Real bachelierPrice,
889 Real discount) {
890
891 Real theta = optionType == Option::Call ? 1.0 : -1.0;
892
893 // compound bechelierPrice, so that effectively discount = 1
894
895 bachelierPrice /= discount;
896
897 // handle case strike = forward
898
899 if (close_enough(strike, forward)) {
900 return bachelierPrice / (std::sqrt(tte) * phi(0.0));
901 }
902
903 // handle case strike != forward
904
905 Real timeValue = bachelierPrice - std::max(theta * (forward - strike), 0.0);
906
907 if (close_enough(timeValue, 0.0))
908 return 0.0;
909
910 QL_REQUIRE(timeValue > 0.0, "bachelierBlackFormulaImpliedVolExact(theta="
911 << theta << ",strike=" << strike << ",forward=" << forward
912 << ",tte=" << tte << ",price=" << bachelierPrice
913 << "): option price implies negative time value ("
914 << timeValue << ")");
915
916 Real PhiTildeStar = -std::abs(timeValue / (strike - forward));
917 Real xstar = inversePhiTilde(PhiTildeStar);
918 Real impliedVol = std::abs((strike - forward) / (xstar * std::sqrt(tte)));
919
920 return impliedVol;
921 }
922
923 Real bachelierBlackFormulaStdDevDerivative(Rate strike,
924 Rate forward,
925 Real stdDev,
926 Real discount)
927 {
928 QL_REQUIRE(stdDev>=0.0,
929 "stdDev (" << stdDev << ") must be non-negative");
930 QL_REQUIRE(discount>0.0,
931 "discount (" << discount << ") must be positive");
932
933 if (stdDev==0.0)
934 return 0.0;
935
936 Real d1 = (forward - strike)/stdDev;
937 return discount *
938 CumulativeNormalDistribution().derivative(d1);
939 }
940
941 Real bachelierBlackFormulaStdDevDerivative(
942 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
943 Real forward,
944 Real stdDev,
945 Real discount) {
946 return bachelierBlackFormulaStdDevDerivative(payoff->strike(), forward,
947 stdDev, discount);
948 }
949
950 Real bachelierBlackFormulaAssetItmProbability(
951 Option::Type optionType,
952 Real strike,
953 Real forward,
954 Real stdDev) {
955 QL_REQUIRE(stdDev>=0.0,
956 "stdDev (" << stdDev << ") must be non-negative");
957 Real d = (forward - strike) * Integer(optionType), h = d / stdDev;
958 if (stdDev==0.0)
959 return std::max(d, 0.0);
960 CumulativeNormalDistribution phi;
961 Real result = phi(h);
962 return result;
963 }
964
965 Real bachelierBlackFormulaAssetItmProbability(
966 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
967 Real forward,
968 Real stdDev) {
969 return bachelierBlackFormulaAssetItmProbability(payoff->optionType(),
970 payoff->strike(), forward, stdDev);
971 }
972}
blackformula.hpp
Black formula.
QuantLib::CumulativeNormalDistribution
Cumulative normal distribution function.
Definition: normaldistribution.hpp:68
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:59
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_ENSURE
#define QL_ENSURE(condition, message)
throw an error if the given post-condition is not verified
Definition: errors.hpp:130
QL_REQUIRE
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
d
Date d
Definition: exchangeratemanager.cpp:32
b
std::function< Real(Real)> b
Definition: extendedornsteinuhlenbeckprocess.cpp:30
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
theta
Real theta
Definition: hestonrndcalculator.cpp:36
payoff
ext::shared_ptr< QuantLib::Payoff > payoff
Definition: integralhestonvarianceoptionengine.cpp:350
functional.hpp
functionals and combinators not included in the STL
M_SQRTPI
#define M_SQRTPI
Definition: mathconstants.hpp:74
M_SQRT2
#define M_SQRT2
Definition: mathconstants.hpp:94
M_PI_2
#define M_PI_2
Definition: mathconstants.hpp:62
M_2_PI
#define M_2_PI
Definition: mathconstants.hpp:82
M_PI
#define M_PI
Definition: mathconstants.hpp:54
QuantLib::exponential_integrals_helper::g
Real g(Real x)
Definition: exponentialintegrals.cpp:57
QuantLib
Definition: any.hpp:37
QuantLib::blackFormulaAssetItmProbability
Real blackFormulaAssetItmProbability(Option::Type optionType, Real strike, Real forward, Real stdDev, Real displacement)
Definition: blackformula.cpp:593
QuantLib::blackFormulaCashItmProbability
Real blackFormulaCashItmProbability(Option::Type optionType, Real strike, Real forward, Real stdDev, Real displacement)
Definition: blackformula.cpp:563
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:483
QuantLib::blackFormulaImpliedStdDevChambers
Real blackFormulaImpliedStdDevChambers(Option::Type optionType, Real strike, Real forward, Real blackPrice, Real blackAtmPrice, Real discount, Real displacement)
Definition: blackformula.cpp:199
QuantLib::bachelierBlackFormulaAssetItmProbability
Real bachelierBlackFormulaAssetItmProbability(Option::Type optionType, Real strike, Real forward, Real stdDev)
Definition: blackformula.cpp:950
QuantLib::bachelierBlackFormulaStdDevDerivative
Real bachelierBlackFormulaStdDevDerivative(Rate strike, Rate forward, Real stdDev, Real discount)
Definition: blackformula.cpp:923
QuantLib::bachelierBlackFormula
Real bachelierBlackFormula(Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount)
Definition: blackformula.cpp:705
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:382
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:109
QuantLib::bachelierBlackFormulaForwardDerivative
Real bachelierBlackFormulaForwardDerivative(Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount)
Definition: blackformula.cpp:738
QuantLib::blackFormulaVolDerivative
Real blackFormulaVolDerivative(Rate strike, Rate forward, Real stdDev, Real expiry, Real discount, Real displacement)
Definition: blackformula.cpp:624
QuantLib::blackFormula
Real blackFormula(Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount, Real displacement)
Definition: blackformula.cpp:59
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:884
QuantLib::bachelierBlackFormulaImpliedVolChoi
Real bachelierBlackFormulaImpliedVolChoi(Option::Type optionType, Real strike, Real forward, Real tte, Real bachelierPrice, Real discount)
Definition: blackformula.cpp:798
QuantLib::blackFormulaStdDevSecondDerivative
Real blackFormulaStdDevSecondDerivative(Rate strike, Rate forward, Real stdDev, Real discount, Real displacement)
Definition: blackformula.cpp:671
QuantLib::blackFormulaImpliedStdDevApproximation
Real blackFormulaImpliedStdDevApproximation(Option::Type optionType, Real strike, Real forward, Real blackPrice, Real discount, Real displacement)
Definition: blackformula.cpp:148
QuantLib::blackFormulaImpliedStdDevApproximationRS
Real blackFormulaImpliedStdDevApproximationRS(Option::Type type, Real K, Real F, Real marketValue, Real df, Real displacement)
Definition: blackformula.cpp:269
QuantLib::blackFormulaStdDevDerivative
Real blackFormulaStdDevDerivative(Rate strike, Rate forward, Real stdDev, Real discount, Real displacement)
Definition: blackformula.cpp:638
newtonsafe.hpp
Safe (bracketed) Newton 1-D solver.
normaldistribution.hpp
normal, cumulative and inverse cumulative distributions
q
ext::shared_ptr< YieldTermStructure > q
Definition: perturbativebarrieroptionengine.cpp:1464
v
ext::shared_ptr< BlackVolTermStructure > v
Definition: perturbativebarrieroptionengine.cpp:1487
beta
Real beta
Definition: sabr.cpp:200
F
Real F
Definition: sabr.cpp:200
nu
Real nu
Definition: sabr.cpp:200
alpha
Real alpha
Definition: sabr.cpp:200
Generated by Doxygen 1.9.5