如何解决HAC函数报错：Error using hac (line 485) - Index exceeds array bounds

Background #

From your description, you're hitting an index out-of-bounds error at line 485 of hac.m (b = getBW(V,weights,model,iFlag);), which traces back to the getBW function failing when ARfit.AR{1} becomes empty after the second loop iteration. This happens whether you call hac(X,y) directly or pass a fitted fitlm object, with your 224-sample dataset containing 17 variables (1 constant + 16 regressors).

Why This Happens #

The crash occurs during the automatic bandwidth selection step of the HAC procedure. The getBW function tries to fit an AR model to estimate the order of serial correlation in the residuals. If your residuals have very weak serial correlation, or the combination of 17 variables and 224 samples leads to insufficient degrees of freedom for the AR fit, the model can fail to produce valid coefficients—resulting in an empty ARfit.AR{1} and the subsequent index error.

Fixes to Try #

1. Manually Specify the Bandwidth

Bypass the problematic automatic bandwidth selection by setting a fixed bandwidth value. You can calculate a reasonable starting point using the Newey-West rule of thumb:

n = size(X,1); % 224 in your case
bw = floor(4*(n/100)^(2/9)); ~5 for n=224

Then call HAC with this bandwidth:

% Direct X/y call
EstCov = hac(X,y,'bw',bw);

% Or with the fitlm object
[EstCov,se,coeff] = hac(OLSModel,'bw',bw,'display','full');

Adjust the bandwidth value (e.g., 4, 6) if needed—this should avoid the AR fitting step entirely.

2. Check Residual Serial Correlation & Switch to White Standard Errors

If your residuals have little to no serial correlation, you might not need HAC at all. First, verify the residual structure:

% Get OLS residuals
[b,~,resid] = regress(y,X);
figure;
subplot(2,1,1)
autocorr(resid); title('Residual Autocorrelation (ACF)');
subplot(2,1,2)
parcorr(resid); title('Residual Partial Autocorrelation (PACF)');

If the ACF/PACF plots show no significant serial correlation, use White's heteroskedasticity-consistent covariance matrix instead:

EstCov = hac(X,y,'type','white');

3. Check for Missing/Outlier Values

Hidden missing values or extreme outliers can disrupt the AR fitting process. Run these checks:

% Check for NaNs
sum(isnan(X))
sum(isnan(y))

% Check for outliers with boxplots
figure; boxplot(y); title('Y Variable Outliers');
figure; boxplot(X); title('Regressor Outliers');

If you find outliers, consider winsorizing them or removing problematic observations before running HAC.

4. Reduce Variable Dimension

17 variables with 224 samples might stretch the degrees of freedom for residual analysis. Try reducing the number of regressors:

• Use stepwise regression to select significant variables:

  OLSModelStep = stepwiselm(DataTable);
  [EstCov,se,coeff] = hac(OLSModelStep,'display','full');
  

• Or use PCA to reduce dimensionality (exclude the constant term when doing PCA):

  [pcaCoeff, pcaScore] = pca(X(:,2:end)); % Skip the constant column
  XReduced = [ones(size(X,1),1) pcaScore(:,1:10)]; % Keep top 10 PCs + constant
  EstCov = hac(XReduced,y);
  

内容的提问来源于stack exchange，提问作者Jerros