**DiffusePrioR » R**, and kindly contributed to R-bloggers)

Specification testing is an important part of econometric practice. However, from what I can see, few researchers perform heteroskedasticity tests after estimating probit/logit models. This is not a trivial point. Heteroskedasticity in these models can represent a major violation of the probit/logit specification, both of which assume homoskedastic errors.

Thankfully, tests for heteroskedasticity in these models exist, and it is also possible to estimate modified binary choice models that are robust to heteroskedastic errors. In this blog post I present an example of how to estimate a heteroskedastic probit in R, and also test for heteroskedasticity.

The standard probit model assumes that the error distribution of the latent model has a unit variance. The heteroskedastic probit model relaxes this assumption, and allows the error variance to depend on some of the predictors in the regression model. Those interested in further details of this model, and the potential implications of this form of model misspecification, should consult these notes.

In the code below, I simulate some data, specify the log-likelihood function for the heteroskedastic probit model, estimate this model via maximum likelihood, and then perform a simple LR test of homoskedasticity. Note the log-likelihood function can be simplified from:

to:

where uses the fact that the PDF of the normal distribution is symmetric to cut out some algebra (and possibly lead to an easier log-likelihood function to maximize). The log-liklihood function in the above may prove difficult to maximize, although there is an LM test for homoskedasticity, that only requires the restricted (standard probit model) to be estimated. An interesting application would be to see if Bayesian methods estimated using a Gibbs sampling (via STAN, Bugs, or JAGS) can help deal with computationally difficult log likelihoods in this regard. A topic, I will hopefully turn to in a future blog post.

# hetprob rm(list=ls()) # clear ws library(maxLik) # use max lik package # dgp n <- 1000 # no. obs x1 <- runif(n,-1,1) # predictor 1 x2 <- runif(n,-1,1) # " 2 e1 <- rnorm(n,0,1) # normal error e2 <- (1 + 0.45*(x1+x2))*e1 # hetero error y <- ifelse(0.5 + 0.5*x1 -0.5*x2 - e2 >0, 1, 0) # outcome # estimate normal probit r1 <- glm(y~x1+x2,family=binomial(link="probit")) # hetprob llik hll <- function(beta){ beta1 <- beta[1:dim(X)[2]] beta2 <- beta[(1+dim(X)[2]):((dim(X)[2]+dim(X)[2])-1)] q <- 2*y-1 # sym tranform xbeta <- pnorm(q*((X %*% beta1)/exp(X[,-1] %*% beta2))) sum(log(xbeta)) } X <- cbind(1,x1,x2) ml1 <- maxLik(hll,start=c(0,0,0,0,0)) # maximize sqrt(diag(-solve(ml1$hessian))) # get standard errors # LR test of homosked # 2*(LR_{unrestriced}-LR_{restriced}) LRtest <- 2*(ml1$maximum+1/2*r1$deviance) # LS is chi^2 distributed with 2 dof 1-pchisq(LRtest,df=2) # REJECT (as expected)

**leave a comment**for the author, please follow the link and comment on his blog:

**DiffusePrioR » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...