File Coverage

lib/Math/Business/BlackScholes/Binaries.pm
Criterion Covered Total %
statement 163 168 97.0
branch 47 52 90.3
condition 16 24 66.6
subroutine 24 24 100.0
pod 16 16 100.0
total 266 284 93.6


line stmt bran cond sub pod time code
1             package Math::Business::BlackScholes::Binaries;
2 10     10   663878 use strict;
  10         17  
  10         402  
3 10     10   51 use warnings;
  10         18  
  10         767  
4              
5             our $VERSION = '1.21';
6              
7             my $SMALLTIME = 1 / ( 60 * 60 * 24 * 365 ); # 1 second in years;
8              
9 10     10   50 use List::Util qw(max);
  10         18  
  10         1194  
10 10     10   5553 use Math::CDF qw(pnorm);
  10         29810  
  10         747  
11 10     10   6765 use Math::Trig;
  10         170669  
  10         1890  
12 10     10   6358 use Machine::Epsilon;
  10         3624  
  10         28159  
13              
14             =head1 NAME
15              
16             Math::Business::BlackScholes::Binaries
17              
18             =head1 VERSION
19              
20             Version 1.21
21              
22             =head1 SYNOPSIS
23              
24             use Math::Business::BlackScholes::Binaries;
25              
26             # price of a Call option
27             my $price_call_option = Math::Business::BlackScholes::Binaries::call(
28             1.35, # stock price
29             1.36, # barrier
30             (7/365), # time
31             0.002, # payout currency interest rate (0.05 = 5%)
32             0.001, # quanto drift adjustment (0.05 = 5%)
33             0.11, # volatility (0.3 = 30%)
34             );
35              
36             =head1 DESCRIPTION
37              
38             Prices options using the GBM model, all closed formulas.
39              
40             Important(a): Basically, onetouch, upordown and doubletouch have two cases of
41             payoff either at end or at hit. We treat them differently. We use parameter
42             $w to differ them.
43              
44             $w = 0: payoff at hit time.
45             $w = 1: payoff at end.
46              
47             Our current contracts pay rebate at hit time, so we set $w = 0 by default.
48              
49             Important(b) :Furthermore, for all contracts, we allow a different
50             payout currency (Quantos).
51              
52             Paying domestic currency (JPY if for USDJPY) = correlation coefficient is ZERO.
53             Paying foreign currency (USD if for USDJPY) = correlation coefficient is ONE.
54             Paying another currency = correlation is between negative ONE and positive ONE.
55              
56             See [3] for Quanto formulas and examples
57              
58             =head1 SUBROUTINES
59              
60             =head2 vanilla_call
61              
62             USAGE
63             my $price = vanilla_call($S, $K, $t, $r_q, $mu, $sigma)
64              
65             DESCRIPTION
66             Price of a Vanilla Call
67              
68             =cut
69              
70             sub vanilla_call {
71 2     2 1 827 my ( $S, $K, $t, $r_q, $mu, $sigma ) = @_;
72              
73 2         9 my $d1 =
74             ( log( $S / $K ) + ( $mu + $sigma * $sigma / 2.0 ) * $t ) /
75             ( $sigma * sqrt($t) );
76 2         8 my $d2 = $d1 - ( $sigma * sqrt($t) );
77              
78             return
79 2         17 exp( -$r_q * $t ) *
80             ( $S * exp( $mu * $t ) * pnorm($d1) - $K * pnorm($d2) );
81             }
82              
83             =head2 vanilla_put
84              
85             USAGE
86             my $price = vanilla_put($S, $K, $t, $r_q, $mu, sigma)
87              
88             DESCRIPTION
89             Price a standard Vanilla Put
90              
91             =cut
92              
93             sub vanilla_put {
94 2     2 1 799 my ( $S, $K, $t, $r_q, $mu, $sigma ) = @_;
95              
96 2         6 my $d1 =
97             ( log( $S / $K ) + ( $mu + $sigma * $sigma / 2.0 ) * $t ) /
98             ( $sigma * sqrt($t) );
99 2         3 my $d2 = $d1 - ( $sigma * sqrt($t) );
100              
101 2         15 return -1 *
102             exp( -$r_q * $t ) *
103             ( $S * exp( $mu * $t ) * pnorm( -$d1 ) - $K * pnorm( -$d2 ) );
104             }
105              
106             =head2 call
107              
108             USAGE
109             my $price = call($S, $K, $t, $r_q, $mu, $sigma)
110              
111             PARAMS
112             $S => stock price
113             $K => barrier
114             $t => time (1 = 1 year)
115             $r_q => payout currency interest rate (0.05 = 5%)
116             $mu => quanto drift adjustment (0.05 = 5%)
117             $sigma => volatility (0.3 = 30%)
118              
119             DESCRIPTION
120             Price a Call and remove the N(d2) part if the time is too small
121              
122             EXPLANATION
123             The definition of the contract is that if S > K, it gives
124             full payout (1). However the formula DC(T,K) = e^(-rT) N(d2) will not be
125             correct when T->0 and K=S. The value of DC(T,K) for this case will be 0.5.
126            
127             The formula is actually "correct" because when T->0 and S=K, the probability
128             should just be 0.5 that the contract wins, moving up or down is equally
129             likely in that very small amount of time left. Thus the only problem is
130             that the math cannot evaluate at T=0, where divide by 0 error occurs. Thus,
131             we need this check that throws away the N(d2) part (N(d2) will evaluate
132             "wrongly" to 0.5 if S=K).
133              
134             NOTE
135             Note that we have call = - dCall/dStrike
136             pair Foreign/Domestic
137              
138             see [3] for $r_q and $mu for quantos
139              
140             =cut
141              
142             sub call {
143 13     13 1 11607 my ( $S, $K, $t, $r_q, $mu, $sigma ) = @_;
144              
145 13 100       48 if ( $t < $SMALLTIME ) {
146 2 100       9 return ( $S > $K ) ? exp( -$r_q * $t ) : 0;
147             }
148              
149 11         59 return exp( -$r_q * $t ) * pnorm( d2( $S, $K, $t, $r_q, $mu, $sigma ) );
150             }
151              
152             =head2 put
153              
154             USAGE
155             my $price = put($S, $K, $t, $r_q, $mu, $sigma)
156              
157             PARAMS
158             $S => stock price
159             $K => barrier
160             $t => time (1 = 1 year)
161             $r_q => payout currency interest rate (0.05 = 5%)
162             $mu => quanto drift adjustment (0.05 = 5%)
163             $sigma => volatility (0.3 = 30%)
164              
165             DESCRIPTION
166             Price a standard Digital Put
167              
168             =cut
169              
170             sub put {
171 13     13 1 3299 my ( $S, $K, $t, $r_q, $mu, $sigma ) = @_;
172              
173 13 100       38 if ( $t < $SMALLTIME ) {
174 2 100       15 return ( $S < $K ) ? exp( -$r_q * $t ) : 0;
175             }
176              
177             return
178 11         29 exp( -$r_q * $t ) * pnorm( -1 * d2( $S, $K, $t, $r_q, $mu, $sigma ) );
179             }
180              
181             =head2 d2
182              
183             returns the DS term common to many BlackScholes formulae.
184              
185             =cut
186              
187             sub d2 {
188 22     22 1 31 my ( $S, $K, $t, $r_q, $mu, $sigma ) = @_;
189              
190 22         268 return ( log( $S / $K ) + ( $mu - $sigma * $sigma / 2.0 ) * $t ) /
191             ( $sigma * sqrt($t) );
192             }
193              
194             =head2 expirymiss
195              
196             USAGE
197             my $price = expirymiss($S, $U, $D, $t, $r_q, $mu, $sigma)
198              
199             PARAMS
200             $S => stock price
201             $t => time (1 = 1 year)
202             $U => barrier
203             $D => barrier
204             $r_q => payout currency interest rate (0.05 = 5%)
205             $mu => quanto drift adjustment (0.05 = 5%)
206             $sigma => volatility (0.3 = 30%)
207              
208             DESCRIPTION
209             Price an expiry miss contract (1 Call + 1 Put)
210              
211             [3] for $r_q and $mu for quantos
212              
213             =cut
214              
215             sub expirymiss {
216 6     6 1 706 my ( $S, $U, $D, $t, $r_q, $mu, $sigma ) = @_;
217              
218 6         10 my ($call_price) = call( $S, $U, $t, $r_q, $mu, $sigma );
219 6         14 my ($put_price) = put( $S, $D, $t, $r_q, $mu, $sigma );
220              
221 6         14 return $call_price + $put_price;
222             }
223              
224             =head2 expiryrange
225              
226             USAGE
227             my $price = expiryrange($S, $U, $D, $t, $r_q, $mu, $sigma)
228              
229             PARAMS
230             $S => stock price
231             $U => barrier
232             $D => barrier
233             $t => time (1 = 1 year)
234             $r_q => payout currency interest rate (0.05 = 5%)
235             $mu => quanto drift adjustment (0.05 = 5%)
236             $sigma => volatility (0.3 = 30%)
237              
238             DESCRIPTION
239             Price an Expiry Range contract as Foreign/Domestic.
240              
241             [3] for $r_q and $mu for quantos
242              
243             =cut
244              
245             sub expiryrange {
246 3     3 1 1024 my ( $S, $U, $D, $t, $r_q, $mu, $sigma ) = @_;
247              
248 3         25 return exp( -$r_q * $t ) - expirymiss( $S, $U, $D, $t, $r_q, $mu, $sigma );
249             }
250              
251             =head2 onetouch
252              
253             PARAMS
254             $S => stock price
255             $U => barrier
256             $t => time (1 = 1 year)
257             $r_q => payout currency interest rate (0.05 = 5%)
258             $mu => quanto drift adjustment (0.05 = 5%)
259             $sigma => volatility (0.3 = 30%)
260              
261             [3] for $r_q and $mu for quantos
262              
263             =cut
264              
265             sub onetouch {
266 36     36 1 1959 my ( $S, $U, $t, $r_q, $mu, $sigma, $w ) = @_;
267              
268             # w = 0, rebate paid at hit (good way to remember is that waiting
269             # time to get paid = 0)
270             # w = 1, rebate paid at end.
271              
272             # When the contract already reached it expiry and not yet reach it
273             # settlement time, it is consider an unexpired contract but will come to
274             # here with t=0 and it will caused the formula to die hence set it to the
275             # SMALLTIME which is 1 second
276 36         80 $t = max( $SMALLTIME, $t );
277              
278 36   100     134 $w ||= 0;
279              
280             # eta = -1, one touch up
281             # eta = 1, one touch down
282 36 100       59 my $eta = ( $S < $U ) ? -1 : 1;
283              
284 36         48 my $sqrt_t = sqrt($t);
285              
286 36         50 my $theta = ( ($mu) / $sigma ) + ( 0.5 * $sigma );
287 36         39 my $theta_ = ( ($mu) / $sigma ) - ( 0.5 * $sigma );
288              
289             # Floor v_ squared at zero in case negative interest rates push it negative.
290             # See: Barrier Options under Negative Rates in Black-Scholes (Le Floc’h and Pruell, 2014)
291 36         84 my $v_ = sqrt( max( 0, ( $theta_ * $theta_ ) + ( 2 * ( 1 - $w ) * $r_q ) ) ) ;
292              
293 36         148 my $e = ( log( $S / $U ) - ( $sigma * $v_ * $t ) ) / ( $sigma * $sqrt_t );
294 36         62 my $e_ = ( -log( $S / $U ) - ( $sigma * $v_ * $t ) ) / ( $sigma * $sqrt_t );
295              
296 36         308 my $price =
297             ( ( $U / $S )**( ( $theta_ + $v_ ) / $sigma ) ) * pnorm( -$eta * $e ) +
298             ( ( $U / $S )**( ( $theta_ - $v_ ) / $sigma ) ) * pnorm( $eta * $e_ );
299              
300 36         98 return exp( -$w * $r_q * $t ) * $price;
301             }
302              
303             =head2 notouch
304              
305             USAGE
306             my $price = notouch($S, $U, $t, $r_q, $mu, $sigma, $w)
307              
308             PARAMS
309             $S => stock price
310             $U => barrier
311             $t => time (1 = 1 year)
312             $r_q => payout currency interest rate (0.05 = 5%)
313             $mu => quanto drift adjustment (0.05 = 5%)
314             $sigma => volatility (0.3 = 30%)
315              
316             DESCRIPTION
317             Price a No touch contract.
318              
319             Payoff with domestic currency
320             Identity:
321             price of notouch = exp(- r t) - price of onetouch(rebate paid at end)
322              
323             [3] for $r_q and $mu for quantos
324              
325             =cut
326              
327             sub notouch {
328 7     7 1 1533 my ( $S, $U, $t, $r_q, $mu, $sigma ) = @_;
329              
330             # No touch contract always pay out at end
331 7         9 my $w = 1;
332              
333 7         27 return exp( -$r_q * $t ) - onetouch( $S, $U, $t, $r_q, $mu, $sigma, $w );
334             }
335              
336             # These variables require 'our' only because they need to be
337             # accessed by a test script.
338             our $MAX_ITERATIONS_UPORDOWN_PELSSER_1997 = 1000;
339             our $MIN_ITERATIONS_UPORDOWN_PELSSER_1997 = 16;
340              
341             #
342             # This variable requires 'our' only because it needs to be
343             # accessed via test script.
344             # Min accuracy. Accurate to 1 dollar for 100,000 notional
345             #
346             our $MIN_ACCURACY_UPORDOWN_PELSSER_1997 = 1.0 / 100000.0;
347             our $SMALL_VALUE_MU = 1e-10;
348              
349             # The smallest (in magnitude) floating-point number which,
350             # when added to the floating-point number 1.0, produces a
351             # floating-point result different from 1.0 is termed the
352             # machine accuracy, e.
353             #
354             # This value is very important for knowing stability to
355             # certain formulas used. e.g. Pelsser formula for UPORDOWN
356             # and RANGE contracts.
357             #
358             my $MACHINE_EPSILON = machine_epsilon();
359              
360             =head2 upordown
361              
362             USAGE
363             my $price = upordown(($S, $U, $D, $t, $r_q, $mu, $sigma, $w))
364              
365             PARAMS
366             $S stock price
367             $U barrier
368             $D barrier
369             $t time (1 = 1 year)
370             $r_q payout currency interest rate (0.05 = 5%)
371             $mu quanto drift adjustment (0.05 = 5%)
372             $sigma volatility (0.3 = 30%)
373              
374             see [3] for $r_q and $mu for quantos
375              
376             DESCRIPTION
377             Price an Up or Down contract
378              
379             =cut
380              
381             sub upordown {
382 15     15 1 5371 my ( $S, $U, $D, $t, $r_q, $mu, $sigma, $w ) = @_;
383              
384             # When the contract already reached it's expiry and not yet reach it
385             # settlement time, it is considered an unexpired contract but will come to
386             # here with t=0 and it will caused the formula to die hence set it to the
387             # SMALLTIME whiich is 1 second
388 15         36 $t = max( $t, $SMALLTIME );
389              
390             # $w = 0, paid at hit
391             # $w = 1, paid at end
392 15 100       40 if ( not defined $w ) { $w = 0; }
  8         15  
393              
394             # spot is outside [$D, $U] --> contract is expired with full payout,
395             # one barrier is already hit (can happen due to shift markup):
396 15 100 100     73 if ( $S >= $U or $S <= $D ) {
397 4 100       14 return $w ? exp( -$t * $r_q ) : 1;
398             }
399              
400             #
401             # SANITY CHECKS
402             #
403             # For extreme cases, the price will be wrong due the values in the
404             # infinite series getting too large or too small, which causes
405             # roundoff errors in the computer. Thus no matter how many iterations
406             # you make, the errors will never go away.
407             #
408             # For example try this:
409             #
410             # my ($S, $U, $D, $t, $r, $q, $vol, $w)
411             # = (100.00, 118.97, 99.00, 30/365, 0.1, 0.02, 0.01, 1);
412             # $up_price = Math::Business::BlackScholes::Binaries::ot_up_ko_down_pelsser_1997(
413             # $S,$U,$D,$t,$r,$q,$vol,$w);
414             # $down_price= Math::Business::BlackScholes::Binaries::ot_down_ko_up_pelsser_1997(
415             # $S,$U,$D,$t,$r,$q,$vol,$w);
416             #
417             # Thus we put a sanity checks here such that
418             #
419             # CONDITION 1: UPORDOWN[U,D] < ONETOUCH[U] + ONETOUCH[D]
420             # CONDITION 2: UPORDOWN[U,D] > ONETOUCH[U]
421             # CONDITION 3: UPORDOWN[U,D] > ONETOUCH[D]
422             # CONDITION 4: ONETOUCH[U] + ONETOUCH[D] >= $MIN_ACCURACY_UPORDOWN_PELSSER_1997
423             #
424 11         27 my $onetouch_up_prob = onetouch( $S, $U, $t, $r_q, $mu, $sigma, $w );
425 11         29 my $onetouch_down_prob = onetouch( $S, $D, $t, $r_q, $mu, $sigma, $w );
426              
427 11         14 my $upordown_prob;
428              
429 11 100 75     79 if ( $onetouch_up_prob + $onetouch_down_prob <
    100          
430             $MIN_ACCURACY_UPORDOWN_PELSSER_1997 )
431             {
432              
433             # CONDITION 4:
434             # The probability is too small for the Pelsser formula to be correct.
435             # Do this check first to avoid PELSSER stability condition to be
436             # triggered.
437             # Here we assume that the ONETOUCH formula is perfect and never give
438             # wrong values (e.g. negative).
439 1         5 return 0;
440             }
441             elsif ( $onetouch_up_prob xor $onetouch_down_prob ) {
442              
443             # One of our ONETOUCH probabilities is 0.
444             # That means our upordown prob is equivalent to the other one.
445             # Pelsser recompute will either be the same or wrong.
446             # Continuing to assume the ONETOUCH is perfect.
447 4         43 $upordown_prob = max( $onetouch_up_prob, $onetouch_down_prob );
448             }
449             else {
450              
451             # THIS IS THE ONLY PLACE IT SHOULD BE!
452 6         12 $upordown_prob =
453             ot_up_ko_down_pelsser_1997( $S, $U, $D, $t, $r_q, $mu, $sigma, $w ) +
454             ot_down_ko_up_pelsser_1997( $S, $U, $D, $t, $r_q, $mu, $sigma, $w );
455             }
456              
457             # CONDITION 4:
458             # Now check on the other end, when the contract is too close to payout.
459             # Not really needed to check for payout at hit, because RANGE is
460             # always at end, and thus the value (DISCOUNT - UPORDOWN) is not
461             # evaluated.
462 10 100       25 if ( $w == 1 ) {
463              
464             # Since the difference is already less than the min accuracy,
465             # the value [payout - upordown], which is the RANGE formula
466             # can become negative.
467 5 50       23 if (
468             abs( exp( -$r_q * $t ) - $upordown_prob ) <
469             $MIN_ACCURACY_UPORDOWN_PELSSER_1997 )
470             {
471 0         0 $upordown_prob = exp( -$r_q * $t );
472             }
473             }
474              
475             # CONDITION 1-3
476             # We use hardcoded small value of $SMALL_TOLERANCE, because if we were to increase
477             # the minimum accuracy, and this small value uses that min accuracy, it is
478             # very hard for the conditions to pass.
479 10         19 my $SMALL_TOLERANCE = 0.00001;
480 10 50 33     65 if (
      33        
481             not( $upordown_prob <
482             $onetouch_up_prob + $onetouch_down_prob + $SMALL_TOLERANCE )
483             or not( $upordown_prob + $SMALL_TOLERANCE > $onetouch_up_prob )
484             or not( $upordown_prob + $SMALL_TOLERANCE > $onetouch_down_prob )
485             )
486             {
487 0         0 die "UPORDOWN price sanity checks failed for S=$S, U=$U, "
488             . "D=$D, t=$t, r_q=$r_q, mu=$mu, sigma=$sigma, w=$w. "
489             . "UPORDOWN PROB=$upordown_prob , "
490             . "ONETOUCH_UP PROB=$onetouch_up_prob , "
491             . "ONETOUCH_DOWN PROB=$onetouch_down_prob";
492             }
493              
494 10         27 return $upordown_prob;
495             }
496              
497             =head2 common_function_pelsser_1997
498              
499             USAGE
500             my $c = common_function_pelsser_1997($S, $U, $D, $t, $r_q, $mu, $sigma, $w, $eta)
501              
502             DESCRIPTION
503             Return the common function from Pelsser's Paper (1997)
504              
505             =cut
506              
507             sub common_function_pelsser_1997 {
508              
509             # h: normalized high barrier, log(U/L)
510             # x: normalized spot, log(S/L)
511 14     14 1 2196 my ( $S, $U, $D, $t, $r_q, $mu, $sigma, $w, $eta ) = @_;
512              
513 14         16 my $pi = Math::Trig::pi;
514              
515 14         32 my $h = log( $U / $D );
516 14         13 my $x = log( $S / $D );
517              
518             # $eta = 1, onetouch up knockout down
519             # $eta = 0, onetouch down knockout up
520             # This variable used to check stability
521 14 100       29 if ( not defined $eta ) {
522 1         21 die "Wrong usage of this function for S=$S, U=$U, D=$D, "
523             . "t=$t, r_q=$r_q, mu=$mu, sigma=$sigma, w=$w, eta not defined.";
524             }
525 13 100       20 if ( $eta == 0 ) { $x = $h - $x; }
  7         9  
526              
527             # $w = 0, paid at hit
528             # $w = 1, paid at end
529              
530 13         13 my $mu_new = $mu - ( 0.5 * $sigma * $sigma );
531 13         32 my $mu_dash = sqrt(
532             ( $mu_new * $mu_new ) + ( 2 * $sigma * $sigma * $r_q * ( 1 - $w ) ) );
533              
534 13         14 my $series_part = 0;
535 13         7 my $hyp_part = 0;
536              
537             # These constants will determine whether or not this contract can be
538             # evaluated to a predefined accuracy. It is VERY IMPORTANT because
539             # if these conditions are not met, the prices can be complete nonsense!!
540 13         24 my $stability_constant =
541             get_stability_constant_pelsser_1997( $S, $U, $D, $t, $r_q, $mu, $sigma,
542             $w, $eta, 1 );
543              
544             # The number of iterations is important when recommending the
545             # range of the upper/lower barriers on our site. If we recommend
546             # a range that is too big and our iteration is too small, the
547             # price will be wrong! We must know the rate of convergence of
548             # the formula used.
549 13         21 my $iterations_required =
550             get_min_iterations_pelsser_1997( $S, $U, $D, $t, $r_q, $mu, $sigma, $w );
551              
552 13         26 for ( my $k = 1 ; $k < $iterations_required ; $k++ ) {
553 195         217 my $lambda_k_dash = (
554             0.5 * (
555             ( $mu_dash * $mu_dash ) / ( $sigma * $sigma ) +
556             ( $k * $k * $pi * $pi * $sigma * $sigma ) / ( $h * $h )
557             )
558             );
559              
560 195         224 my $phi =
561             ( $sigma * $sigma ) /
562             ( $h * $h ) *
563             exp( -$lambda_k_dash * $t ) *
564             $k /
565             $lambda_k_dash;
566              
567 195         209 $series_part += $phi * $pi * sin( $k * $pi * ( $h - $x ) / $h );
568              
569             #
570             # Note that greeks may also call this function, and their
571             # stability constant will differ. However, for simplicity
572             # we will not bother (else the code will get messy), and
573             # just use the price stability constant.
574             #
575 195 50 66     429 if ( $k == 1 and ( not( abs($phi) < $stability_constant ) ) ) {
576 0         0 die
577             "PELSSER VALUATION formula for S=$S, U=$U, D=$D, t=$t, r_q=$r_q, "
578             . "mu=$mu, vol=$sigma, w=$w, eta=$eta, cannot be evaluated because"
579             . "PELSSER VALUATION stability conditions ($phi less than "
580             . "$stability_constant) not met. This could be due to barriers "
581             . "too big, volatilities too low, interest/dividend rates too high, "
582             . "or machine accuracy too low. Machine accuracy is "
583             . $MACHINE_EPSILON . ".";
584             }
585             }
586              
587             #
588             # Some math basics: When A -> 0,
589             #
590             # sinh(A) -> 0.5 * [ (1 + A) - (1 - A) ] = 0.5 * 2A = A
591             # cosh(A) -> 0.5 * [ (1 + A) + (1 - A) ] = 0.5 * 2 = 1
592             #
593             # Thus for sinh(A)/sinh(B) when A & B -> 0, we have
594             #
595             # sinh(A) / sinh(B) -> A / B
596             #
597             # Since the check of the spot == lower/upper barrier has been done in the
598             # _upordown subroutine, we can assume that $x and $h will never be 0.
599             # So we only need to check that $mu_dash is too small. Also note that
600             # $mu_dash is always positive.
601             #
602             # For example, even at 0.0001 the error becomes small enough
603             #
604             # 0.0001 - Math::Trig::sinh(0.0001) = -1.66688941837835e-13
605             #
606             # Since h > x, we only check for (mu_dash * h) / (vol * vol)
607             #
608 13 100       30 if ( abs( $mu_dash * $h / ( $sigma * $sigma ) ) < $SMALL_VALUE_MU ) {
609 1         2 $hyp_part = $x / $h;
610             }
611             else {
612 12         34 $hyp_part =
613             Math::Trig::sinh( $mu_dash * $x / ( $sigma * $sigma ) ) /
614             Math::Trig::sinh( $mu_dash * $h / ( $sigma * $sigma ) );
615             }
616              
617 13         167 return ( $hyp_part - $series_part ) * exp( -$r_q * $t * $w );
618             }
619              
620             =head2 get_stability_constant_pelsser_1997
621              
622             USAGE
623             my $constant = get_stability_constant_pelsser_1997($S, $U, $D, $t, $r_q, $mu, $sigma, $w, $eta, $p)
624              
625             DESCRIPTION
626             Get the stability constant (Pelsser 1997)
627              
628             =cut
629              
630             sub get_stability_constant_pelsser_1997 {
631 15     15 1 3317 my ( $S, $U, $D, $t, $r_q, $mu, $sigma, $w, $eta, $p ) = @_;
632              
633             # $eta = 1, onetouch up knockout down
634             # $eta = 0, onetouch down knockout up
635              
636 15 100       26 if ( not defined $eta ) {
637 1         29 die "Wrong usage of this function for S=$S, U=$U, D=$D, t=$t, "
638             . "r_q=$r_q, mu=$mu, sigma=$sigma, w=$w, Eta not defined.";
639             }
640              
641             # p is the power of pi
642             # p=1 for price/theta/vega/vanna/volga
643             # p=2 for delta
644             # p=3 for gamma
645 14 50 66     38 if ( $p != 1 and $p != 2 and $p != 3 ) {
      66        
646 1         21 die "Wrong usage of this function for S=$S, U=$U, D=$D, t=$t, "
647             . "r_q=$r_q, mu=$mu, sigma=$sigma, w=$w, Power of PI must "
648             . "be 1, 2 or 3. Given $p.";
649             }
650              
651 13         15 my $h = log( $U / $D );
652 13         13 my $x = log( $S / $D );
653 13         12 my $mu_new = $mu - ( 0.5 * $sigma * $sigma );
654 13         33 my $mu_dash = sqrt(
655             ( $mu_new * $mu_new ) + ( 2 * $sigma * $sigma * $r_q * ( 1 - $w ) ) );
656              
657 13         28 my $numerator = $MIN_ACCURACY_UPORDOWN_PELSSER_1997 *
658             exp( 1.0 - $mu_new * ( ( $eta * $h ) - $x ) / ( $sigma * $sigma ) );
659 13         39 my $denominator =
660             ( exp(1) * ( Math::Trig::pi + $p ) ) +
661             ( max( $mu_new * ( ( $eta * $h ) - $x ), 0.0 ) *
662             Math::Trig::pi /
663             ( $sigma**2 ) );
664 13         22 $denominator *= ( Math::Trig::pi**( $p - 1 ) ) * $MACHINE_EPSILON;
665              
666 13         15 my $stability_condition = $numerator / $denominator;
667              
668 13         15 return $stability_condition;
669             }
670              
671             =head2 ot_up_ko_down_pelsser_1997
672              
673             USAGE
674             my $price = ot_up_ko_down_pelsser_1997($S, $U, $D, $t, $r_q, $mu, $sigma, $w)
675              
676             DESCRIPTION
677             This is V_{RAHU} in paper [5], or ONETOUCH-UP-KNOCKOUT-DOWN,
678             a contract that wins if it touches upper barrier, but expires
679             worthless if it touches the lower barrier first.
680              
681             =cut
682              
683             sub ot_up_ko_down_pelsser_1997 {
684 6     6 1 10 my ( $S, $U, $D, $t, $r_q, $mu, $sigma, $w ) = @_;
685              
686 6         8 my $mu_new = $mu - ( 0.5 * $sigma * $sigma );
687 6         5 my $h = log( $U / $D );
688 6         9 my $x = log( $S / $D );
689              
690             return
691 6         13 exp( $mu_new * ( $h - $x ) / ( $sigma * $sigma ) ) *
692             common_function_pelsser_1997( $S, $U, $D, $t, $r_q, $mu, $sigma, $w, 1 );
693             }
694              
695             =head2 ot_down_ko_up_pelsser_1997
696              
697             USAGE
698             my $price = ot_down_ko_up_pelsser_1997($S, $U, $D, $t, $r_q, $mu, $sigma, $w)
699              
700             DESCRIPTION
701             This is V_{RAHL} in paper [5], or ONETOUCH-DOWN-KNOCKOUT-UP,
702             a contract that wins if it touches lower barrier, but expires
703             worthless if it touches the upper barrier first.
704              
705             =cut
706              
707             sub ot_down_ko_up_pelsser_1997 {
708 6     6 1 9 my ( $S, $U, $D, $t, $r_q, $mu, $sigma, $w ) = @_;
709              
710 6         13 my $mu_new = $mu - ( 0.5 * $sigma * $sigma );
711 6         9 my $h = log( $U / $D );
712 6         5 my $x = log( $S / $D );
713              
714             return
715 6         12 exp( -$mu_new * $x / ( $sigma * $sigma ) ) *
716             common_function_pelsser_1997( $S, $U, $D, $t, $r_q, $mu, $sigma, $w, 0 );
717             }
718              
719             =head2 get_min_iterations_pelsser_1997
720              
721             USAGE
722             my $min = get_min_iterations_pelsser_1997($S, $U, $D, $t, $r_q, $mu, $sigma, $w, $accuracy)
723              
724             DESCRIPTION
725             An estimate of the number of iterations required to achieve a certain
726             level of accuracy in the price.
727              
728             =cut
729              
730             sub get_min_iterations_pelsser_1997 {
731 16     16 1 2776 my ( $S, $U, $D, $t, $r_q, $mu, $sigma, $w, $accuracy ) = @_;
732              
733 16 100       30 if ( not defined $accuracy ) {
734 14         13 $accuracy = $MIN_ACCURACY_UPORDOWN_PELSSER_1997;
735             }
736              
737 16 100       41 if ( $accuracy > $MIN_ACCURACY_UPORDOWN_PELSSER_1997 ) {
    100          
738 1         1 $accuracy = $MIN_ACCURACY_UPORDOWN_PELSSER_1997;
739             }
740             elsif ( $accuracy <= 0 ) {
741 1         2 $accuracy = $MIN_ACCURACY_UPORDOWN_PELSSER_1997;
742             }
743              
744 16         31 my $h = log( $U / $D );
745 16         13 my $x = log( $S / $D );
746              
747 16         26 my $it_up =
748             _get_min_iterations_ot_up_ko_down_pelsser_1997( $S, $U, $D, $t, $r_q, $mu,
749             $sigma, $w, $accuracy );
750 16         24 my $it_down =
751             _get_min_iterations_ot_down_ko_up_pelsser_1997( $S, $U, $D, $t, $r_q, $mu,
752             $sigma, $w, $accuracy );
753              
754 16         27 my $min = max( $it_up, $it_down );
755              
756 16         30 return $min;
757             }
758              
759             =head2 _get_min_iterations_ot_up_ko_down_pelsser_1997
760              
761             USAGE
762             my $k_min = _get_min_iterations_ot_up_ko_down_pelsser_1997($S, $U, $D, $t, $r_q, $mu, $sigma, $w, $accuracy)
763              
764             DESCRIPTION
765             An estimate of the number of iterations required to achieve a certain
766             level of accuracy in the price for ONETOUCH-UP-KNOCKOUT-DOWN.
767              
768             =cut
769              
770             sub _get_min_iterations_ot_up_ko_down_pelsser_1997 {
771 35     35   917 my ( $S, $U, $D, $t, $r_q, $mu, $sigma, $w, $accuracy ) = @_;
772              
773 35 100       56 if (!defined $accuracy) {
774 1         7 die "accuracy required";
775             }
776              
777 34         25 my $pi = Math::Trig::pi;
778              
779 34         36 my $h = log( $U / $D );
780 34         22 my $x = log( $S / $D );
781 34         33 my $mu_new = $mu - ( 0.5 * $sigma * $sigma );
782 34         52 my $mu_dash = sqrt(
783             ( $mu_new * $mu_new ) + ( 2 * $sigma * $sigma * $r_q * ( 1 - $w ) ) );
784              
785 34         35 my $A = ( $mu_dash * $mu_dash ) / ( 2 * $sigma * $sigma );
786 34         34 my $B = ( $pi * $pi * $sigma * $sigma ) / ( 2 * $h * $h );
787              
788 34         24 my $delta_dash = $accuracy;
789 34         69 my $delta =
790             $delta_dash *
791             exp( -$mu_new * ( $h - $x ) / ( $sigma * $sigma ) ) *
792             ( ( $h * $h ) / ( $pi * $sigma * $sigma ) );
793              
794             # This can happen when stability condition fails
795 34 100       54 if ( $delta * $B <= 0 ) {
796 1         27 die "(_get_min_iterations_ot_up_ko_down_pelsser_1997) Cannot "
797             . "evaluate minimum iterations because too many iterations "
798             . "required!! delta=$delta, B=$B for input parameters S=$S, "
799             . "U=$U, D=$D, t=$t, r_q=$r_q, mu=$mu, sigma=$sigma, w=$w, "
800             . "accuracy=$accuracy";
801 0         0 return $MAX_ITERATIONS_UPORDOWN_PELSSER_1997;
802             }
803              
804             # Check that condition is satisfied
805 33         69 my $condition = max( exp( -$A * $t ) / ( $B * $delta ), 1 );
806              
807 33         36 my $k_min = log($condition) / ( $B * $t );
808 33         23 $k_min = sqrt($k_min);
809              
810 33 100       48 if ( $k_min < $MIN_ITERATIONS_UPORDOWN_PELSSER_1997 ) {
    50          
811              
812 32         52 return $MIN_ITERATIONS_UPORDOWN_PELSSER_1997;
813             }
814             elsif ( $k_min > $MAX_ITERATIONS_UPORDOWN_PELSSER_1997 ) {
815              
816 1         3 return $MAX_ITERATIONS_UPORDOWN_PELSSER_1997;
817             }
818              
819 0         0 return int($k_min);
820             }
821              
822             =head2 _get_min_iterations_ot_down_ko_up_pelsser_1997
823              
824             USAGE
825              
826             DESCRIPTION
827             An estimate of the number of iterations required to achieve a certain
828             level of accuracy in the price for ONETOUCH-UP-KNOCKOUT-UP.
829              
830             =cut
831              
832             sub _get_min_iterations_ot_down_ko_up_pelsser_1997 {
833 16     16   21 my ( $S, $U, $D, $t, $r_q, $mu, $sigma, $w, $accuracy ) = @_;
834              
835 16         17 my $h = log( $U / $D );
836 16         17 my $x = log( $S / $D );
837 16         16 my $mu_new = $mu - ( 0.5 * $sigma * $sigma );
838              
839 16         22 $accuracy = $accuracy * exp( $mu_new * $h / ( $sigma * $sigma ) );
840              
841 16         17 return _get_min_iterations_ot_up_ko_down_pelsser_1997( $S, $U, $D, $t, $r_q,
842             $mu, $sigma, $w, $accuracy );
843             }
844              
845             =head2 range
846              
847             USAGE
848             my $price = range($S, $U, $D, $t, $r_q, $mu, $sigma, $w)
849              
850             PARAMS
851             $S stock price
852             $t time (1 = 1 year)
853             $U barrier
854             $D barrier
855             $r_q payout currency interest rate (0.05 = 5%)
856             $mu quanto drift adjustment (0.05 = 5%)
857             $sigma volatility (0.3 = 30%)
858              
859             see [3] for $r_q and $mu for quantos
860              
861             DESCRIPTION
862             Price a range contract.
863              
864             =cut
865              
866             sub range {
867              
868             # payout time $w is only a dummy. range contracts always payout at end.
869 7     7 1 2479 my ( $S, $U, $D, $t, $r_q, $mu, $sigma, $w ) = @_;
870              
871             # range always pay out at end
872 7         10 $w = 1;
873              
874             return
875 7         32 exp( -$r_q * $t ) - upordown( $S, $U, $D, $t, $r_q, $mu, $sigma, $w );
876             }
877              
878             =head1 DEPENDENCIES
879              
880             * Math::CDF
881             * Machine::Epsilon
882              
883             =head1 SOURCE CODE
884              
885             https://github.com/binary-com/perl-math-business-blackscholes-binaries
886              
887             =head1 REFERENCES
888              
889             [1] P.G Zhang [1997], "Exotic Options", World Scientific
890             Another good refernce is Mark rubinstein, Eric Reiner [1991], "Binary Options", RISK 4, pp 75-83
891              
892             [2] Anlong Li [1999], "The pricing of double barrier options and their variations".
893             Advances in Futures and Options, 10, 1999. (paper).
894              
895             [3] Uwe Wystup. FX Options and Strutured Products. Wiley Finance, England, 2006. pp 93-96 (Quantos)
896              
897             [4] Antoon Pelsser, "Pricing Double Barrier Options: An Analytical Approach", Jan 15 1997.
898             http://repub.eur.nl/pub/7807/1997-0152.pdf
899              
900             =head1 AUTHOR
901              
902             binary.com, C<< >>
903              
904             =head1 BUGS
905              
906             Please report any bugs or feature requests to
907             C, or through the web
908             interface at
909             L.
910             I will be notified, and then you'll automatically be notified of progress on
911             your bug as I make changes.
912              
913             =head1 SUPPORT
914              
915             You can find documentation for this module with the perldoc command.
916              
917             perldoc Math::Business::BlackScholes::Binaries
918              
919              
920             You can also look for information at:
921              
922             =over 4
923              
924             =item * RT: CPAN's request tracker (report bugs here)
925              
926             L
927              
928             =item * AnnoCPAN: Annotated CPAN documentation
929              
930             L
931              
932             =item * CPAN Ratings
933              
934             L
935              
936             =item * Search CPAN
937              
938             L
939              
940             =back
941              
942              
943             =head1 LICENSE AND COPYRIGHT
944              
945             Copyright 2014 binary.com.
946              
947             This program is free software; you can redistribute it and/or modify it
948             under the terms of the the Artistic License (2.0). You may obtain a
949             copy of the full license at:
950              
951             L
952              
953             Any use, modification, and distribution of the Standard or Modified
954             Versions is governed by this Artistic License. By using, modifying or
955             distributing the Package, you accept this license. Do not use, modify,
956             or distribute the Package, if you do not accept this license.
957              
958             If your Modified Version has been derived from a Modified Version made
959             by someone other than you, you are nevertheless required to ensure that
960             your Modified Version complies with the requirements of this license.
961              
962             This license does not grant you the right to use any trademark, service
963             mark, tradename, or logo of the Copyright Holder.
964              
965             This license includes the non-exclusive, worldwide, free-of-charge
966             patent license to make, have made, use, offer to sell, sell, import and
967             otherwise transfer the Package with respect to any patent claims
968             licensable by the Copyright Holder that are necessarily infringed by the
969             Package. If you institute patent litigation (including a cross-claim or
970             counterclaim) against any party alleging that the Package constitutes
971             direct or contributory patent infringement, then this Artistic License
972             to you shall terminate on the date that such litigation is filed.
973              
974             Disclaimer of Warranty: THE PACKAGE IS PROVIDED BY THE COPYRIGHT HOLDER
975             AND CONTRIBUTORS "AS IS' AND WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES.
976             THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
977             PURPOSE, OR NON-INFRINGEMENT ARE DISCLAIMED TO THE EXTENT PERMITTED BY
978             YOUR LOCAL LAW. UNLESS REQUIRED BY LAW, NO COPYRIGHT HOLDER OR
979             CONTRIBUTOR WILL BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, OR
980             CONSEQUENTIAL DAMAGES ARISING IN ANY WAY OUT OF THE USE OF THE PACKAGE,
981             EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
982              
983              
984             =cut
985              
986             1;
987