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)
