Somers'D and ASE1 - dirty implementation in F#

// uses F# PowerPack

#r "FSharp.PowerPack.dll"

let PD             = [ 0.05; 0.10; 0.50; 1.00; 2.00; 5.00; 25.00 ]
let counterparties = [ 5.; 10.; 20.; 25.; 20.; 15.; 5. ]

let groups = PD.Length // risk groups no.

let div100 x = x / 100.0

let PDprct = [ for x in PD do yield div100 x ]
let CPprct = [ for x in counterparties do yield div100 x ]

let n = CPprct |> Seq.sum

let defaulted    = [ for i in 1..groups do yield CPprct.[i-1]*PDprct.[i-1]/n ]
let nondefaulted = [ for i in 1..groups do yield CPprct.[i-1]*(1.0-PDprct.[i-1])/n ]

let x = matrix [ defaulted; nondefaulted ]

let xr = x.NumRows

let rowSum (x : matrix) (i : int) = Array.sum (RowVector.toArray (x.Row(i-1)))

// WR / DR

let wr =

    let mutable mat_sum = 0.0

    for i in 1..xr do
        let row2  = rowSum x i ** 2.0
        mat_sum   <- mat_sum + row2

    n ** 2.0 - mat_sum

// A / D

let A (x : matrix) i j =

    let xr = x.NumRows
    let xc = x.NumCols

    let rowIdx1 = List.filter (fun x -> x>i) [ 1..xr ]
    let colIdx1 = List.filter (fun x -> x>j) [ 1..xc ]

    let rowIdx2 = List.filter (fun x -> x<i) [ 1..xr ]
    let colIdx2 = List.filter (fun x -> x<j) [ 1..xc ]

    let mutable Asum = 0.0

    for r_i in rowIdx1 do
        for r_j in colIdx1 do
            Asum <- Asum + x.[r_i-1,r_j-1]

    for r_i in rowIdx2 do
        for r_j in colIdx2 do
            Asum <- Asum + x.[r_i-1,r_j-1]

    Asum

let D (x : matrix) i j =

    let xr = x.NumRows
    let xc = x.NumCols

    let rowIdx1 = List.filter (fun x -> x>i) [ 1..xr ]
    let colIdx1 = List.filter (fun x -> x<j) [ 1..xc ]

    let rowIdx2 = List.filter (fun x -> x<i) [ 1..xr ]
    let colIdx2 = List.filter (fun x -> x>j) [ 1..xc ]

    let mutable Dsum = 0.0

    for r_i in rowIdx1 do
        for r_j in colIdx1 do
            Dsum <- Dsum + x.[r_i-1,r_j-1]

    for r_i in rowIdx2 do
        for r_j in colIdx2 do
            Dsum <- Dsum + x.[r_i-1,r_j-1]

    Dsum

let Pf (x : matrix) =

    let xr = x.NumRows
    let xc = x.NumCols

    let mutable Pfsum = 0.0

    for i in 1..xr do
        for j in 1..xc do
            Pfsum <- Pfsum + x.[i-1,j-1] * A x i j

    Pfsum

let Qf (x : matrix) =

    let xr = x.NumRows
    let xc = x.NumCols

    let mutable Qfsum = 0.0

    for i in 1..xr do
        for j in 1..xc do
            Qfsum <- Qfsum + x.[i-1,j-1] * D x i j

    Qfsum

let ASE1 (x : matrix) =

   let xr = x.NumRows
   let xc = x.NumCols

   let mutable sumX = 0.0

   for i in 1..xr do
        for j in 1..xc do
            sumX <- sumX + x.[i-1,j-1]

   let mutable tmp = 0.0

   for i in 1..xr do
        for j in 1..xc do
            let sub_sum = sumX - rowSum x i
            tmp <- tmp + x.[i-1,j-1] * ( wr * (A x i j - D x i j) - (Pf x - Qf x) * sub_sum ) ** 2.0

   2.0 / wr ** 2.0 * sqrt tmp

//

let SomersD = (Pf x - Qf x)/wr

let SomersD_ASE1 = ASE1 x