I want to estimate a hierarchical Bayesian model that has more than one dependent variable, i.e., a hierarchical model for >1 equations, where the equations are tied together via correlation between the error terms. The question is how to program this in Stan.
A twodependent variable model would be:
Sales_Online(it) = Online sales in time t for retailer i.
Sales_Offline(it) = Offline sales in time t for retailer i.

Sales_Online(it) = a0(i) + a1(i)*X(it) + e1(it)

Sales_Offline(it) = b0(i) + b1(i)*X(it) + e2(it)

a0(i) = c01 + c02*W(i) + e3(i)

a1(i) = d01 + d02*W(i) + e4(i)

b0(i) = g01 + g02*W(i) + e5(i)

b1(i) = h01 + h02*W(i) + e6(i)
Equations 1 and 2 are the equations for my 2 dependent variables. They are functions of X(it). X(it) might be price for retailer i in time t.
The function parameters a and b are in turn functions of retailer characteristics W(i) (the hierarchical structure)
A singleequation hierarchical model would just have equations 1, 3, and 4. I have above two equations (1 and 2). Note their errors e1 and e2 might be correlated for a given retailer i. Similarly e3 and e4, and e5 and e6, might be correlated among retailers. That makes the estimation nontrivial.
Question is whether Stan can handle this and how.