# Cumulative Normal Distributions

# Black-Scholes and Error Functions

# Ralph Selfridge

selfrid@cns.ufl.edu

Recent comments in the J Forum about normal probabilities and distributions have led to a revival of interest in Black-Scholes algorithms. Black-Scholes is a set of equations that provide premiums for ‘call’ and ‘put’ options, and there has been an interest in providing appropriate algorithms in J. That, in turn, opens the question of an algorithm that will run in APL. This study was an attempt to create such an algorithm, which may already exist in some APL library but was not available to this author. We use APL2 as the APL medium. A serendipitous result is also a speed-up in computing cumulative normal distributions and an error function for APL.

We start here with Eugene McDonnell’s article
“At Play with J: Beware Scholes” [1],
which provides Black-Scholes equations in reasonable mathematical form
and then considers some algorithms in J. This article
has several algorithms created, among others, by Hu Zhe and
Oleg Kobchenko. McDonnell uses the `cnd`

of Zhe and his own
verb `bs`

for the algorithm of page 141, which we label A1,
followed on page 142 by an algorithm, which we label A2, using a
formula from Ewart Shaw. All algorithms, using essentially the same
verb `bs`

, are below in one of two appendices.

We know that A2 is much faster in execution than A1, but it has the
problem that it uses `H.`

(an included verb in J) not
available in APL2. Thus we must create a replacement for
`H.`

, or go back to A1 and use the `cnd`

verb
from Zhe. This relies only on a verb `normalprob`

from
`statdist`

which can be converted so that we have
Black-Scholes in APL2.

However if we examine `cnd`

in A1 there is the noun
`__`

(negative infinity), which results in a domain error
if we replace `__`

with a large negative number for use in
APL2. That, in turn, forces examination of `normalprob`

to
find, and work around, `__`

. With that covered we have a
working Black-Scholes in APL2. Since fairly simple APL is used these
algorithms should all work in any APL.

The conversion from the verb `cnd`

(J) to the function
`CND`

(APL) allows for considerable simplification, and
even further if `CND`

is only for use in Black-Scholes. We
make these changes, and any other improvements we can think of, and
obtain `CND`

and `BS`

in Appendix 2. Now if we
consider A2 we have a formula for the error function but that in turn
allows us to come back from `CND`

to get an easy APL
function for the error function of a variable (it is entirely possible
that such a function exists in an APL library elsewhere). We add this
function to Appendix 2.

Now if we consider the changes for APL2, which shorten
`CND`

considerably, we can ask if similar changes can be
made in `CND`

of A1. These changes produce A3 in Appendix
1.

S = Current stock price

X = Option strike price

R = Risk-free interest rate

T = time in years until strike date

V = volatility, or standard deviation of asset price

All four algorithms are entered with `BS S,X,T,R,V`

and
return two numbers: the first is the ‘call’ premium, the second the
‘put’ premium, with an example

BS 60 65 0.25 0.08 0.3 2.13337 5.84628

We need to compare execution times for the differing algorithms, both
for BS and `CND/ERF`

. In the case of J we use
`(10000)6!:2'BS 60 65 0.25 0.08 0.3'`

;
for APL we must create a timing function that
runs 10000 executions of the desired function. In the following table
we set w=1.2 and ww a vector of length 500 of numbers in the range 0
to 5 `((i.500)%100)`

. All results are normalized so that A2
values are 100.

BS yc | CND w | CND ww | ERF w | ERF ww | |
---|---|---|---|---|---|

A1 | 288 | 1490 | 1514 | ||

A2 | 100 | 100 | 100 | 100 | 100 |

A3 | 87 | 154 | 4.8 | 289 | 13.9 |

APL | 75 | 130 | 27 | 230 | 69 |

When comparing results numeric output across the algorithms agree to several significant places. Since we must use system clocks, without knowing precisely how they run (nor do we know the resolution of the clocks), we use a repetition of 10000. Even so, timing comparisons should not be taken too seriously.

We find that A2 and A3 are substantially the same for Black-Scholes,
and the benefit of avoiding `H.`

is debatable. It is also interesting
that APL appears faster for Black-Scholes. Both `CND`

