# How calculate variance-covariance matrix of coefficients for multivariate (multiple) linear regression?

How to calculate a variance-covariance matrix of coefficients for multivariate (multiple) linear regression?

Something like (equation below), but for the multivariate case.

Being more specific I'm interested in equations for diagonal terms. The question is different from:

How to derive variance-covariance matrix of coefficients in linear regression

Because of the MULTIVARIATE regression case.

$$V[\hat{\beta}] = V[(X^{T}X)^{-1}X^{T}Y]\\ = (X^{T}X)^{-1}X^{T}V[Y]X(X^{T}X)^{-1}\\ = \sigma^2(X^{T}X)^{-1}$$
This is unfortunately challenging to calculate. Examining $$X^{T}X$$ we see (assuming an intercept is included in the model),
$$X^{T}X = \begin{bmatrix} n & \sum_{i=1}^{n}x_{i1} & \sum_{i=1}^{n}x_{i2} & \ldots& \sum_{i=1}^{n}x_{ik}\\ \sum_{i=1}^{n} x_{i1} & \sum_{i=1}^{n} x_{i1}^2 & \sum_{i=1}^{n} x_{i1}x_{i2} & \ldots & \sum_{i=1}^{n}x_{i1}x_{ik} \\ \vdots & \vdots & \vdots & & \vdots\\ \sum_{i=1}^{n}x_{ik} & \sum_{i=1}^{n} x_{ik}x_{i1} & \sum_{i=1}^{n} x_{ik}x_{i2} & \ldots & \sum_{i=1}^{n} x_{ik}^2 \end{bmatrix}$$