We derive a matrix product formula for symmetric Macdonald polynomials. Our
results are obtained by constructing polynomial solutions of deformed
Knizhnik--Zamolodchikov equations, which arise by considering representations
of the Zamolodchikov--Faddeev and Yang--Baxter algebras in terms of
t-deformed bosonic operators. These solutions form a basis of the ring of
polynomials in n variables, whose elements are indexed by compositions. For
weakly increasing compositions (anti-dominant weights), these basis elements
coincide with non-symmetric Macdonald polynomials. Our formulas imply a natural
combinatorial interpretation in terms of solvable lattice models. They also
imply that normalisations of stationary states of multi-species exclusion
processes are obtained as Macdonald polynomials at q=1.Comment: 27 pages; typos corrected, references added and some better
conventions adopted in v