What distribution do OLS estimators follow when dependent variable is not normally distributed? I understand that when $Y$ is normally distributed, then OLS yields the same estimators as maximum likelihood, which implies that the estimators are sufficient and will be approximately normally distributed in large samples.
However, what distribution do OLS estimators follow when $Y$ is not normal?
As far as I'm concerned, all the regression result tables I've seen in Python and R return $t$ statistics for each estimated coefficient. Does this mean that the estimators always follow the $t$ distribution? And if so, how many degrees of freedom do the $t$ distributions have?
 A: I am going to assume that you are referring to the conditional distribution of $Y$ in the regression (i.e., given the explanatory variables), which follows directly from the underlying error distribution.  So you are really asking what happens when the underlying error terms are not normally distributed.
The distribution of the OLS estimator is quite robust to non-normality of the error terms in the model, so long as you have a reasonable amount of data, and you have non-pathological behaviour in your explanatory variables.  To see this, note that the OLS estimator can be written in terms of the error terms in the model as:
$$\begin{align}
\hat{\boldsymbol{\beta}} 
&= (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \mathbf{Y} \\[6pt]
&= (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} (\mathbf{x} \boldsymbol{\beta} + \boldsymbol{\varepsilon}) \\[6pt]
&= \boldsymbol{\beta} + (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \boldsymbol{\varepsilon} \\[6pt]
&= \boldsymbol{\beta} + \sum_{i=1}^n \varepsilon_i \mathbf{w}_i, \\[6pt]
\end{align}$$
where the vectors $\mathbf{w}_i = [(\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T}]_{\cdot, i}$ are weight vectors that are fully determined by $\mathbf{x}$.  Observe that the deviation of the OLS estimator from the true coefficient vector is a linear function of the error terms.  Now, suppose that the error terms are independent with some distributon that has zero mean and a finite variance  $\sigma^2 < \infty$, but which is not a normal distribution.  Under broad conditions, we can appeal to the multivariate version of the Lyaponov central limit theorem (CLT) to establish that when $n$ is large we have:
$$\sum_{i=1}^n \mathbf{w}_i \varepsilon_i \overset{\text{Approx}}{\sim} \text{N} \Bigg( 0, 
\sigma^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \Bigg).$$
Consequently, for large $n$ you have:
$$\hat{\boldsymbol{\beta}} \overset{\text{Approx}}{\sim} \text{N} \Bigg( \boldsymbol{\beta}, 
\sigma^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \Bigg).$$
Now, the specific conditions required to apply the CLT here are a bit complicated.  Roughly speaking, you need to show that the Lyapunov condition on the weighted sum is satisfied, which requires limiting conditions on the explanatory variables (see e.g., the Grenander conditions discussed here).  However, under non-pathological behaviour for the explanatory variables, and assuming that the error terms are IID with finite variance, this is usually sufficient to allow application of the CLT, which means that the OLS estimator is approximately normally distributed when $n$ is large.  Note that this result applies even if the underlying error distribution is not normal.
By the way, this is one of the big reasons why most of the standard tests in regression analysis are robust to loss of the normality assumption.  All of the coefficient tests and goodness-of-fit tests can be derived using the CLT approximation under broad conditions that do not require the error terms to be normally distributed.  The normality assumption for the error terms is important for prediction purposes, and you can et very bad predictions of new response variables if you apply this assumption without proper scrutiny.  However, so long as you have a reasonable amount of data to fit your model, the normality assumption is not usually important for the internal T-tests and F-tests, and related distributional results for the coefficient estimators and goodness-of-fit statistics.
