A stochastic model for chemical reactions is presented, which represents the population of various species involved in a chemical reaction as the continuous state of a polynonmial stochastic hybrid system (pSHS). pSHSs correspond to stochastic hybrid systems with polynomial continuous vector fields, reset maps, and transition intensities. We show that for pSHSs, the dynamics of the statistical moments of its continuous states, evolves according to infinite-dimensional linear ordinary differential equations (ODEs), which can be approximated by finite-dimensional nonlinear ODEs with arbitrary precision. Based on this result, a procedure to build this types of approximation is provided. This procedure is used to construct approximate stochastic models for a variety of chemical reactions that have appeared in literature. These reactions include a simple bimolecular reaction, for which one can solve the Master equation; a decaying-dimerizing reaction set which exhibits two distinct time scales; a reaction for which the chemical rate equations have a continuum of equilibrium points; and the bistable Schogl reaction. The accuracy of the approximate models is investigated by comparing with Monte Carlo simulations or the solution to the Master equation, when available. Copyright (c) 2005 John Wiley & Sons, Ltd.