2 people like it.

Cumulative bivariate normal distribution

Cumulative bivariate normal distribution in F#. I implemented this using both the code given in Haug's complete guide to option pricing formulae book (which in turn credits his version to Graeme West who in term adapted code from Alan Genz) and using the code available on Alan Genz' website. This code requires the availability of another function called cnd which implements the (univariate) cumulative normal distribution. I haven't included my code there, since that's readily available in .NET from various sources (e.g. Math.NET, alglib, ...)

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
42: 
43: 
44: 
45: 
46: 
47: 
48: 
49: 
50: 
51: 
52: 
53: 
54: 
55: 
56: 
57: 
58: 
59: 
60: 
61: 
62: 
63: 
64: 
65: 
66: 
67: 
68: 
69: 
70: 
71: 
72: 
73: 
74: 
75: 
76: 
77: 
78: 
79: 
80: 
81: 
82: 
83: 
84: 
85: 
86: 
87: 
88: 
89: 
90: 
91: 
92: 

/// Some constants required
let private weights1 = [ 0.17132449237917; 0.17132449237917;0.360761573048138;
                            0.360761573048138;0.46791393457269;0.46791393457269 ]
let private values1 =  [ 0.966234757101576;0.033765242898424;0.830604693233133;
                            0.169395306766867;0.619309593041598;0.380690406958402 ]

let private weights2 = [ 0.0471753363865118; 0.0471753363865118; 0.106939325995318; 0.106939325995318;
                            0.160078328543346;  0.160078328543346;  0.203167426723066; 0.20316742672306; 
                            0.233492536538355;  0.233492536538355;  0.249147045813403; 0.249147045813403]
let private values2 =  [ 0.99078031712336; 0.00921968287664;0.95205862818524; 0.04794137181476;
                            0.88495133709715;0.11504866290285;0.79365897714331;0.20634102285669;
                            0.68391574949909;0.31608425050091;0.56261670425574;0.43738329574427]


let private weights3 = [0.0176140071391521; 0.0176140071391521; 0.0406014298003869;
                        0.0406014298003869; 0.0626720483341091; 0.0626720483341091;
                        0.0832767415767048; 0.0832767415767048; 0.1019301198172400;
                        0.1019301198172400; 0.1181945319615180; 0.1181945319615180;
                        0.1316886384491770; 0.1316886384491770; 0.1420961093183820;
                        0.1420961093183820; 0.1491729864726040; 0.1491729864726040;
                        0.1527533871307260; 0.1527533871307260]
let private values3 = [ 0.99656429959255; 0.00343570040745; 0.98198596363896;
                        0.01801403636104; 0.95611721412566; 0.04388278587434;
                        0.91955848591111; 0.08044151408889; 0.87316595323008;
                        0.12683404676992; 0.81802684036326; 0.18197315963674;
                        0.75543350097541; 0.24456649902459; 0.68685304435771;
                        0.31314695564229; 0.61389292557082; 0.38610707442918;
                        0.53826326056675; 0.46173673943325 ]

let private values4 = [ 0.000047216149159;3.972561612889520;0.001298022024068;
                        3.857185731135710;0.007702795584372;3.656640508589690;
                        0.025883348755653;3.382351236044530;0.064296126805960;
                        3.050028889390040;0.136503050785829;2.658650279846430;
                        0.239251089780572;2.282719097583890;0.392244063312137;
                        1.887068418173820;0.596314691697034;1.507458096263630;
                        0.984753258857668;1.015363867311070 ]


/// Cumulative bivariate normal distribution, based from Alan Genz's matlab code and Haug's VBA code in the complete guide to option pricing formulae (all error's are mine of course)
let public cbnd = fun x y rho ->
    let rho = min 1. (max -1. rho)
    let h = -x
    let k = if rho < -0.925 then y else -y
    let hk = h * k

    match rho with
        | x when x = 0.         ->  cnd(-h) * cnd(-k)
        | x when rho = 1.       ->  cnd(-(max h k))
        | x when rho = -1.      ->  if k > h then cnd(k)-cnd(h) else 0.
        | x when abs(x) < 0.925 ->  let weights, values = match x with
                                                            | r when abs(r) < 0.3    -> weights1,values1
                                                            | r when abs(r) < 0.75   -> weights2,values2
                                                            | _                      -> weights3,values3
                                    let hs = 0.5*(h*h+k*k)
                                    let asr_ = System.Math.Asin(rho)
                                    values  |> List.map(fun x -> let sn = System.Math.Sin(asr_ * x)
                                                                 exp((sn*hk-hs)/(1.-sn*sn)) )
                                            |> List.zip weights
                                            |> List.sumBy (fun (w,v) -> w*v)
                                            |> fun x -> x * asr_/(4.*System.Math.PI) + cnd(-h) *cnd(-k)

        | _                     ->  // 0.925 < abs(x) < 1
                                    let ass = (1.-rho)*(1.+rho)
                                    let A = sqrt(ass)
                                    let bs = (h-k)*(h-k)
                                    let c = (4.-hk) / 8.
                                    let d = (12.-hk)/16.
                                    let asr_ = -0.5*(bs / ass + hk)
                                    let BVN = if asr_ <= -100. then 0.
                                              else A*exp(asr_)*(1.-c*(bs-ass)*(1.-0.2*d*bs) / 3. + 0.2*c* d*ass*ass)
                                    let BVN = if -hk >= 100. then BVN
                                              else
                                                  let B = sqrt(bs)
                                                  BVN - sqrt(2.*System.Math.PI)*exp(-0.5*hk)*cnd(-B/A)*B*(1.-c*bs*(1.-0.2*d*bs)/3.)    

                                    let A = 0.5*A
                                    let A2 = A*A
                                    let values = values4

                                    let weights = weights3

                                    values  |> List.zip weights
                                            |> List.sumBy(fun (w,x) ->  let xs = A2*x
                                                                        let asr_ = -0.5*(bs / xs + hk)
                                                                        if asr_ < -100. then 0.
                                                                        else
                                                                            let rs = sqrt(1.-xs)
                                                                            A*w*exp(asr_)*(exp(-hk*(1.-rs)/(2.*(1.+rs)))/rs - (1.+c*xs*(1.+d*xs))))
                                            |> fun x -> let BVN = -(BVN+x)/(2.*System.Math.PI)
                                                        match rho > 0. with
                                                            | true  ->  BVN + cnd(-(max h k))
                                                            | false ->  if k > h then -BVN + cnd(k)-cnd(h) else -BVN
