I have a multivariate linear model which imposes the specific structure for each vector observation \mathbf{y}_i
\mathbf{y}_i = \mathbf{b}_0 +
\begin{bmatrix}
X_{1,i} & 0 & X_{c,i} \\
0 & X_{2,i} & X_{c,i} \\
X_{1,i} & X_{2,i} & X_{c,i}
\end{bmatrix} \cdot \begin{bmatrix} B_1 & B_2 & B_c \end{bmatrix} + \varepsilon_i
where the B_1, B_2 have size m \times 3 and the B_c is n \times 1, so each \mathbf{y}_i has length 7. \varepsilon_i is a 7-dimensional IID multivariate normal with variance \Sigma.
I am looking for a Julia package that would allow me to estimate \mathbf{b}_0, B_1, B_2, B_c, \Sigma without manually building the design matrix. Or, if that cannot be avoided, allowing me to build a design matrix using a 3-dimensional array.