The Mad Economist

Well, I have finally decided to write some econometrics routines in F#.  I first did this type of stuff in FORTRAN when I was a grad student, so I figured I'd try to see if I could still figure it all out.

But, first things first.

We need routines to Transpose a matrix and to Multiply two matrices.

My transpose method is this:

--------------------------------------------------------------
//
//  Let's make up some data
//

let n = 6
let k = 3

//  Create data array - n x k

let x0 = Array2D.zeroCreate<float> n k

x0.[0,0] <- 1.0
x0.[1,0] <- 1.0
x0.[2,0] <- 1.0
x0.[3,0] <- 1.0
x0.[4,0] <- 1.0
x0.[5,0] <- 1.0

x0.[0,1] <- 22.0
x0.[1,1] <- 4.0
x0.[2,1] <- 13.0
x0.[3,1] <- 3.0
x0.[4,1] <- 9.0
x0.[5,1] <- -3.0

x0.[0,2] <- 12.0
x0.[1,2] <- 14.0
x0.[2,2] <- 19.0
x0.[3,2] <- 9.0
x0.[4,2] <- 13.0
x0.[5,2] <- 23.0

//---------------------------------------------------
// Transpose a Matrix
// Takes matrix x and returns the transpose as xp
// NEWLY UPDATED VERSION 8/12/2009
//---------------------------------------------------

let rec transposeACol(x : float [,] ) (xp : float [,] ) (i : int) (j : int) =
    match j with
    | -1 -> xp
    | _ -> xp.[j, i] <- x.[i, j]
           //printfn "i: %i j: %i x: %10.5f" i j (if i >=0 && j >=0 then x.[i, j] else -1.0)
           //let userresp = Console.ReadLine()
           let xp = transposeACol x xp i (j-1)
           xp
          
let rec transpose (x : float [,] ) (xp : float [,] ) (i : int) (j : int) =
    match i with
    | -1 -> xp
    | _  -> let xp = transposeACol x xp i j 
            let xp = transpose x xp (i-1) j
            xp

// Create a shell array, k x n
let zero_xp = Array2D.zeroCreate<float> k n

//-----------------------------------------
// Create X Prime
//-----------------------------------------
//
// Handle the fact that matrices are zero based: use n-1, k-1 as starting points
//

let xp = transpose x0 zero_xp (n-1) (k-1)

//
//Print xp
//
// Print it out
for i = 0 to k - 1 do
     printfn "%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f" xp.[i,0] xp.[i,1] xp.[i,2] xp.[i,3] xp.[i,4] xp.[i,5]

--------------------------------------------------------------

The method transpose is a recursive method that starts at the top for both n and k and works down to zero, flipping the elements as it goes.  You send it a blank array and it returns the transpose.  I know I could have done this with a couple of for..do loops, but I wanted to practice writing recursive functions in F#.

Note the Array function Array2D.zeroCreate<float> k n.  This function creates a two dimensional array (floating point) with k rows and n columns.  It is prefilled with zeros.  The array is zero-based.  There is also a method Array.create that allows you to specify the initial value (instead of zero).

 So we run the code, and here's what we get:

      1.00       1.00       1.00       1.00       1.00       1.00
     22.00       4.00      13.00       3.00       9.00      -3.00
     12.00      14.00      19.00       9.00      13.00      23.00

The next thing we need to code is a matrix multiplication routine.

If you aren't familiar with matrix multiplication, here is Wikipedia's page: http://en.wikipedia.org/wiki/Matrix_multiplication

The idea here is that you have a matrix that is dimensioned n x k and you want to multiply it bay a matrix of dimension k x m.  The two k's must be the same.  The resulting matrix is dimensioned n x m.

So, for example, if we have two matrices xp (2 x 4) and x (4 x 2), the resulting matrix (xpx) will have be dimensioned 2 x 2 and the i,jth item will be:

xpx[i, j] = (xp[i,0] * x[0,j]) + (xp[i,1] * x[1,j]) + (xp[i,2] * x[2,j]) + (xp[i,3] * x[3,j])

The routine goes as follows:

The function accum, just takes two values, multiplies them, and adds then to an accumulator value (ac) that is passed in.  Obviously the accumulator value would be initialized at 0.0

The function xply iterates over a row (row j0)  in the left matrix and a column (j1) in the right matrix. calling accum to accumulate the result.

Note that in the xply function, this line

    | _  -> xply xp x j0 j1 (i-1) (accum xp.[j0, i] x.[i, j1]  cell_accum)

could be written in two lines as:

    | _  -> let tmp = accum xp.[j0, i] x.[i, j1]  cell_accum
            xply xp x j0 j1 (i-1) tmp

but the former is somewhat more 'elegant.'  The return value from xply is the value for one cell (j0, j1) the resultant matrix from the multiplication operation.

The function xMultiply is a recursive function that iterates over all the cells in the matrix, calling xply each time to fill in a cell.  It has two matching levels, one outer level that iterates over the rows; one inner level that iterates over the rows in a given column.

The top level function  matMult just takes the two input matrices (left & right) and a empty output matrix as inputs.  Note the method GetUpperBound i.  This method returns the upper bound for a matrix for the dimension you pass in.  In a 2D array, The 0th dimension is for the rows, the 1st dimension is for the columns.  This is a very handy function because it means the user does not have to pass in the dimensions of the matrices to be multiplied.

The out of the matrix multiplication is:

  6.00      48.00      90.00
 48.00     768.00     642.00
 90.00     642.00    1480.00