Hi, @brijo:
A Bayesian approach is going to be similar to a frequentist approach when we’re talking about observables. If there is structure in the residuals for a regression, it helps to model it. A typical approach is to use a time-series model for the residuals if they’re autocorrelated (i.e., the error at one time step is related to the error at previous time steps).
Presumably what you mean here is that the errors in Y are autocorrelated? If you want to model that, you need to figure out what the time series is to model the correlation. A vector autoregressive (VAR) prior would work in place of the multi_normal currently specified. Or you could go even further and adopt some kind of stochastic volatility model where the covariance would also vary, but those can be very hard to fit.
I’m not sure what “IV with a joint distribution” means. What’s “IV”?
With the part you document as “a little hacky”, it might be easier to define Y as an L x T matrix and then fill that and use to_vector to compute to a sequence (though you have to be careful of order, because our to_vector is column-major for matrices and row-major for arrays, so use whichever one’s appropriate).
You can code a non-centered parameterization directly now without all the intermediate transformed parameter fiddling. For example, this encodes a non-centered parameterization for beta_E whereas the version without the affine transform implemented a centered parameterization.
real<offset = mu_beta_E, multiplier = sigma_beta_E> beta_E;
...
beta_E ~ normal(mu_beta_E, sigma_beta_E);
It’s much more efficient to use a Cholesky factor formulation of the multivariate normal. You’re using Rho for the correlation matrix and theta for the scale, so you can just declare L_Rho as cholesky_factor_corr type and then just multiply it once by theta rather than using the quadratic form, and then use multi_normal_cholesky. We explain how to do this pretty early on in the User’s Guide chapter on regression.