[R] maximum likelihood accuracy - comparison with Stata

From: Alex Olssen <alex.olssen_at_gmail.com>
Date: Mon, 28 Mar 2011 17:37:52 +1300

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)
sigma2 <- p$par[4]


Kind regards,

Alex Olssen
Motu Research
New Zealand

R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Mon 28 Mar 2011 - 04:40:23 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 28 Mar 2011 - 15:20:25 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive