Gravity models have been utilized by many applied econometricians for analyzing the effect of public policies. The traditional estimation had to log-linearized and estimated by OLS but its results get inconsistent when facing zero-valued bilateral trade flows or heteroskedasticity. Although Poisson pseudo maximum likelihood estimator gives consistent parameters when heteroskedasticity occurs but it has got a problem when zero flows arise. Therefore this paper desires to propose a simple estimation that gives consistent parameters, robust to both different forms of heteroskedasticity and zero flows. This paper used simulation to compare the performance of different estimators. From the result, I suggest that we should estimate gravity equation by two steps estimation method and use Negative binomial pseudo maximum likelihood estimator in the second step when both heteroskedasticity and zero flows occur because it gives consistent parameter, robust to different forms of heteroskedasticity and greatly deal with zero flows.