How does multivariate regression consider the correlation among response when estimating the parameters? I have read a few posts about multivariate regression, where the response is a matrix instead of a vector. However, when I really look into the literatures, there is nothing.
I have read the multivariate books (the books recommended by people answering similar questions), all they talk about is multivariate distribution, Wishart distribution, MANOVA, PCA, factor analysis, etc. None of them discussed how to take the dependency among y into account.
Is there any literature about how to do the 'multivariate regression' instead of talking about how to deal with multivariate data, etc.?
 A: Assume your multivariate regression means multivariate linear model. Multivariate linear model is $$Y=X\beta + \epsilon$$ where the dimensions of these matrices are: $Y: n \times k$, $X: n \times p$, $\beta: p \times k$, $\epsilon: n \times k$. Let $Z_i$ be the row vector of $i$-th row in matrix $Z$. According to the assumptions of multivariate linear model, $$\epsilon_i' \sim N_k(0_k, \Sigma)$$ and $\epsilon_i$ and $\epsilon_j$ are independent for $i \ne j$. Let reshape the matrices as following:
$$Y_R = (Y_1, Y_2,... Y_n)'$$
$$X_R = \left(\begin{matrix} X_1  & 0_k&...&0_k \\0_k &X_1 & ... &0_k\\...&...&...\\0_k&0_k&...&X_1\\
...&...&...&...\\X_n  & 0_k&...&0_k \\0_k &X_n & ... &0_k\\...&...&...\\0_k&0_k&...&X_n
        \end {matrix} \right)$$
$$\epsilon_R = (\epsilon_1, \epsilon_2,... \epsilon_n)'$$
$$\beta_R= (\beta_{11}, \beta_{21}, ..., \beta_{p1},\beta_{12}, \beta_{22}, ..., \beta_{p2},...... \beta_{1k}, \beta_{2k}, ..., \beta_{pk})'$$
We have univaraite linear model:
$$Y_R = X_R\beta_R +\epsilon_R$$
Here 
 $$\epsilon_R \sim N_{kn}(0_{kn}, I_n \times\Sigma)$$ 
The $\times$ here means  Kronecker product, should have a circle. 
It should be easy to fit this model using current computer and software.
A: If you're fitting a linear model, you can try looking at seemingly unrelated regression. 
It turns out it's not all that important to model dependencies among different outcome variables. Estimating a series of models one-by-one using OLS is unbiased and consistent for the true coefficients (though asymptotically less efficient than SUR) and is standard practice in almost all social science research. 
