LCOV - code coverage report
Current view: top level - scaddins/source/pricing - black_scholes.cxx (source / functions) Hit Total Coverage
Test: commit c8344322a7af75b84dd3ca8f78b05543a976dfd5 Lines: 0 413 0.0 %
Date: 2015-06-13 12:38:46 Functions: 0 22 0.0 %
Legend: Lines: hit not hit

          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: */

Generated by: LCOV version 1.11