conditional_logit.R

# Simulations to understand the standard error of conditional logit
# In this particular case, my covariates are alternative specific,
# and do NOT vary across respondents
# (in Cameron and Triverdi, even alternative-specific covariates vary
# across respondents)

negloglik <- function(beta) {
  xb <- xx %*% beta
  - sum(y * (xb - log(1 + exp(xb))) + (1 - y) * (-log(1 + exp(xb))))
}

# Calculate the variance by simulation
alpha <- c(0.2)
num_choices <- 500
num_obs <- 100
S <- 500
res <- matrix(NA, nrow = S, ncol = length(alpha))
for (s in 1:S) {
  ww <- mvtnorm::rmvnorm(num_choices, sigma = diag(length(alpha)))
  mat <- matrix(NA, nrow = num_obs, ncol = num_choices)
  for (i in 1:num_choices) {
    mat[, i] <- sum(alpha * ww[i, ]) + evd::rgumbel(num_obs)
  }
  y <- max.col(mat)
  cl_nllik <- function(alpha) {
    wa <- c(ww %*% alpha)
    lse_wa <- log(sum(exp(wa)))
    - sum( wa[y] - lse_wa )
  }
  res[s, ] <- optimize(cl_nllik, lower = -100, upper = 100)$minimum
}
colMeans(res)
apply(res, 2, var)

# Calculate by hand
alpha <- 0.05
num_choices <- 5
num_obs <- 100
x <- mvtnorm::rmvnorm(num_choices, sigma = diag(length(alpha)))
Uij <- matrix(NA, nrow = num_obs, ncol = num_choices)
for (j in 1:num_choices) {
  Uij[, j] <- (alpha * x)[j] + evd::rgumbel(num_obs)  
}
pij <- exp(Uij) / rowSums(exp(Uij))
xij <- matrix(rep(x, num_obs), byrow = T, nrow = num_obs)
xibar <- rowSums(pij * xij)
1 / sum(pij * (xij - xibar) * (xij - xibar))