F# - GLS - In Progress

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:

s2 * Z; where Z is n x n.  The product is then assumed to be distributed normally with a mean of 0.0 and a constant variance.

In this case the coefficient estimates for beta hat are:

bgls = (X'Z-1X)-1X'Z-1y.  Again Z-1 is just the inverse of the Z matrix above.

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.

That is Var(ebari) = s2 / ni , where ni is the count of observations in the group.

In this case, then, Z looks like:

[ 1/n1    0    0    0 ]
|    0 1/n2    0    0 |
|    0    0 1/n3    0 |
[    0    0    0 1/n4 ]

in a 4 x 4 case.  The ni and the counts of observations in each group

What's the inverse of Z?  In this case it is simple:

[ n1  0  0    0 ]
|    0 n2  0  0 |
|    0  0 n3  0 |
[    0  0  0 n4 ]

But, in the more general case we would just go and calculate the inverse.  This is what I have done in the F# code below.

First, let's make up an array that contains the count 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



Print | posted @ Monday, August 10, 2009 3:40 PM

Comments on this entry:

Gravatar # re: F# - GLS - In Progress
by Web Design London at 3/10/2010 3:31 AM

Useful information shared..Iam very happy to read this article..thanks for giving us nice info.Fantastic walk-through. I appreciate this post.
  
Gravatar # re: F# - GLS - In Progress
by Electronics at 3/10/2010 3:43 AM

Excellent post.I want to thank you for this informative read, I really appreciate sharing this great post. Keep up your work
  
Gravatar # re: F# - GLS - In Progress
by Remote Infrastructure Management at 3/10/2010 4:57 AM

Hi webmaster, commenters and everybody else !!! The blog was absolutely fantastic! Lots of great information and inspiration, both of which we all need!b Keep 'em coming... you all do such a great job at such Concepts... can't tell you how much I, for one appreciate all you do!
  
Gravatar # re: F# - GLS - In Progress
by IT Services Melbourne at 3/10/2010 11:52 PM

This is a really good read for me, Must admit that you are one of the best bloggers I ever saw.Thanks for posting this informative article.
  
Gravatar # re: F# - GLS - In Progress
by Mobile Phone Reviews at 3/11/2010 3:31 AM

This is Mobile Phone Reviews Blog. This is blog about new, used and upcoming Mobile Phones of such world famous brands as Sharp, Samsung, LG, Apple, Nokia and etc. Our blog help you choose right mobile phone for you.
  
Gravatar # re: F# - GLS - In Progress
by Cs Hacks at 3/12/2010 4:43 AM

Hi..it was  a gripping post hi......
  
Gravatar # re: F# - GLS - In Progress
by Email Marketing Software at 4/2/2010 3:25 AM

I have been visiting related blogs and sites lately and i have to admit you have a nice design and content. I have bookmarked your page and hope to mention your post in my potential blog.Good post, but have you thought about F# - GLS - In Progress
  
Gravatar # re: F# - GLS - In Progress
by chatroulette application at 4/16/2010 7:49 AM

This is a very nice article that gives in depth information on the subject. I am looking forward to new ones! Thanks in advance
  
Gravatar # re: F# - GLS - In Progress
by flex programming at 4/23/2010 8:54 AM

There are much more things to think over but this particular one makes me puzzled completely.
  
Gravatar # re: F# - GLS - In Progress
by Kids Games at 4/24/2010 7:04 AM

This site is very informative to those Engineering and Math students. Keep posting informative articles for them. Good work!
  
Gravatar # re: F# - GLS - In Progress
by Ben10 games at 4/27/2010 7:09 AM

Thank you for the samples, very nice site.
  
Gravatar # Your content is really useful
by dress at 7/15/2010 2:24 AM

sale
evening dresses
Prom dresses
wedding dresses
on best wedding dresses for 2009 and 2010. You can find latest collection of woman's dresses and casual dresses on this site
discount Prom dresses
discount wedding dresses
  
Gravatar # re: F# - GLS - In Progress
by christian shoes at 7/16/2010 12:57 AM

Where to buy Cheap and super-fashion [url=http://www.martchristianlouboutin.com/" title="http://www.martchristianlouboutin.com/">http://www.martchristianlouboutin.com/]christian louboutin shoes[/url] ,buy save 70% christian dior,[url=http://www.martchristianlouboutin.com/" title="http://www.martchristianlouboutin.com/">http://www.martchristianlouboutin.com/]christian louboutin pumps[/url] , Free Shipping,7 Days To worldwide!
  
Gravatar # re: F# - GLS - In Progress
by deexu at 7/16/2010 1:45 AM

great range of Ed Hardy products. Ed Hardy Women's Ellerise Lowrise Sneaker · Ed Hardy Women's
ed hardy jeans, ed hardy hoody, ed hardy shirt, ed hardy clothing, ed hardy cap, ed glasses, ed belts,
women fashion shoes, men's clothes. helping .perhaps you will like
Ed Hardy
Ed Hardy shoes
Ed Hardy shirts
Ed Hardy clothes
Ed Hardy clothing
Ed Hardy shoes
Don Ed Hardy is an American tattoo collector raised in Southern California
Ed Hardy Clothing,Christian Audigier,Ed Hardy Shoes,Ed Hardy Swimwear,Ed Hardy Hat,
ED Hardy Caps
Ed Hardy Sunglasses
Ed Hardy Wallets
EdHardy
Gucci outlet store online, numerous cheap Gucci bags, handbags, wallets, purses, totes, shoes on sale,
cheap prices and authentic qualities
gucci handbags
gucci jewelryFDGJTOIGHJJ
  
Gravatar # re: F# - GLS - In Progress
by Marian34GOODMAN at 7/19/2010 2:51 PM

I would like to propose not to hold off until you earn enough money to order different goods! You should just get the business loans or credit loan and feel fine
  
Gravatar # re: F# - GLS - In Progress
by chat software at 7/20/2010 2:12 AM

Well, the info your share here is great and informative to me as I am very new to the subject. But I love reading and getting some more knowledge on it. Thanks
  
Gravatar # shi
by fiwedding at 7/23/2010 3:16 AM

supply in stock and custom lace front wigs, full lace wigs, lace wigs, human hair wigs,
remy lace front wigs, cheap wigs, cheap, buy, celebrity
full lace wigs
lace wigs
lace wigs sale
lace front wigs
synthetic front lace wigs
A Famous Dresses Shop which sell directly Wedding Dresses, Evening Dress, Bridesmaid Dresses,Prom dresses
cheap wedding dresses
cheap evening dresses
cheap prom dresses
cheap evening dresses
cheap prom dresses
Elegant evening dresses are always associated with the brides and their bridesmaids.Shopping for evening dresses and your wedding dress in stylish bridal ..
  
Gravatar # http://www.efox-shop.com
by china handy at 8/6/2010 6:57 PM

efox-shop the best place to buy dual SIM dual standby phone. The efox-shop service is good, and the full range, such as chinesische handy kaufen china handy kaufen Großhandel Handy Grosshandel Handy Großhandel Handys chinesische handy TV Handy Chinesische Handys welcome to purchase http://www.efox-shop.com" title="http://www.efox-shop.com">http://www.efox-shop.com http://www.efox-shop.com" title="http://www.efox-shop.com">http://www.efox-shop.com"target=blank>chinesische handy kaufen china handy tv handy Chinesische Handys
  
Gravatar # re: F# - GLS - In Progress
by pay essays at 8/7/2010 2:21 PM

Don't buy useless stuff and pay your cash for custom research paper service, which would aid you to build your career.
  
Gravatar # re: F# - GLS - In Progress
by essay for sale at 8/7/2010 2:21 PM

Not every university student is perfect. Furthermore, I guess that a masters essay service would assist all the university students with an essay research paper completing.
  
Gravatar # re: F# - GLS - In Progress
by buy essays at 8/9/2010 3:33 AM

Time of party is a good thing but not every student can use that. Generally college students sit at home and research their assignments having no free time for their own affairs. However, it is over at present time because you have a chance to buy essay and enjoy your free time.
  
Gravatar # re: F# - GLS - In Progress
by essay papers at 8/10/2010 8:21 PM

Good post.I appreciate perceive it. Any female all over the world is willing to stay unique, but doesn’t get know the way to do this. But billions of different people serch for writing services.
  
Gravatar # re: F# - GLS - In Progress
by Vibram Five Fingers at 8/13/2010 12:23 AM

I recently came across your blog and have been reading along.
I thought I would leave my first comment. I don’t know what to say except that I have enjoyed reading.Nice blog,I will keep visiting this blog very often.
Christian Louboutin Sale
Cheap Christian Louboutin
Buy Christian Louboutin
Discount Vibram Five Fingers
Cheap Vibram Five Fingers
Vibram Five Fingers Sale

  
Gravatar # re: F# - GLS - In Progress
by ed hardy at 8/16/2010 3:50 AM

Find the right [url=http://www.been2u.com/" title="http://www.been2u.com/">http://www.been2u.com/][b]ed hardy clothing[/b][/url] for yourself, just visit Been.com. There are varieties of ed hardy for your reference. Some of them are real great hot sell in our store such as [url=http://www.been2u.com/" title="http://www.been2u.com/">http://www.been2u.com/men-ed-hardy-jeans][b]ed hardy jeans[/b][/url]. You can enjoy free and fast shipping if you buy [url=http://www.been2u.com/" title="http://www.been2u.com/">http://www.been2u.com/][b]ed hardy[/b][/url] from our store! Good luck and welcome you!
  
Gravatar # re: F# - GLS - In Progress
by paper writing service at 8/16/2010 5:15 PM

When you are confused and don’t know the correct way to perform the research paper writing, you would have an opportunity to buy an essay at the modern essay writing service. That would really save time and money.
  
Gravatar # re:Blu-ray Copy
by Blu-ray Copy at 8/19/2010 8:55 PM

With [url=http://www.bluraycopy.biz]Blu-ray Copy[/url], you canMake and Burn High-Definition Blu-ray Movies with super fast speed. Windows 7 Supported.

More info you can visit: [url=http://bluraycopy.biz" title="http://bluraycopy.biz">http://bluraycopy.biz]http://bluraycopy.biz" title="http://bluraycopy.biz">http://bluraycopy.biz[/url]

More Related Products : * [url=http://bluraycopy.biz/blu-ray-burner.html]blu-ray burning[/url] * [url=www.bluraycopyformac.com/dvd-copy-for-mac.html]blu-ray copy for mac[/url] * [url=bluraycopy.biz/bdburner/avi-to-blu-ray-burner.html]avi to blu-ray[/url]* [url=bluraycopy.biz/.../mpeg-to-blu-ray-burner.html] mpeg to blu-ray[/url] * [url=bluraycopy.biz/bdburner/wmv-to-blu-ray-burner.html] wmv to blu-ray[/url] * [url=bluraycopy.biz/bdburner/mkv-to-blu-ray-burner.html]mkv to blu-ray[/url] * [url=bluraycopy.biz/.../m2ts-to-blu-ray-burner.html]m2ts to blu-ray[/url] *
  
Gravatar # re: F# - GLS - In Progress
by A level coursework at 8/21/2010 12:01 PM

Do not be afraid to buy course work writing service! Doing that you will probably have good results.
  
Gravatar # re: F# - GLS - In Progress
by Sports Betting at 8/26/2010 11:25 AM

Its a great posting.I really like it.

Sports Betting
Online Casino
Web Hosting
Credit Card
Insurance
Biofuel
  
Gravatar # re: F# - GLS - In Progress
by replicas bag at 8/28/2010 7:30 AM

The post is written in very a good manner and it contains many useful information for me.
You have a very impressive writing style.
  
Gravatar # re: F# - GLS - In Progress
by Business Telephone Systems at 9/1/2010 9:59 AM

Great post. How did you come across this yourself?
  
Gravatar # Chanel Handbags
by Chanel Handbags at 9/1/2010 6:33 PM

I will always attempt to find fault with Cameron Diaz because she dates my first no-chance-in-hell-he-will-date-me boyfriend, Justin Timberlake. Chanel Handbags 0202 Since then, I have seriously wondered why I even found him so sexy (ok so he dances from his baby toe to his left ear perfectly) and ever thought I could just maybe meet him someday. Cameron and him make a perfect commitment phobic couple. Cam does not want kids until she gets enough sleep and Justin does not want to settle down. Lately Cameron has been seen with her dark tresses, a seriously sexy pair of heels, skinny jeans, and a hot Chanel Handbags 0145 . Cameron Diaz has been carrying the Balenciaga Work with Giant Hardware Spring/Summer 2007. The giant hardware is a hit or miss for many Monogram Suede Embossed Whisper Gm Black lovers, but there is no doubt this outfit is gorgeous on Cam. Ah well, I let JT slip by me to an older woman with a gorgeous new Chanel Handbags 0175. Call Balenciaga NY for further information,212 206 0872 or contact your local Balenciaga boutique.
  
Gravatar # timberland shoe company
by timberland for you at 9/1/2010 8:19 PM

On a certain cheap timberland boots day at a certain hour, we will pull into the station. Bands will be playing and flags discount timberland boots waving. Once we get there, so many wonderful timberland winter boots dreams will come true and the pieces of our lives will fit together like a completed jigsaw puzzle. How restlessly we pace the aisles, womens timberland boot the minutes for loitering --waiting, waiting, waiting for the station.But uppermost in our timberland shoes store minds is the final destination. On a certain day at a certain hour, we will pull into the station. Bands will be playing and timberland eye boat flags waving. Tucked away in our timberland for you subconscious is an idyllic vision. We see ourselves on a long trip that timberland 6 inch spans the continent. We are traveling by train. Out timberland hiking boots windows, we drink in the passing scene of cars on nearby highways, of children timberland shoe company waving at a crossing, of cattle grazing on a distant timberland boots hillside, of smoke pouring from a power plant, of row upon row of corn and wheat, of flatlands and timberland wheat shoes valleys, of mountains and rolling classic 3 eye timberland boat hillsides, of city skylines and village halls.But uppermost in our black timberland boots minds is the final destination.
  
Gravatar # lace wigs
by lace wigs at 9/2/2010 12:46 AM

what a good and colourfull website. lace wigs
  
Gravatar # re: F# - GLS - In Progress
by nannan at 9/3/2010 12:30 AM

NIUUYH HYJTH Tiffany jewellery
Tiffany
Tiffany & Co
Tiffany Co Bracelets
I came across a story in the WSJ last week (WSJ August 13th) in which some major food companies were warning of a sugar shortage in the US. Excuse me. A sugar shortage? I'm not making this up. Here's the text of the letter
Tiffany Co Charms
Tiffany Co Earrings
Tiffany sale
Tiffany Co Necklaces
Tiffany uk
Tiffany Co Rings
tiffany jewelry
tiffany co
  
Gravatar # re:Acceptance Letter - Sample Letters of Acceptance
by Acceptance Letter - Sample Lette at 9/4/2010 1:40 AM

Thanks to nice post please update latest information....
  
Gravatar # www.favre-jersey.com
by <p>www.favre-jersey.com is a ul at 9/5/2010 10:43 PM

www.favre-jersey.com is a ultimate Minnesota Vikings Brett Favre Jerseys Shop and offer professional Vikings Brett Favre Jersey, specialized in supplying Authentic vikings jersey and Minnesota Vikings jersey Brett Favre jerseys are on sale in different styles and colors with top quality at cheaper prices.In Minnesota Vikings 50 anniversary, we launched the 50th Anniversary purple Jersey . "Love Sports,Love Brett Favre,Love Shopping!" Welcome your visit ! We have earned good reputation in foreign markets, more than 90% customers are satisfied with our products and service, till now our online members are beyond 80,000. In 2010, we will strive to improve the quality of our products and service. We are worth your trust! Buy Authentic favre jerseys and Favre vikings jersey at favre-jersey.com, supporting Minnesota Vikings, supporting Brett Favre,supporting our pro bowl jerseys ,whose success needs your attention and concern, hurry up your action now!

  
Gravatar # re: F# - GLS - In Progress
by baby shower gifts at 9/7/2010 3:28 AM

I did't like this article at all. I think I could dissagree with the main ideas. I won't share it with my friends.. You should think of other ways to express your ideas.
baby shower gifts
  
Gravatar # re: billig tandblekning
by billig tandblekning at 9/7/2010 4:11 AM

Great post,The post is written in very a good manner and it contains many useful information for me.
You have a very impressive writing style.
  
Gravatar # re: F# - GLS - In Progress
by nfl jerseys wholesale at 9/8/2010 7:18 PM


If you like football, you may like the Cheap Jerseys .
We produce the high quality Replica Jerseys .
So you can Wholesale MLB Jerseys from us.
  

Your comment:

Title:
Name:
Email:
Website:
 
Italic Underline Blockquote Hyperlink
 
 
Please add 3 and 7 and type the answer here: