################################################################################ # AEB 6572 - Econometric Methods II # # Charles B. Moss - cbmoss@ufl.edu # # R Regression Program IllinoisProd.R # ################################################################################ # Example Program for Feasible Generalized Least Squares # ################################################################################ # Declare vec command # ################################################################################ v <- cbind(rbind(11,21,31,41),rbind(12,22,32,42)) vec <- function(z) { dimz <- dim(z) zret <- t(z[1,]) for ( i in 2:dimz[1] ) { zret <- cbind(zret,t(z[i,])) } t(zret) } ################################################################################ # Dutch Consumption Data (Modified) # # dta[,1] - Food consumption (aggregated 1 2 3 4 5 6 7 9) # # dta[,2] - Tobacco # # dta[,3] - Clothing (aggregated 10 11) # # dta[,4] - Other Durables # # dta[,5] - Utilities # # dta[,6] - Other Goods # # dta[,7] - Expenditures # # dta[,8] - Price of Food (aggregated by budget share) # # dta[,9] - Price of Tobacco # # dta[,10] - Price of Clothing (aggregated by budget share) # # dta[,11] - Price of Other Durables # # dta[,12] - Price of Utilities # # dta[,13] - Price of Other Goods # ################################################################################ dta <- read.csv("Dutch06.csv") ################################################################################ # Initial Regressions (dropping Other Goods) # ################################################################################ eqD1 <- dta[,1] ~ dta[,7] + dta[,8] + dta[,9] + dta[,10] + dta[,11] + dta[,12] + dta[,13] - 1 eqD2 <- dta[,2] ~ dta[,7] + dta[,8] + dta[,9] + dta[,10] + dta[,11] + dta[,12] + dta[,13] - 1 eqD3 <- dta[,3] ~ dta[,7] + dta[,8] + dta[,9] + dta[,10] + dta[,11] + dta[,12] + dta[,13] - 1 eqD4 <- dta[,4] ~ dta[,7] + dta[,8] + dta[,9] + dta[,10] + dta[,11] + dta[,12] + dta[,13] - 1 eqD5 <- dta[,5] ~ dta[,7] + dta[,8] + dta[,9] + dta[,10] + dta[,11] + dta[,12] + dta[,13] - 1 eq1.ols <- lm(eqD1) eq2.ols <- lm(eqD2) eq3.ols <- lm(eqD3) eq4.ols <- lm(eqD4) eq5.ols <- lm(eqD5) olsEst <- cbind(eq1.ols$coefficients,eq2.ols$coefficients,eq3.ols$coefficients, eq4.ols$coefficients,eq5.ols$coefficients) y <- as.matrix(dta[,1:5]) x <- as.matrix(dta[,7:13]) ii <- diag(1,nrow=5,ncol=5) xx <- x%x%ii yy <- vec(y) bb <- solve(t(xx)%*%xx)%*%(t(xx)%*%yy) ################################################################################ # Restriction Matrix # # 11 21 31 41 51 12 22 32 42 52 13 23 33 43 53 14 24 # # 34 44 54 15 25 35 45 55 16 26 36 46 56 # rr <- rbind(cbind(0,0,0,0,0, 0, 1, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), cbind(0,0,0,0,0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), cbind(0,0,0,0,0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), cbind(0,0,0,0,0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0), cbind(0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), cbind(0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), cbind(0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0), cbind(0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), cbind(0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0), cbind(0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0)) r <- matrix(0,nrow=10,ncol=1) ################################################################################ # Restricted Least Squares # ################################################################################ br <- bb - solve(t(xx)%*%xx)%*%t(rr)%*%solve(rr%*%solve(t(xx)%*%xx)%*%t(rr))%*% (rr%*%bb - r) ################################################################################ # Compute the residual Vector # ################################################################################ err <- yy - xx%*%br ################################################################################ # Reconfigure Error Vector into a Matrix to Compute the variance # ################################################################################ for (i in 1:nrow(dta)) { bi <- (i-1)*5+1 ei <- i*5 if (i==1) serr <- t(err[bi:ei,1]) else serr <- rbind(serr,err[bi:ei,1]) } ################################################################################ # Estimate Restricted Least Squares with Heteroscedasticity # ################################################################################ vhat <- t(serr)%*%serr/nrow(serr) ivhat <- diag(1,nrow=31,ncol=31)%x%solve(vhat) dd <- t(xx)%*%ivhat%*%xx bsur <- solve(dd)%*%(t(xx)%*%ivhat%*%yy)-solve(dd)%*%t(rr)%*%solve(rr%*% solve(dd)%*%t(rr))%*%(rr%*%(solve(dd)%*%(t(xx)%*%ivhat%*%yy))) print(cbind(bb,vec(olsEst),br,bsur))