Statistical uncertainty is bound to occur due to the randomness in material and geometric properties, support conditions, soil variability, etc. in structural engineering problems. An attempt has been made to study the stochastic structural responses, in particular, their mean and variances under such uncertain system parameters. The random parameters are modeled as homogeneous Gaussian stochastic processes and discretized by efficient local averaging method. The discretized Gaussian field is simulated by Cholesky decomposition of the respective covariance matrix. The present paper takes the advantage of Neumann expansion technique in deriving the finite element solution of response variability within the framework of Monte Carlo simulation. Neumann expansion technique needs inversion of only the deterministic part of the stiffness matrix for all sample structures, and thus increases the computational efficiency. Numerical examples are presented to study the advantage of Neumann expansion based simulation method with respect to accuracy and efficiency. The comparison of the results shows that the values approach towards that obtained by direct Monte Carlo simulation as the order of expansion in Neumann series is increased.