In this paper we characterize the moments of stochastic linear systems by means of the solution of a stochastic matrix equation which generalizes the classical Sylvester equation. The solution of the matrix equation is used to define the steady-state response of the system which is then exploited to define a family of stochastic reduced order models. In addition, the notions of stochastic reduced order model in the mean and stochastic reduced order model in the variance are introduced. While the determination of a reduced order model based on the stochastic notion of moment has high computational complexity, stochastic reduced order models in the mean and variance can be determined more easily, yet they preserve some of the stochastic properties of the system to be reduced. The differences between these three families of models are illustrated by means of numerical simulations.