###############################################################################
# AEB 6572 - Econometric Methods II                                           #
# Charles B. Moss - cbmoss@ufl.edu                                            #
# R Regression Program Maximum Likelihood - MaxLike03-20-12.R                 #
###############################################################################

library(numDeriv)

###############################################################################
# Data from Moss, Erickson, Ball, and Mishra                                  #
# dta[,1]  - g1 D[z1] Livestock share of revenue times change in Livestock    #
# dta[,2]  - g2 D[z2] Crop share of revenue times change in Crop              #
# dta[,3]  - f1 D[x1] Purchased Input share of cost times change in Purchased #
# dta[,4]  - f2 D[x2] Labor share of cost times change in Labor               #
# dta[,5]  - f3 D[x3] Capital share of cost times change in Capital           #
# dta[,6]  - Change in Livestock Price                                        #
# dta[,7]  - Change in Crop Price                                             #
# dta[,8]  - Change in Purchased Input Price                                  #
# dta[,9]  - Change in Labor Price                                            #
# dta[,10] - Change in Capital Price                                          #
###############################################################################

dta <- read.csv("AggregDiff.csv")

y1 <- dta[,3]
y2 <- dta[,4]
x1 <- dta[,1]
x2 <- dta[,2] 
x3 <- dta[,8]-dta[,10]
x4 <- dta[,9]-dta[,10]
x5 <- dta[,6]
x6 <- dta[,7]

eqD1 <- y1 ~ x1 + x2 + x3 + x4 - 1
res1 <- lm(eqD1)
print(summary.lm(res1))

eqD2 <- y2 ~ x1 + x2 + x3 + x4 - 1
res2 <- lm(eqD2)
print(summary.lm(res2))

xx <- matrix(0,nrow=2*nrow(dta),ncol=7)
yy <- matrix(0,nrow=2*nrow(dta),ncol=1)
for (i in 1:nrow(dta)) {
  xx[2*(i-1)+1,1] <- x1[i]
  xx[2*(i-1)+1,2] <- x2[i]
  xx[2*(i-1)+1,3] <- x3[i]
  xx[2*(i-1)+1,4] <- x4[i]
  xx[2*(i-1)+2,4] <- x3[i]
  xx[2*(i-1)+2,5] <- x1[i]
  xx[2*(i-1)+2,6] <- x2[i]
  xx[2*(i-1)+2,7] <- x4[i]
  yy[2*(i-1)+1]   <- y1[i]
  yy[2*(i-1)+2]   <- y2[i]
}

bols <- solve(t(xx)%*%xx)%*%t(xx)%*%yy
print(bols)
mse  <- t(yy-xx%*%bols)%*%(yy-xx%*%bols)/nrow(yy)
vols <- mse[1,1]*solve(t(xx)%*%xx)

###############################################################################
# Maximum Likelihood Function with gradients and hessian                      #
###############################################################################

ml <- function(b) {
###############################################################################
# Subroutine for computing the gradient of the objective function             #
###############################################################################
  gr <- function(b) {
    gr <- matrix(0,nrow=1,ncol=ncol(xx))
    for (i in 1:nrow(dta)) {
      gr <- gr - err[i,]%*%solve(sigma2)%*%
            rbind(xx[2*(i-1)+1,],xx[2*(i-1)+2,])}
    }
###############################################################################
# Subroutine for computing the Hessian of the objective function              #
###############################################################################
  h  <- function(b) {
    h <- matrix(0,nrow=ncol(xx),ncol=ncol(xx))
    for (i in 1:nrow(dta)) {
       h <- h + t(rbind(xx[2*(i-1)+1,],xx[2*(i-1)+2,]))%*%
            solve(sigma2)%*%rbind(xx[2*(i-1)+1,],xx[2*(i-1)+2,]) }
    }
###############################################################################
# Common values for the objective function, gradient, and Hessian             #
###############################################################################
  errv   <- yy - xx%*%b
  err    <- matrix(0,nrow=nrow(yy),ncol=2)
  for (i in 1:nrow(dta)) err[i,] <- cbind(errv[2*(i-1)+1],errv[2*(i-1)+2])
  sigma2 <- t(err)%*%err/nrow(err)
  res    <- nrow(yy)/2*log(det(sigma2))
  attr(res, "gradient") <- gr(b)
  attr(res, "hessian")  <- h(b)
  return(res)
  }

