Well, I've been reading about simulating correlated data and I've come across Cholesky decomposition. Everything seemed clear until I found a couple of posts on this site and Cross-Validated that showed a way to alter mean and variance of simulated data. The proposed solution is as follows:
Let $Z$ be a set of uncorrelated random variables normally distributed with mean 0 and variance 1, i.e.
$$Z \sim N(0, I)$$ Then if we make an affine transformation
$$X \equiv A + BZ $$ $X$ will have a distribution $$ X \sim N(A, B{B}^{T}) $$ So given a covariance matrix $\Sigma$ we can find $B$ using the Cholesky decomposition $ \Sigma = B{B}^{T}$.
So, I don't feel like I understand why the application of an affine transformation of form $X \equiv A + BZ $ results in $ X \sim N(A, B{B}^{T}) $ instead of $X \sim N(A, B)$.
Since $X$ is a vector
$$\mathbb{E}X=\mathbb{E}[A+BZ]=A$$
and
$$Var(X)=\mathbb{E}[(X-\mathbb{E}X)(X-\mathbb{E}X)']=\mathbb{E}[BZ(BZ)']$$ $$=\mathbb{E}[BZZ'B']=B\mathbb{E}[ZZ']B'=BIB'=BB'$$