################################
# R code for MLE example       #
################################

logitdata <- read.table("C:\\MLE.txt", header=T)

indvars <- subset(logitdata,select=c(yrseduc))
constant <- 1
X <- as.matrix(cbind(constant,indvars))
y <- as.matrix(subset(logitdata,select=c(voted2004)))

K <- as.numeric(ncol(X))
varnames <- colnames(X)


## start values ##

startv <- c(0,0)

## logit log-likelihood function ##

logit.lf <- function(b, y, X) {

beta <- b[1:K]

exb <- exp(X%*%beta) 
prob1 <- exb/(1+exb) 

logexb <- log(prob1)

y0 <- 1 - y
logexb0 <- log(1 - prob1)

yt <- t(y)
y0t <- t(y0)

logl <- -sum(yt%*%logexb + y0t%*%logexb0)

return(logl)

}

## logit gradient function ##

logit.gr <- function(b, y, X) {

beta <- b[1:K]
grad <- beta*0


exb <- exp(X%*%beta) 
prob1 <- exb/(1+exb) 

for (k in 1:K) { 
  grad[k] <- sum(X[,k] * (y - prob1))
  }

return(-grad)

}

logitmodel <- optim(startv, logit.lf, gr=logit.gr, method="BFGS", 
control=list(trace=TRUE, REPORT=1), hessian=TRUE, y=y, X=X)
inverted <- solve(logitmodel$hessian)
results <- cbind(logitmodel$par, sqrt(diag(inverted)), logitmodel$par/sqrt(diag(inverted)))
colnames(results) <- c("Coefficient", "Std. Err.", "z")
rownames(results) <- varnames
print(results)