library(lfe)
options(lfe.threads = 1, digits = 5, warn = 1)
set.seed(65318)
x <- rnorm(500)
x2 <- rnorm(length(x))

## create individual and firm
id <- factor(sample(10, length(x), replace = TRUE))
firm <- factor(sample(6, length(x), replace = TRUE, prob = c(2, rep(1, 5))))

## effects
id.eff <- rnorm(nlevels(id))
firm.eff <- rnorm(nlevels(firm))

## left hand side
y <- x + 0.25 * x2 + id.eff[id] + firm.eff[firm] + rnorm(length(x))
x[sample(500, 20)] <- NA
y[sample(500, 20)] <- NA
## estimate
summary(est <- felm(y ~ x + x2 | id + firm))
getfe(est)
summary(lm(y ~ x + x2 + id + firm - 1))