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:
Yi = a + B1Xi1 + B2Xi2 + B3Xi3 + B4Xi4 + ei
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 (sigma2) 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) = sigma2*(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 R2 and a printout of the results.