val private weights1 : float list

Full name: Script.weights1


 Some constants required
val private values1 : float list

Full name: Script.values1
val private weights2 : float list

Full name: Script.weights2
val private values2 : float list

Full name: Script.values2
val private weights3 : float list

Full name: Script.weights3
val private values3 : float list

Full name: Script.values3
val private values4 : float list

Full name: Script.values4
val cbnd : x:float -> y:float -> rho:float -> float

Full name: Script.cbnd


 Cumulative bivariate normal distribution, based from Alan Genz's matlab code and Haug's VBA code in the complete guide to option pricing formulae (all error's are mine of course)
val x : float
val y : float
val rho : float
val min : e1:'T -> e2:'T -> 'T (requires comparison)

Full name: Microsoft.FSharp.Core.Operators.min
val max : e1:'T -> e2:'T -> 'T (requires comparison)

Full name: Microsoft.FSharp.Core.Operators.max
val h : float
val k : float
val hk : float
val abs : value:'T -> 'T (requires member Abs)

Full name: Microsoft.FSharp.Core.Operators.abs
val weights : float list
val values : float list
val r : float
val hs : float
val asr_ : float
namespace System
type Math =
  static val PI : float
  static val E : float
  static member Abs : value:sbyte -> sbyte + 6 overloads
  static member Acos : d:float -> float
  static member Asin : d:float -> float
  static member Atan : d:float -> float
  static member Atan2 : y:float * x:float -> float
  static member BigMul : a:int * b:int -> int64
  static member Ceiling : d:decimal -> decimal + 1 overload
  static member Cos : d:float -> float
  ...

Full name: System.Math
System.Math.Asin(d: float) : float
Multiple items
module List

from Microsoft.FSharp.Collections

--------------------
type List<'T> =
  | ( [] )
  | ( :: ) of Head: 'T * Tail: 'T list
  interface IEnumerable
  interface IEnumerable<'T>
  member Head : 'T
  member IsEmpty : bool
  member Item : index:int -> 'T with get
  member Length : int
  member Tail : 'T list
  static member Cons : head:'T * tail:'T list -> 'T list
  static member Empty : 'T list

Full name: Microsoft.FSharp.Collections.List<_>
val map : mapping:('T -> 'U) -> list:'T list -> 'U list

Full name: Microsoft.FSharp.Collections.List.map
val sn : float
System.Math.Sin(a: float) : float
val exp : value:'T -> 'T (requires member Exp)

Full name: Microsoft.FSharp.Core.Operators.exp
val zip : list1:'T1 list -> list2:'T2 list -> ('T1 * 'T2) list

Full name: Microsoft.FSharp.Collections.List.zip
val sumBy : projection:('T -> 'U) -> list:'T list -> 'U (requires member ( + ) and member get_Zero)

Full name: Microsoft.FSharp.Collections.List.sumBy
val w : float
val v : float
field System.Math.PI = 3.14159265359
val ass : float
val A : float
val sqrt : value:'T -> 'U (requires member Sqrt)

Full name: Microsoft.FSharp.Core.Operators.sqrt
val bs : float
val c : float
val d : float
val BVN : float
val B : float
val A2 : float
val xs : float
val rs : float
Raw view Test code New version

More information

Link:http://fssnip.net/jv
Posted:11 years ago
Author:Bram Jochems
Tags: scientific computing

This web site is created using F# and Suave web server. It is hosted on Azure and the source code is on GitHub. Contributions are welcome!

The first version of fssnip.net has been created by @tomaspetricek back in 2010. This web site is a new, open-source and contribution-friendly version.