I have implemented a linear regression model where the error terms follow an ARMA(p, q) process.

But when fitting the model I see divergent transitions after warmup and the AR and MA parameters are very different of the coefficients estimated with `forecast::Arima`

in R, while the regression coefficients \beta_0 and \beta_1 seem reasonable.

So I am wondering if there is an error in the implementation below?

I am new to Stan, so hope someone can help.

This is the model formula:

y_n = \beta_0 + \beta_1 * x_n + \eta_n

\eta_n = \phi_1 * \eta_{n-1} + ... + \phi_{p} * \eta_{n-p} + \theta_1 * \epsilon_{n-1} + ... + \theta_q * \epsilon_{n-q} + \epsilon_n

This is the Stan implementation I came up with:

```
data {
int<lower=0> P;
int<lower=0> Q;
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
real beta_0;
real beta;
real<lower=-1, upper=1> phi[P];
real<lower=-1, upper=1> theta[Q];
real<lower=0> sigma;
}
transformed parameters {
vector[N] eta;
vector[N] mu;
vector[N] epsilon;
for (n in 1:N) {
// Compute mean mu as linear regression of explanatory variables
mu[n] = beta_0 + beta * x[n];
// Compute ARMA error terms eta
eta[n] = y[n] - mu[n];
epsilon[n] = eta[n];
// AR-part
for (p in 1:P) {
if (n > p) {
real ar_p = phi[p] * eta[n - p];
mu[n] = mu[n] + ar_p;
epsilon[n] = epsilon[n] - ar_p;
}
}
// MA-part
for (q in 1:Q) {
if (n > q) {
real ma_q = theta[q] * epsilon[n - q];
mu[n] = mu[n] + ma_q;
epsilon[n] = epsilon[n] - ma_q;
}
}
}
}
model {
y[(max(P, Q) + 1):N] ~ normal(mu[(max(P, Q) + 1):N], sigma);
}
```