QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
bjerksundstenslandengine.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 Ferdinando Ametrano
5 Copyright (C) 2007 StatPro Italia srl
6 Copyright (C) 2023 Klaus Spanderen
7
8 This file is part of QuantLib, a free-software/open-source library
9 for financial quantitative analysts and developers - http://quantlib.org/
10
11 QuantLib is free software: you can redistribute it and/or modify it
12 under the terms of the QuantLib license. You should have received a
13 copy of the license along with this program; if not, please email
14 <quantlib-dev@lists.sf.net>. The license is also available online at
15 <http://quantlib.org/license.shtml>.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the license for more details.
20*/
21
22#include <ql/any.hpp>
23#include <ql/exercise.hpp>
29#include <utility>
30#include <cmath>
31
32namespace QuantLib {
33
34 namespace {
35
36 CumulativeNormalDistribution cumNormalDist;
37 Real phi(Real S, Real gamma, Real H, Real I,
38 Real rT, Real bT, Real variance) {
39
40 Real lambda = (-rT + gamma * bT + 0.5 * gamma * (gamma - 1.0)
41 * variance);
42 Real d = -(std::log(S / H) + (bT + (gamma - 0.5) * variance) )
43 / std::sqrt(variance);
44 Real kappa = 2.0 * bT / variance + (2.0 * 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)));
48 }
49
50 Real phi_S(Real S, Real gamma, Real H, Real I, Real rT, Real bT, Real v) {
51 const Real lsh = std::log(S/H);
52 const Real lis = std::log(I/S);
53 const Real sv = std::sqrt(v);
54
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));
56 }
57
58 Real phi_SS(Real S, Real gamma, Real H, Real I, Real rT, Real bT, Real 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));
64
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));
66 }
67
68 Real phi_gamma(Real S, Real gamma, Real H, Real I, Real rT, Real bT, Real v) {
69 const Real lsh = std::log(S/H);
70 const Real lis = std::log(I/S);
71 const Real sv = std::sqrt(v);
72
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.);
74 }
75
76 Real phi_H(Real S, Real gamma, Real H, Real I, Real rT, Real bT, Real v) {
77 const Real lsh = std::log(S/H);
78
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));
80 }
81
82 Real phi_I(Real S, Real gamma, Real H, Real I, Real rT, Real bT, Real v) {
83 const Real lsh = std::log(S/H);
84 const Real lis = std::log(I/S);
85 const Real sv = std::sqrt(v);
86
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);
88 }
89
90 Real phi_rt(Real S, Real gamma, Real H, Real I, Real rT, Real bT, Real v) {
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);
93 }
94
95 Real phi_bt(Real S, Real gamma, Real H, Real I, Real rT, Real bT, Real v) {
96 const Real lsh = std::log(S/H);
97 const Real lis = std::log(I/S);
98 const Real sv = std::sqrt(v);
99
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);
101 }
102
103 Real phi_v(Real S, Real gamma, Real H, Real I, Real rT, Real bT, Real 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));
108
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.;
110 }
111 }
112
114 ext::shared_ptr<GeneralizedBlackScholesProcess> process)
115 : process_(std::move(process)) {
116 registerWith(process_);
117 }
118
121 Real S, Real X, Real rfD, Real dD, Real variance) const {
122
124
125 const Real forwardPrice = S * dD/rfD;
126 const BlackCalculator black(
127 Option::Call, X, forwardPrice, std::sqrt(variance), rfD);
128
129 results.value = black.value();
130 results.delta = black.delta(S);
131 results.gamma = black.gamma(S);
132
133 const DayCounter rfdc = process_->riskFreeRate()->dayCounter();
134 const DayCounter divdc = process_->dividendYield()->dayCounter();
135 const DayCounter voldc = process_->blackVolatility()->dayCounter();
136 Time t =
137 rfdc.yearFraction(process_->riskFreeRate()->referenceDate(),
138 arguments_.exercise->lastDate());
139 results.rho = black.rho(t);
140
141 t = divdc.yearFraction(process_->dividendYield()->referenceDate(),
142 arguments_.exercise->lastDate());
143 results.dividendRho = black.dividendRho(t);
144
145 t = voldc.yearFraction(process_->blackVolatility()->referenceDate(),
146 arguments_.exercise->lastDate());
147 results.vega = black.vega(t);
148 results.theta = black.theta(S, t);
149 results.thetaPerDay = black.thetaPerDay(S, t);
150
151 results.strikeSensitivity = black.strikeSensitivity();
152 results.additionalResults["strikeGamma"] = Real(results.gamma*squared(S/X));
153 results.additionalResults["exerciseType"] = std::string("European");
154
155 return results;
156 }
157
161 results.value = std::max(0.0, S - X);
162 results.delta = (S >= X)? 1.0 : 0.0;
163 results.gamma = 0.0;
164 results.rho = 0.0;
165 results.dividendRho = 0.0;
166 results.vega = 0.0;
167 results.theta = 0.0;
168 results.thetaPerDay = 0.0;
169
170 results.strikeSensitivity = -results.delta;
171 results.additionalResults["strikeGamma"] = Real(0.0);
172 results.additionalResults["exerciseType"] = std::string("Immediate");
173
174 return results;
175 }
176
179 Real S, Real X, Real rfD, Real dD, Real variance) const {
180
181 const OneAssetOption::results europeanResults
182 = europeanCallResults(S, X, rfD, dD, variance);
183
185
186 const Real bT = std::log(dD/rfD);
187 const Real rT = std::log(1.0/rfD);
188
189 const Real beta = (0.5 - bT/variance) +
190 std::sqrt(squared(bT/variance - 0.5) + 2.0 * rT/variance);
191
192 const Real BInfinity = beta / (beta - 1.0) * X;
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);
195
196 const Real I = B0 + (BInfinity - B0) * (1 - std::exp(ht));
197
198 const Real fwd = S * dD/rfD;
199 const Real q = std::log(I/fwd)/std::sqrt(variance);
200
201 if (S >= I) {
203 }
204 else if (q > 12.5) {
205 // We have a run away exercise boundary. It is numerically
206 // more accurate to use the Greeks of the European engine.
207 results = europeanResults;
208 }
209 else {
210 const Real phi_S_beta_I_I_rT_bT_v
211 = phi(S, beta, I, I, rT, bT, variance);
212 const Real phi_S_1_I_I_rT_bT_v
213 = phi(S, 1.0, I, I, rT, bT, variance);
214 const Real phi_S_1_X_I_rT_bT_V
215 = phi(S, 1.0, X, I, rT, bT, variance);
216 results.value = (I - X) * std::pow(S/I, beta)
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);
222
223 const Real phi_S_S_beta_I_I_rT_bT_v
224 = phi_S(S, beta, I, I, rT, bT, variance);
225 const Real phi_S_S_1_I_I_rT_bT_v
226 = phi_S(S, 1.0, I, I, rT, bT, variance);
227 const Real phi_S_S_1_X_I_rT_bT_v
228 = phi_S(S, 1.0, X, I, rT, bT, variance);
229 results.delta = (I - X) * std::pow(S/I, beta-1)*beta/I
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);
239
240 const Date refDate = process_->riskFreeRate()->referenceDate();
241 const Date exerciseDate = arguments_.exercise->lastDate();
242 const DayCounter qdc = process_->dividendYield()->dayCounter();
243 const Time tq = qdc.yearFraction(refDate, exerciseDate);
244
245 const Real betaDq = tq*(1/variance
246 - 1/(2*std::sqrt(squared(bT/variance - 0.5) + 2.0 * rT/variance))
247 * 2*(bT/variance-0.5)/variance);
248 const Real BInfinityDq = -X/squared(beta-1.0)*betaDq;
249 const Real B0Dq = (dD <= rfD) ? Real(0.0)
250 : Real(X*std::log(rfD)/squared(std::log(dD))*tq);
251
252 const Real htDq = tq * B0 / (BInfinity - B0)
253 - (bT + 2.0*std::sqrt(variance))
254 *(B0Dq*(BInfinity - B0) - B0*(BInfinityDq - B0Dq))
255 /squared(BInfinity - B0);
256 const Real IDq = B0Dq + (BInfinityDq - B0Dq) * (1 - std::exp(ht))
257 - (BInfinity - B0) * std::exp(ht)*htDq;
258
259 const Real phi_H_S_beta_I_I_rT_bT_v
260 = phi_H(S, beta, I, I, rT, bT, variance);
261 const Real phi_I_S_beta_I_I_rT_bT_v
262 = phi_I(S, beta, I, I, rT, bT, variance);
263 const Real phi_gamma_S_beta_I_I_rT_bT_v
264 = phi_gamma(S, beta, I, I, rT, bT, variance);
265 const Real phi_bt_S_beta_I_I_rT_bT_v
266 = phi_bt(S, beta, I, I, rT, bT, variance);
267 const Real phi_H_S_1_I_I_rT_bT_v
268 = phi_H(S, 1.0, I, I, rT, bT, variance);
269 const Real phi_I_S_1_I_I_rT_bT_v
270 = phi_I(S, 1.0, I, I, rT, bT, variance);
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
274 = phi_I(S, 1.0, X, I, rT, bT, variance);
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
278 = phi_H(S, 0.0, I, I, rT, bT, variance);
279 const Real phi_I_S_0_I_I_rT_bT_v
280 = phi_I(S, 0.0, I, I, rT, bT, variance);
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
284 = phi_I(S, 0.0, X, I, rT, bT, variance);
285 const Real phi_bt_S_0_X_I_rT_bT_v
286 = phi_bt(S, 0.0, X, I, rT, bT, variance);
287
288 results.dividendRho =
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);
307
308 const DayCounter rdc = process_->riskFreeRate()->dayCounter();
309 const Time tr = rdc.yearFraction(refDate, exerciseDate);
310
311 const Real betaDr = tr*(-1/variance
312 + 1/(2*std::sqrt(squared(bT/variance - 0.5) + 2.0 * rT/variance))
313 * 2*((bT/variance-0.5)/variance + 1/variance));
314 const Real BInfinityDr = -X/squared(beta-1.0)*betaDr;
315 const Real B0Dr = (dD <= rfD) ? Real(0) : Real(-X*tr/std::log(dD));
316 const Real htDr = -tr * B0 / (BInfinity - B0)
317 - (bT + 2.0*std::sqrt(variance))
318 *(B0Dr*(BInfinity - B0) - B0*(BInfinityDr - B0Dr))
319 /squared(BInfinity - B0);
320 const Real IDr = B0Dr + (BInfinityDr - B0Dr) * (1 - std::exp(ht))
321 - (BInfinity - B0) * std::exp(ht)*htDr;
322
323 results.rho =
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
331 + phi_rt(S, beta, I, I, rT, bT, variance)*tr
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);
347
348 const Real beta = (0.5 - bT/variance) +
349 std::sqrt(squared(bT/variance - 0.5) + 2.0 * rT/variance);
350
351 const DayCounter vdc = process_->blackVolatility()->dayCounter();
352 const Time tv = vdc.yearFraction(refDate, exerciseDate);
353 const Real varianceDv = 2*std::sqrt(variance*tv);
354
355 const Real betaDv = bT/squared(variance)*varianceDv +
356 - 1/(2*std::sqrt(squared(bT/variance - 0.5) + 2.0 * rT/variance))
357 *( 2*(bT/variance - 0.5)*bT*varianceDv/squared(variance)
358 +2*rT/squared(variance)*varianceDv );
359 const Real BInfinityDv = -X/squared(beta-1.0)*betaDv;
360 const Real htDv = -1/std::sqrt(variance)*varianceDv*B0/(BInfinity-B0)
361 + (bT + 2*std::sqrt(variance))*B0/squared(BInfinity-B0)*BInfinityDv;
362
363 const Real IDv = BInfinityDv*(1-std::exp(ht))
364 - (BInfinity-B0)*std::exp(ht)*htDv;
365
366 results.vega =
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
374 + phi_v(S, beta, I, I, rT, bT, variance)*varianceDv)
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);
385
386 results.gamma =
387 (I - X) * std::pow(S/I, beta-2)*beta*(beta-1)/squared(I)
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)
392 * phi_SS(S, beta, I, I, rT, bT, variance)
393
394 + 2*phi_S_S_1_I_I_rT_bT_v
395 + S*phi_SS(S, 1.0, I, I, rT, bT, variance)
396
397 - 2*phi_S_S_1_X_I_rT_bT_v
398 - S*phi_SS(S, 1.0, X, I, rT, bT, variance)
399
400 - X*phi_SS(S, 0.0, I, I, rT, bT, variance)
401 + X*phi_SS(S, 0.0, X, I, rT, bT, variance);
402
403 const Volatility vol = std::sqrt(variance/tv);
404
405 const Date tomorrow = refDate + Period(1, Days);
406 const Time dtq = qdc.yearFraction(refDate, exerciseDate)
407 - qdc.yearFraction(tomorrow, exerciseDate);
408 const Time dtr = rdc.yearFraction(refDate, exerciseDate)
409 - rdc.yearFraction(tomorrow, exerciseDate);
410 const Time dtv = vdc.yearFraction(refDate, exerciseDate)
411 - vdc.yearFraction(tomorrow, exerciseDate);
412
413 results.thetaPerDay = -(0.5*results.vega*vol/tv*dtv
414 + results.rho*rT/(tr*tr)*dtr + results.dividendRho*(rT-bT)/(tq*tq)*dtq);
415 results.theta = 365*results.thetaPerDay;
416
417 results.strikeSensitivity = results.value/X - S/X*results.delta;
418 results.additionalResults["strikeGamma"] = Real(results.gamma*squared(S/X));
419
420 results.additionalResults["exerciseType"] = std::string("American");
421 }
422
423 // check if European engine gives higher NPV
424 if (results.value < europeanResults.value) {
425 results = europeanResults;
426 }
427
428 return results;
429 }
430
431
433
434 QL_REQUIRE(arguments_.exercise->type() == Exercise::American,
435 "not an American Option");
436
437 ext::shared_ptr<AmericanExercise> ex =
438 ext::dynamic_pointer_cast<AmericanExercise>(arguments_.exercise);
439 QL_REQUIRE(ex, "non-American exercise given");
440 QL_REQUIRE(!ex->payoffAtExpiry(),
441 "payoff at expiry not handled");
442
443 ext::shared_ptr<PlainVanillaPayoff> payoff =
444 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
445 QL_REQUIRE(payoff, "non-plain payoff given");
446
447 Real variance =
448 process_->blackVolatility()->blackVariance(ex->lastDate(),
449 payoff->strike());
450 DiscountFactor dividendDiscount =
451 process_->dividendYield()->discount(ex->lastDate());
452 DiscountFactor riskFreeDiscount =
453 process_->riskFreeRate()->discount(ex->lastDate());
454 Real spot = process_->stateVariable()->value();
455 QL_REQUIRE(spot > 0.0, "negative or null underlying given");
456 Real strike = payoff->strike();
457
458 if (payoff->optionType()==Option::Put) {
459 // use put-call symmetry
460 std::swap(spot, strike);
461 std::swap(riskFreeDiscount, dividendDiscount);
462 payoff = ext::make_shared<PlainVanillaPayoff>(
463 Option::Call, strike);
464 }
465
466 if (dividendDiscount > 1.0 && riskFreeDiscount > dividendDiscount)
467 QL_FAIL("double-boundary case r<q<0 for a call given");
468
469
470 if (dividendDiscount >= 1.0 && dividendDiscount >= riskFreeDiscount) {
472 spot, strike, riskFreeDiscount, dividendDiscount, variance);
473 } else {
474 // early exercise can be optimal - use approximation
476 spot, strike, riskFreeDiscount, dividendDiscount, variance);
477 }
478
479 // check if immediate exercise gives higher NPV
480 if (results_.value < (spot - strike)*(1+10*QL_EPSILON) ) {
481 results_ = immediateExercise(spot, strike);
482 }
483
484 if (ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff)
485 ->optionType() == Option::Put) {
486
487 std::swap(results_.delta, results_.strikeSensitivity);
488
489 Real tmp = results_.gamma;
490 results_.gamma =
491 ext::any_cast<Real>(results_.additionalResults["strikeGamma"]);
492 results_.additionalResults["strikeGamma"] = tmp;
493
494 std::swap(results_.rho, results_.dividendRho);
495
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());
502
503 results_.rho *= tr/tq;
504 results_.dividendRho *= tq/tr;
505 }
506 }
507}
Maps any to either the boost or std implementation.
Bjerksund and Stensland approximation engine.
Black-formula calculator class.
const Instrument::results * results_
Definition: cdsoption.cpp:63
BjerksundStenslandApproximationEngine(ext::shared_ptr< GeneralizedBlackScholesProcess >)
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
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
Concrete date class.
Definition: date.hpp:125
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
Results from single-asset option calculation
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
Date d
Option exercise classes and payoff function.
LinearInterpolation variance
#define QL_EPSILON
Definition: qldefines.hpp:178
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
Real Volatility
volatility
Definition: types.hpp:78
Real kappa
ext::shared_ptr< QuantLib::Payoff > payoff
functionals and combinators not included in the STL
#define M_SQRTPI
#define M_SQRT2
#define M_PI
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
STL namespace.
normal, cumulative and inverse cumulative distributions
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< BlackVolTermStructure > v
Real beta
Definition: sabr.cpp:200