Echo　state　networks　(ESNs)　have　been　widely　used　in　the　field　of　time　series　prediction.　In　conventional　ESNs,　the　spectral　radius　of　reservoir　is　always　scaled　to　lower　than　1　to　satisfy　the　necessary　condition　for　echo　state　property　(ESP),　while　the　sufficient　condition　is　unregarded.　Meanwhile,　the　output　weights　are　always　trained　without　considering　the　high　order　statistics　of　input　signals.　To　solve　above　problems,　the　original　ESN　is　extended　to　polynomial　ESNs　(PESNs)　by　employing　the　polynomial　functions　of　input　variables　into　output　weights.　Firstly,　the　reservoir　of　the　PESN　is　built　by　the　singular　value　decomposition　(SVD)　method.　Secondly,　the　prime　PESN　(P-PESN)　is　implemented　and　its　polynomial　output　weights　are　augmented　by　the　high　order　statistics　of　inputs.　Thirdly,　the　simplified　PESN　(S-PESN)　is　constructed　by　decomposing　the　polynomial　output　weights　of　the　P-PESN　into　randomly　generated　polynomial　nodes　and　tuned　output　weights.　Furthermore,　the　regression　matrix　properties　of　P-PESN　and　S-PESN　are　theoretically　analyzed,　respectively.　Finally,　the　simulation　results　show　that　the　proposed　PESNs　obtain　better　performance　than　other　methods　in　terms　of　prediction　accuracy　and　learning　speed.　(c)　2018　Elsevier　B.V.　All　rights　reserved.