Cosmological parameters encoding our current understanding of the expansion history of the Universe can be constrained by the accurate estimation of time delays arising in gravitationally lensed systems. We propose TD-CARMA, a Bayesian method to estimate cosmological time delays by modelling the observed and irregularly sampled light curves as realizations of a Continuous Auto-Regressive Moving Average (CARMA) process. Our model accounts for heteroskedastic measurement errors and microlensing, an additional source of independent extrinsic long-term variability in the source brightness. The CARMA formulation admits a linear state-space representation, that allows for efficient and scalable likelihood computation using the Kalman Filter. We obtain a sample from the joint posterior distribution of the model parameters using a nested sampling approach. This allows for "painless" Bayesian Computation, dealing with the expected multi-modality of the posterior distribution in a straightforward manner and not requiring the specification of starting values or an initial guess for the time delay, unlike existing methods. In addition, the proposed sampling procedure automatically evaluates the Bayesian evidence, allowing us to perform principled Bayesian model selection. TD-CARMA is parsimonious, and typically includes no more than a dozen unknown parameters. We apply TD-CARMA to three doubly lensed quasars HS 2209+1914, SDSS J1001+5027 and SDSS J1206+4332, estimating their time delays as −21.96±1.448 (6.6% precision), 120.93±1.015 (0.8%), and 111.51±1.452 (1.3%), respectively. A python package, TD-CARMA, is publicly available to implement the proposed method.