################################
# R code for probit            #
################################

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

indvars <- subset(probitdata,select=c(age,educ,gender,partyID,limit_imports))
constant <- 1
X <- as.matrix(cbind(constant,indvars))
y <- as.matrix(subset(probitdata,select=c(vote)))

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


## start values ##

startv <- c(0,0,0,0,0,0)

## probit log-likelihood function ##

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

beta <- b[1:K]

cdfn <- pnorm(X%*%beta) 

logcdfn <- log(cdfn)

y0 <- 1 - y
logcdfn0 <- log(1 - cdfn)

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

logl <- -sum(yt%*%logcdfn + y0t%*%logcdfn0)

return(logl)

}

## probit gradient function ##

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

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


cdfn <- pnorm(X%*%beta) 

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

return(-grad)

}

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

####################################
### Hypothetical Individual code ###
####################################

# grab coefficients #

probit_coeffs <- probitmodel$par

# set values of independent variables #
# right now constant=1, all ind vars = 0 #
# modify as necessary #

hypind <- c(1,0,0,0,0,0)

# calculate probability #

hypprob <- pnorm(probit_coeffs%*%hypind)

# display results #

print(hypprob)

#####################################
# simulation to get standard errors #
#####################################

install.packages("MASS")
library(MASS)

# number of draws #

ndraws <- 1000

# grab covariance matrix #

probit_covmat <- inverted

# draw from MVN #

betadraw <- mvrnorm(ndraws, probit_coeffs, probit_covmat)

# multiply to get arguments for CDF #

cdfargs <- betadraw%*%hypind

# calculate probabilities #

hypprobs <- pnorm(cdfargs)

# get mean and standard deviation #

meanprob <- mean(hypprobs)
sdprob <- sd(hypprobs)

# display results #

print(meanprob)
print(sdprob)