QuantLib: a free/open-source library for quantitative finance
Fully annotated sources - version 1.32
Loading...
Searching...
No Matches
bjerksundstenslandengine.cpp
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>
24#include <ql/math/functional.hpp>
25#include <ql/math/comparison.hpp>
26#include <ql/math/distributions/normaldistribution.hpp>
27#include <ql/pricingengines/blackcalculator.hpp>
28#include <ql/pricingengines/vanilla/bjerksundstenslandengine.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
113 BjerksundStenslandApproximationEngine::BjerksundStenslandApproximationEngine(
114 ext::shared_ptr<GeneralizedBlackScholesProcess> process)
115 : process_(std::move(process)) {
116 registerWith(process_);
117 }
118
119 OneAssetOption::results
120 BjerksundStenslandApproximationEngine::europeanCallResults(
121 Real S, Real X, Real rfD, Real dD, Real variance) const {
122
123 OneAssetOption::results results;
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
158 OneAssetOption::results
159 BjerksundStenslandApproximationEngine::immediateExercise(Real S, Real X) const {
160 OneAssetOption::results results;
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
177 OneAssetOption::results
178 BjerksundStenslandApproximationEngine::americanCallApproximation(
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
184 OneAssetOption::results results;
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) {
202 results = immediateExercise(S, X);
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
432 void BjerksundStenslandApproximationEngine::calculate() const {
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) {
471 results_ = europeanCallResults(
472 spot, strike, riskFreeDiscount, dividendDiscount, variance);
473 } else {
474 // early exercise can be optimal - use approximation
475 results_ = americanCallApproximation(
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}
QuantLib::BjerksundStenslandApproximationEngine::BjerksundStenslandApproximationEngine
BjerksundStenslandApproximationEngine(ext::shared_ptr< GeneralizedBlackScholesProcess >)
Definition: bjerksundstenslandengine.cpp:113
QuantLib::BjerksundStenslandApproximationEngine::americanCallApproximation
OneAssetOption::results americanCallApproximation(Real S, Real X, Real rfD, Real dD, Real variance) const
Definition: bjerksundstenslandengine.cpp:178
QuantLib::BjerksundStenslandApproximationEngine::europeanCallResults
OneAssetOption::results europeanCallResults(Real S, Real X, Real rfD, Real dD, Real variance) const
Definition: bjerksundstenslandengine.cpp:120
QuantLib::BjerksundStenslandApproximationEngine::immediateExercise
OneAssetOption::results immediateExercise(Real S, Real X) const
Definition: bjerksundstenslandengine.cpp:159
QuantLib::BjerksundStenslandApproximationEngine::process_
ext::shared_ptr< GeneralizedBlackScholesProcess > process_
Definition: bjerksundstenslandengine.hpp:53
QuantLib::BlackCalculator
Black 1976 calculator class.
Definition: blackcalculator.hpp:37
QuantLib::BlackCalculator::dividendRho
Real dividendRho(Time maturity) const
Definition: blackcalculator.cpp:328
QuantLib::BlackCalculator::delta
virtual Real delta(Real spot) const
Definition: blackcalculator.cpp:202
QuantLib::BlackCalculator::vega
Real vega(Time maturity) const
Definition: blackcalculator.cpp:301
QuantLib::BlackCalculator::gamma
virtual Real gamma(Real spot) const
Definition: blackcalculator.cpp:255
QuantLib::BlackCalculator::thetaPerDay
virtual Real thetaPerDay(Real spot, Time maturity) const
Definition: blackcalculator.hpp:120
QuantLib::BlackCalculator::theta
virtual Real theta(Real spot, Time maturity) const
Definition: blackcalculator.cpp:290
QuantLib::BlackCalculator::rho
Real rho(Time maturity) const
Definition: blackcalculator.cpp:316
QuantLib::Date
Concrete date class.
Definition: date.hpp:125
QuantLib::DayCounter
day counter class
Definition: daycounter.hpp:44
QuantLib::DayCounter::yearFraction
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
QuantLib::OneAssetOption::results
Results from single-asset option calculation
Definition: oneassetoption.hpp:71
QuantLib::Days
@ Days
Definition: timeunit.hpp:37
QL_EPSILON
#define QL_EPSILON
Definition: qldefines.hpp:178
QuantLib::Time
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QuantLib::Real
QL_REAL Real
real number
Definition: types.hpp:50
QuantLib::DiscountFactor
Real DiscountFactor
discount factor between dates
Definition: types.hpp:66
QuantLib::Volatility
Real Volatility
volatility
Definition: types.hpp:78
QuantLib
Definition: any.hpp:35
QuantLib::squared
T squared(T x)
Definition: functional.hpp:37
std
STL namespace.
Generated by Doxygen 1.9.5