Seeing as how I had a similar question earlier and came across this long-unanswered question through a simple web search, I'll take a stab and post what I think is one possible solution to your situation that others may also be encountering.
According to SAS Support, you can take the time-series you have and fit an intercept-only regression model to the series. The estimated intercept for this regression model will be the sample mean of the series. You can then pass this intercept-only regression model through the SAS commands used to retrieve Newey-West standard errors of a regression model.
Here is the link to the SAS Support page:
http://support.sas.com/kb/40/098.html
Look for "Example 2. Newey-West standard error correction for the sample mean of a series"
In your case, simply try the same approach with Matlab.If someone has a better approach, please enlighten us.
That would indeed be surprising, and although @Aksakal is right that without serial correlation, the estimates should not differ too much, it would still be a fluke if they were exactly equal (and I presume you have run your code more than once).
One possible explanation (hence at best a partial answer) could be that your bandwidth (I am not sure how they are precisely chosen) is such that your HAC estimates boil down to HC estimates, which, in turn, boil down to the classical ones when only regressing on a constant:
The robust variance matrix is
$\hat V= n(X'X)^{-1}(X'\hat\Omega X)(X'X)^{-1}$
with $\hat\Omega$ the covariance matrix of the residuals. Here, $X$ is just a column of 1s. Hence, $(X'X)^{-1}=1/n$. Now, if $\hat\Omega$ is a diagonal matrix (possibly because the weights are such that only the off-diagonal elements receive zero weight), we obtain
$$X'\hat\Omega X=\sum_{i=1}^n\hat u_i^2$$
such that
$\hat V= n\frac{1}{n}\sum_{i=1}^n\hat u_i^2\frac{1}{n}=\frac{1}{n}\sum_{i=1}^n\hat u_i^2$
Upon inserting the degrees of freedom correction $n/(n-1)$ used in the default OLS formula, we get
$$\frac{1}{n-1}\sum_{i=1}^n\hat u_i^2.$$ In R (I do not have access to MATLAB):
library(sandwich)
n <- 10
y <- rnorm(n,1)
reg <- lm(y~1)
vcov(reg)
vcovHC(reg,type="HC0")*n/(n-1)
Best Answer
HAC procedures are just about providing consistent estimates of the standard errors. They do not change the estimation of the coefficients. If you have strict exogeneity with serial correlation, your coefficients are unbiased, but the standard errors are incorrect. HAC standard errors address the latter point.
As you allude to, this does not give efficient coefficient estimates. To achieve efficiency, in economics, at least, we typically use a Cochrane-Orcutt/Prais-Whinston procedure. This requires much stronger modeling assumptions to estimate the structure of the serial correlation, however.
They are analogous to Eicker-White heteroskedasticity robust standard errors. This procedure does not alter estimation, it only changes the estimates of the standard errors to ensure that they are consistent in the presence of heteroskedasticity. The efficient fix would be weighted least squares, but this requires modeling the form of the heteroskedasticity.