Abstract Blazars exhibit stochastic flux variability across the electromagnetic spectrum, often exhibiting heavy-tailed flux distributions, commonly modeled as lognormal. However, Tavecchio et al. and Adams et al. found that the high-energy gamma-ray flux distributions of several of the brightest flaring Fermi-LAT flat-spectrum radio quasars (FSRQs) are well modeled by an even heavier-tailed distribution, which we show is the inverse gamma distribution. We propose an autoregressive inverse gamma variability model in which an inverse gamma flux distribution arises as a consequence of a shot-noise process. In this model, discrete bursts are individually unresolved and averaged over within time bins, as in the analysis of Fermi-LAT data. Stochastic variability on timescales longer than the time-bin duration is modeled using first-order autoregressive structure. The flux distribution becomes approximately lognormal in the limiting case of many weak bursts. The fractional variability is predicted to decrease as the time-bin duration increases. Using simulated light curves, we show that the proposed model is consistent with the typical gamma-ray variability properties of FSRQs and BL Lac objects. The model parameters can be physically interpreted as the average burst rate, the burst fluence, and the timescale of long-term stochastic fluctuations.