Being normal

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 x

Lastly, 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.