Polyspectral estimation is a problem of great importance in the analysis of nonlinear time series that has applications in biomedical signal processing, communications, geophysics, image, radar, sonar and speech processing, etc. Higher order spectra (HOS) have been used in unsupervised and supervised clustering in big data scenarios, in testing for Gaussianity, to suppress Gaussian noise, to characterize nonlinearities in time series data, and so on . Any algorithm for computing the $k$th order spectra of a time series of length $n$ needs $Ω(n^{k-1})$ time since the output size will be $Ω(n^{k-1})$ as well. Given that we live in an era of big data, $n$ could be very large. In this case, sequential algorithms might take unacceptable amounts of time. Thus it is essential to develop parallel algorithms. There is also room for improving existing sequential algorithms. In addition, parallel algorithms in the literature are nongeneric. In this paper we offer generic sequential algorithms for computing higher order spectra that are asymptotically faster than any published algorithm for HOS. Further, we offer memory efficient algorithms. We also present optimal parallel implementations of these algorithms on parallel computing models such as the PRAM and the mesh. We provide experimental results on our sequential and parallel algorithms. Our parallel implementation achieves very good speedups.