Line data Source code
1 : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : /*
3 : * This file is part of the LibreOffice project.
4 : *
5 : * This Source Code Form is subject to the terms of the Mozilla Public
6 : * License, v. 2.0. If a copy of the MPL was not distributed with this
7 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 : *
9 : * Copyright (C) 2012 Tino Kluge <tino.kluge@hrz.tu-chemnitz.de>
10 : *
11 : */
12 :
13 : #include <cstdio>
14 : #include <cstdlib>
15 : #include <cmath>
16 : #include <cassert>
17 : #include <algorithm>
18 : #include <rtl/math.hxx>
19 : #include "black_scholes.hxx"
20 :
21 : // options prices and greeks in the Black-Scholes model
22 : // also known as TV (theoretical value)
23 :
24 : // the code is structured as follows:
25 :
26 : // (1) basic assets
27 : // - cash-or-nothing option: bincash()
28 : // - asset-or-nothing option: binasset()
29 :
30 : // (2) derived basic assets, can all be priced based on (1)
31 : // - vanilla put/call: putcall() = +/- ( binasset() - K*bincash() )
32 : // - truncated put/call (barriers active at maturity only)
33 :
34 : // (3) write a wrapper function to include all vanilla pricers
35 : // - this is so we don't duplicate code when pricing barriers
36 : // as this is derived from vanillas
37 :
38 : // (4) single barrier options (knock-out), priced based on truncated vanillas
39 : // - it follows from the reflection principle that the price W(S) of a
40 : // single barrier option is given by
41 : // W(S) = V(S) - (B/S)^a V(B^2/S), a = 2(rd-rf)/vol^2 - 1
42 : // where V(S) is the price of the corresponding truncated vanilla
43 : // option
44 : // - to reduce code duplication and in anticipation of double barrier
45 : // options we write the following function
46 : // barrier_term(S,c) = V(c*S) - (B/S)^a V(c*B^2/S)
47 :
48 : // (5) double barrier options (knock-out)
49 : // - value is an infinite sum over option prices of the corresponding
50 : // truncated vanillas (truncated at both barriers):
51 :
52 : // W(S)=sum (B2/B1)^(i*a) (V(S(B2/B1)^(2i)) - (B1/S)^a V(B1^2/S (B2/B1)^(2i))
53 :
54 : // (6) write routines for put/call barriers and touch options which
55 : // mainly call the general double barrier pricer
56 : // the main routines are touch() and barrier()
57 : // both can price in/out barriers, double/single barriers as well as
58 : // vanillas
59 :
60 :
61 : // the framework allows any barriers to be priced as long as we define
62 : // the value/greek functions for the corresponding truncated vanilla
63 : // and wrap them into internal::vanilla() and internal::vanilla_trunc()
64 :
65 : // disadvantage of that approach is that due to the rules of
66 : // differentiations the formulas for greeks become long and possible
67 : // simplifications in the formulas won't be made
68 :
69 : // other code inefficiency due to multiplication with pm (+/- 1)
70 : // cvtsi2sd: int-->double, 6/3 cycles
71 : // mulsd: double-double multiplication, 5/1 cycles
72 : // with -O3, however, it compiles 2 versions with pm=1, and pm=-1
73 : // which are efficient
74 : // note this is tiny anyway as compared to exp/log (100 cycles),
75 : // pow (200 cycles), erf (70 cycles)
76 :
77 : // this code is not tested for numerical instability, ie overruns,
78 : // underruns, accuracy, etc
79 :
80 :
81 : namespace sca {
82 : namespace pricing {
83 :
84 : namespace bs {
85 :
86 :
87 : // helper functions
88 :
89 0 : inline double sqr(double x) {
90 0 : return x*x;
91 : }
92 : // normal density (see also ScInterpreter::phi)
93 0 : inline double dnorm(double x) {
94 : //return (1.0/sqrt(2.0*M_PI))*exp(-0.5*x*x); // windows may not have M_PI
95 0 : return 0.39894228040143268*exp(-0.5*x*x);
96 : }
97 : // cumulative normal distribution (see also ScInterpreter::integralPhi)
98 0 : inline double pnorm(double x) {
99 : //return 0.5*(erf(sqrt(0.5)*x)+1.0); // windows may not have erf
100 0 : return 0.5 * ::rtl::math::erfc(-x * 0.7071067811865475);
101 : }
102 :
103 : // binary option cash (domestic)
104 : // call - pays 1 if S_T is above strike K
105 : // put - pays 1 if S_T is below strike K
106 0 : double bincash(double S, double vol, double rd, double rf,
107 : double tau, double K,
108 : types::PutCall pc, types::Greeks greeks) {
109 : assert(tau>=0.0);
110 : assert(S>0.0);
111 : assert(vol>0.0);
112 : assert(K>=0.0);
113 :
114 0 : double val=0.0;
115 :
116 0 : if(tau<=0.0) {
117 : // special case tau=0 (expiry)
118 0 : switch(greeks) {
119 : case types::Value:
120 0 : if( (pc==types::Call && S>=K) || (pc==types::Put && S<=K) ) {
121 0 : val = 1.0;
122 : } else {
123 0 : val = 0.0;
124 : }
125 0 : break;
126 : default:
127 0 : val = 0.0;
128 : }
129 0 : } else if(K==0.0) {
130 : // special case with zero strike
131 0 : if(pc==types::Put) {
132 : // up-and-out (put) with K=0
133 0 : val=0.0;
134 : } else {
135 : // down-and-out (call) with K=0 (zero coupon bond)
136 0 : switch(greeks) {
137 : case types::Value:
138 0 : val = 1.0;
139 0 : break;
140 : case types::Theta:
141 0 : val = rd;
142 0 : break;
143 : case types::Rho_d:
144 0 : val = -tau;
145 0 : break;
146 : default:
147 0 : val = 0.0;
148 : }
149 : }
150 : } else {
151 : // standard case with K>0, tau>0
152 0 : double d1 = ( log(S/K)+(rd-rf+0.5*vol*vol)*tau ) / (vol*sqrt(tau));
153 0 : double d2 = d1 - vol*sqrt(tau);
154 0 : int pm = (pc==types::Call) ? 1 : -1;
155 :
156 0 : switch(greeks) {
157 : case types::Value:
158 0 : val = pnorm(pm*d2);
159 0 : break;
160 : case types::Delta:
161 0 : val = pm*dnorm(d2)/(S*vol*sqrt(tau));
162 0 : break;
163 : case types::Gamma:
164 0 : val = -pm*dnorm(d2)*d1/(sqr(S*vol)*tau);
165 0 : break;
166 : case types::Theta:
167 0 : val = rd*pnorm(pm*d2)
168 0 : + pm*dnorm(d2)*(log(S/K)/(vol*sqrt(tau))-0.5*d2)/tau;
169 0 : break;
170 : case types::Vega:
171 0 : val = -pm*dnorm(d2)*d1/vol;
172 0 : break;
173 : case types::Volga:
174 0 : val = pm*dnorm(d2)/(vol*vol)*(-d1*d1*d2+d1+d2);
175 0 : break;
176 : case types::Vanna:
177 0 : val = pm*dnorm(d2)/(S*vol*vol*sqrt(tau))*(d1*d2-1.0);
178 0 : break;
179 : case types::Rho_d:
180 0 : val = -tau*pnorm(pm*d2) + pm*dnorm(d2)*sqrt(tau)/vol;
181 0 : break;
182 : case types::Rho_f:
183 0 : val = -pm*dnorm(d2)*sqrt(tau)/vol;
184 0 : break;
185 : default:
186 0 : printf("bincash: greek %d not implemented\n", greeks );
187 0 : abort();
188 : }
189 : }
190 0 : return exp(-rd*tau)*val;
191 : }
192 :
193 : // binary option asset (foreign)
194 : // call - pays S_T if S_T is above strike K
195 : // put - pays S_T if S_T is below strike K
196 0 : double binasset(double S, double vol, double rd, double rf,
197 : double tau, double K,
198 : types::PutCall pc, types::Greeks greeks) {
199 : assert(tau>=0.0);
200 : assert(S>0.0);
201 : assert(vol>0.0);
202 : assert(K>=0.0);
203 :
204 0 : double val=0.0;
205 0 : if(tau<=0.0) {
206 : // special case tau=0 (expiry)
207 0 : switch(greeks) {
208 : case types::Value:
209 0 : if( (pc==types::Call && S>=K) || (pc==types::Put && S<=K) ) {
210 0 : val = S;
211 : } else {
212 0 : val = 0.0;
213 : }
214 0 : break;
215 : case types::Delta:
216 0 : if( (pc==types::Call && S>=K) || (pc==types::Put && S<=K) ) {
217 0 : val = 1.0;
218 : } else {
219 0 : val = 0.0;
220 : }
221 0 : break;
222 : default:
223 0 : val = 0.0;
224 : }
225 0 : } else if(K==0.0) {
226 : // special case with zero strike (forward with zero strike)
227 0 : if(pc==types::Put) {
228 : // up-and-out (put) with K=0
229 0 : val = 0.0;
230 : } else {
231 : // down-and-out (call) with K=0 (type of forward)
232 0 : switch(greeks) {
233 : case types::Value:
234 0 : val = S;
235 0 : break;
236 : case types::Delta:
237 0 : val = 1.0;
238 0 : break;
239 : case types::Theta:
240 0 : val = rf*S;
241 0 : break;
242 : case types::Rho_f:
243 0 : val = -tau*S;
244 0 : break;
245 : default:
246 0 : val = 0.0;
247 : }
248 : }
249 : } else {
250 : // normal case
251 0 : double d1 = ( log(S/K)+(rd-rf+0.5*vol*vol)*tau ) / (vol*sqrt(tau));
252 0 : double d2 = d1 - vol*sqrt(tau);
253 0 : int pm = (pc==types::Call) ? 1 : -1;
254 :
255 0 : switch(greeks) {
256 : case types::Value:
257 0 : val = S*pnorm(pm*d1);
258 0 : break;
259 : case types::Delta:
260 0 : val = pnorm(pm*d1) + pm*dnorm(d1)/(vol*sqrt(tau));
261 0 : break;
262 : case types::Gamma:
263 0 : val = -pm*dnorm(d1)*d2/(S*sqr(vol)*tau);
264 0 : break;
265 : case types::Theta:
266 0 : val = rf*S*pnorm(pm*d1)
267 0 : + pm*S*dnorm(d1)*(log(S/K)/(vol*sqrt(tau))-0.5*d1)/tau;
268 0 : break;
269 : case types::Vega:
270 0 : val = -pm*S*dnorm(d1)*d2/vol;
271 0 : break;
272 : case types::Volga:
273 0 : val = pm*S*dnorm(d1)/(vol*vol)*(-d1*d2*d2+d1+d2);
274 0 : break;
275 : case types::Vanna:
276 0 : val = pm*dnorm(d1)/(vol*vol*sqrt(tau))*(d2*d2-1.0);
277 0 : break;
278 : case types::Rho_d:
279 0 : val = pm*S*dnorm(d1)*sqrt(tau)/vol;
280 0 : break;
281 : case types::Rho_f:
282 0 : val = -tau*S*pnorm(pm*d1) - pm*S*dnorm(d1)*sqrt(tau)/vol;
283 0 : break;
284 : default:
285 0 : printf("binasset: greek %d not implemented\n", greeks );
286 0 : abort();
287 : }
288 : }
289 0 : return exp(-rf*tau)*val;
290 : }
291 :
292 : // just for convenience we can combine bincash and binasset into
293 : // one function binary
294 : // using bincash() if fd==types::Domestic
295 : // using binasset() if fd==types::Foreign
296 0 : double binary(double S, double vol, double rd, double rf,
297 : double tau, double K,
298 : types::PutCall pc, types::ForDom fd,
299 : types::Greeks greek) {
300 0 : double val=0.0;
301 0 : switch(fd) {
302 : case types::Domestic:
303 0 : val = bincash(S,vol,rd,rf,tau,K,pc,greek);
304 0 : break;
305 : case types::Foreign:
306 0 : val = binasset(S,vol,rd,rf,tau,K,pc,greek);
307 0 : break;
308 : default:
309 : // never get here
310 : assert(false);
311 : }
312 0 : return val;
313 : }
314 :
315 : // further wrapper to combine single/double barrier binary options
316 : // into one function
317 : // B1<=0 - it is assumed lower barrier not set
318 : // B2<=0 - it is assumed upper barrier not set
319 0 : double binary(double S, double vol, double rd, double rf,
320 : double tau, double B1, double B2,
321 : types::ForDom fd, types::Greeks greek) {
322 : assert(tau>=0.0);
323 : assert(S>0.0);
324 : assert(vol>0.0);
325 :
326 0 : double val=0.0;
327 :
328 0 : if(B1<=0.0 && B2<=0.0) {
329 : // no barriers set, payoff 1.0 (domestic) or S_T (foreign)
330 0 : val = binary(S,vol,rd,rf,tau,0.0,types::Call,fd,greek);
331 0 : } else if(B1<=0.0 && B2>0.0) {
332 : // upper barrier (put)
333 0 : val = binary(S,vol,rd,rf,tau,B2,types::Put,fd,greek);
334 0 : } else if(B1>0.0 && B2<=0.0) {
335 : // lower barrier (call)
336 0 : val = binary(S,vol,rd,rf,tau,B1,types::Call,fd,greek);
337 0 : } else if(B1>0.0 && B2>0.0) {
338 : // double barrier
339 0 : if(B2<=B1) {
340 0 : val = 0.0;
341 : } else {
342 0 : val = binary(S,vol,rd,rf,tau,B2,types::Put,fd,greek)
343 0 : - binary(S,vol,rd,rf,tau,B1,types::Put,fd,greek);
344 : }
345 : } else {
346 : // never get here
347 : assert(false);
348 : }
349 :
350 0 : return val;
351 : }
352 :
353 : // vanilla put/call option
354 : // call pays (S_T-K)^+
355 : // put pays (K-S_T)^+
356 : // this is the same as: +/- (binasset - K*bincash)
357 0 : double putcall(double S, double vol, double rd, double rf,
358 : double tau, double K,
359 : types::PutCall putcall, types::Greeks greeks) {
360 :
361 : assert(tau>=0.0);
362 : assert(S>0.0);
363 : assert(vol>0.0);
364 : assert(K>=0.0);
365 :
366 0 : double val = 0.0;
367 0 : int pm = (putcall==types::Call) ? 1 : -1;
368 :
369 0 : if(K==0 || tau==0.0) {
370 : // special cases, simply refer to binasset() and bincash()
371 0 : val = pm * ( binasset(S,vol,rd,rf,tau,K,putcall,greeks)
372 0 : - K*bincash(S,vol,rd,rf,tau,K,putcall,greeks) );
373 : } else {
374 : // general case
375 : // we could just use pm*(binasset-K*bincash), however
376 : // since the formula for delta and gamma simplify we write them
377 : // down here
378 0 : double d1 = ( log(S/K)+(rd-rf+0.5*vol*vol)*tau ) / (vol*sqrt(tau));
379 0 : double d2 = d1 - vol*sqrt(tau);
380 :
381 0 : switch(greeks) {
382 : case types::Value:
383 0 : val = pm * ( exp(-rf*tau)*S*pnorm(pm*d1)-exp(-rd*tau)*K*pnorm(pm*d2) );
384 0 : break;
385 : case types::Delta:
386 0 : val = pm*exp(-rf*tau)*pnorm(pm*d1);
387 0 : break;
388 : case types::Gamma:
389 0 : val = exp(-rf*tau)*dnorm(d1)/(S*vol*sqrt(tau));
390 0 : break;
391 : default:
392 : // too lazy for the other greeks, so simply refer to binasset/bincash
393 0 : val = pm * ( binasset(S,vol,rd,rf,tau,K,putcall,greeks)
394 0 : - K*bincash(S,vol,rd,rf,tau,K,putcall,greeks) );
395 : }
396 : }
397 0 : return val;
398 : }
399 :
400 : // truncated put/call option, single barrier
401 : // need to specify whether it's down-and-out or up-and-out
402 : // regular (keeps monotonicity): down-and-out for call, up-and-out for put
403 : // reverse (destroys monoton): up-and-out for call, down-and-out for put
404 : // call pays (S_T-K)^+
405 : // put pays (K-S_T)^+
406 0 : double putcalltrunc(double S, double vol, double rd, double rf,
407 : double tau, double K, double B,
408 : types::PutCall pc, types::KOType kotype,
409 : types::Greeks greeks) {
410 :
411 : assert(tau>=0.0);
412 : assert(S>0.0);
413 : assert(vol>0.0);
414 : assert(K>=0.0);
415 : assert(B>=0.0);
416 :
417 0 : int pm = (pc==types::Call) ? 1 : -1;
418 0 : double val = 0.0;
419 :
420 0 : switch(kotype) {
421 : case types::Regular:
422 0 : if( (pc==types::Call && B<=K) || (pc==types::Put && B>=K) ) {
423 : // option degenerates to standard plain vanilla call/put
424 0 : val = putcall(S,vol,rd,rf,tau,K,pc,greeks);
425 : } else {
426 : // normal case with truncation
427 0 : val = pm * ( binasset(S,vol,rd,rf,tau,B,pc,greeks)
428 0 : - K*bincash(S,vol,rd,rf,tau,B,pc,greeks) );
429 : }
430 0 : break;
431 : case types::Reverse:
432 0 : if( (pc==types::Call && B<=K) || (pc==types::Put && B>=K) ) {
433 : // option degenerates to zero payoff
434 0 : val = 0.0;
435 : } else {
436 : // normal case with truncation
437 0 : val = binasset(S,vol,rd,rf,tau,K,types::Call,greeks)
438 0 : - binasset(S,vol,rd,rf,tau,B,types::Call,greeks)
439 0 : - K * ( bincash(S,vol,rd,rf,tau,K,types::Call,greeks)
440 0 : - bincash(S,vol,rd,rf,tau,B,types::Call,greeks) );
441 : }
442 0 : break;
443 : default:
444 : assert(false);
445 : }
446 0 : return val;
447 : }
448 :
449 : // wrapper function for put/call option which combines
450 : // double/single/no truncation barrier
451 : // B1<=0 - assume no lower barrier
452 : // B2<=0 - assume no upper barrier
453 0 : double putcalltrunc(double S, double vol, double rd, double rf,
454 : double tau, double K, double B1, double B2,
455 : types::PutCall pc, types::Greeks greek) {
456 :
457 : assert(tau>=0.0);
458 : assert(S>0.0);
459 : assert(vol>0.0);
460 : assert(K>=0.0);
461 :
462 0 : double val=0.0;
463 :
464 0 : if(B1<=0.0 && B2<=0.0) {
465 : // no barriers set, plain vanilla
466 0 : val = putcall(S,vol,rd,rf,tau,K,pc,greek);
467 0 : } else if(B1<=0.0 && B2>0.0) {
468 : // upper barrier: reverse barrier for call, regular barrier for put
469 0 : if(pc==types::Call) {
470 0 : val = putcalltrunc(S,vol,rd,rf,tau,K,B2,pc,types::Reverse,greek);
471 : } else {
472 0 : val = putcalltrunc(S,vol,rd,rf,tau,K,B2,pc,types::Regular,greek);
473 : }
474 0 : } else if(B1>0.0 && B2<=0.0) {
475 : // lower barrier: regular barrier for call, reverse barrier for put
476 0 : if(pc==types::Call) {
477 0 : val = putcalltrunc(S,vol,rd,rf,tau,K,B1,pc,types::Regular,greek);
478 : } else {
479 0 : val = putcalltrunc(S,vol,rd,rf,tau,K,B1,pc,types::Reverse,greek);
480 : }
481 0 : } else if(B1>0.0 && B2>0.0) {
482 : // double barrier
483 0 : if(B2<=B1) {
484 0 : val = 0.0;
485 : } else {
486 0 : int pm = (pc==types::Call) ? 1 : -1;
487 0 : val = pm * (
488 0 : putcalltrunc(S,vol,rd,rf,tau,K,B1,pc,types::Regular,greek)
489 0 : - putcalltrunc(S,vol,rd,rf,tau,K,B2,pc,types::Regular,greek)
490 0 : );
491 : }
492 : } else {
493 : // never get here
494 : assert(false);
495 : }
496 0 : return val;
497 : }
498 :
499 : namespace internal {
500 :
501 : // wrapper function for all non-path dependent options
502 : // this is only an internal function, used to avoid code duplication when
503 : // going to path-dependent barrier options,
504 : // K<0 - assume binary option
505 : // K>=0 - assume put/call option
506 0 : double vanilla(double S, double vol, double rd, double rf,
507 : double tau, double K, double B1, double B2,
508 : types::PutCall pc, types::ForDom fd,
509 : types::Greeks greek) {
510 0 : double val = 0.0;
511 0 : if(K<0.0) {
512 : // binary option if K<0
513 0 : val = binary(S,vol,rd,rf,tau,B1,B2,fd,greek);
514 : } else {
515 0 : val = putcall(S,vol,rd,rf,tau,K,pc,greek);
516 : }
517 0 : return val;
518 : }
519 0 : double vanilla_trunc(double S, double vol, double rd, double rf,
520 : double tau, double K, double B1, double B2,
521 : types::PutCall pc, types::ForDom fd,
522 : types::Greeks greek) {
523 0 : double val = 0.0;
524 0 : if(K<0.0) {
525 : // binary option if K<0
526 : // truncated is actually the same as the vanilla binary
527 0 : val = binary(S,vol,rd,rf,tau,B1,B2,fd,greek);
528 : } else {
529 0 : val = putcalltrunc(S,vol,rd,rf,tau,K,B1,B2,pc,greek);
530 : }
531 0 : return val;
532 : }
533 :
534 : } // namespace internal
535 :
536 : // path dependent options
537 :
538 :
539 : namespace internal {
540 :
541 : // helper term for any type of options with continuously monitored barriers,
542 : // internal, should not be called from outside
543 : // calculates value and greeks based on
544 : // V(S) = V1(sc*S) - (B/S)^a V1(sc*B^2/S)
545 : // (a=2 mu/vol^2, mu drift in logspace, ie. mu=(rd-rf-1/2vol^2))
546 : // with sc=1 and V1() being the price of the respective truncated
547 : // vanilla option, V() would be the price of the respective barrier
548 : // option if only one barrier is present
549 0 : double barrier_term(double S, double vol, double rd, double rf,
550 : double tau, double K, double B1, double B2, double sc,
551 : types::PutCall pc, types::ForDom fd,
552 : types::Greeks greek) {
553 :
554 : assert(tau>=0.0);
555 : assert(S>0.0);
556 : assert(vol>0.0);
557 :
558 : // V(S) = V1(sc*S) - (B/S)^a V1(sc*B^2/S)
559 0 : double val = 0.0;
560 0 : double B = (B1>0.0) ? B1 : B2;
561 0 : double a = 2.0*(rd-rf)/(vol*vol)-1.0; // helper variable
562 0 : double b = 4.0*(rd-rf)/(vol*vol*vol); // helper variable -da/dvol
563 0 : double c = 12.0*(rd-rf)/(vol*vol*vol*vol); // helper -db/dvol
564 0 : switch(greek) {
565 : case types::Value:
566 0 : val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
567 0 : - pow(B/S,a)*
568 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
569 0 : break;
570 : case types::Delta:
571 0 : val = sc*vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
572 0 : + pow(B/S,a) * (
573 0 : a/S*
574 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
575 0 : + sqr(B/S)*sc*
576 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
577 0 : );
578 0 : break;
579 : case types::Gamma:
580 0 : val = sc*sc*vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
581 0 : - pow(B/S,a) * (
582 0 : a*(a+1.0)/(S*S)*
583 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
584 0 : + (2.0*a+2.0)*B*B/(S*S*S)*sc*
585 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Delta)
586 0 : + sqr(sqr(B/S))*sc*sc*
587 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Gamma)
588 0 : );
589 0 : break;
590 : case types::Theta:
591 0 : val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
592 0 : - pow(B/S,a)*
593 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
594 0 : break;
595 : case types::Vega:
596 0 : val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
597 0 : - pow(B/S,a) * (
598 0 : - b*log(B/S)*
599 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
600 0 : + 1.0*
601 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
602 0 : );
603 0 : break;
604 : case types::Volga:
605 0 : val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
606 0 : - pow(B/S,a) * (
607 0 : log(B/S)*(b*b*log(B/S)+c)*
608 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
609 0 : - 2.0*b*log(B/S)*
610 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Vega)
611 0 : + 1.0*
612 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Volga)
613 0 : );
614 0 : break;
615 : case types::Vanna:
616 0 : val = sc*vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
617 0 : - pow(B/S,a) * (
618 0 : b/S*(log(B/S)*a+1.0)*
619 0 : vanilla_trunc(B*B/S*sc,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
620 0 : + b*log(B/S)*sqr(B/S)*sc*
621 0 : vanilla_trunc(B*B/S*sc,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Delta)
622 0 : - a/S*
623 0 : vanilla_trunc(B*B/S*sc,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Vega)
624 0 : - sqr(B/S)*sc*
625 0 : vanilla_trunc(B*B/S*sc,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Vanna)
626 0 : );
627 0 : break;
628 : case types::Rho_d:
629 0 : val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
630 0 : - pow(B/S,a) * (
631 0 : 2.0*log(B/S)/(vol*vol)*
632 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
633 0 : + 1.0*
634 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
635 0 : );
636 0 : break;
637 : case types::Rho_f:
638 0 : val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
639 0 : - pow(B/S,a) * (
640 0 : - 2.0*log(B/S)/(vol*vol)*
641 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
642 0 : + 1.0*
643 0 : vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
644 0 : );
645 0 : break;
646 : default:
647 0 : printf("barrier_term: greek %d not implemented\n", greek );
648 0 : abort();
649 : }
650 0 : return val;
651 : }
652 :
653 : // one term of the infinite sum for the valuation of double barriers
654 0 : double barrier_double_term( double S, double vol, double rd, double rf,
655 : double tau, double K, double B1, double B2,
656 : double fac, double sc, int i,
657 : types::PutCall pc, types::ForDom fd, types::Greeks greek) {
658 :
659 0 : double val = 0.0;
660 0 : double b = 4.0*i*(rd-rf)/(vol*vol*vol); // helper variable -da/dvol
661 0 : double c = 12.0*i*(rd-rf)/(vol*vol*vol*vol); // helper -db/dvol
662 0 : switch(greek) {
663 : case types::Value:
664 0 : val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek);
665 0 : break;
666 : case types::Delta:
667 0 : val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek);
668 0 : break;
669 : case types::Gamma:
670 0 : val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek);
671 0 : break;
672 : case types::Theta:
673 0 : val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek);
674 0 : break;
675 : case types::Vega:
676 0 : val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek)
677 0 : - b*log(B2/B1)*fac *
678 0 : barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Value);
679 0 : break;
680 : case types::Volga:
681 0 : val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek)
682 0 : - 2.0*b*log(B2/B1)*fac *
683 0 : barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Vega)
684 0 : + log(B2/B1)*fac*(c+b*b*log(B2/B1)) *
685 0 : barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Value);
686 0 : break;
687 : case types::Vanna:
688 0 : val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek)
689 0 : - b*log(B2/B1)*fac *
690 0 : barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Delta);
691 0 : break;
692 : case types::Rho_d:
693 0 : val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek)
694 0 : + 2.0*i/(vol*vol)*log(B2/B1)*fac *
695 0 : barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Value);
696 0 : break;
697 : case types::Rho_f:
698 0 : val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek)
699 0 : - 2.0*i/(vol*vol)*log(B2/B1)*fac *
700 0 : barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Value);
701 0 : break;
702 : default:
703 0 : printf("barrier_double_term: greek %d not implemented\n", greek );
704 0 : abort();
705 : }
706 0 : return val;
707 : }
708 :
709 : // general pricer for any type of options with continuously monitored barriers
710 : // allows two, one or zero barriers, only knock-out style
711 : // payoff profiles allowed based on vanilla_trunc()
712 0 : double barrier_ko(double S, double vol, double rd, double rf,
713 : double tau, double K, double B1, double B2,
714 : types::PutCall pc, types::ForDom fd,
715 : types::Greeks greek) {
716 :
717 : assert(tau>=0.0);
718 : assert(S>0.0);
719 : assert(vol>0.0);
720 :
721 0 : double val = 0.0;
722 :
723 0 : if(B1<=0.0 && B2<=0.0) {
724 : // no barriers --> vanilla case
725 0 : val = vanilla(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
726 0 : } else if(B1>0.0 && B2<=0.0) {
727 : // lower barrier
728 0 : if(S<=B1) {
729 0 : val = 0.0; // knocked out
730 : } else {
731 0 : val = barrier_term(S,vol,rd,rf,tau,K,B1,B2,1.0,pc,fd,greek);
732 : }
733 0 : } else if(B1<=0.0 && B2>0.0) {
734 : // upper barrier
735 0 : if(S>=B2) {
736 0 : val = 0.0; // knocked out
737 : } else {
738 0 : val = barrier_term(S,vol,rd,rf,tau,K,B1,B2,1.0,pc,fd,greek);
739 : }
740 0 : } else if(B1>0.0 && B2>0.0) {
741 : // double barrier
742 0 : if(S<=B1 || S>=B2) {
743 0 : val = 0.0; // knocked out (always true if wrong input B1>B2)
744 : } else {
745 : // more complex calculation as we have to evaluate an infinite
746 : // sum
747 : // to reduce very costly pow() calls we define some variables
748 0 : double a = 2.0*(rd-rf)/(vol*vol)-1.0; // 2 (mu-1/2vol^2)/sigma^2
749 0 : double BB2=sqr(B2/B1);
750 0 : double BBa=pow(B2/B1,a);
751 0 : double BB2inv=1.0/BB2;
752 0 : double BBainv=1.0/BBa;
753 0 : double fac=1.0;
754 0 : double facinv=1.0;
755 0 : double sc=1.0;
756 0 : double scinv=1.0;
757 :
758 : // initial term i=0
759 0 : val=barrier_double_term(S,vol,rd,rf,tau,K,B1,B2,fac,sc,0,pc,fd,greek);
760 : // infinite loop, 10 should be plenty, normal would be 2
761 0 : for(int i=1; i<10; i++) {
762 0 : fac*=BBa;
763 0 : facinv*=BBainv;
764 0 : sc*=BB2;
765 0 : scinv*=BB2inv;
766 : double add =
767 0 : barrier_double_term(S,vol,rd,rf,tau,K,B1,B2,fac,sc,i,pc,fd,greek) +
768 0 : barrier_double_term(S,vol,rd,rf,tau,K,B1,B2,facinv,scinv,-i,pc,fd,greek);
769 0 : val += add;
770 : //printf("%i: val=%e (add=%e)\n",i,val,add);
771 0 : if(fabs(add) <= 1e-12*fabs(val)) {
772 0 : break;
773 : }
774 : }
775 : // not knocked-out double barrier end
776 : }
777 : // double barrier end
778 : } else {
779 : // no such barrier combination exists
780 : assert(false);
781 : }
782 :
783 0 : return val;
784 : }
785 :
786 : // knock-in style barrier
787 0 : double barrier_ki(double S, double vol, double rd, double rf,
788 : double tau, double K, double B1, double B2,
789 : types::PutCall pc, types::ForDom fd,
790 : types::Greeks greek) {
791 0 : return vanilla(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
792 0 : -barrier_ko(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
793 : }
794 :
795 : // general barrier
796 0 : double barrier(double S, double vol, double rd, double rf,
797 : double tau, double K, double B1, double B2,
798 : types::PutCall pc, types::ForDom fd,
799 : types::BarrierKIO kio, types::BarrierActive bcont,
800 : types::Greeks greek) {
801 :
802 0 : double val = 0.0;
803 0 : if( kio==types::KnockOut && bcont==types::Maturity ) {
804 : // truncated vanilla option
805 0 : val = vanilla_trunc(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
806 0 : } else if ( kio==types::KnockOut && bcont==types::Continuous ) {
807 : // standard knock-out barrier
808 0 : val = barrier_ko(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
809 0 : } else if ( kio==types::KnockIn && bcont==types::Maturity ) {
810 : // inverse truncated vanilla
811 0 : val = vanilla(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
812 0 : - vanilla_trunc(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
813 0 : } else if ( kio==types::KnockIn && bcont==types::Continuous ) {
814 : // standard knock-in barrier
815 0 : val = barrier_ki(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
816 : } else {
817 : // never get here
818 : assert(false);
819 : }
820 0 : return val;
821 : }
822 :
823 : } // namespace internal
824 :
825 :
826 : // touch/no-touch options (cash/asset or nothing payoff profile)
827 0 : double touch(double S, double vol, double rd, double rf,
828 : double tau, double B1, double B2, types::ForDom fd,
829 : types::BarrierKIO kio, types::BarrierActive bcont,
830 : types::Greeks greek) {
831 :
832 0 : double K=-1.0; // dummy
833 0 : types::PutCall pc = types::Call; // dummy
834 0 : double val = 0.0;
835 0 : if( kio==types::KnockOut && bcont==types::Maturity ) {
836 : // truncated vanilla option
837 0 : val = internal::vanilla_trunc(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
838 0 : } else if ( kio==types::KnockOut && bcont==types::Continuous ) {
839 : // standard knock-out barrier
840 0 : val = internal::barrier_ko(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
841 0 : } else if ( kio==types::KnockIn && bcont==types::Maturity ) {
842 : // inverse truncated vanilla
843 0 : val = internal::vanilla(S,vol,rd,rf,tau,K,-1.0,-1.0,pc,fd,greek)
844 0 : - internal::vanilla_trunc(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
845 0 : } else if ( kio==types::KnockIn && bcont==types::Continuous ) {
846 : // standard knock-in barrier
847 0 : val = internal::vanilla(S,vol,rd,rf,tau,K,-1.0,-1.0,pc,fd,greek)
848 0 : - internal::barrier_ko(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
849 : } else {
850 : // never get here
851 : assert(false);
852 : }
853 0 : return val;
854 : }
855 :
856 : // barrier option (put/call payoff profile)
857 0 : double barrier(double S, double vol, double rd, double rf,
858 : double tau, double K, double B1, double B2,
859 : double rebate,
860 : types::PutCall pc, types::BarrierKIO kio,
861 : types::BarrierActive bcont,
862 : types::Greeks greek) {
863 : assert(tau>=0.0);
864 : assert(S>0.0);
865 : assert(vol>0.0);
866 : assert(K>=0.0);
867 0 : types::ForDom fd = types::Domestic;
868 0 : double val=internal::barrier(S,vol,rd,rf,tau,K,B1,B2,pc,fd,kio,bcont,greek);
869 0 : if(rebate!=0.0) {
870 : // opposite of barrier knock-in/out type
871 : types::BarrierKIO kio2 = (kio==types::KnockIn) ? types::KnockOut
872 0 : : types::KnockIn;
873 0 : val += rebate*touch(S,vol,rd,rf,tau,B1,B2,fd,kio2,bcont,greek);
874 : }
875 0 : return val;
876 : }
877 :
878 : // probability of hitting a barrier
879 : // this is almost the same as the price of a touch option (domestic)
880 : // as it pays one if a barrier is hit; we only have to offset the
881 : // discounting and we get the probability
882 0 : double prob_hit(double S, double vol, double mu,
883 : double tau, double B1, double B2) {
884 0 : double rd=0.0;
885 0 : double rf=-mu;
886 : return 1.0 - touch(S,vol,rd,rf,tau,B1,B2,types::Domestic,types::KnockOut,
887 0 : types::Continuous, types::Value);
888 : }
889 :
890 : // probability of being in-the-money, ie payoff is greater zero,
891 : // assuming payoff(S_T) > 0 iff S_T in [B1, B2]
892 : // this the same as the price of a cash or nothing option
893 : // with no discounting
894 0 : double prob_in_money(double S, double vol, double mu,
895 : double tau, double B1, double B2) {
896 : assert(S>0.0);
897 : assert(vol>0.0);
898 : assert(tau>=0.0);
899 0 : double val = 0.0;
900 0 : if( B1<B2 || B1<=0.0 || B2<=0.0 ) {
901 0 : val = binary(S,vol,0.0,-mu,tau,B1,B2,types::Domestic,types::Value);
902 : }
903 0 : return val;
904 : }
905 0 : double prob_in_money(double S, double vol, double mu,
906 : double tau, double K, double B1, double B2,
907 : types::PutCall pc) {
908 : assert(S>0.0);
909 : assert(vol>0.0);
910 : assert(tau>=0.0);
911 :
912 : // if K<0 we assume a binary option is given
913 0 : if(K<0.0) {
914 0 : return prob_in_money(S,vol,mu,tau,B1,B2);
915 : }
916 :
917 0 : double val = 0.0;
918 : double BM1, BM2; // range of in the money [BM1, BM2]
919 : // non-sense parameters with no positive payoff
920 0 : if( (B1>B2 && B1>0.0 && B2>0.0) ||
921 0 : (K>=B2 && B2>0.0 && pc==types::Call) ||
922 0 : (K<=B1 && pc==types::Put) ) {
923 0 : val = 0.0;
924 : // need to figure out between what barriers payoff is greater 0
925 0 : } else if(pc==types::Call) {
926 0 : BM1=std::max(B1, K);
927 0 : BM2=B2;
928 0 : val = prob_in_money(S,vol,mu,tau,BM1,BM2);
929 0 : } else if (pc==types::Put) {
930 0 : BM1=B1;
931 0 : BM2= (B2>0.0) ? std::min(B2,K) : K;
932 0 : val = prob_in_money(S,vol,mu,tau,BM1,BM2);
933 : } else {
934 : // don't get here
935 : assert(false);
936 : }
937 0 : return val;
938 : }
939 :
940 : } // namespace bs
941 :
942 : } // namespace pricing
943 : } // namespace sca
944 :
945 :
946 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
