Having the following development, how to explain the transformation that the squared expected value undergoes to become a sum of covariances. I can't find anywhere in internet which is the relation between both that is made to obtain this.
$\qquad\begin{align}&\quad \mathsf E((X -2Y + Z + 1)^2) \\ &= \mathsf V [X -2Y + Z + 1] + (\mathsf E [X -2Y + Z + 1])^2 \\ &= \mathsf V [X] + (-2)^2\mathsf V[Y ] + \mathsf V [Z] +Σ\mathsf {Cov}+ (\mathsf E [X] -2\mathsf E [Y ] + \mathsf E[Z] + 1)^2\end{align}$
The definition of variance is:
$$\mathsf{V}(U)~=~\mathsf E(U^2)-\mathsf E(U)^2$$
Likewise covariance is defined as: $$\mathsf{Cov}(U,V)=\mathsf E(UV)-\mathsf E(U)\,\mathsf E(V)$$
The variance of a sum of two random variables is thus:
$$\begin{align}\mathsf V(aU+bV) &=\mathsf E((aU+bV)^2)-\mathsf E(aU+bV)^2\\&=\mathsf E(a^2U^2+2abUV+b^2V)-(a\mathsf E(U)+b\mathsf E(V))^2\\&=a^2\mathsf E(U^2)+2ab\mathsf E(UV)+b^2\mathsf E(V^2)-a^2\mathsf E(U)^2-2ab\mathsf E(U)\mathsf E(V)-b^2\mathsf E(V)^2\\&=a^2\mathsf{V}(U)+2ab\mathsf{Cov}(U,V)+b^2\mathsf{V}(V) \end{align}$$
And you can easily extend this for the sum of multiple random variables.
Thus we obtain the rule of Bilinearity of Covariance, for constants $(a_i)$ and random variables $(X_i)$:
$$\mathsf V(\sum_i a_iX_i)~=~\sum_ia_i^2\mathsf V(X_i)+2\sum\limits_{i,j~:~i < j}a_ia_j\mathsf {Cov}(X_i,X_j)$$
So, this is how we have:
$$\begin{align}&\quad~\mathsf V(X-2Y+Z+1)\\&=\mathsf V(X)+(-2)^2\mathsf V(Y)+\mathsf V(Z)+2\big((-2)\mathsf {Cov}(X,Y)+\mathsf {Cov}(X,Z)+(-2)\mathsf {Cov}(Y,Z)\big)\end{align}$$
And since $$\mathsf E((X-2Y+Z+1)^2) = \mathsf V(X-2Y+Z+1) +(\mathsf E(X-2Y+Z+1))^2$$
... then...