and `ERF`

are faster
in A2 for a scalar argument but considerably slower as the argument
length grows; the fact that `BS`

has an internal 2 by 2 is probably the
explanation for speed-ups in A3 and APL. In the case of APL, avoiding
`H.`

becomes crucial, unless an `H`

function is written either in APL
(very slow) or some better choice, say Fortran and a DLL. There is
little choice in APL, but it would appear A3 is better in J, unless
the arguments for `cnd`

and `erf`

are scalars.

If we create a Fortran subroutine to replace `H.`

specifically for
Black-Scholes, and hook across to APL2 by way of a DLL the
resulting timing for `BS`

is slightly slower than the given APL2 version
(though much faster than using a Fortran version of `H.`

with the
attendant use of `H"0`

).

### References

- Eugene McDonnell, “At Play with J: Beware Scholes” Vector 19.3, pp137-142

## Appendix 1

A1 load 'statdist' cnd=: 3 : 'normalprob 0 1 __, y' "0 NB. Hu Zhe BS=: 3 : 0 NB. Eugene McDonnell 'S X T R V'=.y d=. ((^.S%X) + T*R(+,−)−:*:V)%V*%:T │(S,X*^−R*T)−/ .*cnd d*/1 _1 ) A2 erf=:1 H. 1.5@*: * 2p_0.5 &* %^@:*: NB. Ewart Shaw cnd=: −:@ >:@erf %&(%:2) BS=: 3 : 0 'S X T R V'=.y d=. ((^.S%X) + T*R(+,−)−:*:V)%V*%:T │(S,X*^−R*T)−/ .*cnd d*/1 _1 ) A3 (modified A1) cnd=: 3 : 0 b=. 0 0.31938153 _0.356563782 1.781477937 _1.821255978 1.330274429 1−│(0.398942* ^_0.5**:y)* b p. %>:y*0.2316419 ) BS=: 3 : 0 'S X T R V'=.y d=. ((^.S%X) + T*R(+,−)−:*:V)%V*%:T │(S,X*^−R*T)−/ .*cnd d*/1 _1 )

## Appendix 2

`CND`

provides for cumulative normal distributions, but
only for use in Black-Scholes. It can be modified for any shape
argument by simple changes in line 2 (made, of necessity, for timing
comparisons). `ERF`

returns the error function value for
its argument, and must also be modified for timing comparisons (again
in line 2) or multiple arguments.

∇ Z←CND Y [1] Z←6 1⍴1.330274 ¯1.821256 1.78148 ¯0.356565 0.31938153 0 [2] Z←1-|(0.398942×*¯0.5×Y*2)×2 2⍴(4 1⍴÷1+Y×0.2316419)⊥Z ∇ ∇ Z←BS Y;S;X;T;R;V [1] (S X T R V)←Y [2] Z←|(S,X×*-R×T)-.×CND(((⍟S÷X)+T×R+Z,-Z←0.5×V*2)÷V×T*0.5)∘.×1 ¯1 ∇ ∇ Z←ERF Y [1] Z←6 1⍴1.330274 ¯1.821256 1.78148 ¯0.356565 0.31938153 0 [2] Z←1-2×|(0.398942×*-Y*2)×(÷1+Y×0.3275911166)⊥Z ∇

Readers in Scotland, noted for their economy with code, might prefer
the following slightly shorter version in Dyalog, using the
direct-definition form implemented by John Scholes. We understand John
is widely referred to north of the Border as “The Black Scholes”.
*Ed.*

CONSTANTS←6 1⍴1.330274 ¯1.821256 1.78148 ¯0.356565 0.31938153 0 CND←CONSTANTS∘{1-|(0.398942×*¯0.5×⍵*2)×2 2⍴(4 1⍴÷1+⍵×0.2316419)⊥⍺} ERF←CONSTANTS∘{1-2×|(0.398942×*-Y*2)×(÷1+⍵×0.3275911166)⊥⍺} ∇ Z←BS (S X T R V) [1] Z←|(S,X×*-R×T)-.×CND(((⍟S÷X)+T×R+1 ¯1×0.5×V*2)÷V×T*0.5)∘.×1 ¯1 ∇