################################################################################
# AEB 6572 - Econometric Methods II                                            #
# Charles B. Moss - cbmoss@ufl.edu                                             #
# R Regression Program IllinoisProd.R                                          #
################################################################################
# Example Program for Feasible Generalized Least Squares                       #
################################################################################
# Illinois Corn/Input data from the USDA Chemical Use Survey                   #
# dta[,1] - Nitrogen (lbs per acre)                                            #
# dta[,2] - Phosphorous (lbs per acre)                                         #
# dta[,3] - Potash (lbs per arcre)                                             #
# dta[,4] - Corn Yield (bu. per acre)                                          #
################################################################################

dta <- read.table("illinois.dta")
z   <- as.matrix(dta)

################################################################################
# Setup the quadratic production function                                      #
################################################################################

y   <- z[,4]
x   <- cbind(matrix(1,nrow=nrow(dta),ncol=1),z[,1:3],0.5*z[,1]*z[,1],
           z[,1]*z[,2],z[,1]*z[,3],0.5*z[,2]*z[,2],z[,2]*z[,3],
           0.5*z[,3]*z[,3])

################################################################################
# Estimate ordinary least squares production function                          #
################################################################################

b   <- solve(t(x)%*%x)%*%(t(x)%*%y)
mse <- t(y - x%*%b)%*%(y-x%*%b)/(nrow(x)-ncol(x))
vb  <- mse[1,1]*solve(t(x)%*%x)

res1 <- cbind(b,sqrt(diag(vb)),b/sqrt(diag(vb)),2*(1-pnorm(abs(b/sqrt(diag(vb))))))
print(res1)

################################################################################
# Estimate Residual Squared                                                    #
################################################################################

yh <- x%*%b
errsq <- (y - yh)^2

plot(cbind(z[,2],errsq))

################################################################################
# Regress error squared on inputs                                              #
################################################################################

err1 <- lm( errsq ~ z[,1] + z[,2] + z[,3])
print(summary.lm(err1))

plot(cbind(z[,2],errsq))
pred1 <- err1$fitted.values

################################################################################
# Regress square root of error squared on inputs                               #
################################################################################

errsqrt <- sqrt(errsq)
err2 <- lm(errsqrt ~ z[,1] + z[,2] + z[,3])
print(summary.lm(err2))
plot(cbind(z[,2],errsqrt))
pred2 <- err2$fitted.values

################################################################################
# Use estimated values from the square root of error regression to construct   #
#  Omega hat                                                                   #
################################################################################

p <- diag(pred2)
ohatinv <- solve(p%*%p)

################################################################################
# Apply Feasible Generalized Least Squares                                     #
################################################################################

bgls <- solve(t(x)%*%ohatinv%*%x)%*%(t(x)%*%ohatinv%*%y)
msegls <- t(y-x%*%bgls)%*%(y-x%*%bgls)/(nrow(x)-ncol(x))
vgls <- msegls[1,1]*solve(t(x)%*%ohatinv%*%x)

print(cbind(b,bgls))