Lag-1 autocovariance incorrectly normalized in autocorrelation coefficient calculation

[cmicek1]

It looks like there is some legacy code included in the package (ar1nv.m) responsible for estimating the AR(1) model parameters that are used to make the red noise for significance testing. However, the calculation for the autocorrelation coefficient (g) that is used incorrectly normalizes the lag-1 autocovariance. The existing implementation divides by N-1, which makes sense in isolation because there is one fewer element in the numerator, but when divided by the lag-0 autocovariance this can in some cases result in a value of g that is greater than 1. (Say the dot products in the numerators of both c0 and c1 are very close, then g = c0/c1 ≈ N/(N - 1), which will be greater than 1.) This leads to the "Subscript indices must either be real positive integers or logicals" error in rednoise.m that some users have observed.
A more correct approach I think would be to either normalize c1 by N instead of N - 1, or trim the numerator of c0 to be length N - 1, and normalize c0 using N - 1 as well. (Or better yet, remove this function altogether and substitute with an existing MATLAB implementation instead.)

[grinsted]

The "nv" in ar1nv is intended to stand for "naive". There is a better estimator in the AR1.m file. But really you can use any estimator you like. I'm reluctant to change the code because:

1. I very rarely work in matlab myself anymore.
2. it is adopted as provided from the source given in the comments.
3. In practice I think that it would make virtually no difference for any time series of a length where you would consider applying wavelet analysis.

[cmicek1]

That's fair! It definitely has made a difference for some of the data I was analyzing, though. In some cases the code will error when using the naive estimator, but using other more accurate estimators does not cause this issue.

[LucyRah]

Hey,
as I run into the same issue of erroring out.. do you have a recommendation which AR1 estimation would be best or which MATLAB implementation to use instead?
Thankful for any help and all the best

[cmicek1]

@LucyRah Since that bit of code is only for hypothesis testing, and I used a permutation-based approach (similar to the random pair analysis described here) for that instead, it actually never ended up mattering for me.

But MATLAB has plenty of functions you can substitute for it instead, e.g. https://www.mathworks.com/help/signal/ref/arcov.html. When I was testing them out back in November they all seemed to perform similarly.