The stochastic solution with Gaussian stationary increments is establihsed for the symmetric space-time fractional diffusion equation when $0 < β< α\le 2$, where $0 < β\le 1$ and $0 < α\le 2$ are the fractional derivation orders in time and space, respectively. This solution is provided by imposing the identity between two probability density functions resulting (i) from a new integral representation formula of the fundamental solution of the symmetric space-time fractional diffusion equation and (ii) from the product of two independent random variables. This is an alternative method with respect to previous approaches such as the scaling limit of the continuos time random walk, the parametric subordination and the subordinated Langevin equation. A new integral representation formula for the fundamental solution of the space-time fractional diffusion equation is firstly derived. It is then shown that, in the symmetric case, a stochastic solution can be obtained by a Gaussian process with stationary increments and with a random wideness scale variable distributed according to an arrangement of two extremal Lévy stable densities. This stochastic solution is self-similar with stationary increments and uniquely defined in a statistical sense by the mean and the covariance structure. Numerical simulations are carried out by choosing as Gaussian process the fractional Brownian motion. Sample paths and probability densities functions are shown to be in agreement with the fundamental solution of the symmetric space-time fractional diffusion equation.