Considere el siguiente modelo
$$y_i = \sigma_{c(i)} + \mathbf x_i^\top\beta + u^y_i $$ $$\sigma_{c} = z_c\lambda + \eta_c$$
donde para todos $i$
$$\mathbb E[u^y_i \lvert x_i] = 0$$
Los datos se refieren a una muestra aleatoria $\{y_i,\mathbf x_i,z_{c(i)}\}_{i=1}^N$ dejando $u^y_i, \sigma_c, \eta_c,\lambda$ y $\beta$ no se ha observado.
Intuitivamente, el modelo puede interpretarse como un modelo de dos niveles en el que $i$ es un trabajador observado que recibe un salario $y_i$ en la ciudad $c$ y $z_c$ son covariables observadas a nivel de ciudad mientras que $\eta_c$ son factores no observados específicos de la ciudad que afectan al salario de forma aditiva a través de $\sigma_c$ . La función $c(i)$ simplemente connota la ciudad donde el individuo $i$ funciona.
Es evidente que la estimación de la primera ecuación puede llevarse a cabo utilizando variables ficticias específicas de la ciudad, lo que daría lugar a una estimación $\hat \sigma_c$ para cada ciudad (tengo muchas observaciones para cada ciudad/grupo, así que supongo que esto está bien). Entonces, para estimar $\lambda$ la segunda etapa de regresión
$$\hat \sigma_c = z_c \lambda + \eta_c$$
se realiza. ¿Puede justificarse este enfoque (dar una estimación coherente de $\lambda$ ) cuando
$$\mathbb E[\eta_c \lvert z_c] = 0$$
pero
$$\mathbb E[\eta_c \lvert \mathbf x_i] \not = 0,$$
quizás considerando el DAG del modelo que creo que podría ser algo así 
que debe ser implementado en el siguiente código, que creo que muestra que el enfoque funciona. Pero no estoy seguro de cómo demostrarlo utilizando, por ejemplo, argumentos de la autoría de Pearl en los DAG o cualquier otro argumento dados los supuestos.
library(data.table)
library(lfe)
N <- 100000
C <- 300
# Make index over what cities individual worker are in
city_index <- sample(1:C,N,replace=TRUE)
# Make unobserved city productivity effect eta and observed z
eta <- 6*runif(C)
z <- 2*runif(C)
# Calculate city level effect
a <- 1
c_i <- z[city_index]*a + eta[city_index]
# Simulate worker specific skill x
u_x <- rnorm(N)
x <- u_x + c_i
b <- 2
u_y <- rnorm(N)
# Simulate wages
y <- c_i + x*b + u_y
mydata <- data.table(wage=y,city=city_index,skill=x,city_chr=z[city_index])
model_1 <- felm(wage ~ skill + city_chr,data=mydata)
model_2 <- felm(wage ~ skill - 1|city,data=mydata)
model_1
model_2
city_data <- data.table(getfe(model_2))[,.(idx,effect)]
city_data$city_chr <- z
lm(effect ~ city_chr,data=city_data)
plot(city_data$effect[city_index],c_i)
