Estimating standard error in a probit: econometrics or programming problem?

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
set.seed(5849)
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 -> 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)
-llik
}

### 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.

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:

\$par
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")).