What is the benefit of regression with student-t residuals over OLS regression? Sometimes I see advice to fit regressions with student-t residuals rather than using OLS (which is equivalent to assuming normally distributed residuals) if the distribution of the residuals is heavy-tailed. However, since the OLS estimator is BLUE (by Gauss-Markov), it should have lower variance (and therefore MSE) than a regression that assumes student-t residuals fit via maximum likelihood. This is true even if the residuals truly are t-distributed.
In a simple simulation (see bottom for R code), OLS and t-regression are essentially equivalent with respect to out-of-sample MSE and both achieve correct coverage for CIs for the $\beta$ coefficient even though the true residuals follow a student-t distribution. If OLS performs equally for standard tasks such as prediction (as judged by MSE) or inference (as judged by coverage of CIs for model parameters), what are the advantages of fitting a regression assuming student-t errors? Shouldn't I benefit by estimating parameters assuming a correctly specified likelihood rather than a misspecified likelihood? Or is my simulation misleading?
An answer to this post suggests that prediction intervals will be wrong if one fails to use t-distributed errors. Is that really the only benefit?
Here's the simulation code:
library(hett) # for fitting t-based regressions

get_regression_function <- function(beta) {
  reg_function <- function(x){
    beta[1] + beta[2] * x
  }
  return(reg_function)
}
get_t_prediction <- function(df, model) {
  t_beta_hat <- model$loc.fit$coefficients
  design_mat <- df$x
  intercept <- rep(1, length(df$x))
  design_mat <- cbind(intercept, design_mat)                 
  pred <- design_mat %*% as.matrix(t_beta_hat)
  return(pred)
}

get_ols_ci_coverage <- function(ols_fit, par_name, true_theta) {
  CI_mat <- confint(ols_fit)
  lw <- CI_mat[par_name, '2.5 %']
  up <- CI_mat[par_name, '97.5 %']
  cover_beta_1_ols <- (lw <= true_theta) & (true_theta <= up)
  return(cover_beta_1_ols)
}

get_t_ci_coverage <- function(t_fit, par_name, true_theta) {
  t_summary <- as.data.frame(summary(t_fit)$loc.summary$coefficients)
  t_pt_est <- t_summary[par_name, 'Estimate']
  t_se <- t_summary[par_name, 'Std. Error']
  t_ci <- t_pt_est + 1.96 * t_se * c(-1, 1)
  cover_t <- (t_ci[1] <= true_theta) &
    (true_theta <= t_ci[2])
  return(cover_t)
}

## simulation
true_beta_0 = 0.5
true_beta_1 = 1.5
reg_function <- get_regression_function(
  beta = c(true_beta_0, true_beta_1))

ssr_ols <- c()
ssr_t <- c()
ols_ci_coverage <- c()
t_ci_coverage <- c()

num_sims <- 10000
set.seed(1)
for (i in 1:num_sims){
  # Generate a dataset of size N
  N = 10000
  x <- rnorm(n = N, mean = 0, sd = 1)
  errors <- rt(N, df = 3)
  y_mean <- reg_function(x = x)
  y <- y_mean + errors
  df <- data.frame(y = y, x = x)
  df_train <- df[1:(N/2), ]
  df_test <- df[(N/2 + 1):N, ]
  
  # fit linear models
  OLS_fit <- lm(y ~ x, data = df_train)
  t_fit <- tlm(y ~ x,
               data = df_train)
  
  # generate test set predictions
  lm_pred <- predict(OLS_fit, newdata = df_test)
  t_pred <- get_t_prediction(
    df = df_test, model = t_fit)
  
  # test residuals
  lm_res <- lm_pred - df_test$y
  t_res <- t_pred - df_test$y
  
  # MSE -- test set
  ssr_ols[i] <- mean(lm_res^2)
  ssr_t[i] <- mean(t_res^2)
  
  # CI coverage for beta1
  ols_ci_coverage[i] <- get_ols_ci_coverage(
    ols_fit = OLS_fit, par_name = 'x', true_theta = true_beta_1)
  t_ci_coverage[i] <- get_t_ci_coverage(
    t_fit = t_fit, par_name = 'x', true_theta = true_beta_1)
}

mean(ssr_ols)
mean(ssr_t)
mean(ols_ci_coverage)
mean(t_ci_coverage)


Student-t Distribution is the distribution of small sample; as opposed: one of the important assumptions of using OLS - is Normal Distribution, that needs great amount of independent variables in the sample... and in order to cope with outliers - you'd better to increase the size of your sample-data to reach its normality (though, of course, you can just remove your outliers from your data analysis - the other way to cope with robust vars)... to check the data for normality you can use Shapiro-Wilk normality test or Kolmogorov–Smirnov test...
Anyway, you cannot compare OLS-regression & testing residuals - because in order to trust R^2 score of your regression curve (you're getting with OLS instrument) - you have to check your residuals for Normality as well after OLS...
So, OLS is having rather strict requirements about the structure & size of the sample - Normality (of data & residuals as well)... but if you cannot prove the normality you'd better use MLE (maximum likelyhood estimation) instead of OLS - to find out more appropriate distribution to your emperical data... Also speaking about OLS - it also requires the f(x) to be unambiguously defined in apriory analysis of your data dependences, otherwise MLE is advised to find the theoretical curve that mostly suits your emperical distribution...
But to prove that your regression coeficients for each X are meaningful - yes, you see t-Student criteria for each of coef's of regression model found, - additionaly to the F-criteria for the whole regression model found...
import statsmodels.api as sm
...
model= sm.OLS(y,x)
results= model.fit()
print(results.summary())

be carefull with terms: t-Student tests can be used in different stat.investigations with different sense & t-Student Distribution can describe a "small sample"...
