/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
/*
* This file is part of the LibreOffice project.
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* Copyright (C) 2012 Tino Kluge <tino.kluge@hrz.tu-chemnitz.de>
*
*/
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <rtl/math.hxx>
#include "black_scholes.hxx"
// options prices and greeks in the Black-Scholes model
// also known as TV (theoretical value)
// the code is structured as follows:
// (1) basic assets
// - cash-or-nothing option: bincash()
// - asset-or-nothing option: binasset()
// (2) derived basic assets, can all be priced based on (1)
// - vanilla put/call: putcall() = +/- ( binasset() - K*bincash() )
// - truncated put/call (barriers active at maturity only)
// (3) write a wrapper function to include all vanilla prices
// - this is so we don't duplicate code when pricing barriers
// as this is derived from vanillas
// (4) single barrier options (knock-out), priced based on truncated vanillas
// - it follows from the reflection principle that the price W(S) of a
// single barrier option is given by
// W(S) = V(S) - (B/S)^a V(B^2/S), a = 2(rd-rf)/vol^2 - 1
// where V(S) is the price of the corresponding truncated vanilla
// option
// - to reduce code duplication and in anticipation of double barrier
// options we write the following function
// barrier_term(S,c) = V(c*S) - (B/S)^a V(c*B^2/S)
// (5) double barrier options (knock-out)
// - value is an infinite sum over option prices of the corresponding
// truncated vanillas (truncated at both barriers):
// W(S)=sum (B2/B1)^(i*a) (V(S(B2/B1)^(2i)) - (B1/S)^a V(B1^2/S (B2/B1)^(2i))
// (6) write routines for put/call barriers and touch options which
// mainly call the general double barrier pricer
// the main routines are touch() and barrier()
// both can price in/out barriers, double/single barriers as well as
// vanillas
// the framework allows any barriers to be priced as long as we define
// the value/greek functions for the corresponding truncated vanilla
// and wrap them into internal::vanilla() and internal::vanilla_trunc()
// disadvantage of that approach is that due to the rules of
// differentiations the formulas for greeks become long and possible
// simplifications in the formulas won't be made
// other code inefficiency due to multiplication with pm (+/- 1)
// cvtsi2sd: int-->double, 6/3 cycles
// mulsd: double-double multiplication, 5/1 cycles
// with -O3, however, it compiles 2 versions with pm=1, and pm=-1
// which are efficient
// note this is tiny anyway as compared to exp/log (100 cycles),
// pow (200 cycles), erf (70 cycles)
// this code is not tested for numerical instability, ie overruns,
// underruns, accuracy, etc
namespace sca::pricing::bs {
// helper functions
static double sqr(
double x) {
return x*x;
}
// normal density (see also ScInterpreter::phi)
static double dnorm(
double x) {
//return (1.0/sqrt(2.0*M_PI))*exp(-0.5*x*x); // windows may not have M_PI
return 0.
39894228040143268*exp(-
0.
5*x*x);
}
// cumulative normal distribution (see also ScInterpreter::integralPhi)
static double pnorm(
double x) {
return 0.
5 * std::erfc(-x * M_SQRT1_2);
}
// binary option cash (domestic)
// call - pays 1 if S_T is above strike K
// put - pays 1 if S_T is below strike K
double bincash(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
types::PutCall pc, types::Greeks greeks) {
assert(tau>=
0.
0);
assert(S>
0.
0);
assert(vol>
0.
0);
assert(K>=
0.
0);
double val=
0.
0;
if(tau<=
0.
0) {
// special case tau=0 (expiry)
switch(greeks) {
case types::Value:
if( (pc==types::Call && S>=K) || (pc==types::Put && S<=K) ) {
val =
1.
0;
}
else {
val =
0.
0;
}
break;
default:
val =
0.
0;
}
}
else if(K==
0.
0) {
// special case with zero strike
if(pc==types::Put) {
// up-and-out (put) with K=0
val=
0.
0;
}
else {
// down-and-out (call) with K=0 (zero coupon bond)
switch(greeks) {
case types::Value:
val =
1.
0;
break;
case types::Theta:
val = rd;
break;
case types::Rho_d:
val = -tau;
break;
default:
val =
0.
0;
}
}
}
else {
// standard case with K>0, tau>0
double d1 = ( log(S/K)+(rd-rf+
0.
5*vol*vol)*tau ) / (vol*sqrt(tau));
double d2 = d1 - vol*sqrt(tau);
int pm = (pc==types::Call) ?
1 : -
1;
switch(greeks) {
case types::Value:
val = pnorm(pm*d2);
break;
case types::Delta:
val = pm*dnorm(d2)/(S*vol*sqrt(tau));
break;
case types::Gamma:
val = -pm*dnorm(d2)*d1/(sqr(S*vol)*tau);
break;
case types::Theta:
val = rd*pnorm(pm*d2)
+ pm*dnorm(d2)*(log(S/K)/(vol*sqrt(tau))-
0.
5*d2)/tau;
break;
case types::Vega:
val = -pm*dnorm(d2)*d1/vol;
break;
case types::Volga:
val = pm*dnorm(d2)/(vol*vol)*(-d1*d1*d2+d1+d2);
break;
case types::Vanna:
val = pm*dnorm(d2)/(S*vol*vol*sqrt(tau))*(d1*d2-
1.
0);
break;
case types::Rho_d:
val = -tau*pnorm(pm*d2) + pm*dnorm(d2)*sqrt(tau)/vol;
break;
case types::Rho_f:
val = -pm*dnorm(d2)*sqrt(tau)/vol;
break;
default:
printf(
"bincash: greek %d not implemented\n", greeks );
abort();
}
}
return exp(-rd*tau)*val;
}
// binary option asset (foreign)
// call - pays S_T if S_T is above strike K
// put - pays S_T if S_T is below strike K
double binasset(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
types::PutCall pc, types::Greeks greeks) {
assert(tau>=
0.
0);
assert(S>
0.
0);
assert(vol>
0.
0);
assert(K>=
0.
0);
double val=
0.
0;
if(tau<=
0.
0) {
// special case tau=0 (expiry)
switch(greeks) {
case types::Value:
if( (pc==types::Call && S>=K) || (pc==types::Put && S<=K) ) {
val = S;
}
else {
val =
0.
0;
}
break;
case types::Delta:
if( (pc==types::Call && S>=K) || (pc==types::Put && S<=K) ) {
val =
1.
0;
}
else {
val =
0.
0;
}
break;
default:
val =
0.
0;
}
}
else if(K==
0.
0) {
// special case with zero strike (forward with zero strike)
if(pc==types::Put) {
// up-and-out (put) with K=0
val =
0.
0;
}
else {
// down-and-out (call) with K=0 (type of forward)
switch(greeks) {
case types::Value:
val = S;
break;
case types::Delta:
val =
1.
0;
break;
case types::Theta:
val = rf*S;
break;
case types::Rho_f:
val = -tau*S;
break;
default:
val =
0.
0;
}
}
}
else {
// normal case
double d1 = ( log(S/K)+(rd-rf+
0.
5*vol*vol)*tau ) / (vol*sqrt(tau));
double d2 = d1 - vol*sqrt(tau);
int pm = (pc==types::Call) ?
1 : -
1;
switch(greeks) {
case types::Value:
val = S*pnorm(pm*d1);
break;
case types::Delta:
val = pnorm(pm*d1) + pm*dnorm(d1)/(vol*sqrt(tau));
break;
case types::Gamma:
val = -pm*dnorm(d1)*d2/(S*sqr(vol)*tau);
break;
case types::Theta:
val = rf*S*pnorm(pm*d1)
+ pm*S*dnorm(d1)*(log(S/K)/(vol*sqrt(tau))-
0.
5*d1)/tau;
break;
case types::Vega:
val = -pm*S*dnorm(d1)*d2/vol;
break;
case types::Volga:
val = pm*S*dnorm(d1)/(vol*vol)*(-d1*d2*d2+d1+d2);
break;
case types::Vanna:
val = pm*dnorm(d1)/(vol*vol*sqrt(tau))*(d2*d2-
1.
0);
break;
case types::Rho_d:
val = pm*S*dnorm(d1)*sqrt(tau)/vol;
break;
case types::Rho_f:
val = -tau*S*pnorm(pm*d1) - pm*S*dnorm(d1)*sqrt(tau)/vol;
break;
default:
printf(
"binasset: greek %d not implemented\n", greeks );
abort();
}
}
return exp(-rf*tau)*val;
}
// just for convenience we can combine bincash and binasset into
// one function binary
// using bincash() if fd==types::Domestic
// using binasset() if fd==types::Foreign
static double binary(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
types::PutCall pc, types::ForDom fd,
types::Greeks greek) {
double val=
0.
0;
switch(fd) {
case types::Domestic:
val = bincash(S,vol,rd,rf,tau,K,pc,greek);
break;
case types::Foreign:
val = binasset(S,vol,rd,rf,tau,K,pc,greek);
break;
default:
// never get here
assert(
false);
}
return val;
}
// further wrapper to combine single/double barrier binary options
// into one function
// B1<=0 - it is assumed lower barrier not set
// B2<=0 - it is assumed upper barrier not set
static double binary(
double S,
double vol,
double rd,
double rf,
double tau,
double B1,
double B2,
types::ForDom fd, types::Greeks greek) {
assert(tau>=
0.
0);
assert(S>
0.
0);
assert(vol>
0.
0);
double val=
0.
0;
if(B1<=
0.
0 && B2<=
0.
0) {
// no barriers set, payoff 1.0 (domestic) or S_T (foreign)
val = binary(S,vol,rd,rf,tau,
0.
0,types::Call,fd,greek);
}
else if(B1<=
0.
0 && B2>
0.
0) {
// upper barrier (put)
val = binary(S,vol,rd,rf,tau,B2,types::Put,fd,greek);
}
else if(B1>
0.
0 && B2<=
0.
0) {
// lower barrier (call)
val = binary(S,vol,rd,rf,tau,B1,types::Call,fd,greek);
}
else if(B1>
0.
0 && B2>
0.
0) {
// double barrier
if(B2<=B1) {
val =
0.
0;
}
else {
val = binary(S,vol,rd,rf,tau,B2,types::Put,fd,greek)
- binary(S,vol,rd,rf,tau,B1,types::Put,fd,greek);
}
}
else {
// never get here
assert(
false);
}
return val;
}
// vanilla put/call option
// call pays (S_T-K)^+
// put pays (K-S_T)^+
// this is the same as: +/- (binasset - K*bincash)
double putcall(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
types::PutCall putcall, types::Greeks greeks) {
assert(tau>=
0.
0);
assert(S>
0.
0);
assert(vol>
0.
0);
assert(K>=
0.
0);
double val =
0.
0;
int pm = (putcall==types::Call) ?
1 : -
1;
if(K==
0 || tau==
0.
0) {
// special cases, simply refer to binasset() and bincash()
val = pm * ( binasset(S,vol,rd,rf,tau,K,putcall,greeks)
- K*bincash(S,vol,rd,rf,tau,K,putcall,greeks) );
}
else {
// general case
// we could just use pm*(binasset-K*bincash), however
// since the formula for delta and gamma simplify we write them
// down here
double d1 = ( log(S/K)+(rd-rf+
0.
5*vol*vol)*tau ) / (vol*sqrt(tau));
double d2 = d1 - vol*sqrt(tau);
switch(greeks) {
case types::Value:
val = pm * ( exp(-rf*tau)*S*pnorm(pm*d1)-exp(-rd*tau)*K*pnorm(pm*d2) );
break;
case types::Delta:
val = pm*exp(-rf*tau)*pnorm(pm*d1);
break;
case types::Gamma:
val = exp(-rf*tau)*dnorm(d1)/(S*vol*sqrt(tau));
break;
default:
// too lazy for the other greeks, so simply refer to binasset/bincash
val = pm * ( binasset(S,vol,rd,rf,tau,K,putcall,greeks)
- K*bincash(S,vol,rd,rf,tau,K,putcall,greeks) );
}
}
return val;
}
// truncated put/call option, single barrier
// need to specify whether it's down-and-out or up-and-out
// regular (keeps monotonicity): down-and-out for call, up-and-out for put
// reverse (destroys monoton): up-and-out for call, down-and-out for put
// call pays (S_T-K)^+
// put pays (K-S_T)^+
double putcalltrunc(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
double B,
types::PutCall pc, types::KOType kotype,
types::Greeks greeks) {
assert(tau>=
0.
0);
assert(S>
0.
0);
assert(vol>
0.
0);
assert(K>=
0.
0);
assert(B>=
0.
0);
int pm = (pc==types::Call) ?
1 : -
1;
double val =
0.
0;
switch(kotype) {
case types::Regular:
if( (pc==types::Call && B<=K) || (pc==types::Put && B>=K) ) {
// option degenerates to standard plain vanilla call/put
val = putcall(S,vol,rd,rf,tau,K,pc,greeks);
}
else {
// normal case with truncation
val = pm * ( binasset(S,vol,rd,rf,tau,B,pc,greeks)
- K*bincash(S,vol,rd,rf,tau,B,pc,greeks) );
}
break;
case types::Reverse:
if( (pc==types::Call && B<=K) || (pc==types::Put && B>=K) ) {
// option degenerates to zero payoff
val =
0.
0;
}
else {
// normal case with truncation
val = binasset(S,vol,rd,rf,tau,K,types::Call,greeks)
- binasset(S,vol,rd,rf,tau,B,types::Call,greeks)
- K * ( bincash(S,vol,rd,rf,tau,K,types::Call,greeks)
- bincash(S,vol,rd,rf,tau,B,types::Call,greeks) );
}
break;
default:
assert(
false);
}
return val;
}
// wrapper function for put/call option which combines
// double/single/no truncation barrier
// B1<=0 - assume no lower barrier
// B2<=0 - assume no upper barrier
double putcalltrunc(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
double B1,
double B2,
types::PutCall pc, types::Greeks greek) {
assert(tau>=
0.
0);
assert(S>
0.
0);
assert(vol>
0.
0);
assert(K>=
0.
0);
double val=
0.
0;
if(B1<=
0.
0 && B2<=
0.
0) {
// no barriers set, plain vanilla
val = putcall(S,vol,rd,rf,tau,K,pc,greek);
}
else if(B1<=
0.
0 && B2>
0.
0) {
// upper barrier: reverse barrier for call, regular barrier for put
if(pc==types::Call) {
val = putcalltrunc(S,vol,rd,rf,tau,K,B2,pc,types::Reverse,greek);
}
else {
val = putcalltrunc(S,vol,rd,rf,tau,K,B2,pc,types::Regular,greek);
}
}
else if(B1>
0.
0 && B2<=
0.
0) {
// lower barrier: regular barrier for call, reverse barrier for put
if(pc==types::Call) {
val = putcalltrunc(S,vol,rd,rf,tau,K,B1,pc,types::Regular,greek);
}
else {
val = putcalltrunc(S,vol,rd,rf,tau,K,B1,pc,types::Reverse,greek);
}
}
else if(B1>
0.
0 && B2>
0.
0) {
// double barrier
if(B2<=B1) {
val =
0.
0;
}
else {
int pm = (pc==types::Call) ?
1 : -
1;
val = pm * (
putcalltrunc(S,vol,rd,rf,tau,K,B1,pc,types::Regular,greek)
- putcalltrunc(S,vol,rd,rf,tau,K,B2,pc,types::Regular,greek)
);
}
}
else {
// never get here
assert(
false);
}
return val;
}
namespace internal {
// wrapper function for all non-path dependent options
// this is only an internal function, used to avoid code duplication when
// going to path-dependent barrier options,
// K<0 - assume binary option
// K>=0 - assume put/call option
static double vanilla(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
double B1,
double B2,
types::PutCall pc, types::ForDom fd,
types::Greeks greek) {
double val =
0.
0;
if(K<
0.
0) {
// binary option if K<0
val = binary(S,vol,rd,rf,tau,B1,B2,fd,greek);
}
else {
val = putcall(S,vol,rd,rf,tau,K,pc,greek);
}
return val;
}
static double vanilla_trunc(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
double B1,
double B2,
types::PutCall pc, types::ForDom fd,
types::Greeks greek) {
double val =
0.
0;
if(K<
0.
0) {
// binary option if K<0
// truncated is actually the same as the vanilla binary
val = binary(S,vol,rd,rf,tau,B1,B2,fd,greek);
}
else {
val = putcalltrunc(S,vol,rd,rf,tau,K,B1,B2,pc,greek);
}
return val;
}
}
// namespace internal
// path dependent options
namespace internal {
// helper term for any type of options with continuously monitored barriers,
// internal, should not be called from outside
// calculates value and greeks based on
// V(S) = V1(sc*S) - (B/S)^a V1(sc*B^2/S)
// (a=2 mu/vol^2, mu drift in logspace, ie. mu=(rd-rf-1/2vol^2))
// with sc=1 and V1() being the price of the respective truncated
// vanilla option, V() would be the price of the respective barrier
// option if only one barrier is present
static double barrier_term(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
double B1,
double B2,
double sc,
types::PutCall pc, types::ForDom fd,
types::Greeks greek) {
assert(tau>=
0.
0);
assert(S>
0.
0);
assert(vol>
0.
0);
// V(S) = V1(sc*S) - (B/S)^a V1(sc*B^2/S)
double val =
0.
0;
double B = (B1>
0.
0) ? B1 : B2;
double a =
2.
0*(rd-rf)/(vol*vol)-
1.
0;
// helper variable
double b =
4.
0*(rd-rf)/(vol*vol*vol);
// helper variable -da/dvol
double c =
12.
0*(rd-rf)/(vol*vol*vol*vol);
// helper -db/dvol
switch(greek) {
case types::Value:
case types::Theta:
val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
- pow(B/S,a)*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
break;
case types::Delta:
val = sc*vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
+ pow(B/S,a) * (
a/S*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
+ sqr(B/S)*sc*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
);
break;
case types::Gamma:
val = sc*sc*vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
- pow(B/S,a) * (
a*(a+
1.
0)/(S*S)*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
+ (
2.
0*a+
2.
0)*B*B/(S*S*S)*sc*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Delta)
+ sqr(sqr(B/S))*sc*sc*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Gamma)
);
break;
case types::Vega:
val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
- pow(B/S,a) * (
- b*log(B/S)*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
+
1.
0*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
);
break;
case types::Volga:
val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
- pow(B/S,a) * (
log(B/S)*(b*b*log(B/S)+c)*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
-
2.
0*b*log(B/S)*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Vega)
+
1.
0*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Volga)
);
break;
case types::Vanna:
val = sc*vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
- pow(B/S,a) * (
b/S*(log(B/S)*a+
1.
0)*
vanilla_trunc(B*B/S*sc,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
+ b*log(B/S)*sqr(B/S)*sc*
vanilla_trunc(B*B/S*sc,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Delta)
- a/S*
vanilla_trunc(B*B/S*sc,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Vega)
- sqr(B/S)*sc*
vanilla_trunc(B*B/S*sc,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Vanna)
);
break;
case types::Rho_d:
val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
- pow(B/S,a) * (
2.
0*log(B/S)/(vol*vol)*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
+
1.
0*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
);
break;
case types::Rho_f:
val = vanilla_trunc(sc*S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
- pow(B/S,a) * (
-
2.
0*log(B/S)/(vol*vol)*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,types::Value)
+
1.
0*
vanilla_trunc(sc*B*B/S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
);
break;
default:
printf(
"barrier_term: greek %d not implemented\n", greek );
abort();
}
return val;
}
// one term of the infinite sum for the valuation of double barriers
static double barrier_double_term(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
double B1,
double B2,
double fac,
double sc,
int i,
types::PutCall pc, types::ForDom fd, types::Greeks greek) {
double val =
0.
0;
double b =
4.
0*i*(rd-rf)/(vol*vol*vol);
// helper variable -da/dvol
double c =
12.
0*i*(rd-rf)/(vol*vol*vol*vol);
// helper -db/dvol
switch(greek) {
case types::Value:
val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek);
break;
case types::Delta:
val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek);
break;
case types::Gamma:
val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek);
break;
case types::Theta:
val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek);
break;
case types::Vega:
val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek)
- b*log(B2/B1)*fac *
barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Value);
break;
case types::Volga:
val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek)
-
2.
0*b*log(B2/B1)*fac *
barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Vega)
+ log(B2/B1)*fac*(c+b*b*log(B2/B1)) *
barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Value);
break;
case types::Vanna:
val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek)
- b*log(B2/B1)*fac *
barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Delta);
break;
case types::Rho_d:
val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek)
+
2.
0*i/(vol*vol)*log(B2/B1)*fac *
barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Value);
break;
case types::Rho_f:
val = fac*barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,greek)
-
2.
0*i/(vol*vol)*log(B2/B1)*fac *
barrier_term(S,vol,rd,rf,tau,K,B1,B2,sc,pc,fd,types::Value);
break;
default:
printf(
"barrier_double_term: greek %d not implemented\n", greek );
abort();
}
return val;
}
// general pricer for any type of options with continuously monitored barriers
// allows two, one or zero barriers, only knock-out style
// payoff profiles allowed based on vanilla_trunc()
static double barrier_ko(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
double B1,
double B2,
types::PutCall pc, types::ForDom fd,
types::Greeks greek) {
assert(tau>=
0.
0);
assert(S>
0.
0);
assert(vol>
0.
0);
double val =
0.
0;
if(B1<=
0.
0 && B2<=
0.
0) {
// no barriers --> vanilla case
val = vanilla(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
}
else if(B1>
0.
0 && B2<=
0.
0) {
// lower barrier
if(S<=B1) {
val =
0.
0;
// knocked out
}
else {
val = barrier_term(S,vol,rd,rf,tau,K,B1,B2,
1.
0,pc,fd,greek);
}
}
else if(B1<=
0.
0 && B2>
0.
0) {
// upper barrier
if(S>=B2) {
val =
0.
0;
// knocked out
}
else {
val = barrier_term(S,vol,rd,rf,tau,K,B1,B2,
1.
0,pc,fd,greek);
}
}
else if(B1>
0.
0 && B2>
0.
0) {
// double barrier
if(S<=B1 || S>=B2) {
val =
0.
0;
// knocked out (always true if wrong input B1>B2)
}
else {
// more complex calculation as we have to evaluate an infinite
// sum
// to reduce very costly pow() calls we define some variables
double a =
2.
0*(rd-rf)/(vol*vol)-
1.
0;
// 2 (mu-1/2vol^2)/sigma^2
double BB2=sqr(B2/B1);
double BBa=pow(B2/B1,a);
double BB2inv=
1.
0/BB2;
double BBainv=
1.
0/BBa;
double fac=
1.
0;
double facinv=
1.
0;
double sc=
1.
0;
double scinv=
1.
0;
// initial term i=0
val=barrier_double_term(S,vol,rd,rf,tau,K,B1,B2,fac,sc,
0,pc,fd,greek);
// infinite loop, 10 should be plenty, normal would be 2
for(
int i=
1; i<
10; i++) {
fac*=BBa;
facinv*=BBainv;
sc*=BB2;
scinv*=BB2inv;
double add =
barrier_double_term(S,vol,rd,rf,tau,K,B1,B2,fac,sc,i,pc,fd,greek) +
barrier_double_term(S,vol,rd,rf,tau,K,B1,B2,facinv,scinv,-i,pc,fd,greek);
val += add;
//printf("%i: val=%e (add=%e)\n",i,val,add);
if(fabs(add) <=
1e-
12*fabs(val)) {
break;
}
}
// not knocked-out double barrier end
}
// double barrier end
}
else {
// no such barrier combination exists
assert(
false);
}
return val;
}
// knock-in style barrier
static double barrier_ki(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
double B1,
double B2,
types::PutCall pc, types::ForDom fd,
types::Greeks greek) {
return vanilla(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
-barrier_ko(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
}
// general barrier
static double barrier(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
double B1,
double B2,
types::PutCall pc, types::ForDom fd,
types::BarrierKIO kio, types::BarrierActive bcont,
types::Greeks greek) {
double val =
0.
0;
if( kio==types::KnockOut && bcont==types::Maturity ) {
// truncated vanilla option
val = vanilla_trunc(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
}
else if ( kio==types::KnockOut && bcont==types::Continuous ) {
// standard knock-out barrier
val = barrier_ko(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
}
else if ( kio==types::KnockIn && bcont==types::Maturity ) {
// inverse truncated vanilla
val = vanilla(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek)
- vanilla_trunc(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
}
else if ( kio==types::KnockIn && bcont==types::Continuous ) {
// standard knock-in barrier
val = barrier_ki(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
}
else {
// never get here
assert(
false);
}
return val;
}
}
// namespace internal
// touch/no-touch options (cash/asset or nothing payoff profile)
double touch(
double S,
double vol,
double rd,
double rf,
double tau,
double B1,
double B2, types::ForDom fd,
types::BarrierKIO kio, types::BarrierActive bcont,
types::Greeks greek) {
double K=-
1.
0;
// dummy
types::PutCall pc = types::Call;
// dummy
double val =
0.
0;
if( kio==types::KnockOut && bcont==types::Maturity ) {
// truncated vanilla option
val = internal::vanilla_trunc(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
}
else if ( kio==types::KnockOut && bcont==types::Continuous ) {
// standard knock-out barrier
val = internal::barrier_ko(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
}
else if ( kio==types::KnockIn && bcont==types::Maturity ) {
// inverse truncated vanilla
val = internal::vanilla(S,vol,rd,rf,tau,K,-
1.
0,-
1.
0,pc,fd,greek)
- internal::vanilla_trunc(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
}
else if ( kio==types::KnockIn && bcont==types::Continuous ) {
// standard knock-in barrier
val = internal::vanilla(S,vol,rd,rf,tau,K,-
1.
0,-
1.
0,pc,fd,greek)
- internal::barrier_ko(S,vol,rd,rf,tau,K,B1,B2,pc,fd,greek);
}
else {
// never get here
assert(
false);
}
return val;
}
// barrier option (put/call payoff profile)
double barrier(
double S,
double vol,
double rd,
double rf,
double tau,
double K,
double B1,
double B2,
double rebate,
types::PutCall pc, types::BarrierKIO kio,
types::BarrierActive bcont,
types::Greeks greek) {
assert(tau>=
0.
0);
assert(S>
0.
0);
assert(vol>
0.
0);
assert(K>=
0.
0);
types::ForDom fd = types::Domestic;
double val=internal::barrier(S,vol,rd,rf,tau,K,B1,B2,pc,fd,kio,bcont,greek);
if(rebate!=
0.
0) {
// opposite of barrier knock-in/out type
types::BarrierKIO kio2 = (kio==types::KnockIn) ? types::KnockOut
: types::KnockIn;
val += rebate*touch(S,vol,rd,rf,tau,B1,B2,fd,kio2,bcont,greek);
}
return val;
}
// probability of hitting a barrier
// this is almost the same as the price of a touch option (domestic)
// as it pays one if a barrier is hit; we only have to offset the
// discounting and we get the probability
double prob_hit(
double S,
double vol,
double mu,
double tau,
double B1,
double B2) {
double const rd=
0.
0;
double rf=-mu;
return 1.
0 - touch(S,vol,rd,rf,tau,B1,B2,types::Domestic,types::KnockOut,
types::Continuous, types::Value);
}
// probability of being in-the-money, ie payoff is greater zero,
// assuming payoff(S_T) > 0 iff S_T in [B1, B2]
// this the same as the price of a cash or nothing option
// with no discounting
double prob_in_money(
double S,
double vol,
double mu,
double tau,
double B1,
double B2) {
assert(S>
0.
0);
assert(vol>
0.
0);
assert(tau>=
0.
0);
double val =
0.
0;
if( B1<B2 || B1<=
0.
0 || B2<=
0.
0 ) {
val = binary(S,vol,
0.
0,-mu,tau,B1,B2,types::Domestic,types::Value);
}
return val;
}
double prob_in_money(
double S,
double vol,
double mu,
double tau,
double K,
double B1,
double B2,
types::PutCall pc) {
assert(S>
0.
0);
assert(vol>
0.
0);
assert(tau>=
0.
0);
// if K<0 we assume a binary option is given
if(K<
0.
0) {
return prob_in_money(S,vol,mu,tau,B1,B2);
}
double val =
0.
0;
double BM1, BM2;
// range of in the money [BM1, BM2]
// non-sense parameters with no positive payoff
if( (B1>B2 && B1>
0.
0 && B2>
0.
0) ||
(K>=B2 && B2>
0.
0 && pc==types::Call) ||
(K<=B1 && pc==types::Put) ) {
val =
0.
0;
// need to figure out between what barriers payoff is greater 0
}
else if(pc==types::Call) {
BM1=std::max(B1, K);
BM2=B2;
val = prob_in_money(S,vol,mu,tau,BM1,BM2);
}
else if (pc==types::Put) {
BM1=B1;
BM2= (B2>
0.
0) ? std::min(B2,K) : K;
val = prob_in_money(S,vol,mu,tau,BM1,BM2);
}
else {
// don't get here
assert(
false);
}
return val;
}
}
// namespace sca
/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
