F# - Ordinary Least Squares (OLS)

<< F# Matrix Inversion (2) - The Inverse | Home | F# - OLS (R-Square) - Array2D.init >>

F# - Ordinary Least Squares (OLS)

So we're finally ready to run a regression.

Let's first make up some data.

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

#light

open System
open System.Collections.Generic
open System.IO
open Xunit
open MathMod
open Microsoft.FSharp.Math

//---------------------------------------------------------
//Regression
//---------------------------------------------------------

let ni = 10
let k = 4
let x0 = Array2D.zeroCreate<float> ni k
let y = Array2D.zeroCreate<float> ni 1

y.[0,0] <- 290.68
y.[1,0] <- 236.44
y.[2,0] <- 397.78
y.[3,0] <- 209.13
y.[4,0] <- 266.98
y.[5,0] <- 359.83
y.[6,0] <- 56.08
y.[7,0] <- 299.26
y.[8,0] <- 254.52
y.[9,0] <- 496.51

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.[6,0] <- 1.0
x0.[7,0] <- 1.0
x0.[8,0] <- 1.0
x0.[9,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] <- -13.0
x0.[6,1] <- -12.0
x0.[7,1] <- 7.0
x0.[8,1] <- 19.0
x0.[9,1] <- 21.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
x0.[6,2] <- 7.0
x0.[7,2] <- 8.0
x0.[8,2] <- 3.0
x0.[9,2] <- 15.0

x0.[0,3] <- 5.0
x0.[1,3] <- 7.0
x0.[2,3] <- 9.0
x0.[3,3] <- 11.0
x0.[4,3] <- 7.0
x0.[5,3] <- 17.0
x0.[6,3] <- 2.0
x0.[7,3] <- 14.0
x0.[8,3] <- 5.0
x0.[9,3] <- 22.0

I have created a small data set of 10 observations and 4 coefficients.

In scaler notation, these observations represent the following relationship:

Y_i = a + B₁X_i1 + B₂X_i2 + B₃X_i3 + B4Xi4 + e_i

In matrix notation, this equation is represented as:

y = XB + e

where Y is an n x 1 vector, X is n x k matrix, B is a k x 1 vector, and e is an n x 1 vector.

The OLS problem is to solve for (estimate) B-hat using the following equation:

B-hat = (X'X)^-1X'y

This is the standard estimator for the OLS found in every econometrics text book.

Solving this equation is simple given all the F# functions we have written up until now.

First off, I want to refine my transpose and matrix multiplication functions so that I don't have to pass in an empty array. These are below. They just do the work that I was already doing in code, and make life a little bit easier.

//
//NEWLY UPDATED 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

let rec transpose2 (x : float [,] ) =
    let n = x.GetUpperBound 0
    let k = x.GetUpperBound 1
    let xp = Array2D.zeroCreate<float> (k+1) (n+1)
    transpose x xp n k

let rec matMult2 (xp : float [,] ) (x : float [,] ) =
    let n = xp.GetUpperBound 0
    let k = xp.GetUpperBound 1
    let m = x.GetUpperBound 1
    let xpx = Array2D.zeroCreate<float> (n+1) (m+1)
    matMult xp x xpx

Next my main program takes my data array x0 and creates its transpose in xp. I then multiply xp and x0 to create xpx. This result is a square k x k matrix that I want to invert.

let xp = transpose2 x0
let xpx = matMult2 xp x0

I then assign xpx to X for no good reason other than to simplify my subsequent code. I then pass this matrix to my matInverse function, returning XI as the inverse

let X = xpx
let XI, ret = matInverse X

Below is my matInverse function that was derived in my prior post:

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

let matInverse (X : float [,]) =

let ret, L, U = liDecomp X

let k = X.GetUpperBound 0

let zeroX = Array2D.zeroCreate<float> (k+1) (k+1)
let LU = matMult L U zeroX

    //--------------------------------
    //Now Invert L
    //--------------------------------
    let LI, ret_LI = invertL L

    //-----------------------------------------
    //Now Invert U
    //First Transpose U - UP is U Transpose
    //-----------------------------------------
    let Uzero = Array2D.zeroCreate<float> (k+1) (k+1)
    let UP = transpose U Uzero k k

    //-----------------------------------------
    //Now Invert UP
    //Return UPI as Inverse of UP
    //-----------------------------------------
    let UPI, ret_UI = invertL UP

    //-----------------------------------------
    //Now Transpose UPI
    //UI is now U inverse
    //-----------------------------------------
    let UI = transpose UPI Uzero k k

    //-----------------------------------------
    //Multiply UI and LI to calculate XI
    //XI is now X inverse
    //-----------------------------------------
    let XI = Array2D.zeroCreate<float> (k+1) (k+1)
    let XI = matMult UI LI XI

    //-----------------------------------------
    //return the inverse
    //-----------------------------------------
    XI, ret

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

So now calculating B-hat is simply one line of code:

let b = matMult2 (matMult2 XI xp) y

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

The stuff in the parentheses: (matMult2 XI xp) is the product (X'X)^-1X. This result is then multiplied by y, yielding the vector of B-Hat coefficient estimates.

Piece of cake.

Now, let's work on the Standard Errors.

    //---------------------------------------
    // Calc Std Errors
    //---------------------------------------

    // Calc Residual Vector = Y - Yhat
    // First Calc Y Hat = X BetaHat

    let YHat = matMult2 x0 b

    let e = Array2D.zeroCreate<float> ni 1
    for i = 0 to (ni-1) do
        e.[i,0] <- YHat.[i,0] - y.[i,0]

    for i = 0 to (ni-1) do
         printfn "Y, Hat, e: %i %10.5f %10.5f %10.5f " i y.[i,0] YHat.[i,0] e.[i,0]

    //--------------------------------------------------
    //Calc SSE -> Sum of Squared Errors e'e
    //--------------------------------------------------
    let ssV = matMult2 (transpose2 e) e

    let sse = ssV.[0,0]
    printfn "SSE: %10.2f" sse

    let sigma2 = sse / float (ni - (k+1))
    printfn "SE: %10.2f" sigma2

    //--------------------------------------------------
    //Calc Variance/Covariance Matrix of B Hat
    //--------------------------------------------------
    //---------------------------------------
    // Calc T-Stats
    //---------------------------------------

    let bVar = Array2D.zeroCreate<float> (k+1) 1
    let bSE = Array2D.zeroCreate<float> (k+1) 1
    let TStat = Array2D.zeroCreate<float> (k+1) 1

    for i = 0 to k do
        bVar.[i,0] <- sigma2 * XI.[i,i]
        bSE.[i,0] <- Math.Sqrt bVar.[i,0]
        TStat.[i,0] <- b.[i,0] / bSE.[i,0]
        printfn "B, Var, SE, TStat: %10.5f %10.5f %10.5f %10.5f" b.[i,0] bVar.[i,0] bSE.[i,0] TStat.[i,0]
I first calculate the estimated values for Y: YHat. YHat equals is simply the product of X & b: YHat = Xb.

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

The vector of residuals is just YHat - Y. The vector e stores the residuals.

The product e'e is the sum of squared errors. The result of this product is a 1 x 1 matrix. I assign this value to the variable sse.

The variance of the regression (sigma²) is the sum of squared errors divided by n - k.

Now to calculate the standard errors for each of the regression coefficients, we first need to calculate the variance of each estimate:

Var(BHat) = sigma²*(X'X)^-1. That is, the variance of BHat is just the variance of the regression multiplied by (X'X)^-1. The standard error for BHat (bSE) is just the square root of the variance.

And, of course, a coefficient's T-Statistic (TStat) is just the ratio of the coefficient to it's standard error.

See my next post for the calculation of the R² and a printout of the results.

Print | posted @ Thursday, August 06, 2009 9:28 AM