Stationarity is a cornerstone property that facilitates the analysis and processing of random signals in the time domain. Although time-varying signals are abundant in nature, in many practical scenarios the information of interest resides in more irregular graph domains. This lack of regularity hampers the generalization of the classical notion of stationarity to graph signals. The contribution in this paper is twofold. Firstly, we propose a definition of weak stationarity for random graph signals that takes into account the structure of the graph where the random process takes place, while inheriting many of the meaningful properties of the classical definition in the time domain. Our definition requires that stationary graph processes can be modeled as the output of a linear graph filter applied to a white input. We will show that this is equivalent to requiring the correlation matrix to be diagonalized by the graph Fourier transform. Secondly, we analyze the properties of the power spectral density and propose a number of methods to estimate it. We start with nonparametric approaches, including periodograms, window-based average periodograms, and filter banks. We then shift the focus to parametric approaches, discussing the estimation of moving-average (MA), autoregressive (AR) and ARMA processes. Finally, we illustrate the power spectral density estimation in synthetic and real-world graphs.