Fitting insurance loss data can be challenging because of their non-negativity, asymmetry, skewness, and possible multi-modality. Though many parametric models have been used in the actuarial literature, these difficulties call for more flexible models for actuarial applications. In this paper, we propose a new class of gamma kernel density estimators (GKDEs) based on the gamma density. We prove that the density of the proposed model converges to that of any loss random variable which is non-negative and continuous, and establish its rate of convergence, under some technical conditions. The proposed model has several advantages over the existing gamma kernel class by Chen (2000) in that it is a valid density for any finite sample and has standard distributional quantities, such as the moments, the conditional tail moments, and the compound distribution with GKDE claim amounts, in analytic form. The model is also a competing model of the Erlang mixture by Lee and Lin (2010) in its flexibility, but with a straightforward implementation and optimization. As numerical examples, we fit the gamma kernel density estimator to actual insurance data and find that the proposed model gives adequate results compared to the Erlang mixture and the Phase-type models. (C) 2013 Elsevier B.V. All rights reserved.