The field of topology optimization has progressed substantially in recent years, with applications varying in terms of the type of structures, boundary conditions, loadings, and materials. Nevertheless, topology optimization of stochastically excited structures has received relatively little attention. Most current approaches replace the dynamic loads with either equivalent static or harmonic loads. In this study, a direct approach to problem is pursued, where the excitation is modeled as a stationary zero-mean filtered white noise. The excitation model is combined with the structural model to form an augmented representation, and the stationary covariances of the structural responses of interest are obtained by solving a Lyapunov equation. An objective function of the optimization scheme is then defined in terms of these stationary covariances. A fast large-scale solver of the Lyapunov equation is implemented for sparse matrices, and an efficient adjoint method is proposed to obtain the sensitivities of the objective function. The proposed topology optimization framework is illustrated for four examples: (i) minimization of the displacement of a mass at the free end of a cantilever beam subjected to a stochastic dynamic base excitation, (ii) minimization of tip displacement of a cantilever beam subjected to a stochastic dynamic tip load, (iii) minimization of tip displacement and acceleration of a cantilever beam subjected to a stochastic dynamic tip load, and (iv) minimization of a plate subjected to multiple stochastic dynamic loads. The results presented herein demonstrate the efficacy of the proposed approach for efficient multi-objective topology optimization of stochastically excited structures, as well as multiple input-multiple output systems.