The multi-factor model is a multivariate gaussian distribution?

#1
I'm sure that the following question is trivial, but I can't found a solid demonstration.

Can we state that multi-factor model is a multivariate gaussian distribution?

multi-factor model:
X_ij = sigma_i * Z_i + sqrt(1-sigma_i^2)* E_ij
where
i=1,...,k
j=1,...,n_i
0 < sigma_i < 1
Z_i ~ N(0,1) univariate Normal
Z ~ N(0,M) multivariate Normal with correlation matrix M
E_ij ~ N(0,1) iid for all i,j
Z independent from E_ij for all i,j

X_ij are marginals distributed as N(0,1), but this don't grant that the multi-factor distribution is the multivariate normal. The following numeric example using Mathematica suggest that the statement is true:

(* PDF conditional to z_i *)
F[x1_, x2_, x3_, z1_, z2_, s1_, s2_, s12_] :=
PDF[MultinormalDistribution[{0, 0}, {{1, s12}, {s12, 1}}],{z1,z2}]*
PDF[NormalDistribution[0, Sqrt[1 - s1^2]], (x1 - s1*z1)]*
PDF[NormalDistribution[0, Sqrt[1 - s1^2]], (x2 - s1*z1)]*
PDF[NormalDistribution[0, Sqrt[1 - s2^2]], (x3 - s2*z2)]

(* multi-factor PDF *)
Integrate[F[x1, x2, x3, z1, z2, s1, s2, s12], {z1, -Infinity, +Infinity},
{z2, -Infinity, +Infinity}]

(* multi-factor PDF as analytical expression (see previous output) *)
Ffactor[x1_, x2_, x3_, s1_, s2_, s12_] := -E^((((-1 + s1^2 s12^2 s2^2) x1^2 +
(-1 + s1^2 s12^2 s2^2) x2^2 -
2 s1 (-1 + s1^2) s12 s2 x2 x3 + (-1 + s1^4) x3^2 -
2 s1 x1 (s1 (-1 + s12^2 s2^2) x2 - s12 s2 x3 + s1^2 s12 s2 x3))/(
2 (-1 + s1) (1 + s1) (-1 + s1^2 (-1 + 2 s12^2 s2^2)))))/(2 Sqrt[
2] \[Pi]^(3/2) (-1 + s1^2) Sqrt[1 - s12^2] Sqrt[1 - s2^2] Sqrt[(
1 - s12^2 s2^2)/((-1 + s12^2) (-1 + s2^2))] Sqrt[(
1 + s1^2 (1 - 2 s12^2 s2^2))/((-1 + s1^2) (-1 + s12^2 s2^2))])

(* we check that it is a valid PDF -> integral = 1 *)
Integrate[
Ffactor[x1, x2, x3, s1, s2, s12], {x1, -Infinity, +Infinity},
{x2, -Infinity, +Infinity}, {x3, -Infinity, +Infinity}]

(* previous output distinct than 1, it is an expression involving s1,s2,s12 *)
(* fixed s1,s2,s12, we check that integral = 1 *)
% /. {s1->0.1, s2->0.2, s12->0.3}

(* PDF multivariate Normal *)
Fgaussian[x1_, x2_, x3_, s1_, s2_, s12_] :=
PDF[MultinormalDistribution[{0, 0, 0}, {{1, s1^2, s1*s2*s12},
{s1^2, 1, s1*s2*s12}, {s1*s2*s12, s1*s2*s12, 1}}], {x1, x2, x3}]

(* we define the difference *)
Fdiff[x1_, x2_, x3_, s1_, s2_, s12_] := Fgaussian[x1, x2, x3, s1, s2, s12] -
Ffactor[x1, x2, x3, s1, s2, s12]

(* fixed a point, we check that difference is near to 0 *)
Fdiff[1, 1, 1, 0.2, 0.3, 0.4]
 
#2
Josep Fortiana suggests to use the characterization of the multivariate Normal (see 'An Introduction to Multivariate Statistical Analysis' by T.W. Anderson, theorem 2.6.2) that states that if every linear combination of the components of a vector Y is normally distributed, then Y is normally distributed.
Thanks Josep.