This paper describes an original approach for computing the stationary response of linear periodic time variant MDOF systems subjected to stationary stochastic external excitation. The proposed method is derived in the frequency domain, is purely numerical, and provides the explicit power spectral density (PSD) of the response. Its implementation first requires expressing the PSD response as a function of the bilinear Fourier transform of the so-called bitemporal impulse response. Then, the spectral method is used to compute the bispectrum function. The efficiency of this spectral process is demonstrated by comparison with Monte Carlo simulations on three parametrical systems. The computational time required and the accuracy are very satisfactory.