Support for autocorrelation in response variables

[eb8680]

brms supports autoregressive correlation structure on response variables via the cor argument to its fit function and the cor_arma object. It would be nice to support at least a subset of this functionality, say only AR(1) terms with grouping.

In conjunction with support for multiple response variables #42 and unobserved response variables #43 #44 and a syntax extension prev for indicating AR(1) autocorrelation, this would allow brms to provide general and highly scalable support for most multivariate structural time series models with an intuitive syntax:

x ~ 1 + prev(x) + (prev(x) | g)
y ~ 1 + (1 + x | g)

Semantically, formula language support for autocorrelation and latent response variables #44 should be backed by code generation targeting the parallel-scan algorithms for Gaussian and discrete state-space models implemented in Pyro and Funsor, so that we could work with very long time series in O(log(T)) time on parallel hardware.

That might mean implementing a vectorized pyro.markov statement upstream (analogous to the vectorized pyro.plate) and, when autocorrelation terms are included in a model, generating code where the observations are wrapped in a pyro.markov context instead of a pyro.plate. We might also need to slightly modify Pyro's HMM distributions to generate latent Markov chains and implement an elbo variant capable of exploiting this structure:

...
with pyro.markov(dim=-2, window=1):
    x = pyro.sample("x", dist.DiscreteMarkovChain(...))
    y = pyro.sample("y", dist.Normal(loc_y[..., x], scale_y), obs=y_obs)
...

cc @fritzo