/* Stata */
clear

set obs 1000

gen x=invnorm(uniform())

gen e=invnorm(uniform())

gen y=0.5+0.4*x+e

replace y=0 if y<=0

replace y=1 if y>0



capture program drop myprobit

program myprobit /* to utilize maximum likelihood method to estimate the probit model */

version 10

args lnf theta1

quietly replace `lnf’ = ln(normal(`theta1′)) if $ML_y1==1

quietly replace `lnf’ = ln(normal(-`theta1′)) if $ML_y1==0

end
ml model lf myprobit (y = x) /* to estimate the basic probit model using method lf */

ml check /* to verify that the log-likelihood evaluator you have written seems to work */

ml search /* to find a better initial value so that the convergence of log-likelihood takes less time*/

ml maximize

ml graph /* to graph the log-likelihood values against

the iteration number to verify the process of convergence */
# R

n<-1000;

e<-rnorm(n);

x<-rnorm(n);

a<-0.5

b<-0.4

y<-a+b*x+e;

y<-ifelse(y<=0,0,1)
#Define Likelihood

# Probit LL

llik.probit <- function(par, X, y){

Y <- as.matrix(y)

X <- cbind(1,X)

phi <- pnorm(ifelse(Y == 0, -1, 1) * X%*%par, log.p = TRUE)

sum(phi)

}
# another probit function

#probit.ll <- function(beta,X,y){

#X<-cbind(1,X)

#    phi <- pnorm(X %*% beta)

#    out <- sum(y * log(phi) + (1-y) * log(1-phi))

#}
# Maximize LL

out <- optim(par = c(0,1), llik.probit , y = y, X = x, method=”BFGS”,

control=list(fnscale = -1), hessian = T)
# Calculate standard errors

se <- sqrt(diag(solve(-optim.out$hessian)))

VCV <- solve(-optim.out$hessian)

Zv<-out$par/se

pv<-2*(1-pnorm(Zv))
result<-cbind(out$par,se,Zv,pv)

colnames(result)<-c(”Coef”,”Std.err”,”Zvalue”,”P-value”)

rownames(result)<-c(”Intercept”,”x”)

round(result,4)
# check with canned function; an estimated identical

out <- glm(y ~ x, family=binomial(link=”probit”))

summary(out)

# library(car);lfp<-recode(lfp,” ‘yes’=1;’no’=0 “, as.factor=F)