### Being normal

Given the sophistication of today's financial markets, it's hard to imagine that until 1973, despite active trading since at least the 17th century, there was no known way of estimating the price of European-style options!

Black, Scholes and Merton changed all that with the development of a mathematical framework of a financial market in their seminal paper "The Pricing of Options and Corporate Liabilities". From that model, one can deduce the celebrated **Black-Scholes formula** : $p = \theta \cdot \left[ F \cdot \Phi\left( \theta \cdot
\left[ \frac{\ln \left( F/K \right)}{\sigma} + \frac{\sigma}{2}\right] \right) - K \cdot \Phi\left( \theta \cdot
\left[ \frac{\ln \left( F/K \right)}{\sigma} - \frac{\sigma}{2} \right] \right) \right]$ where $\theta = 1$ for call options and $\theta = -1$ for put options,
$F := Se^{\left( r-d \right)T}$, $S$ is the current spot, $K$ is the
strike, $r$ is a flat continuously compounded interest rate to maturity, $d$ is a flat continuously compounded dividend rate, $\sigma = \hat{\sigma} \cdot \sqrt{T}$, $\hat{\sigma}$ is the root-mean-square lognormal volatility, $T$ is the time to expiry and $\Phi(\cdot)$ is the standard cumulative normal distribution function.

Now, the only real problem with evaluating the above equation is in estimating the the values of the **cumulative distribution function of the normal distribution** ($\Phi$). This is a function that is defined by an integral that cannot be expressed in terms of elementary functions and is thus said to be one of a family of *special* functions.

One approach to obtaining a high quality, cross platform implementation of $\Phi$ in OCaml is to leverage the C++ Boost Math Toolkit library. First the C wrapper (let's write this in a file 'special_functions_c.cpp' (say)).

#include <boost/math/distributions/normal.hpp> #include <caml/mlvalues.h> #include <caml/alloc.h> extern "C" value caml_norm_cdf (value v) { using boost::math::cdf; using boost::math::normal; return caml_copy_double (cdf (normal (0.0, 1.0), Double_val (v))); }Next, the forwarding OCaml function ('special_functions.ml').

external caml_norm_cdf : float -> float = "caml_norm_cdf" let norm_cdf x = caml_norm_cdf xLastly, the module interface ('special_functions_sig.mli').

module type S = sig val norm_cdf : float -> float end

Et voilĂ . With this in hand, the Black-Scholes equation follows easily!

let black_scholes ~spot ~strike ~r ~vol ~t ~cp = let open Special_functions in let d1 = (log (spot/.strike) +. (r +. 0.5*.(vol*.vol))*.t)/.(vol*.(sqrt t)) in let d2 = d1 -. vol*.(sqrt t) in cp*.spot*.(norm_cdf (cp*.d1)) -. cp*.strike*.(exp ((~-.r)*.t))*.(norm_cdf (cp*.d2))Building the OCaml special functions library above follows the same procedure as covered in this blog entry. To wrap things up (and marvel at our handiwork,) here's a little test.

let _ = Printf.printf "%.8f\n" (black_scholes ~spot:42. ~strike:40. ~r:0.1 ~vol:0.2 ~t:0.5 ~cp:1.0) ; Printf.printf "%.8f\n" (black_scholes ~spot:42. ~strike:40. ~r:0.1 ~vol:0.2 ~t:0.5 ~cp:(~-.1.0))I observe the values

*4.75942239*for the call and,

*0.80859937*for the put, respectively.