[R] maximum likelihood accuracy - comparison with Stata

[Alex Olssen]

Hi everyone,

I am looking to do some manual maximum likelihood estimation in R.  I
have done a lot of work in Stata and so I have been using output
comparisons to get a handle on what is happening.

I estimated a simple linear model in R with   lm()   and also my own
maximum likelihood program.  I then compared the output with Stata.
Two things jumped out at me.

Firstly, in Stata my coefficient estimates are identical in both the
OLS estimation by   -reg-  and the maximum likelihood estimation using
Stata's   ml lf  command.
In R my coefficient estimates differ slightly between   lm()   and my
own maximum likelihood estimation.

Secondly, the estimates for   sigma2   are very different between R
and Stata.  3.14 in R compared to 1.78 in Stata.

I have copied my maximum likelihood program below.  It is copied from
a great intro to MLE in R by Macro Steenbergen
http://artsci.wustl.edu/~jmonogan/computing/r/MLE_in_R.pdf

Any comments are welcome.  In particular I would like to know why the
estimate of   sigma2   is so different.  I would also like to know
about the accuracy of the coefficient estimates.

## ols
ols <- lm(Kmenta$consump ~ Kmenta$price + Kmenta$income)
coef(summary(ols))

## mle
y <- matrix(Kmenta$consump)
x <- cbind(1, Kmenta$price, Kmenta$income)
ols.lf <- function(theta, y, x) {
  N <- nrow(y)
  K <- ncol(x)
  beta <- theta[1:K]
  sigma2 <- theta[K+1]
  e <- y - x%*%beta
  logl <- -0.5*N*log(2*pi)-0.5*N*log(sigma2)-((t(e)%*%e)/(2*sigma2))
  return(-logl)
}
p <- optim(c(0,0,0,2), ols.lf, method="BFGS", hessian=T, y=y, x=x)
OI <- solve(p$hessian)
se <- sqrt(diag(OI))
se <- se[1:3]
beta <- p$par[1:3]
results <- cbind(beta, se)
results
sigma2 <- p$par[4]
sigma2

Kind regards,

Alex Olssen
Motu Research
Wellington
New Zealand