In this paper, we develop a computational method for model reduction of polynomial dynamical systems. This is achieved using sum of squares relaxations on certain Lyapunov inequalities, which are the nonlinear counterparts of the Lyapunov controllability and observability linear matrix inequalities for linear systems. In our model reduction procedure, we use notions of balanced realization and balanced truncation for a polynomial model. In addition, we derive an a-priori error bound on the approximation error for balanced truncation.