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]