# Log-linear autoregressive Poisson model of order p and q

**URL:** <https://discourse.mc-stan.org/t/log-linear-autoregressive-poisson-model-of-order-p-and-q/13613>\
**Category:** Modeling\
**Tags:** techniques, fitting-issues\
**Created:** [March 9, 2020, 7:46pm UTC](https://discourse.mc-stan.org/t/log-linear-autoregressive-poisson-model-of-order-p-and-q/13613 "2020-03-09T19:46:25Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![hve1964](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hve1964/32/7344_2.png) [@hve1964](https://discourse.mc-stan.org/u/hve1964)\
**Post date:** [March 9, 2020, 7:46pm UTC](https://discourse.mc-stan.org/t/log-linear-autoregressive-poisson-model-of-order-p-and-q/13613/1 "2020-03-09T19:46:25Z")

</div>

Hi, everyone,

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](http://dx.doi.org/10.18637/jss.v082.i05). This is how far I got:

```stan
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?

Many thanks for any feedback provided.

Cheers

@Ara_Winter, @asael_am .

---

<div class="post-metadata">

**Author:** ![asael\_am](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/asael_am/32/5377_2.png) [@asael\_am](https://discourse.mc-stan.org/u/asael_am)\
**Post date:** [March 9, 2020, 8:14pm UTC](https://discourse.mc-stan.org/t/log-linear-autoregressive-poisson-model-of-order-p-and-q/13613/2 "2020-03-09T20:14:19Z")

</div>

Did you ran the code I send you?

Maybe change this:

```
  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];
   }

```

To this:

```
for(i in 1:N){
    // ar iteration
     if(P > 0) for ( j in 1:P )if(i > j) nu[i] += log(1+y[i-j]) * beta[j];
     // ma iteration
     if(Q > 0) for ( j in 1:Q )if(i > j) nu[i] += nu[i-j] * alpha[j];
}

```

And save your self one extra loop :) :).

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [March 11, 2020, 3:38pm UTC](https://discourse.mc-stan.org/t/log-linear-autoregressive-poisson-model-of-order-p-and-q/13613/3 "2020-03-11T15:38:10Z")

</div>

I couple things from some folks smarter then me you might want to consider:

- A gamma distribution might be better than the normal for ARMA-y type models.

- This article [https://personal.utdallas.edu/~pbrandt/pests/pewma.pdf](https://personal.utdallas.edu/~pbrandt/pests/pewma.pdf) has a bunch of really useful advice and references around this topic.  
Article reference  
Brandt PT, Williams JT, Fordham BO, Pollins B. Dynamic modeling for persistent event-count time series. American Journal of Political Science. 2000 Oct 1:823-43.

---

<div class="post-metadata">

**Author:** ![hve1964](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hve1964/32/7344_2.png) [@hve1964](https://discourse.mc-stan.org/u/hve1964)\
**Post date:** [March 22, 2020, 6:27pm UTC](https://discourse.mc-stan.org/t/log-linear-autoregressive-poisson-model-of-order-p-and-q/13613/4 "2020-03-22T18:27:05Z")

</div>

Many thanks for this suggestion, Asael, much appreciated. I have now integrated this section into my code and it is doing well. Cheers.

---

<div class="post-metadata">

**Author:** ![asael\_am](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/asael_am/32/5377_2.png) [@asael\_am](https://discourse.mc-stan.org/u/asael_am)\
**Post date:** [March 22, 2020, 8:05pm UTC](https://discourse.mc-stan.org/t/log-linear-autoregressive-poisson-model-of-order-p-and-q/13613/5 "2020-03-22T20:05:51Z")

</div>

Glad to hear! :) :)