search

crflynn / stochastic

Generate realizations of stochastic processes in python.
http://stochastic.readthedocs.io/en/stable/
MIT License
456 stars 82 forks source link

Wrong implementation of _fgn_autocovariance #48

Closed tsianosk closed 3 years ago

tsianosk commented 3 years ago

Apologies if this is not the right to communicate with the developers. By looking at the code for generating fractional gaussian noise using the Davies and Harte method, I see an issue with the calculation of the circulant matrix first row. As shown in the appendix of the 1987 paper, the first row contains the autocovariances c0, c1,...cn, cn-1, ...c1. However, the current procedure returns a sequence of cs which just gives different values for each i in [0, n]. What I am saying is also corroborated by the description in Wood and Chan (1994 and 1997) whose method in 1 dimension is identical to Davies and Harte.

If you do not agree this is a problem I would love to hear the explanation for my own understanding. Otherwise, an implementation such as 

H2 = 2 * hurst
k = np.hstack([np.arange(0, n // 2 + 1), np.arange(n // 2 - 1, 0, -1)])
return 0.5 * (abs(k - 1) ** H2 - 2 * k ** H2 + (k + 1) ** H2)

Does I believe exhibit the correct behavior.

crflynn commented 3 years ago

This was done in https://github.com/crflynn/stochastic/pull/27 by @anntzer . I'll have another look to double check.

tsianosk commented 3 years ago

If you plot or inspect the values that come out of the function you can easily see they do not repeat. Maybe this is ok. I am just pointing out what I understood from the papers.

anntzer commented 3 years ago

IIRC (it's been a while...) the trick here is that I feed this autocovariance vector into irfft rather than ifft (the DH paper calls for an (i)fft), basically taking advantage of the fact that

irfft([a0, ..., an]) == ifft([a0, a1, ..., an, an-1, ..., a1])

(modulo off-by-1's that I believe I took care of properly in the original PR, but double checking these wouldn't hurt). I agree this could have been better documented, though.

tsianosk commented 3 years ago

Ah thank you for clarifying the hidden step :) I was wondering why irfft or ifft in the first place. The wood and chan papers and also coeurjolly all use fft directly so I was trying to get everything to work in that way.

So no bug then and we close the issue. Thanks again.