Problem in the documentation for the autoregressive moving average model

[Sghizzus]

I'm sorry if i'm wrong, but it looks like an error to me. This is the code shown in the documentation which i link here

data {
  int<lower=1> T;            // num observations
  array[T] real y;                 // observed outputs
}
parameters {
  real mu;                   // mean coeff
  real phi;                  // autoregression coeff
  real theta;                // moving avg coeff
  real<lower=0> sigma;       // noise scale
}
model {
  vector[T] nu;              // prediction for time t
  vector[T] err;             // error for time t
  nu[1] = mu + phi * mu;     // assume err[0] == 0
  err[1] = y[1] - nu[1];
  for (t in 2:T) {
    nu[t] = mu + phi * y[t - 1] + theta * err[t - 1];
    err[t] = y[t] - nu[t];
  }
  mu ~ normal(0, 10);        // priors
  phi ~ normal(0, 2);
  theta ~ normal(0, 2);
  sigma ~ cauchy(0, 5);
  err ~ normal(0, sigma);    // error model
}

in the code nu_1 is set to mu + phi * mu. This is already cryptic to me. If we assume for semplicity that y_0 and err_0 begin at their respective mean, nu[1] should be equal simply to mu. For the same reason the assignment nu[t] = mu + phi * y[t - 1] + theta * err[t - 1]; feels wrong to me. it should be
nu[t] = mu + phi * (y[t - 1] - mu) + theta * err[t - 1]: in the ARMA(1, 1) process and in general in every ARMA(p, q) it should be the sequenze y_t = x_t - mu that follows y_t = phi * y_t-1 + z_t + theta z_t-1.

Beside the above problem that i find most important, i found some inconsistency in the documentation. In the section for the autoregressive model it's suggested to leave unconstrained the parameter phi in order to assess stationarity. But in the ARMA section it's suggested that if there's no suspect of a non stationary time series to constrain the parameters phi. In the documentation it's not clear by the way what should be the constraint for higer order ARMA since constraining them in (-1, 1) it's probably not sufficient, and i don't even know if there's a constrain on the parameters in closed form

[bob-carpenter]

Hi, @Sghizzus and thanks for submitting a carefully thought-through issue. I think I wrote that example code a long time ago and I'm not an ARMA expert, so I can believe there are issues. I may have been looking at something assuming a zero-mean and messed up the translation. Looking at the ARMA description on the moving average Wikipedia page, I'm confused about what happened to mu when moving from MA to ARMA models.

in the code nu[1] is set to mu + phi * mu.

I agree that it is natural to assume nu[0] = mu if the process is stationary. But then what do we do with nu[1]? If set set nu[1] = mu + nu[0] * phi + epsilon, and assume epsilon = 0, I think we get the result in the example code. What feels most natural to me looking at it now would be to set nu[1] = mu + epsilon[1].

For the same reason the assignment nu[t] = mu + phi * y[t - 1] + theta * err[t - 1] feels wrong to me. it should be nu[t] = mu + phi * (y[t - 1] - mu) + theta * err[t - 1].

I can totally believe I slipped a difference term here. I'm having trouble relating all this back to the Wikipedia article without a board and a bit of time.

i found some inconsistency in the documentation

There's more than some as our thinking has evolved over time. Also, it's very hard to keep something like this consistent with so many authors and reviewers.

In the section for the autoregressive model it's suggested to leave unconstrained the parameter phi in order to assess stationarity. But in the ARMA section it's suggested that if there's no suspect of a non stationary time series to constrain the parameters phi.

This may be partly due to harder models being harder to fit without constraints. In general, with 1 step of history, you can constrain, but in higher dimensions, the constraints involve roots of complex equations---to code solutions in Stan would requiring doing that. The most progress I've seen on that is from Sarah Heaps, who created a VAR model constrained to be stationary. See:

• Sarah Heaps. 2024. Enforcing Stationarity through the Prior in Vector Autoregressions. JCGS.

I'm not sure if this can be transferred to ARMA models.

I'm guessing @bgoodri or @SteveBronder might know the answer.

Do you know how to submit a pull request to suggest a fix? We can find a reviewer who's better than me at time series and update the models.

Thanks again for pointing this out.

[Sghizzus]

Sorry, right now i'm in holyday, so i don't have the time to write a weel reasoned pull request. If you're not in a hurry i will propose one in a couple of weeks. For the formula part i think that the most natural development should be err[0] = 0, and nu[1] = mu. At this point we have that nu[i] = mu + phi * (y[i - 1] - mu) + theta * err[i-1] and err[i] = y[i] - nu[i]. The reasoning behind that is this.

$$ \nu_t = \mathbb{E}[X_t|X_{t-1}, Z_{t-1}] = \mu + \phi * (x_{t-1} - \mu) + \mathbb{E}[Z_t|X_{t-1}, Z_{t-1}] + \theta Z_{t-1} $$

The reamining mean fore definition equals 0. And the difference between y_t and nu_t gives the approximation of Z_t that is needed for the next iteration. The choices of the initial values feels just natural to me

[bob-carpenter]

We're not in a hurry. We do incremental doc releases, but we won't have another Stan release for a few months. Let me know if you'd rather have us fix it and I can give it a go myself.