This paper presents a constructive approach to neural network modeling of polynomial harmonic functions. This is an approach to growing higher-order networks like these build by the multilayer GMDH algorithm using activation polynomials. Two contributions for enhancement of the neural network learning are offered: (1) extending the expressive power of the network representation with another compositional scheme for combining polynomial terms and harmonics obtained analytically from the data; (2) space improving the higher-order network performance with a backpropagation algorithm for further gradient descent learning of the weights, initialized by least squares fitting during the growing phase. Empirical results show that the polynomial harmonic version phGMDH outperforms the previous GMDH, a Neurofuzzy GMDH and traditional MLP neural networks on time series modeling tasks. Applying next backpropagation training helps to achieve superior polynomial network performances.