################################################################################
# AEB 6572 - Econometric Methods II                                            #
# Charles B. Moss - cbmoss@ufl.edu                                             #
# R System Derivatives                                                         #
################################################################################
# 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[,3:4])          # y is a matrix where the first column is   #
                                   #   log change in capital and the second    #
                                   #   is the log change in labor              #
x <- cbind(matrix(1,nrow=nrow(y),ncol=1),as.matrix(dta[,7:11]))
                                   # x is the independent variables with       #
                                   #   constant                                #

################################################################################
# Two seperate regressions - b1 is the regression for log change in capital and#
#  b2 is the regression for the log change in labor                            #
################################################################################

b1 <- solve(t(x)%*%x)%*%(t(x)%*%y[,1])
b2 <- solve(t(x)%*%x)%*%(t(x)%*%y[,2])

################################################################################
# bind the coefficients into a beta matrix                                     #
################################################################################

bt <- cbind(b1,b2)

################################################################################
# Compute yhat (the estimated values of the y matrix) and the errors           #
################################################################################

yhat <- x %*% bt
err  <- y-yhat

################################################################################
# Compute the sample variance matrix                                           #
################################################################################

vhat <- t(err)%*%err/nrow(err)

################################################################################
# Compute the gradient                                                         #
#   I want to start with a vector of zeros                                     #
################################################################################

glike <- matrix(0,nrow=1,ncol=12)

################################################################################
# Next setup the loop to do the sum                                            #
################################################################################

for (i in 1:nrow(err)) {
  gerr <- cbind(rbind(x[i,],matrix(0,nrow=1,ncol=6)),
                   rbind(matrix(0,nrow=1,ncol=6),x[i,]))
  glike <- glike - err[i,]%*%solve(vhat)%*%gerr
}

print(glike)