Implications of $X \sim \mathcal{N}(\mu, \Sigma)$ Suppose I have a random vector $X
\equiv (X_1,\ldots, X_M)$ of dimension $M\times 1$ distributed as $
\mathcal{N}(\mu, \Sigma)
$, 
where
$$
\mu\equiv(0,0,\ldots,0)
$$
and
$$
\Sigma\equiv I_M\times \rho
$$
with $\rho\in (-1,1)$ and $I_M$ is the identity matrix of dimension $M\times M$.
Can we derive the distribution of the vector 
$$
\begin{pmatrix}
X_1-X_2\\
X_1-X_3\\
\vdots\\
X_1-X_M
\end{pmatrix} \text{ ?}
$$
When $M=2$ the distribution is $\mathcal{N}(0, 2(1-\rho))$ as explained here. Can we generalise to any $M$?
 A: If the vector $\mathbf X$ has a multivariate Gaussian distribution $\mathcal N(\mathbf \mu, \mathbf \Sigma)$ and the random variable $\mathbf Y$ is defined as $\mathbf Y = \mathbf M \mathbf X$, where $\mathbf M$ is a constant matrix, then $\mathbf Y$ has a multivariate Gaussian distribution $\mathcal N(\mathbf M \mathbf \mu, \mathbf M \mathbf \Sigma \mathbf M^T)$. (See this Wikipedia page.)
In your situation, $\mathbf Y = (X_1 - X_2, \dots, X_1 - X_M)$ is related to $\mathbf X = (X_1, X_2, \dots, X_M)$ by the equation $\mathbf Y = \mathbf M \mathbf X$, where 
$$ \mathbf M = \begin{bmatrix} 1 & - 1 & 0 & \dots & 0 \\ 1& 0 & - 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots& \vdots\\1 & 0 & 0 & \dots & -1\end{bmatrix}$$
Judging from the answer you linked, I suspect you want
$$ \mathbf \mu = \begin{bmatrix} \mu \\ \mu \\ \vdots \\ \mu \end{bmatrix}, \ \ \ \ \ \mathbf \Sigma = \begin{bmatrix} 1 & \rho &  \dots & \rho \\ \rho & 1 & \dots & \rho\\ \vdots & \vdots & \ddots & \vdots \\ \rho & \rho & \dots & 1\end{bmatrix}$$
So by the result above, $\mathbf Y$ has a multivariate normal distribution, with mean
$$ \mathbf M \mathbf \mu = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0\end{bmatrix}$$
and covariance
$$ \mathbf M \mathbf \Sigma \mathbf M^T = \begin{bmatrix} 2(1- \rho) & 1 - \rho & \dots & 1 - \rho \\ 1 - \rho & 2(1 - \rho) & \dots & 1 - \rho \\ \vdots & \vdots & \ddots & \vdots \\ 1 - \rho & 1 - \rho & \dots & 2(1 - \rho)\end{bmatrix} $$
[And by the way, in the answer you linked for the $M = 2$ case,  $\Sigma$ is the matrix $\begin{bmatrix} 1 & \rho \\ \rho  & 1\end{bmatrix}$; it is not $I \times \rho$. The constant $\rho$ can be interpreted as the correlation between $X_1$ and $X_2$.]
A: $\newcommand{\c}{\operatorname{cov}}\newcommand{\v}{\operatorname{var}}$It is clear that each component of the vector in question has expectation $0.$
\begin{align}
& \c(X_1-X_i, X_1-X_j) \\[10pt]
= {} & \c(X_1,X_1) - \c(X_1,X_j) - \c(X_i,X_1) + \c(X_i,X_j) \\[10pt]
= {} & \rho - 0 - 0 + 0. \\[10pt]
& \v(X_1-X_i) = \v(X_1) + (-1)^2\v(X_i) = 2\rho.
\end{align}
Beyond the fact that this is multivariate normal one need only know the expectations and covariances.
