We derive matrix expressions in closed form for the spectral and bispectral densities of Markov switching bilinear models. Under suitable assumptions, we prove that the sample estimators of the spectral and bispectral density matrices are consistent and asymptotically normally distributed. Simulations and empirical applications confirm the validity of the asymptotic properties and the suitability of the proposed methods for the analysis of time series in the frequency domain.