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)


Here's one reason:
If you fit the parameters using maximum likelihood, with an assumption of t-distributed errors the fitted line is less impacted by points that are far away from the bulk of the data (in the y-direction).
If you have heavy tailed errors, it's easy to get points that pull the line up (away from the bulk of the data) if you fit by least squares. Using a heavy-tailed error distribution down-weights those points so they don't pull the line about so much.
However, it's important to note that this robustness to heavy tails in the y-values doesn't offer protection against observations with high influence. If that's likely to be an issue, you need something that's robust against influential outliers. [With a designed experiment where you control the predictor-values this may not be an issue, but you don't always have designed experiments.]
(There are other reasons you might want to use a suitable model for errors, like more efficient estimates of parameters and the possibility of more suitable small-sample inference.)

The point about least-squares being BLUE is a bit misleading. Yes, it is best among linear unbiased estimators, but if your conditional distribution is very heavy tailed all linear estimators may be very poor, even arbitrarily so, and having the best of a very poor collection of estimators is no consolation. If you want a reasonable level of efficiency in the presence of heavy tails, linear estimators are not a good choice.
