# Maximum Likelihood Estimation for Linear Regression Model
#
# ADM 6/5/2006

library(mvtnorm)

regression.loglikelihood <- function(parameters, Y, X) {

   # inputs:
   #
   #   parameters -- beta followed by ln(sigma2)
   #   Y -- (n x 1) matrix Y
   #   X -- (n x k) matrix X (with or without a constant)
   #
   # output:
   #
   #   the evaluated log-likelihood
   
   Y <- as.matrix(Y)
   X <- as.matrix(X) 

   # extract parameters
   k <- ncol(X)
   n <- nrow(X)
   beta <- as.matrix(parameters[1:k])
   lnsigma2 <- as.double(parameters[k+1])
   
   # compute the log-likelihood by evaluating the multivariate
   # Normal density once
   mu <- X %*% beta
   sigma <- exp(lnsigma2) * diag(n)
   output <- dmvnorm(t(Y), mu, sigma, log=TRUE)

   return(output)
}

# Data [BUGS Lines Data]
X <- cbind(1,c(1,2,3,4,5))
Y <- as.matrix(c(1,3,3,3,5))

# Fit Model Using OLS
ols <- lm(Y~X-1)
print(summary(ols))

# Fit Model Using ML
mle <- optim(c(4,3,2),
   regression.loglikelihood,
   hessian=TRUE,
   control=list(fnscale=-1),
   method="BFGS",
   Y=Y, X=X)
cat("MLE Betas\n", mle$par[1:2], "\n")
cat("MLE Sigma\n", sqrt(exp(mle$par[3])), "\n\n")
cat("Hessian\n")
print(mle$hessian)
cat("\nMLE Standard Errors \n",
   sqrt(diag(solve(-1 * mle$hessian))), "\n\n")