For direct statistical analysis of 3‐D structures against lateral earthquake shaking, stochastic simulation of fully nonstationary bi‐directional ground motions is essential. The recorded bi‐directional ground motions (or their projections) along the two orthogonal axes of a structure are usually correlated and their correlation influences the overall structural response. Further, the correlation is, in general, frequency‐dependent, that is, the correlation coefficients are different for different narrow frequency bands of the (composite) pair of motions. So it is important to have a method for stochastic simulation of bi‐directional motions that can maintain a target frequency‐dependent correlation. This type of simulated ensemble of bi‐directional motions for a recording process provides a sound basis for the statistical analysis of structural response vis‐à‐vis the observed/recorded correlation. In the present paper, a Priestley process based framework is proposed to simulate bi‐directional ground motions that can emulate the frequency‐dependent correlation structure of the parent recording process while maintaining the time–frequency characteristics of the recorded bi‐directional motion. Since principal directions are frequency‐dependent, it is imperative to develop a modal combination rule for multi‐directional ground motions which can address the frequency‐dependence of principal directions. For this purpose, a new ground motion frequency‐dependent modal combination rule is also proposed. It is found that the proposed frequency‐dependent combination rule performs better than the conventional frequency‐independent one, especially when the response is dominated by some frequency bands having different individual principal directions from each other or from those of the composite pair of motions.