In the general linear regression model, one of the assumptions is that the variance is constant across all observation points. Often times this assumption is not valid and variance is not uniform across observations. For example, if we were looking at how much people spend on housing based on their income, it probably is the case that the variance of the errors will be larger for high income folks than for poor folks. In this case, we call the errors "heteroskedastic" as opposed to being "homoskedastic". I just love trying to spell those words on an exam.
So, a more general version of the linear model posits that the variance of the errors is:
So, let's work through an example.
Suppose that our data from my earlier example were not actual individual observations, but grouped data, where each point represents the average values for a group of people. This situation happens a lot in cross section data where the data are grouped into different strata in order to mask the values for any specific person.
We a priori in this case that the variances will not be constant -- the averaging across observations tends to reduce the variance of the grouped observation. As a matter of fact, it reduces it proportionately to the number of observations in each group.
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/ni.
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 bhatGLS:
let b = matMult2 (matMult2 (matMult2 XI xp) ZI) y
Again, we have chained our matMult2 function to include ZI.
Recall, bhatGLS equals: bgls = (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]
-----------------------------------------------------------------------
Next, we need to calcula the Std. Errors.
This is pretty much the same as under OLS, except that the estimator for the regression variance is:
sigma2 = ((y'Z-1y) - (bhat'GLSX'Z-1y) ) / (n - k)
In the code below, I break this big ugly into to parts and then subtract one sum from the other.
-----------------------------------------------------------------------
//---------------------------------------
// Calc Std Errors
//---------------------------------------
//--------------------------------------------------
//Calc SSE -> Sum of Squared Errors e'e
//--------------------------------------------------
let yp = transpose2 y
let S1 = matMult2 (matMult2 yp ZI) y
let S2 = matMult2 (matMult2 (matMult2 (transpose2 b) xp) ZI) y
let sse = S1.[0,0] - S2.[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]
When I run it, this is my output:
Beta Hat: 0 25.85775
Beta Hat: 1 5.70154
Beta Hat: 2 8.54613
Beta Hat: 3 11.30829
SSE: 53125.03
SE: 8854.17
B, Var, SE, TStat: 25.85775 1244.21405 35.27342 0.73307
B, Var, SE, TStat: 5.70154 1.29402 1.13755 5.01212
B, Var, SE, TStat: 8.54613 4.00752 2.00188 4.26906
B, Var, SE, TStat: 11.30829 5.85907 2.42055 4.67178