Are the estimated linear regression coefficients still normally distributed if one of the predictors is a factor variable? The linear regression model in matrix form is
$$
\mathbf{y} = \mathbf{X}\mathbf{\beta} + \mathbf{\varepsilon},
$$
where there errors have finite variance $\sigma^2$. The least squares estimated solution is
$$
\hat \beta = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}.
$$
If the errors are normally distributed then
$$
\hat \beta|X \sim \mathcal{N}(\beta,\sigma^2 (\mathbf{X}^T\mathbf{X})^{-1}).
$$
Now suppose one of the predictors is a factor variable, for example, gender. Are the estimated regression coefficients $\hat \beta$ still normally distributed in this case?
 A: I like to address this sort of questions with simulations, which of course don't prove anything but at least they give a sense of how things may be.
So, let's regress $y$ on the binary variable $x$ N times and see if the resulting coefficients are normally distributed:
n <- 10
N <- 5000
b <- rep(NA, N)
for(i in 1:N) {
    set.seed(i)
    y <- c(rnorm(n= n, mean= 10), rnorm(n= n, mean= 12))
    x <- rep(c('M', 'F'), each= n)
    b[i] <- lm(y ~ x)$coefficients[2]
}
hist(b)

Yes, they look normal:

A: The only assumptions you’ve made about $\mathbf{X}$ is that it contains numbers and has full rank. If some of those numbers form a column of all $0$s and $1$s for a factor variable, so be it.
A: Yes, one can use dummy variables, but there is a possible concern on the validity of test results, based on small underlying sample size and/or large variance situations, per some sources.
For example, per Research Gate, to quote a comment by Hensen:

My econometrics book "introduction to econometrics" writes in regard to hypothesis tests and confidence intervals that
"The theoretical underpinnings of this procedure are that the OLS estimator has a large-sample normal distribution that, under the null hypothesis, has as its mean the hypothesized true value and that the variance of this distribution can be estimated consistently ".
What is meant by that the OLS estimator has a large-sample normal distribution? and especially when the variable is a dummy?

Also, same thread comment by Kelvyn Jones, to quote:

When I operate in Bayesian estimation mode, I find that normal priors and hence posteriors are fine for regression fixed effects estimate (that is the estimation of a mean) under 1 above (hence symmetric ci’s) and hence ok for a beta associated with a dummy, but variances need a gamma distribution which can differ from a normal by being positively skewed, no values below zero, and hence asymmetric ci’s. Unless the variance is large, and far away from zero, when the normal approximation is fine. That is I am making explicit distributional assumptions about the parameters as well as the residuals. Indeed I make assumptions about the variance of the residuals (the distribution of the parameter) and not just about the individual residuals.

[EDIT] So bottom line with dummy variables, unless one is dealing with small numbers (say under 10) for cell counts and/or do not have apparent volatility issues, statistical test results are still likely accurate.
