Hi, I’ve got a model that is not converging well and, if the problem is what I think it is, I’m not sure the best way to go about fixing it.
I have a hidden time-varying Poisson rate and an observed discrete number of events per interval. I am trying to use an ARMA model on the rate… estimating its parameters in Stan. This means estimating the AR§ and MA(1) parameters and all the interval perturbations (x_t = \mu + \err_t + \phi_1\x_{t-1} + \phi_2\x_{t-2} … + \theta_1\err_{t-1} ). The interval perturbations’ prior is zero-mean Gaussian noise with a small stddev.
The observations seem like they would have a lot of redundant uncertainty between the Poisson distribution and the perturbation in the rate parameter. And while the posterior predictive check with a 95% interval looks OK against the data, the chains kind of drift aimlessly. It makes some sense that this would be a problem; increasing a perturbation might bring a rate closer to an observation, increasing the observation’s likelihood at the cost of the perturbation’s likelihood going down, and causing everything else in the time series to shift as well. Lots of correlation seems to imply slow exploration and convergence.
I appreciate any suggestions on how to improve this systems convergence properties, or what other better time series model I should be considering for a hidden rate like this. Hopefully lessons learned for me on this will be useful to someone else down the road as well.
It’s hard to reply in much detail without seeing the Stan model. If things are drifting, that can often be from non-identifiabilty. You didn’t say anything about priors, which can be a problem. You can check out the problematic posteriors chapter of the manual for some general advice on identifiability.
Also, you may want non-centered parameterizations.
Thanks. I will look at the priors, and see if I can reproduce the problem with just the ARMA rate and Poisson observations so I can post a concise bit of model code (and if I can’t reproduce it that way maybe that’s not the problem anyway).
As a follow-up I do believe the problem was nonidentifiability. I teased the model out into two portions, using the first model to estimate some general rate related features, and then taking MAP estimates of those features as givens in the ARMA model. Still a work in progress but it seems to be going better. Thanks.
Dear all, the statistical problem I am currently trying to tackle with a decent model parametrised in Stan code appears to be exactly of this nature: a small-number values count variable (per 15 min observation time) in dependence of day-of-week and time-of-day. However, even though I spent some time finding an adequate (Bayesian time series) model code in a number of forums, I have not yet come across anything satisfactory. The exception is the R package tscount which provides a log-linear model of order p and q based on a frequentist framework. Nevertheless, besides conceptual issues, it has the disadvandage of being intransparent as regards underlying code building. Moreover, I don’t see how one could easily subject a fitted model in this framework to posterior predictive checks and a leave-future-out cross-validation. Yes, the Stan Reference Manual has a sample code for an ARMA(1,1) representing a Gaußian data generating process. But I guess in my context I should aim for an ARMA(p,q) model and a Poisson likelihood with a log-link function. Would you have any suggestions as to where to look in order to make progress? Comments are much appreciated. Cheers.
Y ~ gamma(shape_mu + exp(b), 1 + shape_mu* exp(- mu));
shape_mu ~ gamma(0.01, 0.01); // adapted from BRMS
mu ~ normal(0, 2); // mu is your intercept.
then you have to replace exp(b) by your ARMA(p, q) structure. But I would start with an AR(1)
process first. Keep it simple first. The enhance AR(1) to ARMA(p, q).
exp(b) should be replaced by something with Y[i-1] * alpha…
Here I am stuck with just comes directly out of my brain, I would have to study the Brandt papers again.
I have been trying to code in Stan a Bayes-Laplace version of the log-linear autoregressive Poisson model of order p and q according to Eqs. (3) and (4) in Liboschik et al (2017) 10.18637/jss.v082.i05. This is how far I got:
data {
int<lower=1> N; // length of time series
int<lower=1> M; // number of independent variables plus one
int<lower=1> K; // number of groups
matrix[N,M] X; // design matrix
int<lower=0> y[N]; // time series data
int<lower=1,upper=N> P; // number of lagged observations
int<lower=1,upper=N> Q; // number of lagged conditional means
}
parameters {
vector[M] eta; // intercept and slopes
vector[Q] alpha; // slopes for lagged conditional means
vector[P] beta; // slopes for lagged observations
}
transformed parameters {
// Initialise linear predictor
vector[N] nu;
nu = X * eta;
// Add lagged conditional means
for ( i in (Q+1):N ) {
for ( k in 1:Q ) {
nu[i] += nu[i-k] * alpha[k];
}
}
// Add lagged observations
for ( i in (P+1):N ) {
for ( j in 1:P ) {
nu[i] += log(1+y[i-j]) * beta[j];
}
}
}
model {
// regularising fixed priors for intercept and slopes
target += normal_lpdf( eta | 0 , 1 );
target += normal_lpdf( alpha | 0 , 1 );
target += normal_lpdf( beta | 0 , 1 );
// likelihood w/ log-link
target += poisson_log_lpmf( y | nu );
}
I have the following queries:
(i) Could anyone please comment on the (in-)correctness of this code?
(ii) Given it was fine, is there a way to vectorise the contributions by the lagged conditional means and the lagged observations?
(iii) Would it be possible, ultimately, to have the log-linear autoregressive Poisson model of order p and q included in an updated version of the Stan reference manual?
I declare the constant beta0 outside of the matrix, Cause I have a feeling that it might be problematic (Still dont know why yet, so I just put it apart ).
(ii) Given it was fine, is there a way to vectorize the contributions by the lagged conditional means and the lagged observations?
In theory it could be vectorized but could be hard to perfom! My solution to make your code a bit faster is to just use one loop to move in the time series data, and then just see if you can do an ar( p) and ma(q) update (check pac.stan).
(iii) Would it be possible, ultimately, to have the log-linear autoregressive Poisson model of order p and q included in an updated version of the Stan reference manual?
I think you can ask to make a contribution to the stan manual. But @Bob_Carpenter can give you better details of how to do it!
Soon, I will try to add these models to my varstan package for time series :) :).
Sorry, Ara_Winter, how do a generate a “new thread”: delete this one altogether and post it again as a new entry?
Go to the principal page, and select the + new topic bottom! If you do it, please tag me! So I can try to help you more :)