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

library(numDeriv)

###############################################################################
# This program demonstrates the maximum likelihood solution for a single      #
#  equation.                                                                  #
###############################################################################

load("DummyCost.RData")

cost <- z[,1]
x1   <- z[,2]
x2   <- z[,3]
x3   <- z[,4]
w1   <- z[,5]
w2   <- z[,6]
w3   <- z[,7]
y1   <- z[,8]
y2   <- z[,9]

xx <- cbind(matrix(1,nrow=nrow(z),ncol=1),w1,w2,w3,0.5*w1*w1,w1*w2,w1*w3,
            0.5*w2*w2,w2*w3,0.5*w3*w3,y1,y2,0.5*y1*y1,y1*y2,0.5*y2*y2,w1*y1,
            w1*y2,w2*y1,w2*y2,w3*y1,w3*y2)

bols <- solve(t(xx)%*%xx)%*%(t(xx)%*%cost)
mse  <- t(cost - xx%*%bols)%*%(cost - xx%*%bols)/nrow(xx)
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(z)) gr <- gr - err[i]*xx[i,]/sigma2[1,1]
    }
###############################################################################
# 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(z)) h <- h + xx[i,]%*%t(xx[i,])/(2*sigma2[1,1])
    }
###############################################################################
# Common values for the objective function, gradient, and Hessian             #
###############################################################################
  err    <- cost - xx%*%b
  sigma2 <- t(err)%*%err/nrow(err)
  res    <- nrow(z)/2*log(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                                  #
###############################################################################

print(cbind(bols,res$estimate))

numGrad <- grad(ml,res$estimate)
numHess <- hessian(ml,res$estimate)

vml <- solve(res$hessian)

print(cbind(diag(vols),diag(vml)))