<p>Factorization of the spectral density into the product of a causal and an anti-causal components plays a critical role in the frequency domain analysis of a second order time series, including linear forecasting and signal extraction. This paper generalizes the approach to higher order polyspectral densities with a view towards nonlinear filtering, and develops a general theory of polyspectral factorization, providing new mathematical results for polyspectral densities. New bijections between a restricted space of higher-dimensional cepstral coefficients (where the restrictions are induced by the symmetries of the polyspectra) and the autocumulants are derived. Applications to modeling nonlinear time series are developed; in particular, it is shown that semi-parametric nonlinear time series modeling can be accomplished by approximation of the cepstral representation of polyspectra.</p>