This question has two parts, as I do not understand whether my problem is theoretical (identification of the parameters) or practical (insufficient R skills).

  • Econometrics

Most "probit" style models are identified through a normalization of the standard error to one. In my case, I would argue that this is not necessary as the order of magnitude is already set by fixing one coefficient to 1. More specifically, for each observation there is a dummy variable equal to one if the latent index is above a (observation-specific and observed) threshold :

$$d_i=1 \text{ iff } K_i > X_i \beta + \epsilon_i$$

The error term is assumed $\epsilon_i \thicksim \mathcal{N}(0, \sigma)$. In my naive understanding, the likelihood of this problem should somehow look like (where $\Phi$ is the standard normal cdf):

$$L= \prod_{i=1}^N \Phi \left(\frac{K_i - X_i \beta}{\sigma} \right)^{d_i} \Phi \left(\frac{X_i \beta -K_i}{\sigma} \right)^{1-d_i} $$

Is it possible to estimate $\sigma$ without further normalization?

  • Estimation

If the answer to previous part is "yes" -- then why does my R implementation not work?

### simulate data
N <- 2000
b.cons  <- 8 
b.x         <- 10
sig <- 2
x <- cbind(rep(1, N), runif(N)) #"observed variables"
e <- rnorm(N, sig) # "unobserved error"
k <- runif(N)*10+8 # threshold: something random, but high enough to guarantee some variation in i
t <- x%*%c(b.cons, b.x)+e
i <- 1*(k>t) #participation dummy

### likelihood function
probit.sim <- function(params, I, K, X) {
        params[1:2] -> b
        params[3] -> s
        z= (K-X%*%b)/s

    pr.1       = pnorm(z)
    pr.1[pr.1==0] <- 0.001  #seems somehow weird to me, but how is this problem usually treated??
    pr.1[pr.1==1] <- 0.999
    pr.0       = 1-pr.1

    llik = t(I)%*%log(pr.1) + t(1-I)%*%log(pr.0)

### maximize likelihood
optim(c(1,1,1), probit.sim, I = i, K = k, X = x) #using a random starting vector
st <- coef(lm(k*(1-i) ~ x-1)) #searching for better standard values
optim(c(st, 1), probit.sim, I = i, K = k, X = x)

The estimated parameters are clearly not c(8, 10, 2) as they should be.

I asked a related question on already on stack overflow, and the answer was "take better starting values", but this does not seem to do the trick here. Or maybe I don't know how to do it right.

Any ideas?

  • Alternative approach

My alternative was to use standard statistical software and estimate a probit (needs a bit of twisting but should be possible to make it equivalent). This estimates a coefficient for K, which should be equal to $-1/\sigma$; how about taking this $\sigma$ and computing the "non-normalized"/"true" values of the other coefficients?

Many thanks in advance for any suggestion on any of these 3 parts.


1 Answer 1


Right now your error variable e is drawn from a normal with mean 2 and standard deviation 1, so the result you're getting is perfectly expected. Try:

e <- rnorm(N, sd = sig)

This gives me the following estimates:

       x1        x2           
 7.681889 10.507317  1.993629 

Which is in line with what you're expecting.

As a side note, unless you have a good reason, you probably shouldn't be coding your own maximum likelihood optimization. Look up tutorials on using data.frames and glm(..., family = binomial(link = "probit")).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.