The Monte-Carlo method is applied to the Polyakov-loop extended
Nambu--Jona-Lasinio (PNJL) model. This leads beyond the saddle-point
approximation in a mean-field calculation and introduces fluctuations around
the mean fields. We study the impact of fluctuations on the thermodynamics of
the model, both in the case of pure gauge theory and including two quark
flavors. In the two-flavor case, we calculate the second-order Taylor expansion
coefficients of the thermodynamic grand canonical partition function with
respect to the quark chemical potential and present a comparison with
extrapolations from lattice QCD. We show that the introduction of fluctuations
produces only small changes in the behavior of the order parameters for chiral
symmetry restoration and the deconfinement transition. On the other hand, we
find that fluctuations are necessary in order to reproduce lattice data for the
flavor non-diagonal quark susceptibilities. Of particular importance are pion
fields, the contribution of which is strictly zero in the saddle point
approximation