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