Why do my bootstrapped coefficient standard errors differ from statsmodels standard errors

Question

In stats models, I run a regression and get the following SE from the regression output:

SE(const) = 0.028

SE(AveBedrms) = .017

SE(HouseAge) = .0001

From the sample dataset, I randomly sample with replacement and run a regression 50000 times. I then compute the standard deviation of the coefficients from these regressions and get:

std const coeffs = 0.049039

std AveBedrms coeffs = 0.038270

std HouseAge coeffs = 0.000665

These are almost double what statsmodels computes. Am I doing something wrong?

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from sklearn.datasets import fetch_california_housing

# Loading  dataset
cali = fetch_california_housing(as_frame=True)
features = cali.feature_names
sub_features = ["AveBedrms", "HouseAge"] # features I will use to keep simple
df = cali.frame


X = sm.add_constant(df[sub_features])
model = sm.OLS(df.iloc[:, -1], X)
results = model.fit()
print(results.summary())

# SE(const) = 0.028
# SE(AveBedrms) = .017
# SE(HouseAge) = .0001


bootstrap_values = []
N = 50_000
for i in range(N):
  df_sampled = df.sample(len(df), replace=True) # sample all values with replacement
  # fit regression
  X = sm.add_constant(df_sampled[sub_features])
  results = sm.OLS(df_sampled.iloc[:, -1], X).fit()
  # get coeff results and save
  bootstrap_values.append(results.params)
# data frame of coeffs
df_bs = pd.concat(bootstrap_values, axis=1).T
print(df_bs.std(axis=0, ddof=2))

Answer 1

With a large dataset like yours (20,640 observations), the asymptotic assumptions should hold reasonably well, making the bootstrap and OLS standard errors similar unless there are issues:

The bootstrap method captures the variability in the data by sampling with replacement, effectively creating new samples that could have different properties from the original sample. This includes the variability of the predictors, which is not accounted for in the standard errors from a classical regression model, where the predictors are assumed to be fixed and only the variability in the response variable is considered.

The standard errors reported by statsmodels are based on asymptotic theory, which assumes that as the sample size goes to infinity, the distribution of the coefficient estimates will approach a normal distribution with a certain standard deviation. These asymptotic standard errors might not capture the actual finite-sample variability, especially if there are departures from the OLS assumptions.

If there is heteroskedasticity in the data (non-constant variance of the errors), the standard errors from the OLS regression can be biased. This is why @josef suggested using model.fit(cov_type="hc3"), which adjusts the standard errors for heteroskedasticity.

Outliers or influential observations can have a significant impact on the regression coefficients and their estimated variances. Bootstrap resamples may include these points multiple times or not at all, which can lead to greater variability in the bootstrap estimates.

If the distribution of the predictors is skewed or has heavy tails, this can lead to a larger variability in the bootstrapped coefficients compared to the OLS standard errors, which assume that the predictors are normally distributed.

The HC3 standard errors you obtained from model.fit(cov_type="hc3") are closer to your bootstrap estimates, suggesting that heteroskedasticity might be present in the data. It's often a good practice to compare standard errors from different robust estimators and bootstrap methods, especially when you suspect violations of OLS assumptions.

To address these discrepancies, you can:

• Use robust standard errors (cov_type="hc3") in your OLS model. Check for and possibly mitigate the influence of outliers or influential observations.
• Examine the distribution of the predictors and consider transformations if necessary.
• Consider whether other model specifications might better capture the relationships in your data.
• Use other forms of bootstrapping (e.g., residual, wild) that might be more appropriate given potential issues in the data.