We present the results of the Fermi -LAT 10-years-long light curves (LCs) modeling of selected blazars: six flat spectrum radio quasars (FSRQs) and five BL Lacertae (BL Lacs), examined in 7-, 10-, and 14-day binning. The LCs and power spectral densities (PSDs) were investigated with various methods: Fourier transform, Lomb-Scargle periodogram (LSP), wavelet scalogram, autoregressive moving average (ARMA) process, continuous-time ARMA (CARMA), Hurst exponent (H), and the A − T plane. First, with extensive simulations we showed that parametric modeling returns unreliable parameters, with a high dispersion for different realizations of the same stochastic model. Hence any such analysis should be supported with Monte Carlo simulations. For our blazar sample, we find that the power law indices β calculated from the Fourier and LSP modeling mostly fall in the range 1 β 2. Using the wavelet scalograms, we confirm a quasi-periodic oscillation (QPO) in PKS 2155−304 at a 3σ significance level, but do not detect any QPOs in other objects. The ARMA fits reached higher orders for 7-day binned LCs and lower orders for 10-and 14-days binned LCs for the majority of blazars, suggesting there might exist a characteristic timescale for the perturbations in the jet and/or accretion disk to die out. ARMA and CARMA modeling revealed breaks in their PSDs at timescales of few hundred days. The estimation of H was performed with several methods. We find that most blazars exhibit H > 0.5, indicating long-term memory. Finally, the FSRQ and BL Lac subclasses are clearly separated in the A − T plane.