# Comparision: simple AR(1) model in stan vs. pymc

**URL:** <https://discourse.mc-stan.org/t/comparision-simple-ar-1-model-in-stan-vs-pymc/37529>\
**Category:** Modeling\
**Tags:** python, cmdstanpy, prior-predictive, autoregressive-model\
**Created:** [December 6, 2024, 11:00am UTC](https://discourse.mc-stan.org/t/comparision-simple-ar-1-model-in-stan-vs-pymc/37529 "2024-12-06T11:00:49Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![teagum](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/teagum/32/19496_2.png) [@teagum](https://discourse.mc-stan.org/u/teagum)\
**Post date:** [December 6, 2024, 11:00am UTC](https://discourse.mc-stan.org/t/comparision-simple-ar-1-model-in-stan-vs-pymc/37529/1 "2024-12-06T11:00:50Z")

</div>

I’am trying to replicate an AR model I once code in PyMC. Th Stan version gives yield some different results and would like to know why this is the case.

This is the PyMC code:

```python
with pm.Model() as pm_model:
    rho = pm.Normal("alpha", mu=[0, 0], sigma=[10, 1])
    sigma = pm.HalfNormal("sigma", 1)
    
    ar = pm.AR("ar", rho=rho, sigma=sigma, constant=True, init_dist=pm.Normal.dist(4, 1), steps=n_samples-1)
    y_hat = pm.Normal("y_hat", mu=ar, sigma=sigma, observed=data) 
                      
    pm_data = pm.sample()
    pm_data.extend(pm.sample_posterior_predictive(pm_data))

```

And this is the result of translating the above to Stan

```stan
data {
  int<lower=1> N;
  vector[N] y;
}

parameters {
  real init_dist;
  real rho1;
  real<lower=-1, upper=1> rho2;
  real<lower=0> sigma;
}

model {
  init_dist ~ normal(4, 1);
  rho1 ~ normal(0, 10);
  rho2 ~ normal(0, 1);
  sigma ~ normal(0, 1);

  y[1] ~ normal(rho1 + rho2 * init_dist, sigma);
  y[2:] ~ normal(rho1 + rho2 * y[1:(N-1)], sigma);
}

generated quantities {
  vector[N] y_hat;
  y_hat[] = to_vector(normal_rng(rho1 + rho2 * y[], sigma));
}

```

Event though the estimated parameters are similar, the posterior predictive samples are way more smoothed in the result from Stan. See this image below for illustration:

 ![ar1](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/a/e/ae0b7aa1223f070b5f3cdf6abd71ccbd7964d9f7.png)

Can anyone help me to understand what’s happening here?

---

<div class="post-metadata">

**Author:** ![spinkney](https://avatars.discourse-cdn.com/v4/letter/s/dec6dc/32.png) [@spinkney](https://discourse.mc-stan.org/u/spinkney)\
**Post date:** [December 6, 2024, 12:05pm UTC](https://discourse.mc-stan.org/t/comparision-simple-ar-1-model-in-stan-vs-pymc/37529/2 "2024-12-06T12:05:32Z")

</div>

Probably the lower and upper bounds on rho2. I don’t how pymc would know that.

---

<div class="post-metadata">

**Author:** ![teagum](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/teagum/32/19496_2.png) [@teagum](https://discourse.mc-stan.org/u/teagum)\
**Post date:** [December 6, 2024, 1:18pm UTC](https://discourse.mc-stan.org/t/comparision-simple-ar-1-model-in-stan-vs-pymc/37529/3 "2024-12-06T13:18:59Z")

</div>

Thanks for the hint. Interestingly it’s not the bounds. But I found that, I can scale the smoothness, by forcing rho1 in a certain direction. So, for example, setting rho1 ~ normal(2, 0.5) the result is nearly identical to pymc.

However, this poses some further questions. The data was generated with rho1 = 4. rho1 ~ normal(2, 0.5) produces an adequate result without any numerical issues. However, using rho1 ~ normal(4, 0.1) produces a worse result, including numerical issues. Is suspect, this is because the first prior fits this very sample better.

But why is the model in Stan so susceptible to the prior? I’d guess because of the small sample size. However, in pymc I can give a very broad range of prior parameters, and it’s always converging to the expected result without any numerical issues.  
Any further thoughts on that?

---

<div class="post-metadata">

**Author:** ![hhau](https://avatars.discourse-cdn.com/v4/letter/h/2bfe46/32.png) [@hhau](https://discourse.mc-stan.org/u/hhau)\
**Post date:** [December 6, 2024, 1:37pm UTC](https://discourse.mc-stan.org/t/comparision-simple-ar-1-model-in-stan-vs-pymc/37529/4 "2024-12-06T13:37:35Z")

</div>

Can imagine this being a difference in parameterisations – pymc [seems to do it via the initial value + the innovations](https://github.com/pymc-devs/pymc/blob/main/pymc/distributions/timeseries.py#L677-L706). The latter are probably more likely to be zero-mean / more like a non-centred parameterisation and perhaps “easier” to sample from?

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [December 10, 2024, 9:10pm UTC](https://discourse.mc-stan.org/t/comparision-simple-ar-1-model-in-stan-vs-pymc/37529/5 "2024-12-10T21:10:17Z")

</div>

> [@teagum](#):
>
> Thanks for the hint. Interestingly it’s not the bounds.

> [@teagum](#):
>
> I can scale the smoothness, by forcing rho1 in a certain direction

I’m curious how you determined that, because adding bounds will force `rho1` away from the boundaries and can bias estimates compared to the unconstrained solution. I’m pretty sure the PyMC code as written does not put bounds on the `rho` values. It is not typical to put bounds on these values in an AR(1) model.

I would suggest simplifying your model so that `y[1] ~ normal(4, 1);` if you think `normal(4, 1)` is a good distribution for the initial state. Alternatively, you can push all the normals through and keep the same behavior as you have now without the `init_dist` parameter. What you get is additional variance from `rho2 * init_dist`. I think it works out to the following, but you should check my algebra, taking

```no-highlight
y[1] = rho2 * normal(4, 1) + normal(rho1, sigma)
    = normal(rho2 * 4, rho2) + normal(rho1, sigma)
    = rho2 * 4 + rho1 + normal(0, rho2) + normal(0, sigma)
    = rho2 * 4 + rho1 + normal(0, 1 / sqrt(1 / rho2^2 + 1 / sigma^2)))
    = normal(rho1 + rho2 * 4, 1 / sqrt(1 / rho2^2 + 1 / sigma^2))

```

so that

```stan
y[1] ~ normal(rho1 + rho2 * 4, inv(sqrt(inv(square(sigma)) + inv(square(rho2)))));

```

In any case, I wouldn’t recommend calling a scalar value `init_dist` as it’s not a distribution.