# [R] maximum likelihood accuracy - comparison with Stata

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

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

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.