Repository navigation
Problem in the documentation for the autoregressive moving average model #884
Description
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
Activity
bob-carpenter commentedon Jul 6, 2025
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
muwhen moving from MA to ARMA models.I agree that it is natural to assume
nu[0] = muif the process is stationary. But then what do we do withnu[1]? If set setnu[1] = mu + nu[0] * phi + epsilon, and assumeepsilon = 0, I think we get the result in the example code. What feels most natural to me looking at it now would be to setnu[1] = mu + epsilon[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.
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.
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:
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 commentedon Jul 7, 2025
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.
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 commentedon Jul 9, 2025
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.