SummaryA system of nonlinear stochastic equations excited by a Gaussian random term leading to a statistically stationary solution is considered. The Carleman linearization is used to handle the nonlinearity and the statistical characterization of the solution is formulated in terms of a sequence of correlations of increasing order. Truncated in a finite but arbitrary order, the problem leads to a linear system of equations yielding directly the correlations. It is demonstrated in two examples that, if the excitation is not too strong, the solution converges as a function of the truncation order and provides an alternative to the Monte Carlo approach consisting in ensemble averaging a large number of time-dependent random solutions. For a low-dimensional system, it is shown that replacing the tensor indexing by a numbering accounting for redundancies makes it possible to keep the total problem dimension within reasonable limits even if relatively high-order correlations are accounted for. A model reduction based on the singular value decomposition of the second-order correlation matrix is tested with success for a case of a partial differential equation (Burgers’ equation), showing that the method can be potentially applied even to high-dimensional systems originating in partial differential equations.