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

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;

for ( i in (Q+1):N ) {
for ( k in 1:Q ) {
nu[i] += nu[i-k] * alpha[k];
}
}

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

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

1 Like

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

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.