Poisson time series analysis with multi-component random walks

[martinmodrak]

Hi,
sorry for not getting to you earlier. This looks like a tough problem, although I don’t full grasp a lot of the details as I don’t have physics/signal processing background and can’t translate sentences like "random walk to have, say, a power-law power spectrum " into math.

@j13r

Q1 : My first Stan question here is whether it is preferable to make these intermediate random variables Stan parameters instead. My suspicion is yes, because otherwise the derivatives are very noisy, and HMC will complain.

This is correct. Also note that you can build linear state space models (i.e. the stuff Kalman filters are built for) with arbitrary likelihood, but you then need to make the latent states explicit parameters and can’t integrate them out (as Kalman filter does). Not sure if that’s good approach for you, but I’ve built a model with this and it worked OK-ish.

@j13r

( Q2 ) how do I best place a powerlaw slope prior on them.

So if I understand it correctly, PSD is something you can compute from the ARMA parameters. So a direct approach would be to have PSD as an explicit parameter in your model (so you can directly put prior on it), and remove one of the ARMA parameters and instead compute it from the PSD and the rest of the parameters. If there is no analytic solution, you can use algebra_solver to solve the equation numerically (but performance will suffer). Not sure if this is feasible for your case.

@j13r

Q3 : I believe one can place a powerlaw PSD prior on the covariance kernel, but how does one do this exactly? Specifically, the covariance decline with |r| is not the same as the PSD slope, if I understand correctly, but I couldn’t find the right formulae. I would appreciate pointers.

This is unfortunately completely out of my expertise (while the points above are only somewhat out of my expertise :-). I know that there is “Spectral approximation to Gaussian processes” (e.g Approximate GPs with Spectral Stuff and https://discourse.mc-stan.org/t/hermite-polynomials-and-their-inverse-matrix-continued-from-approximate-gps-with-spectral-stuff/37930), but I am not sure if they use “spectrum” in the sense you use it. So just a wild guess that this might be relevant. (once again, I don’t understand the spectral stuff neither in GPs nor in your application :-/)

@j13r

The convergence on these states is poor, but I guess that is just how it is, because they must be highly correlated but are poorly constrained. The situation is different for the times with observations, where the states are well constrained. I tried to fill in these Npremodel points by running the AR(1) backwards, but it is not performing much better. Q4 : Is there a smarter way to go about this?

One thing that sometimes helps with random walk models is that instead of working with the states directly, you take lognorm_raw ~ normal(0, 1) and compute lognorm by scaling and summing lognorm_raw appropriately. This could remove some of the correlation.

Also we usually recommend enforcing constraints on the parameters that are not physical as priors rather than bounds. Actually you kind of do this, so when you have drw_logsigma ~ normal(-1, 0.3);, you should not need real<lower=-3,upper=max_time_logsigma> drw_logsigma; (assuming max_time_logsigma > 0) - those configurations are already quite strongly excluded by the prior. If your chains are moving such far from the prior that the bounds have an effect, it is likely there is an issue with the model (a bug or simply a prior-data conflict).

Other than that the model is too complex for me to understand well, but hope this helps at least a bit.

Good luck!