Hello everybody,

I’m wondering how to implement a bivariate model with an AR error model in which the error terms are correlated. To give an idea of what I mean, I’m attaching an image from *Beyond the Cross-Lagged Panel Model: Next-generation statistical tools for analyzing interdependencies across the life course*:

The key point is to implement a cross lagged model that includes random intercepts. The cross lagged terms are displayed as c1 and c2 in the figure above.

I tried something like

```
formula_model <- (
bf(mvbind(value1, value2) | mi() ~ poly(week, 2) + (1|p|subject))
)
fit_nonlinear <- brm(
formula_model,
autocor = ~ ar(p=1),
data=results,
chains=1,
iter=2000)
```

but this only yields two independent AR processes.

- Operating System: Windows
- brms Version: 2.15.0