################################################################################
# AEB 6572 - Econometric Methods II                                            #
# Charles B. Moss - cbmoss@ufl.edu                                             #
# R Regression Program Prog02-01-12.R                                          #
################################################################################
# Loading the nlme - Nonlinear Mixed - Effects Models Package                  #
#   We want to use the fdHess() Command                                        #
################################################################################

library(nlme);

################################################################################
# Data Set Jorgensen's KLEM dta for Food and Kindred Products                  #
# dta[,1]  - Year                                                              #
# dta[,2]  - log change in output                                              #
# dta[,3]  - log change in quantity of capital                                 #
# dta[,4]  - log change in quantity of labor                                   #
# dta[,5]  - log change in quantity of energy                                  #
# dta[,6]  - log change in quantity of materials                               #
# dta[,7]  - log change in output price                                        #
# dta[,8]  - log change in capital price                                       #
# dta[,9]  - log change in labor price                                         #
# dta[,10] - log change in energy price                                        #
# dta[,11] - log change in materials price                                     #
################################################################################

dta <- read.csv("KLEM-Food01-25-12.csv")

################################################################################
# Set up a basic regression log change in output is a function of log change   #
#  in output price, and the log change in input prices.                        #
################################################################################

y <- as.matrix(dta[,2])            # y is the dependent variable log change in #
                                   #   output                                  #
x <- cbind(matrix(1,nrow=nrow(y),ncol=1),as.matrix(dta[,7:11]))
                                   # x is the independent variables with       #
                                   #   constant                                #

b <- solve(t(x)%*%x)%*%(t(x)%*%y)  # Ordinary Least Squares Estimator          #

print(b)

err <- y - x %*% b                 # Compute the errors                        #
sse <- t(err)%*% err               # Sum of Square Errors                      #
mse <- sse /(nrow(x) - ncol(x))    # Mean Squared Errors                       #
vb  <- mse[1,1]*solve(t(x)%*%x)    # Variance matrix for betas                 #

print(mse)
print(sqrt(mse))

qt <- diag(1,nrow(x),nrow(x)) - 
            (matrix(1,nrow(x),1)%*%matrix(1,1,nrow(x)))/nrow(x)

ssr <- t(b)%*%t(x)%*%qt%*%y       # Sum Squared Regression                     #
sst <- t(y)%*%qt%*%y              # Sum Squared Total                          #
r2 <- ssr/sst                     # r^2 (multiple correlation coefficient)     #
print(r2)
ff <- (ssr/(ncol(x)-1))/(sse/(nrow(x)-ncol(x)))
                                  # Compute the F-test for the model parameters#
print(ff)

################################################################################
# Use the "Canned" Regression command in R [lm() - Linear Model]               #
################################################################################

b1 <- lm(y ~ x[,2] + x[,3] + x[,4] + x[,5] + x[,6] )
print(summary.lm(b1))

################################################################################
# Print our computed variance of beta matrix for comparison                    #
################################################################################

print(sqrt(diag(vb)))

################################################################################
# Estimate the linear model by Maximum Likelihood assuming normality           #
#  (1) Define the likelihood function (In this case the concentrated likelihood#
#      function).                                                              #
################################################################################

ml <- function(b) {
    sighat <- t(y-b[1]-b[2]*x[,2]-b[3]*x[,3]-b[4]*x[,4]-b[5]*x[,5]-b[6]*x[,6])%*%
               (y-b[1]-b[2]*x[,2]-b[3]*x[,3]-b[4]*x[,4]-b[5]*x[,5]-b[6]*x[,6])/nrow(y)
   nrow(y)/2*log(sighat) }

################################################################################
#  (2) Call the function which generates the numerical optimal (by gradient    #
#      search).                                                                #
################################################################################

res <- optim(matrix(0,6,1),ml)

################################################################################
#  (3) Compute the Hessian matrix as the Cramer-Rao lower bound estimate of    #
#      of the variance matrix.                                                 #
################################################################################

resh <- fdHess(res$par,ml)
vest <- solve(resh$Hessian)

print(sqrt(diag(vest)))

print(sqrt(diag(vest)*nrow(y)/(nrow(y)-1)))