QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
baroneadesiwhaleyengine.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) 2003, 2006 Ferdinando Ametrano
5 Copyright (C) 2007 StatPro Italia srl
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#include <ql/exercise.hpp>
27#include <utility>
28
29namespace QuantLib {
30
32 ext::shared_ptr<GeneralizedBlackScholesProcess> process)
33 : process_(std::move(process)) {
34 registerWith(process_);
35 }
36
37
38 // critical commodity price
40 const ext::shared_ptr<StrikedTypePayoff>& payoff,
41 DiscountFactor riskFreeDiscount,
42 DiscountFactor dividendDiscount,
43 Real variance, Real tolerance) {
44
45 // Calculation of seed value, Si
46 Real n= 2.0*std::log(dividendDiscount/riskFreeDiscount)/(variance);
47 Real m=-2.0*std::log(riskFreeDiscount)/(variance);
48 Real bT = std::log(dividendDiscount/riskFreeDiscount);
49
50 Real qu, Su, h, Si;
51 switch (payoff->optionType()) {
52 case Option::Call:
53 qu = (-(n-1.0) + std::sqrt(((n-1.0)*(n-1.0)) + 4.0*m))/2.0;
54 Su = payoff->strike() / (1.0 - 1.0/qu);
55 h = -(bT + 2.0*std::sqrt(variance)) * payoff->strike() /
56 (Su - payoff->strike());
57 Si = payoff->strike() + (Su - payoff->strike()) *
58 (1.0 - std::exp(h));
59 break;
60 case Option::Put:
61 qu = (-(n-1.0) - std::sqrt(((n-1.0)*(n-1.0)) + 4.0*m))/2.0;
62 Su = payoff->strike() / (1.0 - 1.0/qu);
63 h = (bT - 2.0*std::sqrt(variance)) * payoff->strike() /
64 (payoff->strike() - Su);
65 Si = Su + (payoff->strike() - Su) * std::exp(h);
66 break;
67 default:
68 QL_FAIL("unknown option type");
69 }
70
71
72 // Newton Raphson algorithm for finding critical price Si
73 Real Q, LHS, RHS, bi;
74 Real forwardSi = Si * dividendDiscount / riskFreeDiscount;
75 Real d1 = (std::log(forwardSi/payoff->strike()) + 0.5*variance) /
76 std::sqrt(variance);
77 CumulativeNormalDistribution cumNormalDist;
78 Real K = (!close(riskFreeDiscount, 1.0, 1000))
79 ? Real(-2.0*std::log(riskFreeDiscount)
80 / (variance*(1.0-riskFreeDiscount)))
81 : 2.0/variance;
82 Real temp = blackFormula(payoff->optionType(), payoff->strike(),
83 forwardSi, std::sqrt(variance))*riskFreeDiscount;
84 switch (payoff->optionType()) {
85 case Option::Call:
86 Q = (-(n-1.0) + std::sqrt(((n-1.0)*(n-1.0)) + 4 * K)) / 2;
87 LHS = Si - payoff->strike();
88 RHS = temp + (1 - dividendDiscount * cumNormalDist(d1)) * Si / Q;
89 bi = dividendDiscount * cumNormalDist(d1) * (1 - 1/Q) +
90 (1 - dividendDiscount *
91 cumNormalDist.derivative(d1) / std::sqrt(variance)) / Q;
92 while (std::fabs(LHS - RHS)/payoff->strike() > tolerance) {
93 Si = (payoff->strike() + RHS - bi * Si) / (1 - bi);
94 forwardSi = Si * dividendDiscount / riskFreeDiscount;
95 d1 = (std::log(forwardSi/payoff->strike())+0.5*variance)
96 /std::sqrt(variance);
97 LHS = Si - payoff->strike();
98 Real temp2 = blackFormula(payoff->optionType(), payoff->strike(),
99 forwardSi, std::sqrt(variance))*riskFreeDiscount;
100 RHS = temp2 + (1 - dividendDiscount * cumNormalDist(d1)) * Si / Q;
101 bi = dividendDiscount * cumNormalDist(d1) * (1 - 1 / Q)
102 + (1 - dividendDiscount *
103 cumNormalDist.derivative(d1) / std::sqrt(variance))
104 / Q;
105 }
106 break;
107 case Option::Put:
108 Q = (-(n-1.0) - std::sqrt(((n-1.0)*(n-1.0)) + 4 * K)) / 2;
109 LHS = payoff->strike() - Si;
110 RHS = temp - (1 - dividendDiscount * cumNormalDist(-d1)) * Si / Q;
111 bi = -dividendDiscount * cumNormalDist(-d1) * (1 - 1/Q)
112 - (1 + dividendDiscount * cumNormalDist.derivative(-d1)
113 / std::sqrt(variance)) / Q;
114 while (std::fabs(LHS - RHS)/payoff->strike() > tolerance) {
115 Si = (payoff->strike() - RHS + bi * Si) / (1 + bi);
116 forwardSi = Si * dividendDiscount / riskFreeDiscount;
117 d1 = (std::log(forwardSi/payoff->strike())+0.5*variance)
118 /std::sqrt(variance);
119 LHS = payoff->strike() - Si;
120 Real temp2 = blackFormula(payoff->optionType(), payoff->strike(),
121 forwardSi, std::sqrt(variance))*riskFreeDiscount;
122 RHS = temp2 - (1 - dividendDiscount * cumNormalDist(-d1)) * Si / Q;
123 bi = -dividendDiscount * cumNormalDist(-d1) * (1 - 1 / Q)
124 - (1 + dividendDiscount * cumNormalDist.derivative(-d1)
125 / std::sqrt(variance)) / Q;
126 }
127 break;
128 default:
129 QL_FAIL("unknown option type");
130 }
131
132 return Si;
133 }
134
136
137 QL_REQUIRE(arguments_.exercise->type() == Exercise::American,
138 "not an American Option");
139
140 ext::shared_ptr<AmericanExercise> ex =
141 ext::dynamic_pointer_cast<AmericanExercise>(arguments_.exercise);
142 QL_REQUIRE(ex, "non-American exercise given");
143 QL_REQUIRE(!ex->payoffAtExpiry(),
144 "payoff at expiry not handled");
145
146 ext::shared_ptr<StrikedTypePayoff> payoff =
147 ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
148 QL_REQUIRE(payoff, "non-striked payoff given");
149
150 Real variance = process_->blackVolatility()->blackVariance(
151 ex->lastDate(), payoff->strike());
152 DiscountFactor dividendDiscount = process_->dividendYield()->discount(
153 ex->lastDate());
154 DiscountFactor riskFreeDiscount = process_->riskFreeRate()->discount(
155 ex->lastDate());
156 Real spot = process_->stateVariable()->value();
157 QL_REQUIRE(spot > 0.0, "negative or null underlying given");
158 Real forwardPrice = spot * dividendDiscount / riskFreeDiscount;
159 BlackCalculator black(payoff, forwardPrice, std::sqrt(variance),
160 riskFreeDiscount);
161
162 if (dividendDiscount>=1.0 && payoff->optionType()==Option::Call) {
163 // early exercise never optimal
164 results_.value = black.value();
165 results_.delta = black.delta(spot);
166 results_.deltaForward = black.deltaForward();
167 results_.elasticity = black.elasticity(spot);
168 results_.gamma = black.gamma(spot);
169
170 DayCounter rfdc = process_->riskFreeRate()->dayCounter();
171 DayCounter divdc = process_->dividendYield()->dayCounter();
172 DayCounter voldc = process_->blackVolatility()->dayCounter();
173 Time t =
174 rfdc.yearFraction(process_->riskFreeRate()->referenceDate(),
175 arguments_.exercise->lastDate());
176 results_.rho = black.rho(t);
177
178 t = divdc.yearFraction(process_->dividendYield()->referenceDate(),
179 arguments_.exercise->lastDate());
180 results_.dividendRho = black.dividendRho(t);
181
182 t = voldc.yearFraction(process_->blackVolatility()->referenceDate(),
183 arguments_.exercise->lastDate());
184 results_.vega = black.vega(t);
185 results_.theta = black.theta(spot, t);
186 results_.thetaPerDay = black.thetaPerDay(spot, t);
187
188 results_.strikeSensitivity = black.strikeSensitivity();
189 results_.itmCashProbability = black.itmCashProbability();
190 } else {
191 // early exercise can be optimal
192 CumulativeNormalDistribution cumNormalDist;
193 Real tolerance = 1e-6;
194 Real Sk = criticalPrice(payoff, riskFreeDiscount,
195 dividendDiscount, variance, tolerance);
196 Real forwardSk = Sk * dividendDiscount / riskFreeDiscount;
197 Real d1 = (std::log(forwardSk/payoff->strike()) + 0.5*variance)
198 /std::sqrt(variance);
199 Real n = 2.0*std::log(dividendDiscount/riskFreeDiscount)/variance;
200 Real K = (!close(riskFreeDiscount, 1.0, 1000))
201 ? Real(-2.0*std::log(riskFreeDiscount)
202 / (variance*(1.0-riskFreeDiscount)))
203 : 2.0/variance;
204 Real Q, a;
205 switch (payoff->optionType()) {
206 case Option::Call:
207 Q = (-(n-1.0) + std::sqrt(((n-1.0)*(n-1.0))+4.0*K))/2.0;
208 a = (Sk/Q) * (1.0 - dividendDiscount * cumNormalDist(d1));
209 if (spot<Sk) {
210 results_.value = black.value() +
211 a * std::pow((spot/Sk), Q);
212 } else {
213 results_.value = spot - payoff->strike();
214 }
215 break;
216 case Option::Put:
217 Q = (-(n-1.0) - std::sqrt(((n-1.0)*(n-1.0))+4.0*K))/2.0;
218 a = -(Sk/Q) *
219 (1.0 - dividendDiscount * cumNormalDist(-d1));
220 if (spot>Sk) {
221 results_.value = black.value() +
222 a * std::pow((spot/Sk), Q);
223 } else {
224 results_.value = payoff->strike() - spot;
225 }
226 break;
227 default:
228 QL_FAIL("unknown option type");
229 }
230 } // end of "early exercise can be optimal"
231 }
232
233}
Barone-Adesi and Whaley approximation engine.
Black-formula calculator class.
Black formula.
const Instrument::results * results_
Definition: cdsoption.cpp:63
BaroneAdesiWhaleyApproximationEngine(ext::shared_ptr< GeneralizedBlackScholesProcess >)
static Real criticalPrice(const ext::shared_ptr< StrikedTypePayoff > &payoff, DiscountFactor riskFreeDiscount, DiscountFactor dividendDiscount, Real variance, Real tolerance=1e-6)
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
virtual Real gamma(Real spot) const
virtual Real elasticity(Real spot) const
virtual Real thetaPerDay(Real spot, Time maturity) const
virtual Real theta(Real spot, Time maturity) const
Real rho(Time maturity) const
Cumulative normal distribution function.
day counter class
Definition: daycounter.hpp:44
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.
Definition: daycounter.hpp:128
floating-point comparisons
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Option exercise classes and payoff function.
LinearInterpolation variance
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
Real DiscountFactor
discount factor between dates
Definition: types.hpp:66
ext::shared_ptr< QuantLib::Payoff > payoff
Definition: any.hpp:35
bool close(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:163
Real blackFormula(Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount, Real displacement)
STL namespace.
normal, cumulative and inverse cumulative distributions