Poisson Vector Autoregression Model

Modeling
1

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.

2

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)).

1 Like
3

Thank you!

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[1] 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.

4

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.

1 Like
5

@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!

6

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.
Best, Cord

7

@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)

3 Likes
8

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.

1 Like
9

Just refreshing this to briefly state that I’ve incorporated this functionality in my {mvgam} R package. The purpose of the package is to set up and fit complex State-Space models, including time series models in which the latent state can evolve under a variety of dynamic processes such as Vector Autoregressions. There is also functionality to do a lot of the downstream analysis that is commonly used when working with VARs, such as impulse responses or forecast error variance decompositions. Here is the blog post, hopefully it can be useful to others that find this conversation in future.

2 Likes
10

Thanks for sharing, @nicholasjclark.

Super cool that you incorporated Sarah Heaps’s parameterizations for stationary VAR processes.

1 Like