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)



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.

The estimator that assumes student-t residuals is not a linear estimator. OLS is the Best among Linear Unbiased Estimators. So, that does not include the student-t residuals estimator. The variance of the estimator that uses student-t residuals can be smaller than the variance of the OLS estimator.
The student-t estimator has the advantage that it is more robust to violations of the assumption that the error is Gaussian distributed. In cases of relatively more higher values it can have a lower variance than the OLS estimator.
An extreme and simple example is the estimation of the location parameter of a shifted t-distribution with 1 degree of freedom (which is the same as the Cauchy distribution). In that case the OLS estimator is the mean of the sample and will have undefined variance. On the other hand, the MLE estimate of the Cauchy distribution although difficult to compute, has lower variance

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?

The linear estimator, a weighted sum of all the $y$, has the distribution of a sum of the distributions of the $y$. Except in extreme cases, this will approach a normal distribution quickly and it wouldn't be so wrong to for prediction intervals.
