At high angular frequencies, beyond the damping tail of the primary power spectrum, the dominant contribution to the power spectrum of cosmic microwave background (CMB) temperature fluctuations is the thermal Sunyaev-Zeldovich (tSZ) effect. We investigate various important statistical properties of the Sunyaev-Zeldovich maps, using well-motivated models for dark matter clustering to construct statistical descriptions of the tSZ effect to all orders enabling us to determine the entire probability distribution function (PDF). Any generic deterministic biasing scheme can be incorporated in our analysis and the effects of projection, biasing and the underlying density distribution can be analysed separately and transparently in this approach. We introduce the cumulant correlators as tools to analyse tSZ catalogs and relate them to corresponding statistical descriptors of the underlying density distribution. The statistics of hot spots in frequency-cleaned tSZ maps are also developed in a self-consistent way to an arbitrary order, to obtain results complementary to those found using the halo model. We also consider different beam sizes, to check the extent to which the PDF can be extracted from various observational configurations. The formalism is presented with two specific models for underlying matter clustering: (1) the hierarchical ansatz; and (2) the lognormal distribution. We find both models to be in very good agreement with the simulation results, though the lognormal model has an edge over the hierarchical model.