I would like to reproduce a Poisson Vector Autoregression Model by Brandt and Sandler with stan (they used BUGS).
I found a basis for a Vector Autoregression (VAR) model at the depricated google mailing list: here and here.
Here is the respective code from Ryan Batt (thanks!):
VAR_p_cov.stan (1000 Bytes)
Essentially, I would like to insert a hierarchical layer, where the dependent variables are poisson distributed and the latent lambdas are modeled by the VAR process.
I have problems getting the latent lambdas into a matrix form. Forgive me if this is trivial; I am new to stan and any help is greatly appreciated.
That’s still back in our old syntax. Use
= instead of
<- for assignment. And you can also vectorize a lot more now and use
a += b instead of
a = a + b. And you can declare-define and slice arrays, so the transformed data is now a one-liner (I think I got the indexing right, but you should double check):
vector[M] Y_obs[T - P] = y[P + 1: ];
It’s definitely not trivial working with all these structures, especially when they’re rolled out into more efficient array of vector formats.
I don’t have time to go through the paper—is the idea just that you replace the continous
Y with count data with an underlying continous log mean? If so, then you can declare a
lambda with the same structure as
Y for the latent parameter, and use that where you would have had
Y before. Then take the observed
Y and use
Y ~ poisson_log(lambda);
It’s important it be the Poisson on the log scale—that’s equivalent to
Y ~ poisson(exp(lambda)).
I updated the syntax.
The declare-define resulted in errors when I have multipe declarations. As far as I understand, all declarations need to be set first per block.
The idea is exactly that. Currently I am stuck with a sampling error:
Error evaluating the log probability at the initial value.
Exception: multi_normal_cholesky_lpdf: Random variable is nan, but must not be nan! (in ‘…’ at line 46)
I created an R file that generates some data. PVAR-R_2018-03-08.R (2.8 KB)
My current stan model
PVAR-stan_2018-03-08.stan (1.8 KB)
Again, thank you so much: your suppport is greatly appreciated.
Yes, all declarations need to come first. But any or all of those declarations may come with definitions. So what probably resulted in errors was a mismatch between the right-hand and left-hand side types.
The problem here is that you don’t define
lambda before using. It’s declared as a local variable in the model, but never defined. Every variable has to be defined before it can be used. Data variables and parameter variables are defined externally—every other variable is up to the user to initialize. The sampling statement with
~ just increments the log density—it doesn’t literally sample and assign to the left-hand side—that wouldn’t work with Stan’s langauge semantics, which involves computing the log density so that we can differentiate it and use Hamiltonian Monte Carlo.
@cord-otten I ran across this while trying to figure out an implementation of Brandt and Sandler’s 2012 model. Did you get the above code working, and if so, would you mind sharing it?
Thanks for your time!
Hi, sorry but I did not get it to work. That’s very likely due to my limited ability with Stan, though. If you manage to set it up I’d still be curious.
@jkalley I wrote a VAR(1) model with a poisson link in stan that seems to work on simulated data. Disclaimer: I’m an applied social scientist, not a statistician. Feedback is welcome. I’m presently not highly confident in my prior choices. I’m fitting this on activity levels in social media groups of varying scales so the log link seems appropriate. I might extend this to negative binomial models if it seems useful.
Generalizing to VAR(p) models should be straightforward.
var_1_pois.stan (1.7 KB)
Great- thanks for sharing! Wound up not needing it for our project, but I’m glad you put that up in case it’s of interest to someone else.