# Linear Regression with ARMA errors

**URL:** <https://discourse.mc-stan.org/t/linear-regression-with-arma-errors/18183>\
**Category:** Modeling\
**Created:** [September 24, 2020, 1:30pm UTC](https://discourse.mc-stan.org/t/linear-regression-with-arma-errors/18183 "2020-09-24T13:30:18Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![markusdumke](https://avatars.discourse-cdn.com/v4/letter/m/9f8e36/32.png) [@markusdumke](https://discourse.mc-stan.org/u/markusdumke)\
**Post date:** [September 24, 2020, 1:30pm UTC](https://discourse.mc-stan.org/t/linear-regression-with-arma-errors/18183/1 "2020-09-24T13:30:19Z")

</div>

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:

```stan
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);
}

```

---

<div class="post-metadata">

**Author:** ![asael\_am](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/asael_am/32/5377_2.png) [@asael\_am](https://discourse.mc-stan.org/u/asael_am)\
**Post date:** [September 24, 2020, 9:37pm UTC](https://discourse.mc-stan.org/t/linear-regression-with-arma-errors/18183/2 "2020-09-24T21:37:23Z")

</div>

Hola,

The code seems reasonable.

Some priors will help a lot in the performance, be aware of the chosen priors to phi and theta, cause you constrain every paramater to [-1,1] and somehow guarantee stationarity. But a normal prior with those hard constains will throw alot of errors .

Try to check your model’s performance with some simulated data. If it gives reasonable estimates (probably it will) then dont be afraid if the estimates in real data differ from a frequentist approach.

Hope it works :)

---

<div class="post-metadata">

**Author:** ![markusdumke](https://avatars.discourse-cdn.com/v4/letter/m/9f8e36/32.png) [@markusdumke](https://discourse.mc-stan.org/u/markusdumke)\
**Post date:** [September 25, 2020, 12:02pm UTC](https://discourse.mc-stan.org/t/linear-regression-with-arma-errors/18183/3 "2020-09-25T12:02:58Z")

</div>

Thank you for your feedback :)

Yes, it works for simulated data, so I think the implementation is ok. I observed that it does not work well on real data when I have few data points and when the priors on the beta coefficients are very specific, e.g. normal(-0.01, 3.05e-06). But the same model runs fine on the same data if I use normal(-0.01, 0.1) prior. I was expecting that it would be more robust in terms of prior choice. And if it diverges, I cannot generate useful forecasts as the coefficients are far off.

I wasn’t sure if the constraints on theta and phi are needed. I put them in because this is advised in the Stan user guide in the time series section.

---

<div class="post-metadata">

**Author:** ![asael\_am](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/asael_am/32/5377_2.png) [@asael\_am](https://discourse.mc-stan.org/u/asael_am)\
**Post date:** [September 25, 2020, 12:20pm UTC](https://discourse.mc-stan.org/t/linear-regression-with-arma-errors/18183/4 "2020-09-25T12:20:45Z")

</div>

Ok,

1. Those recommendations are if your data behaves stationary. Have you checked that?

The theory says that a process AR(1) is stationary if the phi parameter is between [-1,1], then for an AR§ check if every phi is between [-1,1] is not enough but helps.

1. For guarantee [-1,1] a prior like N(0,0.01) is too strong, try N(0,0.5) maybe I think I use that as default.

2. ARMA processes are really bad models for forecast, you can use State-space models for a better forecast, or use prophet for an automatic forecast.

3. If is you have some time check my package varstan, it can be installed from github:  
[https://github.com/asael697/varstan](https://github.com/asael697/varstan)

The package works similar to the forecast package, I even put the auto.sarima() function that fits the “best” ARIMA model to your data, and you can do dynamic regression as well ;)

Hopefully, I can retake the work of my package to update some other forecast’s functionalities and put it on CRAN. Hope this work for now for more stuff contact me I will be glad to help. :)

---

<div class="post-metadata">

**Author:** ![markusdumke](https://avatars.discourse-cdn.com/v4/letter/m/9f8e36/32.png) [@markusdumke](https://discourse.mc-stan.org/u/markusdumke)\
**Post date:** [September 25, 2020, 1:14pm UTC](https://discourse.mc-stan.org/t/linear-regression-with-arma-errors/18183/5 "2020-09-25T13:14:52Z")

</div>

1. Ah, I missed that these restrictions only make sense when the time series is differenced before. That might explain why it does not work well in some situations.

2. I did not set any priors for the ARMA coefficients, only for the regression coefficients \beta

3. ARMA models are State-space models as far as I know and I think it depends on the predictive power of the regressors. The ARMA term just makes sure that autocorrelation left in the residuals of a classical linear model is used.

4. Your package looks like really good work! Will definitely try this out in the future. Is it possible to use `auto.sarima()` without seasonality like `forecast::auto.arima`? As we are currently handling seasonality with fourier terms in the xreg matrix and not with a seasonal arma.

---

<div class="post-metadata">

**Author:** ![asael\_am](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/asael_am/32/5377_2.png) [@asael\_am](https://discourse.mc-stan.org/u/asael_am)\
**Post date:** [September 25, 2020, 1:34pm UTC](https://discourse.mc-stan.org/t/linear-regression-with-arma-errors/18183/6 "2020-09-25T13:34:59Z")

</div>

> I did not set any priors for the ARMA coefficients, only for the regression coefficients  
> β

If you put some priors it will work a bit better, or at least Stan will complain less (he complains to much XD).

> ARMA models are State-space models as far as I know and I think it depends on the predictive power of the regressors. The ARMA term just makes sure that autocorrelation left in the residuals of a classical linear model is used.

Well yeah SSMs are really a general class of models, but in this particular case, its predictive power is really bad.

> Your package looks like really good work! Will definitely try this out in the future. Is it possible to use `auto.sarima()` without seasonality like `forecast::auto.arima` ? As we are currently handling seasonality with fourier terms in the xreg matrix and not with a seasonal arma.

No, saddly I haven’t put that it ignores the seasonal component. But I copy Rob’s Fourier function to varstan, so It’s available when you load varstan, and you can put **fourier()** results in **reg** argument so **auto.sarima** will look the best model ARIMA with the regression residuals. (don’t tell Rob XD).