First off, I post this question in Physics since it stems from a physics problem, but it may be more pertinent to signal processing; sorry if it's the wrong place.
TLDR: what are the correct normalization factors to compare an analitycal power spectrum with a stochastic one?
Lengthy explanation: I'm studying a simple Lotka-Volterra system, following quite precisely this reference (since I want to understand the procedure in order to apply it to a model I'm doing research on I'm using the same parameter values and everything). I'm solving both the deterministic version, using the function scipy.integrate.solve_ivp in python, and the stochastic version integrating the associated Langevin equation.
A first incongruence I'm meeting is that my solutions seem to be faster than those in fig. 1 of the article, which shouldn't be the case since, as I mentioned, I'm (hopefully) doing everything in the same way as in the paper; the following are my solutions:
Instead, in the solutions shown in the article, the (MF) time series stop oscillating around t=400. Anyway, the good thing is the two solutions agree, in the sense that the stochastic ones oscillate around the MF ones at stationarity (which is exactly the main finding of the article). Given this phenomenon of stochastic amplification, the article goes on evaluating the power spectrum of the trajectories, both analitically from the Langevin equation and numerically from the stochastic time series. In figure 2 of the article, the two power spectra are shown to perfectly agree, meaning they are basically superimposable. The analytical power spectrum is the following (eq. 7 in the article):
The numerical one is obtained as the average of the power spectra of anumber of stochastic simulations; each ps is calculated exploiting the scipy.signal.welch method; following is the line of code calculating the spectra:
freqw, psw = scsig.welch(z[0,ips,ii],fs=1/dtlang,scaling='spectrum',nperseg=2**12,nfft=2**15)
The first two arguments are simply my time series and the sampling frequency, which is the inverse of the timestep used in the stochastic simulations. The following is the resulting power spectrum:
Comparing the two spectra, it seems that the resonant frequencies agree (although they are different from those of the article, but I guess it's the same problem I noticed with the speed of convergence of the time series), but there is a discrepancy of several orders of magnitude in their value, while in the article they perfectly agree. Moreover, both amplitudes are different from those of the spectra shown in the paper.
I tried putting many multiplicative factors to get them to be compatible, including some factors with N=3200 used in the article for individual-based Monte Carlo simulations, but without any success; plus, I'd like to find meaningful factors a priori insted that putting random ones to get the expected result and justify them a posteriori.
So, my question is: what factors should I put to get compatibility between the two spectra? What am I missing?
Thanks!
