Context. Monthly binned γ-ray light curves of 236 bright γ-ray sources, particularly blazars, selected from a sample of 2278 high-galactic latitude objects observed with Fermi-LAT show flux variability characterized by power spectral densities consisting of a single power-law component, ranging from Brownian to white noise. Aims. The main goal here is to assess the Ornstein-Uhlenbeck (OU) model by studying the range of its three parameters that reproduces these statistical properties. Methods. We develop procedures for extracting values of the three OU model parameters (mean flux, correlation length, and random amplitude) from time series data and apply them to compare numerical integrations of the OU process with the Fermi-LAT data. Results. The OU process fully describes the statistical properties of the flux variations of the 236 blazars. The distributions of the extracted OU parameters are narrowly peaked around well-defined values (σ, μ, θ) = (0.2, −8.4, 0.5) with variances (0.004, 0.07, 0.13). The distributions of rise and the decay time scales of flares in the numerical simulations, meaning major flux variations fulfilling pre-defined criteria, are in agreement with the observed ones. The power spectral densities of the synthetic light curves are statistically indistinguishable from those of the measured light curves. Conclusions. The long-term γ-ray flux variability of blazars on monthly time scales is well described by a stochastic model that involves only three parameters. The methods described here are powerful tools for studying randomness in light curves and thereby for constraining the physical mechanisms responsible for the observed flux variations.