hac results vary bewteen MATLAB R2015b and MATLAB R2016b

Hi, my problem was related to the matlab function "hac". I used the same code last year (with R2015b version) and obtained the heteroscedasticity and serial correlation consistent standard errors (StdErr* 7.9709 0.0321 3.7480 0.0404 2.9265 10.7310 )and pvalues (pValue* 1.89E-01 1.49E-02 8.13E-21 3.30E-02 0.0033 0.0414), shown as below. Under R2016b version, the function now reports error message: "Error using hac (line 559) Matrix must be positive definite." Can you help me to figure out what happened? thanks a lot. Note: df_ttest = 27

[Walter Roberson]

Note: df_ttest is not defined in your .mat or your code.

However, that line of your code is not reached until after hac() has worked, so it is not a primary problem.

Answer 1

[Walter Roberson]

A model is constructed fitting the data. The 5th row of the data is predicted exactly, residual 0. That leads to a u matrix that contains a 0, which leads to a scores matrix, V, which contains a row of 0, which leads to a covariance matrix PhiHat which is not full rank, which requires there to be an eigenvalue of 0. The matrix is not positive definite, only positive semi-definite, so chol() fails. Due to round-off error if you eig(PhiHat) you might see a slightly negative eigenvalue (around 1E-18) but if you compute the eigenvalue symbolically you can find it is exactly 0.

I do not know why it happens to be the 5th row which is predicted exactly. Hang Qian's discussion of the third row might be the key point.

[Walter Roberson]

There is no relevant change to the code of hac between R2015b and present. There are some changes to handle the string datatype and there is a change to how the data and missing data is extracted from the model that is created (the older version did not make a relevant property public so there was a hack before), but there is no change to the algorithms of hac.

There are changes in eig() -- when I save the R2015b-produced matrix and ask to eig() it under a newer version, the eigenvalue that should be 0 comes at slightly negative under the new version and slightly positive under R2015b. In both versions chol() of the matrix works, so I will need to go back further to determine where the difference in matrices is coming from.

[Walter Roberson]

The values predicted by the model differ by essentially round-off error. In R2015b the value predicted for the 5th row is on the order of 1E-14 instead of 0 for the newer version. This is just enough for chol() to be happy with the derived values, which would not have a row of 0 as is the case in later versions.

I am still chasing to determine the step at which the predictions vary.

[Walter Roberson]

This issue traces back to different outputs from qr() for the same inputs. I do not see a Release Note related to qr or anything numeric; the closest I see is the R2016b release note,

Performance is improved for MATLAB code that is dominated by tight loops, straightforward indexing, and simple math.

which does not apply for qr()

[Walter Roberson]

I created a service request about the undocumented change to qr() as of R2016a.

Answer 2

[Hang Qian]

Hi Lucia,

Thank you for reporting this edge case.

The codes work on my Windows computer both in R2016a and R2016b. The results are the same as what you reported. Your hardware and software platforms might be different than mine. Occasionally, the results could be different in edge cases. It appears that the error message comes from the Cholesky decomposition of a matrix (named PhiHat in line 557). When I run the codes, I saw that matrix is ill-conditioned, and thus it is possible that Cholesky decomposition fails on your computer.

The predictor matrix X has some problems in this case. The values of the third variable are all zeros except for a single observation being one. Meanwhile, a constant term is included as the first predictor. The data have strong collinearity. A regression might not work well for that dataset.

Also, try to adjust the scale of the regressors, either by the data scale or the OLS estimator scale. It may help improve conditioning and overcome the edge case.

Regards,

Hang Qian

[Lucia Ryan]

Hi Hang, I just realized that my previous version was R2015b, instead of R2016a.

I run the same data under R2015b and R2016b. no issues reported under the older version, but show error message under R2016b.

With respect with variable X, the first column represents intercept, but while running function hac(), only X(:,2:end) was used.

I need to understand why this happens, as I need to document the difference in outputs. Thanks a lot.