## Sample Code for MLE in R
# Create Data
rm(list = ls())
a<-rep(1,4)
b<-rep(2,5)
c<-rep(3,6)
d<-rep(4,5)
e<-rep(5,4)
f<-rep(6,1)
q1.data<-c(a,b,c,d,e,f)

# Geometric Log-Likelihood Fuction
geometric.loglikelihood <- function(pi, y, n) {
    loglikelihood <- n*log(pi)-n*log(1-pi)+log(1-pi)*sum(y)
    return(loglikelihood)
}

# optimize the log-likelihood
test <- optim(c(.5),            # starting value for pi
   geometric.loglikelihood,     # the log-likelihood function
   method="BFGS",               # optimization method
   hessian=TRUE,                # return numerical Hessian
   control=list(fnscale=-1),    # maximize function instead of minimize
   y=q1.data, n=25)                   # the data
print(test)

# standard error of estimator
OI<-solve(-1*test$hessian)
se<-sqrt(diag(OI))
se

## Sample Code for Graphs in R
# plot the log-likelihood from MLE example
ruler <- seq(0,1,0.01)
loglikelihood <- geometric.loglikelihood(pi=ruler, y=q1.data, n=25)

#postscript("C:/CLASSES/mle/pi_loglik.ps", horizontal=FALSE, onefile=FALSE, pointsize=12)
plot(ruler, loglikelihood, type="l", lwd=2, lty=1, col="black",
   xlab="pi", ylab="l(pi)", ylim=c(-100,-40))
abline(v=.3205)
#dev.off()

# plot the likelihood
geometric.likelihood<-function(pi, y, n){
    likelihood<-(pi^n)*((1-pi)^(sum(y)-n))
    return(likelihood)
}
likelihood<-geometric.likelihood(pi=ruler, y=q1.data, n=25)

#postscript("C:/CLASSES/mle/pi_lik.ps", horizontal=FALSE, onefile=FALSE, pointsize=12)
plot(ruler, likelihood, type="l", lwd=2, lty=2, col="black",
    xlab="pi", ylab="L(pi)")
abline(v=.3205)
#dev.off()

# plot predicted probabilities of a logit model
rm(list=ls())

month <- seq(0,12,0.01)

month.util <- function(x,z) {
    util <- -.023*x-.249*z-.795*.2496426+.657*.5273876-1.194*.3684639+.278*.3504426+.008*5.230634-.004*.0069339+.004*9.810553-.114*1.253197-.003*29.00538-.689
    return(util)
}

#postscript("C:/Mandate/month.ps", horizontal=FALSE, onefile=FALSE, width=5, height=5, pointsize=12)
plot(month, plogis(month.util(x=month,z=0)), type="l", lwd=2, lty=1, col="black",
   xlab="Month Congress Convenes", ylab="Predicted Probability", xlim=c(0,12),ylim=c(0,1))
lines(month,plogis(month.util(x=month,z=1)),lwd=2,lty=2,col="black")
legend(x=0,y=.8,legend=c('Majority President','Minority President'),lty=1:2, lwd=2)
box()
#dev.off()