Formula to calculate beta matrix in multivariate analysis I have to implement a multivariate analysis on  $n$ random variables with a sample of $m$ data points. I would like to get a matrix with the $\beta$ (as in $n$ $\beta$ vectors put together). 
Is there an analytical formula for that or do I have to calculate each vector individually and put them back together ?
Also is there a formula that would give me directly a vector of the estimation error for each variable?
 A: I think I understand what you're asking, but correct me if I'm wrong.  The analytical formula for $\beta$ is the same for the multivariate case as the univariate case:
$$
\hat \beta = (X'X)^{-1}X'Y
$$
You find this the same way as for the univariate case, by taking the first derivative of residual sum of squares.  It is relatively straightforward to calculate using matrix calculus (which is covered in the matrix cookbook linked to by queenbee).  You can test whether this solution works in R:
y <- cbind(rnorm(10), rnorm(10), rnorm(10))

x <- cbind(1, rnorm(10), rnorm(10), rnorm(10),
       rnorm(10), rnorm(10), rnorm(10))
colnames(x) <- paste("x", 1:6, sep = "")
colnames(y) <- paste("y", 1:3, sep = "")

fit <- lm(y ~ x - 1)
summary(fit)

anaSol <- solve((t(x) %*% x)) %*% t(x) %*% y
anaSol

coef(fit) - anaSol

Here's another reference, specifically related to multivariate analysis:
http://socserv.mcmaster.ca/jfox/Books/Companion/appendix/Appendix-Multivariate-Linear-Models.pdf
A: If you have $q$ equations and $p$ independent variables (including a constant) that appear in every equation, the parameter estimates are given by the $p \times q$ matrix:
$$M=(X'IX)^{-1}X'IY$$
where


*

*$Y$ is $n \times q$ matrix of dependent variables

*X is $n \times p$ matrix of covariates

*I is the identity matrix

A: Other answers nicely cover how to derive the $\beta$ coefficients. I'm not sure what you mean by $n, \beta$ "put together." But, if it means that you'd like to use the coefficients to derive the model's predicted values using the coefficients, it's simply the product $XB$, where $X$ is an $m \times n$ matrix of $m$ observations, each with $n$ independent variables, and $B$ a $n \times p$ matrix of regression coefficients. ($p$ here is the number of dependent variables.)
To your question about the standard error, for a single independent variable and single coefficient, the formula is:
$s.e.(\beta_j) = \sqrt{s^2 (X'X)^{-1}_{jj} }$
where $s^2$ is the sum of squared residuals, given by $\sum_i y_i -\hat y_i $, over $m - n$. (More here.) To broaden the formula to return a vector of standard errors corresponding to each coefficient:
$s.e.(\beta) = \sqrt{s^2 diag(X'X)^{-1} }$
where $diag$ returns the diagonal of the matrix. To further broaden that formula to return a matrix of standard error corresponding to each $\beta$ coefficient, replace $s^2$ with $S^2$, a row vector containing the $s^2$ for each independent variable:
$s.e.(B) = \sqrt{diag(X'X)^{-1} S^2 }$
Note that the vector returned by $diag$ should be $n \times 1$, and $S^2$ $1 \times p$, making their product $n \times p$, standard errors corresponding to the coefficients in $B$. (Square root again applied element-wise.)
