Line intensity mapping (LIM) is an emerging technique with a unique potential to probe a wide range of scales and redshifts. Realizing the full potential of LIM, however, relies on accurate modeling of the signal. We introduce an extended halo model for the power spectrum of intensity fluctuations of CO rotational lines and [CII] fine transition line in real space, modeling nonlinearities in matter fluctuations and biasing relation between the line intensity fluctuations and the underlying dark matter distribution. We also compute the stochastic contributions beyond the Poisson approximation using the halo model framework. To establish the accuracy of the model, we create the first cosmological-scale simulations of CO and [CII] intensity maps, MithraLIMSims, at redshifts 0.5 ≤ z ≤ 6, using halo catalogs from Hidden-Valley simulations, and painting halos according to mass-redshift-luminosity relations for each line. We show that at z = 1 on scales k max 0.8 Mpc −1 h, the model predictions of clustering power (with only two free parameters) are in agreement with the measured power spectrum at better than 5%. At higher redshift of z = 4.5, this remarkable agreement extends to smaller scale of k max 2 Mpc −1 h. Furthermore, we show that on large scales, the stochastic contributions to CO and CII power spectra are non-Poissonian, with amplitudes reproduced reasonably well by the halo model prescription. Lastly, we assess the performance of the theoretical model of the baryon acoustic oscillations (BAO) and show that hypothetical LIM surveys probing CO lines at z = 1, that can be deployed within this decade, will be able to make a high significance measurement of the BAO. On a longer time scale, a space-based mission probing [CII] line can uniquely measure the BAO on a wide range of redshifts at an unprecedented precision.