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 ______________________________________________ R-help@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.
