Gravity Recovery and Climate Experiment and its Follow On (GRACE (-FO)) missions have resulted in a paradigm shift in understanding the temporal changes in the Earth's gravity field and its drivers. To provide continuous observations to the user community, missing monthly solutions within and between GRACE (-FO) missions (33 solutions) need to be imputed. Here, we modeled GRACE (-FO) data (196 solutions) between 04/2002-04/2021 to infer missing solutions and derive uncertainties in the existing and missing observations using Bayesian inference. First, we parametrized the GRACE (-FO) time series using an additive generative model comprising long-term variability (secular trend + interannual to decadal variations), annual, and semi-annual cycles. Informative priors for each component were used and Markov Chain Monte Carlo (MCMC) was applied to generate 2,000 samples for each component to quantify the posterior distributions. Second, we reconstructed the new data (229 solutions) by joining medians of posterior distributions of all components and adding back the residuals to secure the variability of the original data. Results show that the reconstructed solutions explain 99% of the variability of the original data at the basin scale and 78% at the one-degree grid scale. The results outperform other reconstructed data in terms of accuracy relative to land surface modeling. Our data-driven approach relies only on GRACE (-FO) observations and provides a total uncertainty over GRACE (-FO) data from the data-generation process perspective. Moreover, the predictive posterior distribution can be potentially used for "nowcasting" in GRACE (-FO) near-real-time applications (e.g., data assimilations), which minimize the current mission data latency (40-60 days).