Log-linear autoregressive Poisson model of order p and q

Modeling
,
1

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

Many thanks for any feedback provided.

Cheers

@Ara_Winter, @asael_am .

1 Like
2

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

1 Like
3

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

2 Likes
4

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

1 Like
5

Glad to hear! :) :)