Very slow sampling for time series with non-linear trend

Modeling
1

I have some issues in appropriately modelling a time series with a non-linear trend.

My model is in fact a negative binomial model for the evolution of mortality and the “period/time” effect is usually modelled with a AR(1) process with drift:
\kappa_t=c + \rho \kappa_{t-1}+\epsilon_t .

If you fit this model to death rates of UK, you would find:

kappa
kappa756×375 33.7 KB

Hence, you see that the mortality is decreasing over time but since 2010, there is a deceleration.
To allow for non-linearity, I have considered the following model:
\kappa_t=c +\delta t +\rho \kappa_{t-1}+\epsilon_t .
Such a model would allow for non-linear trend if \rho=1 and \delta \neq 0.
However, when I run such a model in Stan, the sampling becomes very slow (10 times slower compared to the random walk with drift specification). After analyzing the convergence, it seems that such a model has an identifiability problem which, I think, is the cause of the problem.
Indeed, if the time-series is linear, it can be specified either via
\kappa_t=c +\kappa_{t-1}+\epsilon_t or \kappa_t=c +\delta t +\epsilon_t .
This is also confirmed if we plot \delta against \rho which are perfectly correlated (plot from ShinyStan):

kappa2
kappa22788×952 108 KB

My question then is: How can I generalize the AR(1) model to allow for non-linear trend while still having identifiable parameters?

Thanks a lot in advance!

2

Hola,
AR process are linear by definition so, maybe that is a problem there.

Why dont you fit a non linear model with simple regression and then you model the time series dependency with an AR process or a filter.

Other thing is that for a good inference you need stationarity. Another way to get a stationary ts is applying 1 difference to your data. That should make the sampling faster and more efficient.

If doesnt help post the code and data so I can give a look.

Hope it works.

2 Likes
3

Thank you @asael_am for your reply !

I agree with you that with a random walk with drift or AR(1), you can only model linear trend. That is the reason why I proposed above to add a linear trend to the random walk:
\kappa_t- \kappa_{t-1}= c + \delta t

Such model allows for quadratic trend but has some identifiability problems (as explained above) which slow down the sampling considerably. In fact, this time-series is part of a more complicated Negative-Binomial model used to project mortality rates, which has the following specifications:

modelspecifications
modelspecifications1630×1406 132 KB

You can find my stan code enclosed.

CBDmodel.stan (2.3 KB)

P.S: Stupid question: I can’t upload Rdata (list object) in Stan forum to share my data. What is the best way to share it then?

Thanks in advance!

4

Hola,
I check your model… hard to understand and have so many questions in a lot of parts, I send you my stan code version of the described model (or at least what I understand), check if this is a bit faster.

cbd_model.stan (2.1 KB)

A way to make more efficient your code is to estimate sigma using Stan’s recommendations for estimating covariance matrices. If there are still struggles we can keep talking about it. Or send me a Github link where you store the data and I can try too.

1 Like
5

Thank you very much @asael_am for writing your own version of the model!

I was able to run your model on some mortality dataset but still it takes 1755 sec to sample (approximately 30min). Moreover, when I diagnose the convergence, it still didn’t converge…
My opinion is that Stan is unable to fit such a model in less than 30 minutes while it would only take 5 seconds to fit in the frequentist setting…

For your info, I also join a simplified version of your code (there were some errors in the previous version). You can try such a model on mortality data and you will notice that the sampling completely fails. I think the reason is that Stan can’t handle overparametrized model, that is also the reason why I implemented the bivariate normal via

target += normal_lpdf(first_variable | first_means, first_sigmas);
target += normal_lpdf(second_variable | second_means + 
  rhos .* first_sigmas ./ second_sigmas .* (first_variable - first_means),
  second_sigmas .* sqrt(1 - square(rhos)));

following the answer of Ben Goodrich from this post. This is much more efficient than looping over standard normal epsilon.

cbdmodel2.stan (1.9 KB)