The Uncarved Block

next>> [Last update: 2007-07-09 20:49:53 GMT] [37 articles in total] [show 5 10 per page]

Towards a more idiomatic FP pricer

Taking advice from a real expert, the pricer becomes more idiomatic

In my quest for a good ocaml book, I found and bought this one: "OCaml for Scientists" by Jon Harrop. I really can't recommend it too highly- it's an excellent book. Furthermore, Jon Harrop looked at this website and gave me some advice about making the code more idiomatic by using record types. He even wrote the first one for me to show how!

The basic idea is that you don't need classes to encapsulate data (as I have been using them), you can just use bindings which become enclosed within functions in a record type. This is very nice indeed.

(** parameter for models *)
type param = {
    integral :       float -> float -> float;
    integral_sq :    float -> float -> float;
    mean :           float -> float -> float;
    rms :            float -> float -> float;
}
(** Constructor function that makes a parameter record given an
 integral and an integral_sq func *)
let make_param integral integral_sq =
    {
        integral=integral;
        integral_sq=integral_sq;
        mean=(fun t1 t2->(integral t1 t2) /. (t2 -. t1));
        rms=(fun t1 t2->(integral_sq t1 t2) /. (t2 -. t1));
    };;

(** build a constant parameter *)
let param_const x =
    let integral t1 t2 = (t2-.t1) *. x in
    let x_sq = x *. x in
    let integral_sq t1 t2 = (t2-.t1) *. x_sq in
make_param integral integral_sq;;

(** Parameter that models a linear function f(x) -> ax + b *)
let param_linear a b =
    let integral t1 t2 =
      let integrate_to x = 0.5 *. x *. a *. x +. x *. b in
        integrate_to t2 -. integrate_to t1
    in
    let integral_sq t1 t2 =
        a ** 2. *. (t2 ** 3. -. t1 ** 3.) +.
        a *. b *. (t2 ** 2. -. t1 ** 1.) +.
        b ** 2. *. (t2 -. t1)
    in
    make_param integral integral_sq;;

So our classes have become a simple record of four functions and some functions which return this type. Notice how make_param takes two functions and makes two more from them. The code is otherwise similar to the class-based examples, but is a little more straightforward. For piecewise parameters, we do the same sort of thing:

(** For piecewise constants and piecewise linear functions a "piece" is a
  parameter that applies between low and high limits. *)
type piece = { low: float; high: float; p: param; };;

(** Make piecewise parameters.  Takes a list of pieces which each
  specify the limits and the linear or const parameter in each region.

  At present, does no checking that the pieces are continuous and not
  overlapping *)
