  The stochastic solution with Gaussian stationary increments is establihsed
for the symmetric space-time fractional diffusion equation when $0 < \beta <
\alpha \le 2$, where $0 < \beta \le 1$ and $0 < \alpha \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\'evy 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.
