The Mad Economist

[ n₁  0  0    0 ]
|    0 n₂  0  0 |
|    0  0 n₃  0 |
[    0  0  0 n₄ ]

let nObs = Array2D.zeroCreate<float> ni 1  // ni is the record count

nObs.[0,0] <- 10.0
nObs.[1,0] <- 20.0
nObs.[2,0] <- 5.0
nObs.[3,0] <- 15.0
nObs.[4,0] <- 22.0
nObs.[5,0] <- 11.0
nObs.[6,0] <- 4.0
nObs.[7,0] <- 17.0
nObs.[8,0] <- 11.0
nObs.[9,0] <- 1.0

Next we need to make up a matrix Z that is n x n, where the diagonal is filled with the values 1/n_i.

Again we can use our trick method Array2D.init:

    let Z = Array2D.init<float> ni ni (fun i j -> if i=j then 1.0/nObs.[i,0] else 0.0)

This method creates an ni x ni square matrix.  The statement:

fun i j -> if i=j then 1.0/nObs.[i,0] else 0.0

is an anonymous function that says when the row index i equals the column index j, then fill the value with one divided by the ith value of our nObs vector; otherwise fill it with zero.  Cool enough, eh?  That's Canadian, you know?

So, once we have our Z matrix we need to invert it:

    let ZI, ret = matInverse Z

and then recreate xpx = X'Z^-1X:

    let xpx = matMult2 (matMult2 xp ZI) x0

Note: we're just chaining our matrix multiplication function:  first multiplying X' by Z^-1, and then multiplying the result from that operation by X

Now, we need to invert our new xpx

    let X = xpx
    let XI, ret = matInverse X

The result XI is now equal to: (X'Z^-1X)^-1

Finally, we need to calculate bhat_GLS:

    let b = matMult2 (matMult2 (matMult2 XI xp) ZI) y

Again, we have chained our matMult2 function to include ZI.

Recall, bhat_GLS equals: b_gls = (X'Z^-1X)^-1X'Z^-1y. 

So here's our full code:

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

let regressGLS =

    let Z = Array2D.init<float> ni ni (fun i j -> if i=j then 1.0/nObs.[i,0] else 0.0)

    let ZI, ret = matInverse Z

    //--------------------------------
    // Now create X'ZI X
    //--------------------------------

    let xpx = matMult2 (matMult2 xp ZI) x0

    let X = xpx
    let XI, ret = matInverse X

    let k = X.GetUpperBound 0
   
    for i = 0 to k do
         printfn "XI Matrix: %10.5f %10.5f %10.5f %10.5f " XI.[i,0] XI.[i,1] XI.[i,2] XI.[i,3]

    //-----------------------------------------
    //Now Test
    //Multiply XI and X to calculate Identity Matrix
    //-----------------------------------------
    let XTest = Array2D.zeroCreate<float> (k+1) (k+1)
    let XTest = matMult XI X XTest

    for i = 0 to k do
         printfn "XTest Matrix: %10.5f %10.5f %10.5f %10.5f " XTest.[i,0] XTest.[i,1] XTest.[i,2] XTest.[i,3]
   
    //let a0 = Array2D.zeroCreate<float> (k+1) ni
    //let a1 = Array2D.zeroCreate<float> ni 1
    let b = matMult2 (matMult2 (matMult2 XI xp) ZI) y

    for i = 0 to k do
         printfn "Beta Hat: %i %10.5f" i b.[i,0]