###############################################################################
# Set initial guess at maximizing values                                      #
###############################################################################

b0 <- matrix(0.01,nrow=nrow(bols),ncol=ncol(bols))

###############################################################################
# Call nonlinear minimizing routine                                           #
###############################################################################

res <- nlm(ml,b0,hessian = TRUE,iterlim = 500)

###############################################################################
# Compare with ordinary least squares values                                  #
###############################################################################

vml <- solve(res$hessian)

print(cbind(bols,diag(vols),bols/diag(vols),res$estimate,diag(vml),res$estimate/diag(vml)))

ml2 <- function(b) {
###############################################################################
# Subroutine for computing the gradient of the objective function             #
###############################################################################
  gr2 <- function(b) {
    gr2 <- matrix(0,nrow=1,ncol=9)
    for (i in 1:nrow(dta)) {
      dx <- rbind(cbind( -x1[i], -x2[i], -x3[i], -x4[i], 0, 0, 0, 0, 0),
                  cbind( 0, 0, 0, -x3[i], -x1[i], -x2[i], -x4[i], 0, 0),
                  cbind(b[8]*x3[i],b[9]*x3[i], 0, 0, b[8]*x4[i], b[9]*x4[i], 0,
                         -x5[i]+b[1]*x3[i]+b[2]*x4[i], x6[i]+b[5]*x3[i]+b[6]*x4[i]))
      gr2 <- gr2 + err[i,]%*%solve(sigma2)%*%dx}
    }
###############################################################################
# Subroutine for computing the Hessian of the objective function              #
###############################################################################
  h2  <- function(b) {
    h2 <- matrix(0,nrow=9,ncol=9)
    for (i in 1:nrow(dta)) {
      dx <- rbind(cbind( -x1[i], -x2[i], -x3[i], -x4[i], 0, 0, 0, 0, 0),
                  cbind( 0, 0, 0, -x3[i], -x1[i], -x2[i], -x4[i], 0, 0),
                  cbind(b[8]*x3[i],b[9]*x3[i], 0, 0, b[8]*x4[i], b[9]*x4[i], 0,
                         -x5[i]+b[1]*x3[i]+b[2]*x4[i], x6[i]+b[5]*x3[i]+b[6]*x4[i]))
      h2 <- h2 + t(dx)%*%solve(sigma2)%*%dx }
    }
###############################################################################
# Common values for the objective function, gradient, and Hessian             #
###############################################################################
  err     <- matrix(0,nrow=nrow(dta),ncol=3)
  err[,1] <- y1 - (b[1]*x1 + b[2]*x2 + b[3]*x3 + b[4]*x4)
  err[,2] <- y2 - (b[5]*x1 + b[6]*x2 + b[4]*x3 + b[7]*x4)
  err[,3] <- x1 - (b[8]*(x5 - b[1]*x3 - b[5]*x4) +
                   b[9]*(x6 - b[2]*x3 - b[6]*x4))
  sigma2  <- t(err)%*%err/nrow(err)
  res     <- nrow(dta)/2*log(det(sigma2))
  attr(res, "gradient") <- gr2(b)
  attr(res, "hessian")  <- h2(b)
  return(res)
  }

btemp <- cbind(t(as.matrix(res$estimate)),0,0)

print(btemp)

gradml2 <- grad(ml2,btemp)
hessml2 <- hessian(ml2,btemp)

lm <- t(gradml2) %*% solve(hessml2) %*% gradml2
print(lm)