let param_piecewise pieces =
  (* there's probably a better way of doing this, but anyway...
     Apply a function to each piece, with the parameters being the part of
     the interval from t1 to t2 that is inside the interval of the piece *)
  let visit_pieces pieces fn t1 t2 =
      let visit_piece piece low high =
        if (low > t2 || high < t1) then
          0.
        else
            let range_start = max t1 low in
            let range_end = min t2 high in
            (fn piece) range_start range_end in
      let rec visit_list so_far lis =
        match lis with
            [] -> so_far (** we're done *)
          | {low=low; high=high; p=p} :: rest -> visit_list (so_far +. (visit_piece p low high)) rest
      in
        visit_list 0. pieces
  in
  let integral t1 t2 = visit_pieces pieces (fun x->x.integral) t1 t2 in
  let integral_sq t1 t2 = visit_pieces pieces (fun x->x.integral_sq) t1 t2
  in
  make_param integral integral_sq;;

(** helper func to make parts of a piecewise const function *)
let make_const_piece low high x =
  assert (low<high);
  {low=low; high=high; p=(param_const x)};;

(** helper func to make parts of a piecewise linear function *)
let make_linear_piece low high a b =
  assert (low<high);
  {low=low; high=high; p=(param_linear a b)};;

This really is a great deal nicer than what we had before. Flushed by our success, we make a simple function to turn a list of x,y tuples into a piecewise constant function:

let make_const_param_curve lis =
  let rec build_curve (so_far : piece list) x_old lis =
    match lis with
      [] -> List.rev so_far
    | (x, y)::res when x>x_old ->
        let this_piece = make_const_piece x_old x y in
        build_curve (this_piece::so_far) x res
    | (x, y)::res -> invalid_arg "make_const_param_curve: expects sorted x in list"
  in
    build_curve [] 0. lis;;

Now we can build a piecewise const function by doing:

make_const_param_curve [ 0.1, 0.5; 0.5, 0.25; 1., 0.2 ];;

It would be a simple matter to make a function that built a piecewise linear function up in the same manner.

A Simple Ocaml test harness

It's nice to be able to write tests for code as you go along, so you need a test harness...

I usually like to be able to write tests as I write code, and to have them run every time the code builds to make sure I haven't broken anything. To do this, you need a test harness so that adding tests is as painless as possible. So I wrote this one:

(** test.ml - a simple test harness in ocaml

  This is demonstration code only.  You are free to use it under the
  terms of the Creative Commons Attribution 2.5 license.

  @author Sean Hunter <sean\@uncarved.com>

 *)

(** Test harness class.  Tracks numbers of tests run and how many 
 succeed. *)
class harness =
object(self)
    (** The number of tests so far *)
    val mutable n = 0

    (** The number of tests so far which have succeeded *)
    val mutable succeeded = 0

    (** A name for the group of tests currenly running *)
    val mutable group = ""

    (** Runs a predicate function and fails if it throws or
     returns false.  Otherwise it succeeds *)
    method pass_if desc pred =
      n <- n + 1;
      let dots = String.make (50-(min 50 (String.length desc))) '.' in
      Printf.printf "%5.5d: %s ....%s" n desc dots;

      try
          if pred () then
            begin
                Printf.printf "ok\n" ;
                succeeded <- succeeded + 1;
                true
            end
          else
            begin
              Printf.printf "not ok\n";
              false
            end
      with
      _ ->Printf.printf "not ok (threw exception)\n"; false

    (** Runs a predicate function and fails if it throws or
     returns true.  Otherwise it succeeds *)
    method fail_if desc pred = self#pass_if desc (fun () -> not pred)

    (** Takes a bool and marks the test as succeeded if it is true *)
    method ok desc x = self#pass_if desc (fun () -> x)

    (** Takes a bool and marks the test as failed if it is true *)
    method not_ok desc x = self#ok desc (not x)

    (** Evaluate a predicate, ignoring its result. Succeed if it throws
     an exception, fail if not *)
    method pass_if_throws desc (pred : unit -> bool) =
        try
            (ignore (pred ()));
            self#not_ok desc true
        with
            exn -> self#ok desc true

    (** Evaluate a predicate, ignoring its result. Fail if it throws
     an exception, succeed if not *)
    method fail_if_throws desc (pred : unit -> bool) =
        try
            (ignore (pred ()));
            self#ok desc true
        with
            exn -> self#not_ok desc true

    (** Check that two floats are within a certain tolerance *)
    method pass_if_close ?(eps=1e-9) desc x y = 
        let diff = abs_float (x-.y) in
        self#ok desc (diff <= eps)

    (** mark the beginning of a group of tests *)
    method start_group name =
        n <- 0 ;
        succeeded <- 0;
        group <- name ;
        Printf.printf "Begin test group %s\n" group

    (** mark the end of a group of tests, printing out success count *)
    method end_group  =
        Printf.printf "End test group %s: %d of %d tests passed\n"
            group succeeded n
end;;

(* vim: set syn=ocaml sw=4 ts=4 expandtab: *)

Now it's very easy for me to write tests. For example, tests for all my payoff functions might look like this:

(** payoff_tests.ml - tests for payoff classes

  Written by Sean Hunter <sean@uncarved.com>

  This is demonstration code only.  You are free to use it under the
  terms of the Creative Commons Attribution 2.5 license, but don't
  expect it to accurately price real options.

 *)

let tests = new Test.harness;;

tests#start_group "payoff.ml";;

(* tests for vanilla put and call payoffs *)
tests#ok "Vanilla ITM call" ((Payoff.call 100.0 110.0) = 10.0);;
tests#ok "Vanilla ATM call" ((Payoff.call 100.0 100.0) = 0.0);;
tests#ok "Vanilla OTM call" ((Payoff.call 100.0 90.0) = 0.0);;
tests#ok "Vanilla ITM put" ((Payoff.put 100.0 10.0) = 90.0);;
tests#ok "Vanilla ATM put" ((Payoff.put 100.0 100.0) = 0.0);;
tests#ok "Vanilla OTM put" ((Payoff.put 100.0 190.0) = 0.0);;

(* tests for digital payoffs *)
tests#ok "Digital ITM call" ((Payoff.digital_call 100.0 110.0) = 1.0);;
tests#ok "Digital ATM call" ((Payoff.digital_call 100.0 100.0) = 0.0);;
tests#ok "Digital OTM call" ((Payoff.digital_call 100.0 90.0) = 0.0);;
tests#ok "Digital ITM put" ((Payoff.digital_put 100.0 10.0) = 1.0);;
tests#ok "Digital ATM put" ((Payoff.digital_put 100.0 100.0) = 0.0);;
tests#ok "Digital OTM put" ((Payoff.digital_put 100.0 190.0) = 0.0);;
let dd = Payoff.double_digital ~low:90.0 ~high:110.0;;
tests#ok "Double Digital below the low barrier" ((dd 89.0) = 0.0);;
tests#ok "Double Digital at the low barrier" ((dd 90.0) = 1.0);;
tests#ok "Double Digital inside the low barrier" ((dd 90.1) = 1.0);;
tests#ok "Double Digital inside the high barrier" ((dd 109.9) = 1.0);;
tests#ok "Double Digital at the high barrier" ((dd 110.0) = 1.0);;
tests#ok "Double Digital above the high barrier" ((dd 110.1) = 0.0);;

tests#end_group;;

(* vim: set syn=ocaml sw=4 ts=4 expandtab: *)

...and when I run the tests the output looks like this:

    Begin test group payoff.ml
    00001: Vanilla ITM call ......................................ok
    00002: Vanilla ATM call ......................................ok
    00003: Vanilla OTM call ......................................ok
    00004: Vanilla ITM put .......................................ok
    00005: Vanilla ATM put .......................................ok
    00006: Vanilla OTM put .......................................ok
    00007: Digital ITM call ......................................ok
    00008: Digital ATM call ......................................ok
    00009: Digital OTM call ......................................ok
    00010: Digital ITM put .......................................ok
    00011: Digital ATM put .......................................ok
    00012: Digital OTM put .......................................ok
    00013: Double Digital below the low barrier ..................ok
    00014: Double Digital at the low barrier .....................ok
    00015: Double Digital inside the low barrier .................ok
    00016: Double Digital inside the high barrier ................ok
    00017: Double Digital at the high barrier ....................ok
    00018: Double Digital above the high barrier .................ok
    End test group payoff.ml: 18 of 18 tests passed

Ocaml objects part 2

Parameter objects make more sense when we pass piecewise functions to the pricer

Unfortunately the class for constant parameters that we developed before is just not that useful. We would be better off just using a float if that was all we wanted. The advantage is, however, that we have defined the interface we need for these functions, and what we need is something that can return us the integral of a function and the integral of the square of a function in an interval:

   method integral :       float -> float -> float
   method integral_sq :    float -> float -> float

If we can do that, we can use our function in the mc pricer, so let's define a 'parameters' class for linear functions. This will take two coefficients a and b and our parameter will model the function f(x) = ax + b . The integral of f(x) dx is ax^2 + bx because we're trying to find the area under a triangle for the first term and a rectangle for the second. I needed a little help with f(x)^2 dx (thanks Iain!), but here it is:

(** Parameter class that models a linear function f(x) -> ax + b *)
class parameter_linear' a b =
object(self)
    inherit parameter
    method integral t1 t2 =
      let integrate_to x = 0.5 *. x *. a *. x +. x *. b in
        integrate_to t2 -. integrate_to t1
    method integral_sq t1 t2 =
        a ** 2.0 *. (t2 ** 3.0 -. t1 ** 3.0) +.
        a *. b *. (t2 ** 2.0 -. t1 ** 1.0) +.
        b ** 2.0 *. (t2 -. t1)
end;;

class parameter_linear a b = (parameter_linear' a b : parameter_type);;

This works fine, but what we really want is piecewise functions and piecewise constants. These give us the ability to have parameters that we can calibrate to observable data:

(** For piecewise constants and piecewise linear functions a "piece" is a
* parameter that applies between low and high limits. *)
type piece = Piece of float * float * parameter_type;;

This is what we are going to use to make the pieces of our piecewise functions- a tuple of low and high limits, and a parameter object that can give us the integrals in these limits. Integrating the piecewise function then becomes walking the list of pieces and integrating everything that is in the region we are interested in.

(** Class for piecewise parameters.  Takes a list of pieces which each
* specify the limits and the linear or const parameter in each region. 
*
* At present, does no checking that the pieces are continuous and not
* overlapping *)
class parameter_piecewise' pieces =
  (* there's probably a better way of doing this, but anyway...
   * Apply a function to each piece, with the parameters being the part of
   * the interval from t1 to t2 that is inside the interval of the piece *)
  let visit_pieces pieces fn t1 t2 = 
      let visit_piece piece low high =
        if (low > t2 || high < t1) then
          0.0
        else 
            let range_start = max t1 low in
            let range_end = min t2 high in
            (fn piece) range_start range_end in
      let rec visit_list so_far lis =
        match lis with
            [] -> so_far (** we're done *)
          | Piece(low, high, p) :: rest -> visit_list (so_far +. (visit_piece p low high)) rest
      in  
        visit_list 0.0 pieces
  in         
object(self)
    inherit parameter
    method integral t1 t2 = visit_pieces pieces (fun x->x#integral) t1 t2
    method integral_sq t1 t2 = visit_pieces pieces (fun x->x#integral_sq) t1 t2
end;;                          

(** Piecewise parameters class with implementation hidden *)
class parameter_piecewise x = (parameter_piecewise' x : parameter_type);;

As you can see, I have a local visitor function which walks the list of pieces and calls the right method on each one. I do this using a nifty little tail-recursive pattern-matcher. This works but I am a bit dissatisfied with the way I have to trick visit_pieces into calling the integral or integral_sq methods. There is a better way I'm sure- I just don't know what it is yet.

I could leave this implementation and it would be functionally-complete, but it would be a bit of a pain to build the piecewise functions, so I make a couple of little helper functions:

(** helper funcs to make parts of a piecewise function *)
let make_const_piece low high x = Piece(low, high, new parameter_const x);;
let make_linear_piece low high a b = Piece(low, high, new parameter_linear a b);;

They're still not ideal- I would like to give a list of 2-tuples specifying the upper bound and a point for each piece, and have the function return the piecewise function or constant that "join the dots", but anyway this gives us the ability to make lists of pieces by doing this sort of thing:

(* price one option and print the result  *)
let print_mc ?(num_paths=100000) label payoff =
    let pieces = [
        Param.make_const_piece 0.0 0.25 0.35;
        Param.make_const_piece 0.25 1.0 0.25;
    ] in
    let vol=new Param.parameter_piecewise pieces in
    let mc payoff =
        Mc.sim
            ~payoff:payoff
            ~expiry:0.25
            ~spot:1.613
            ~vol:vol 
            ~r:(new Param.parameter_const 0.055)
            ~num_paths:num_paths
         in
    Printf.printf "%s: %f\n" label (mc payoff)
;;

Note that our piecewise const parameter doesn't require the function to be continuous. Given this vol function, we can price as before.

Ocaml objects part 1

Writing some objects in ocaml, we begin to make the pricer more flexible

The next in the somewhat haphazard series that began here in which I learn a bit about ocaml objects, and in so doing, also learn why Joshi proposes a parameter's class and how to calculate simple integrals.

The purpose of this exercise was to learn how to do object-oriented programming in ocaml, so I pick up where Joshi suggests modifying the Monte Carlo simulator to take "parameters" objects instead of floats for vol and rates. The purpose of this modification is that ultimately we will then be able to use functions as parameters to our models, rather than just constants. I begin by defining the interface to these "parameter" classes. This will define what these parameter objects can do.

(** interface for all parameter classes. Users of parameter classes expect
* this interface *)
class type parameter_type =
object
    method integral :       float -> float -> float
    method integral_sq :    float -> float -> float
    method mean :           float -> float -> float
    method root_mean_sq :   float -> float -> float
end;;

This is a direct translation of Joshi's interface. I'm not 100% sure we really need this, but it does mean that differences between various parameters classes and their private members can be hidden behind this interface.

So given that, we can adapt the mc:

(* Price an option with a flexible payoff using Monte Carlo. 
*
* Use 'parameters' objects for vol and rates to allow passing flexible
* parameters *)
let sim ~payoff ~expiry ~spot ~vol ~r ~num_paths =
    let variance = vol#integral_sq 0.0 expiry in
    let root_variance = sqrt variance in
    let ito_correction = -0.5 *. variance in
    let integral_r = r#integral 0.0 expiry in
    let moved_spot = spot *. exp (integral_r +. ito_correction) in
    let rec do_path i running_sum =
        if i < num_paths then begin
            let this_gaussian = Gaussian.get_one_gaussian () in
            let this_spot = moved_spot *. (exp (root_variance *. this_gaussian)) in
            let this_payoff = payoff ~spot:this_spot in
            do_path (i+1) (running_sum +. this_payoff)
        end
        else (running_sum /. (float_of_int num_paths)) *. (exp (-1.0 *. integral_r))
    in
    do_path 0 0.0
    ;;

as you can see, the difference here is that we are calling the "integral_sq" method on the vol parameter, and the "integral" method on the rates parameter. If we put this into the toploop the type has changed to:

val sim :
  payoff:(spot:float -> float) ->
  expiry:'a ->
  spot:float ->
  vol:< integral_sq : float -> 'a -> float; .. > ->
  r:< integral : float -> 'a -> float; .. > -> num_paths:int -> float = <fun>

So vol and r are params which have an integral and integral_sq method, which have two arguments, and return a float in each case. We could get rid of the polymorphic argument by explicitly specifying the type of the expiry argument, but I'm not sure how to do that for named args. ~expiry:float gives a syntax error anyway- I suppose I could write an explicit sig for the function but it works without it.

Now we can write the base parameters class. This is an abstract class- derived classes will specify how to take the integral and the integral of the square. In this mc model we don't use the concrete methods of this class but presumably they come in handy later.

(** base class specifying the mean and root_mean_sq methods in terms of
  * the integral and integral_sq, but don't define those.  This allows us
  * to define them for any type of parameter *)
class virtual parameter =
object(self)
    method virtual integral :       float -> float -> float
    method virtual integral_sq :    float -> float -> float
    method mean t1 t2 = (self#integral t1 t2) /. (t2 -. t1)
    method root_mean_sq t1 t2 = (self#integral_sq t1 t2) /. (t2 -. t1)
end;;

Now it's simplicity itself to write our first class, which models a constant parameter:

(** class for constant parameters.  This class will be hidden behind the
* parameter_type interface *)
class parameter_const' x =
object(self)
    inherit parameter
    val x_sq = x*.x
    method integral t1 t2 = (t2-.t1) *. x;
    method integral_sq t1 t2 = (t2-.t1) *. x_sq;
end;;

I add an alias to this (without the tick), which hides the implementation behind the parameter_type interface:

(** Constant parameters class with implementation hidden *)
class parameter_const x = (parameter_const' x : parameter_type);;

The type of parameter_const' is

class parameter_const' :
  float ->
  object
    val x_sq : float
    method integral : float -> float -> float
    method integral_sq : float -> float -> float
    method mean : float -> float -> float
    method root_mean_sq : float -> float -> float
  end

whereas the type of parameter_const is

class parameter_const : float -> parameter_type

...which seems to express the idea of the class much better.

At this point we can run our Monte-Carlo using constant parameters by doing this sort of thing:

    let mc payoff =
        Mc.sim
            ~payoff:payoff
            ~expiry:0.25
            ~spot:1.613
            ~vol:(new Param.parameter_const 0.35)
            ~r:(new Param.parameter_const 0.055)
            ~num_paths:num_paths
;;

If you're thinking "so what? That's no better than what we had before except probably a bit slower and certainly harder to read", then I am sympathetic to this point of view. Next we look at other parameters classes and suddenly it starts to become useful.

Testing for put/call parity

Now let's check the prices from the ocaml Monte-Carlo pricer against the most obvious arbitrage relationship

Put call parity is an important relationship in options pricing. Let's test our prices to see if they obey this relationship.

To do so I add a simple helper function to payoff.ml to get the npv of some future cash:

(* get the net present value of some amount *)
let npv ~amt ~rate ~time = amt *. exp(-1.0 *. rate *. time);;

...and the tests themselves. Note that the put/call parity relationship for digitals is slightly different:

let assert_nearly ?(tolerance = 0.001) label a b =
    Printf.printf "Test %s - a: %f b: %f\n" label a b;
    assert( abs_float(a-.b) <tolerance );;

let spot = 1.613 in
let expiry = 0.25 in
let r = 0.055 in
let strike = 1.600 in
let discount x = Payoff.npv ~rate:r ~time:expiry ~amt:x in
let mc payoff =
    Mc1c.sim
        ~payoff:payoff
        ~expiry:expiry
        ~spot:spot
        ~vol:0.35
        ~r:r
        ~num_paths:100000 in
let price payoff = mc (payoff ~strike:strike) in
assert_nearly
    "Put call parity"
    ((price Payoff.call) +. (discount strike))
    ((price Payoff.put) +. spot);
assert_nearly
    ~tolerance:0.01
    "Digital put/digital call parity"
    ((price Payoff.digital_put) +. (price Payoff.digital_call))
    (discount 1.0)
;;

Here's the output:

Test Put call parity - a: 1.707748 b: 1.708199
Test Digital put/digital call parity - a: 0.982665 b: 0.986344

Trebles all round! The prices are ok. It's also interesting that ocaml is actually even faster than I thought. I was running bytecode and if you use ocamlopt to create native code my examples run about five times faster.

Using Ocaml in Practice

Starting to learn how to make ocaml programs that are not just toys, and extending the derivatives pricer further

My experience learning ocaml has been pretty enjoyable so far. A friend told me about Jason Hickey's pdf book Introduction to the Objective Caml Programming Language which is commendably brief and has really helped as I go on to try ocaml further. The second thing that has helped is that I have discovered rlwrap, so the toploop is no longer such an unfriendly place to be.

Furthermore, I have split the pricer that I began developing yesterday into a number of files and build them into a single executable using ocamlc. First, gaussian.ml, containing the random number functions:

open Random;;

(* initialize the random number generator *)
Random.self_init;;

(* get a random gaussian using a Box-Muller transform, described
 * here http://en.wikipedia.org/wiki/Box-Muller_transform *)
let rec get_one_gaussian_by_box_muller () =
    (* Generate two uniform numbers from -1 to 1 *)
    let x = Random.float 2.0 -. 1.0 in
    let y = Random.float 2.0 -. 1.0 in
    let s = x*.x +. y*.y in
    if s > 1.0 then get_one_gaussian_by_box_muller ()
    else x *. sqrt (-2.0 *. (log s) /. s)
    ;;

(* get a gaussian through oversampling and subtraction *)
let get_one_gaussian_by_summation () =
    let rec add_one limit count so_far =
        if count==limit then so_far
        else add_one limit (count+1) (so_far +. (Random.float 1.0)) in
    (add_one 12 0 0.0) -. 6.0
    ;;

let get_one_gaussian = get_one_gaussian_by_box_muller

I added the summation method because when I first tried the pricer on real data the numbers were hopeless (now they are just somewhat out of line with market observables), and I suspected a bug in my random numbers. I was correct, I did have incorrect random numbers.

Then I have payoff.ml, containing my payoff functions. I have added a few more simple payoffs, and moved to named function arguments:

(** a vanilla option pays off the difference between the spot price
 ** and the strike, or expires worthless *)
let call ~strike ~spot = max (spot -. strike) 0.0;;
let put ~strike ~spot = max (strike -. spot) 0.0;;

let digital payoff = if payoff> 0.0 then 1.0 else 0.0;;
let digital_call ~strike ~spot = digital (call ~strike:strike ~spot:spot);;
let digital_put ~strike ~spot = digital (put ~strike:strike ~spot:spot);;

(** A double digital pays 1 if spot is between two barriers, zero
 ** otherwise *)
let double_digital ~low ~high ~spot =
    assert (low < high);
    if (low <= spot && spot <= high) then 1.0
    else 0.0;;

mc1c.ml contains the actual Monte Carlo simulator, and it is unchanged except to use named function arguments, and to qualify the name of the get_one_gaussian function, which is now in a seperate file:

(* Price an option with a flexible payoff using Monte Carlo. *)
let sim ~payoff ~expiry ~spot ~vol ~r ~num_paths =
    let variance = vol *. vol *. expiry in
    let root_variance = sqrt variance in
    let ito_correction = -0.5 *. variance in
    let moved_spot = spot *. exp (r *. expiry +. ito_correction) in
    let rec do_path i running_sum =
        if i < num_paths then begin
            let this_gaussian = Gaussian.get_one_gaussian () in
            let this_spot = moved_spot *. (exp (root_variance *. this_gaussian)) in
            let this_payoff = payoff ~spot:this_spot in
            do_path (i+1) (running_sum +. this_payoff)
        end
        else (running_sum /. (float_of_int num_paths)) *. (exp (-1.0 *. r *. expiry))
    in
    do_path 0 0.0
    ;;

Finally I have my test file which runs the test cases. Now that I use named function args I can partially-apply function args in any order, so I make a test harness that sets up a particular marketdata scenario and runs the pricer:

let print_mc ?(num_paths=100000) label payoff =
    let mc payoff =
        Mc1c.sim
            ~payoff:payoff
            ~expiry:0.2
            ~spot:161.3
            ~vol:0.35
            ~r:0.045
            ~num_paths:num_paths
         in
    Printf.printf "%s: %f\n" label (mc payoff)
;;

print_mc "call" (Payoff.call ~strike:160.0);;
print_mc "digital call" (Payoff.digital_call ~strike:160.0);;
print_mc "put" (Payoff.put ~strike:170.0);;
print_mc "digital put" (Payoff.digital_put ~strike:170.0);;
print_mc "double digital" (Payoff.double_digital ~low:160.0 ~high:170.0) ~num_paths:250000;;

(* price one option, test the payoff against a target price and 
 * print the result  *)
let test_mc ?(num_paths=1000) ?(expiry=1.0) ?(r=0.0) label payoff target =
    let mc payoff =
        Mc1c.sim
            ~payoff:payoff
            ~expiry:expiry
            ~spot:161.3
            ~vol:0.35
            ~r:r
            ~num_paths:num_paths
         in
    let price = mc payoff in
    let tolerance = 0.00001 in
    Printf.printf "Test %s - price: %f target: %f\n" label price target;
    assert( abs_float(price-.target) < tolerance )
;;

test_mc "Spot to one" (fun ~spot->1.0) 1.0;;
test_mc "Spot to zero" (fun ~spot->0.0) 0.0;;
let r = 0.05 in
let npv =  exp(-1.0 *. r) in
    test_mc "Spot to one with rates" (fun ~spot->1.0) npv ~r:r;;

(* vim: set sw=4 ts=4 expandtab: *)

This is pretty cool. As you can see, I made the number of mc paths an optional parameter. I have added a few assertions that trivial payoff functions price correctly and that discounting works.

Now I build an object from each ml file using ocamlc and then compile them all to a built object at the end. You need to be careful to do them all in the correct order where you have functions defined in one file, called in another, as we have here.

Partial Function Application is not Currying

These two related concepts are often confused. Until yesterday I thought they were the same thing myself.

Often you will see in computer books and articles a pattern where a function is applied to some but not all of it's required arguments, resulting in a function of fewer arguments. In python, this looks like this (from PEP 309):

def papply(fn, *cargs, **ckwargs):
    def call_fn(*fargs, **fkwargs):
        d = ckwargs.copy()
        d.update(fkwargs)
        return fn(*(cargs + fargs), **d)
    return call_fn

This is called "partial function application" - it returns a function which acts like the function you pass in but which takes fewer arguments, the others having been "bound in". The author of this code, however, had the (very common) misconception that this is currying, and called his function "curry" as a result. I shared this misconception for some time, and thought that currying and partial application were the same thing. In fact they are to certain extent opposites.

Where partial application takes a function and from it builds a function which takes fewer arguments, currying builds functions which take multiple arguments by composition of functions which each take a single argument. Thus we curry in python like this:

def addN(n):
    return lambda x: x + n

def plus(a, b):
    addA=addN(a)
    return addA(b)

Now why would we ever want to do that? Well, in some pure functional languages this is exactly how functions with multiple arguments are built up. In ocaml, a function which takes two ints and returns a float is actually a function which takes an int and returns a function which takes an int and returns a float. In this world, partial application just happens without any extra code:

% rlwrap ocaml
    Objective Caml version 3.09.3

# let add a b=a+b;;
val add : int -> int -> int = <fun>

So the type of add is a function which takes an int and returns a function which takes an int and returns an int.

# let add2=add 2;; 
val add2 : int -> int = <fun>

add is a curried function, so here we can partially apply by just calling with a single arg- it returns the function that takes the other arg and returns the result.

# add2 34;;
- : int = 36

...and we can call add2 with a single argument as you would expect. Because ocaml curries add for us, the function has been partially applied. It's interesting to note that in ocaml if you label your function arguments, they can be partially applied in any order.

Derivatives pricing in ocaml Part 2

Extending the basic mc pricer to handle different payoffs, we see how partial function application works

This is the third article in a series on using functional programming for financial applications which started here.

So given our first Monte Carlo simulator, and the payoff functions that we started with, it's easy to see how we extend the pricer to handle Joshi's first question at the end of chapter 1 (to price puts):

(* mc1b.ml - A rudimentary option pricer to answer the exercises at the
 * end of chapter 1 in Joshi.
 *
 * Written by Sean Hunter <sean@uncarved.com>
 *
 * This is demonstration code only.  You are free to use it under the
 * terms of the Creative Commons Attribution 2.5 license, but don't
 * expect it to accurately price real options.
 *)
open Random;;
open Printf;;

(* initialize the random number generator *)
Random.self_init;;

(* get a random gaussian using a Box-Muller transform, described
 * here http://en.wikipedia.org/wiki/Box-Muller_transform *)
let rec get_one_gaussian_by_box_muller () =
    (* Generate two uniform numbers from -1 to 1 *)
    let x = Random.float 2.0 -. 1.0 in
    let y = Random.float 2.0 -. 1.0 in
    let s = x*.x +. y*.y in
    if s > 1.0 then get_one_gaussian_by_box_muller ()
    else x *. sqrt (-2.0 *. (log s) /. s)
    ;;

 (* a vanilla option pays off the difference between the spot price
 * and the strike, or expires worthless *)
let put_payoff strike spot =
    max ( strike -. spot ) 0.0;;

let call_payoff strike spot =
    max (spot -. strike ) 0.0;;

(* Price an option with a flexible payoff using Monte Carlo. *)
let simple_monte_carlo_1a payoff expiry strike spot vol r num_paths =
    let variance = vol *. vol *. expiry in
    let root_variance = sqrt variance in
    let ito_correction = -0.5 *. variance in
    let moved_spot = spot *. exp (r *. expiry +. ito_correction) in
    let rec do_path i running_sum =
        if i < num_paths then begin
            let this_gaussian = get_one_gaussian_by_box_muller () in
            let this_spot = moved_spot *. (exp (root_variance *. this_gaussian)) in
            let this_payoff = payoff strike this_spot in
            do_path (i+1) (running_sum +. this_payoff)
        end
        else (running_sum /. (float_of_int num_paths)) *. (exp (-1.0 *. r *. expiry))
    in
    do_path 0 0.0
    ;;

(* price one put and one call option struck at the money *)    
printf "%f\n" (simple_monte_carlo_1a call_payoff 0.025 195.5 195.5 0.20 0.045 100000);;
printf "%f\n" (simple_monte_carlo_1a put_payoff 0.025 195.5 195.5 0.20 0.045 100000);;

So we pass a payoff function into the Monte Carlo simulator and call that function for each path. This isn't quite general enough to handle double digitals though. A double digital is an option that pays 1 if the spot is between two barrier prices at expiry and zero otherwise. However, the payoff function here has to take the strike and the spot. This is where one of the great features of functional programming comes in: partial function application. In ocaml, if you call a function that takes 2 parameters, but give it only one, the return type is a function that takes one parameter. Amazingly partial function application is just there by default and there's no need for tedious binders like there are in the STL in C++ for example. This will explain:

% ocaml                                                                                                                          
    Objective Caml version 3.09.3

# let myadd x y=x+y;;
val myadd : int -> int -> int = <fun>
# let add2=myadd 2;;
val add2 : int -> int = <fun>
# add2 5;;
- : int = 7
# add2 16;;
- : int = 18

So myadd adds two numbers, and by calling it with just one (a 2) we create a function that takes one argument and adds two to it. This is exactly what we need for our flexible payoff function. You can think of the partially applied function args as being all the things that are constant on the termsheet of the trade. Our payoff functions remain the same, except we add one for double digitals:

let double_digital_payoff low high spot =
        if (low <= spot && spot <= high) then 1.0
        else 0.0;;

...and we change the mc pricer to just call the payoff func with the spot on the current mc path. We will partially-apply any other arguments the payoff functions need when we invoke the pricer:

(* Price an option with a flexible payoff using Monte Carlo. *)
let simple_monte_carlo_1b payoff expiry spot vol r num_paths =
    let variance = vol *. vol *. expiry in
    let root_variance = sqrt variance in
    let ito_correction = -0.5 *. variance in
    let moved_spot = spot *. exp (r *. expiry +. ito_correction) in
    let rec do_path i running_sum =
        if i < num_paths then begin
            let this_gaussian = get_one_gaussian_by_box_muller () in
            let this_spot = moved_spot *. (exp (root_variance *. this_gaussian)) in
            let this_payoff = payoff this_spot in
            do_path (i+1) (running_sum +. this_payoff)
        end
        else (running_sum /. (float_of_int num_paths)) *. (exp (-1.0 *. r *. expiry))
    in
    do_path 0 0.0
    ;;

Now see how we call this pricer. It's simplicity itself:

printf "%f\n" (simple_monte_carlo_1b (call_payoff 160.0) 0.2 161.3 0.35 0.045 250000);;
printf "%f\n" (simple_monte_carlo_1b (put_payoff 170.0) 0.2 161.3 0.31 0.045 250000);;
printf "%f\n" (simple_monte_carlo_1b (double_digital_payoff 160.0 170.0) 0.2 161.3 0.29 0.045 250000);;

I find this immensely pleasing. I have hardly started learning the language, and yet generalising this code was simplicity itself due to features in the language. Partial function application is available in imperative languages (C++ notably uses it to provide predicates and adaptable functions in the STL), but there is a lot of nasty additional syntax. The way it works so seemlessly in ocaml is terrific.

Its worth coming back to earth a little bit with the observation that the prices that I am seeing from this thing right now don't match those I can observe in the market so I am sure there is some debugging yet to do. I would also like to make the code into a few modules, but I am not sure how you do that in ocaml yet.

Derivatives pricing in ocaml Part 1

Following on from the "Learning Functional Programming" article, we write our first call pricer

So from the first article we move on to write the first Monte Carlo pricer. This is a very close translation of Joshi's implementation in part 1.3 (listings 1.1 to 1.3) except that I don't read the values in from stdin, I just hardcode them.

Writing this was fun, but I am extremely short of sleep so it was actually trickier than I expected. It was made harder by the cryptic error messages that ocaml gives you when things go wrong. These were almost always This expression has type int but is here used with type float which means you are passing int arguments to a float operator or function, rather than the other way around.

% ocaml                                                           
    Objective Caml version 3.09.3

# 1 /. 2;;
This expression has type int but is here used with type float
# 1.0 / 2.0;;
This expression has type float but is here used with type int

I hope that makes it a little clearer. I find it confusing anyway. It's like the error message has been written by someone for whom English is a second language. Perhaps that's the case.

Here's the listing without further ado. On my computer it does a million paths in 13.2 seconds which is pretty decent I think. The ocamlc compiler is fast as lightning too, compiling this to an executable in 0.029 secs on my machine. If I go on learning ocaml I won't be able to go and get coffee while waiting for my code to compile.

   (* mc1.ml - A rudimentary European call option pricer designed to mimic 
     * listings 1.1 to 1.3 in Joshi.
     *
     * Written by Sean Hunter <sean@uncarved.com>
     *
     * This is demonstration code only.  You are free to use it under the
     * Creative Commons Attribution 2.5 license, but don't expect it to
     * accurately price real options.
     *)
    open Random;;
    open Printf;;

    (* initialize the random number generator *)
    Random.self_init;;

    (* get a random gaussian using a Box-Muller transform, described
     * here http://en.wikipedia.org/wiki/Box-Muller_transform *)
    let rec get_one_gaussian_by_box_muller () =
        (* Generate two uniform numbers from -1 to 1 *)
        let x = Random.float 2.0 -. 1.0 in
        let y = Random.float 2.0 -. 1.0 in
        let s = x*.x +. y*.y in
        if s > 1.0 then get_one_gaussian_by_box_muller ()
        else x *. sqrt (-2.0 *. (log s) /. s)
        ;;

    (* Price a European call using Monte Carlo. 
     * 
     * We pre-compute as much as possible before the simulation, then the 
     * actual mc paths are done as a nested recursive function.  This seems 
     * more idiomatically functional even though ocaml has for loops.
     *
     * It's also good because I couldn't get the other way to work.*)    
    let simple_monte_carlo1 expiry strike spot vol r num_paths =
        let variance = vol *. vol *. expiry in
        let root_variance = sqrt variance in
        let ito_correction = -0.5 *. variance in
        let moved_spot = spot *. exp (r *. expiry +. ito_correction) in
        let rec do_path i running_sum =
            if i < num_paths then begin
                let this_gaussian = get_one_gaussian_by_box_muller () in
                let this_spot = moved_spot *. (exp (root_variance *. this_gaussian)) in
                let this_payoff = max (this_spot -. strike) 0.0 in
                do_path (i+1) (running_sum +. this_payoff)
            end
            else (running_sum /. (float_of_int num_paths)) *. (exp (-1.0 *. r *. expiry))
        in
        do_path 0 0.0
        ;;

    (* price a three-month call option near to the money.
     * we are using 35% vol and 4.5% interest rates*)    
    printf "%f\n" (simple_monte_carlo1 0.2 160.0 161.3 0.35 0.045 250000);;

This, however, is crying out for more flexibility. The first thing to do is to see if we can get it to price other simple payoffs. That's the subject of the next article.

Learning Functional Programming

I decided to teach myself ocaml by writing a derivatives pricer

I have been interested in functional programming for some time, but finally decided to bite the bullet and learn properly, and that (for me anyway) means writing some code to accomplish a practical task. My idea is to reimplement most of the C++ code in Mark Joshi's excellent book C++ Design Patterns and Derivatives Pricing, but in ocaml. Now I'm a newcomer to functional programming and to ocaml, so what I write won't be pretty or idiomatic, especially at first. To begin with I'm learning from the main tutorial. If I find and use another I'll post that too.

I will try to explain some of the financial stuff that's going on as I do so, but for the full lowdown on why derivatives price the way they do, you need to learn some financial maths. You could do a lot worse that picking up Joshi's other book, The Concepts and Practice of Mathematical Finance. A lot of people get Hull, but I prefer Joshi for a really practical introduction with great explanations of concepts.

So without further ado, here is my first ocaml program, which defines payoffs for vanilla European put and call options. In fact Joshi starts off straight away with a Monte Carlo pricer, but my copy is downstairs so I'm straying off-piste here. It's my intention to follow Joshi step by step, and write up each one here as I go.

    open Printf

    (* a vanilla option pays off the difference between the spot price
     * and the strike, or expires worthless *)
    let put_payoff strike spot=
            max ( strike -. spot ) 0.0;;

    let call_payoff strike spot=
            max (spot -. strike ) 0.0;;

    let print_payoff payoff strike spot=
            let outcome=payoff strike spot in
            printf "%f\n" outcome;;

    print_payoff call_payoff 195.0 190.0;;
    print_payoff call_payoff 195.0 200.0;;
    print_payoff put_payoff 195.0 190.0;;
    print_payoff put_payoff 190.0 195.0;;

Now I'm running and writing this on fedora Linux, and my ocaml is 3.09.3. When I run this I see:

% ocaml tmp/payoff.ml
5.000000
0.000000
0.000000
5.000000

Which is what I would expect. Now this is very cheesy at present, but it's a start and we'll improve it in the next article. It's worth a couple of observations about this because already this demonstrates a few things that strike me as being interesting about ocaml. For one thing, there isn't any default type promotion or operator overloading, so we need to explicitly qualify our constants with .0 to get them to be floats. Secondly, we need to use -. to subtract them. The max function can operate on any type so it works with floats or ints.