probit ml in stata as well as R

peacelovekenzie — Tue, 16 Sep 2008 19:33:41 -0400

/* 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)

Chang ,JX

probit ml in stata as well as R