TD-CARMA: Painless, Accurate, and Scalable Estimates of Gravitational Lens Time Delays with Flexible CARMA Processes

The data for a doubly lensed quasar observed for n epochs can be written as ${\boldsymbol{D}}=\{{t}_{i},{x}_{i},{\delta }_{i}^{x},{y}_{i},{\delta }_{i}^{y}\}{}_{i=1}^{n}$, where t_i denotes the measurement time (in days) of observation epoch i, x_i and y_i denote the observed brightness of the light curves of the two lensed images of the same quasar at time t_i (typically measured in apparent magnitude), and ${\delta }_{i}^{x}$ and ${\delta }_{i}^{y}$ respectively denote the standard deviation of the measurement error on x_i and y_i (thus encoding the heteroskedasticity of measurement errors). Table 1 summarizes the symbols used in this work. Figure 1 shows the data for the doubly lensed quasar SDSS J1001+5027, where light curve B is delayed from light curve A, for Δ ≈ 119.3 days according to previous analyses (a figure reproduced from Kumar et al. 2013).

Zoom In Zoom Out Reset image size

Table 1. Glossary of the Notation Used in This Work

Data
Symbol	Description
t	Observation times
x	Brightness measurements for light curve A
y	Brightness measurements for light curve B
z	Composite light curve
δ x	Standard deviations for light curve A
δ y	Standard deviations for light curve B
Model Parameters
Symbol	Description
Δ	Time delay (days)
θ	Microlensing polynomial regression coefficients
α	CARMA auto-regressive coefficients
β	CARMA moving average coefficients
τ	DRW mean-reversion timescale
σ2	CARMA and DRW white noise variance
μ	CARMA and DRW long-term mean
Ω	Full set of either DRW or CARMA parameters
Model Hyperparameters a
Symbol	Description
p	Auto-regressive order
q	Moving average order
m	Microlensing polynomial regression order
Other Symbols
Symbol	Description
${ \mathcal Z }$	Bayesian evidence
${{ \mathcal Z }}_{i}$	Bayesian evidence for Model i or Mode i
πi	Relative probability of Model i or Mode i
ϕ	Generic parameter vector
Φ	Range of possible values for generic parameter ϕ

Note.

^a The hyperparameters index competing models. Values for the hyperparameters are selected via the Bayesian evidence or the relative probabilities of the models considered; see Section 3.1.

Download table as:  ASCIITypeset image

The general idea of our time delay estimation framework is to reconstruct the common intrinsic quasar light curve from the measurements of its two lensed counterparts, under the assumption that they are time- and magnitude-shifted versions of one another, up to a microlensing term, where the time shift corresponds to the time delay Δ.

The magnitude measurements ${\boldsymbol{x}}={\{{x}_{i}\}}_{i=1}^{n}$ and ${\boldsymbol{y}}={\{{y}_{i}\}}_{i=1}^{n}$ are assumed to be discrete realizations of unobserved continuous light curves X(t) and Y(t), which represent the true source magnitudes at time $t\in {\mathbb{R}}$. As a result of the strong gravitational lensing phenomenon, we assume that Y(t) is a time-shifted version of X(t), which yields

$$\begin{eqnarray}&&Y(t)=X(t-{\rm{\Delta }}),\end{eqnarray} \tag{ 1 }$$

where Δ is the time delay in days. In addition to the intrinsic brightness fluctuations of the source, two lensing effects induce extra variability in the observed magnitudes. First, strong gravitational lensing magnifies the brightness of each light curve to different degrees. Therefore x and y can exhibit different average magnitudes. Tak et al. (2017) take this effect into account by introducing an intercept term θ₀ as follows

$$\begin{eqnarray}&&Y(t)=X(t-{\rm{\Delta }})+{\theta }_{0}.\end{eqnarray} \tag{ 2 }$$

Second, additional extrinsic long-term variability might be caused by microlensing, a lensing effect occurring when a light ray passes close to a moving object near the lensing galaxy. This effect is typically modeled by a polynomial regression of order m on time (Kochanek et al. 2006; Courbin et al. 2011; Morgan et al. 2012; Tak et al. 2017). That is

$$\begin{eqnarray}&&Y(t)=X(t-{\rm{\Delta }})+{{\boldsymbol{w}}}_{m}(t-{\rm{\Delta }}){\boldsymbol{\theta }},\end{eqnarray} \tag{ 3 }$$

where ${{\boldsymbol{w}}}_{m}(t-{\rm{\Delta }})=\{1,t-{\rm{\Delta }},...,\,{\left(t-{\rm{\Delta }}\right)}^{m}\}$ is the vector containing the polynomial time variables and θ = {θ₀,..., θ_m } is the microlensing (regression) coefficients. We treat m as a hyperparameter. That is, it indexes the choice of model and its value is set by comparing several values as part of the statistical model selection, see Section 3.1.

Given the parameters Δ and θ , we can construct a composite light curve that combines two lensed light curves (as if they were a single light curve) by adjusting one light curve according to the model in Equation (3). The combined observation times are composed of the original observation times and the time-delay-shifted observation times, ${{\boldsymbol{t}}}^{{\rm{\Delta }}}={\{{t}_{i}\}}_{i=1}^{n}\cup {\{{t}_{i}-{\rm{\Delta }}\}}_{i=1}^{n}$. We denote the magnitudes of the composite light curve by ${\boldsymbol{z}}={\{{z}_{j}\}}_{j=1}^{2n}$, where each z_j is defined as

$$\begin{eqnarray}\begin{array}{c}{z}_{j}=\left\{\begin{array}{ll}{x}_{i} & {\rm{f}}{\rm{o}}{\rm{r}}\,{\rm{s}}{\rm{o}}{\rm{m}}{\rm{e}}\,i\,{\rm{i}}{\rm{f}}\,{t}_{j}^{{\rm{\Delta }}}\,{\rm{i}}{\rm{s}}\,{\rm{i}}{\rm{n}}\,{\boldsymbol{t}},\\ {y}_{i}-{{\boldsymbol{w}}}_{m}({t}_{j}-{\rm{\Delta }}){\boldsymbol{\theta }} & {\rm{f}}{\rm{o}}{\rm{r}}\,{\rm{s}}{\rm{o}}{\rm{m}}{\rm{e}}\,i\,{\rm{i}}{\rm{f}}\,{t}_{j}^{{\rm{\Delta }}}\,{\rm{i}}{\rm{s}}\,{\rm{i}}{\rm{n}}\,{\boldsymbol{t}}-{\rm{\Delta }}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 4 }$$

We can similarly define the vector of measurement error standard deviations as ${{\boldsymbol{\delta }}}^{z}={\left\{{\delta }_{i}^{z}\right\}}_{i=1}^{2n}$ with

$$\begin{eqnarray}{\delta }_{j}^{z}=\left\{\begin{array}{ll}{\delta }_{i}^{x} & \mathrm{for}\ \mathrm{some}\ i\ \mathrm{if}\ {t}_{j}^{{\rm{\Delta }}}\ \mathrm{is}\ \mathrm{in}\ {\boldsymbol{t}},\\ {\delta }_{i}^{y} & \mathrm{for}\ \mathrm{some}\ i\ \mathrm{if}\ {t}_{j}^{{\rm{\Delta }}}\ \mathrm{is}\ \mathrm{in}\ {\boldsymbol{t}}-{\rm{\Delta }}.\end{array}\right.\end{eqnarray} \tag{ 5 }$$

We assume that the vector of measurements for the composite light curve z is the discrete realization of an unobserved continuous-time light curve Z(t), which represents the true source magnitude of the AGN at time $t\in {\mathbb{R}}$. Therefore, given Δ and θ , we only need to model the stochastic variability of the single latent light curve Z(t). The DRW and CARMA processes have been widely used to model the stochastic nature of the intrinsic variability in AGN light curves since Kelly et al. (2009, 2014) introduced them to the astronomical community.

The DRW process has already been integrated into time delay estimation methodology by Tak et al. (2017), and we propose an extension of this method by modeling the intrinsic light curve Z(t) with the more general CARMA processes.

The DRW process is defined as the solution of the stochastic differential equation

$$\begin{eqnarray}&&{dZ}(t)=-\displaystyle \frac{1}{\tau }(Z(t)-\mu ){dt}+\epsilon (t),\end{eqnarray} \tag{ 6 }$$

where μ is the long-term mean of the process, τ is the timescale for the process to revert to its long-term mean, and (t) is a white noise process with variance σ² (that is, (t) ∼ N(0, σ²)) governing the stochastic fluctuations in the process (Brockwell & Davis 2002).

A CARMA process is indexed by two hyperparameters: its auto-regressive order p and moving average order q. The process is denoted CARMA(p, q) and is defined as the solution of the stochastic differential equation

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{d}^{p}Z(t)}{{{dt}}^{p}}+{\alpha }_{p-1}\displaystyle \frac{{d}^{p-1}Z(t)}{{{dt}}^{p-1}}+...+{\alpha }_{0}Z(t)\\ \,\,=\,\mu {dt}+{\beta }_{q}\displaystyle \frac{{d}^{q}\epsilon (t)}{{{dt}}^{q}}+{\beta }_{q-1}\displaystyle \frac{{d}^{q-1}\epsilon (t)}{{{dt}}^{q-1}}+...+\epsilon (t),\end{array}\end{eqnarray} \tag{ 7 }$$

where α = {α₀,..., α_p−1} are the auto-regressive parameters, β = {β₁,..., β_q } are the moving average coefficients, and (t) is a white noise process with variance σ². We also define α_p = 1 and β₀ = 1, and note that the DRW is equivalent to the CARMA(1,0) process with τ = 1/α₀.

DRW has a limited applicability in describing the stochastic variability in AGN light curves, since for instance it cannot account for quasiperiodic variability features and multiple auto-correlation timescales. CARMA processes, however, can account for such features and are thus suitable to characterize the variability in a wide range of AGN light curves. The reader is referred to Appendix A, Kelly et al. (2009), and Kelly et al. (2014) for technical reviews of the mathematics governing the dynamics of the DRW and CARMA processes.

DRW and CARMA processes that are stochastically driven by white noise processes are Gaussian (Kelly et al. 2009, 2014) ⁶ and their likelihood functions can be evaluated via a GP. Let Ω be the parameter vector for the intrinsic light curve Z(t), i.e., Ω = (μ, σ², τ) for DRW and Ω = (μ, σ², α , β ) for CARMA. GPs are formulated in terms of a mean function M(t_i ) (parametrized by μ, θ , and Δ in our case), and a covariance or kernel function k(t_i , t_j ) (parametrized by σ² and τ for DRW and by α and β for CARMA). Following the notation in Foreman-Mackey et al. (2017), we can write the log-likelihood function of the model parameters given observed data as

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{ln}L({\rm{\Delta }},{\boldsymbol{\theta }},{\rm{\Omega }}| {\boldsymbol{z}})=\mathrm{ln}p({\boldsymbol{z}}| {\rm{\Delta }},{\boldsymbol{\theta }},{\rm{\Omega }})\\ \quad =\,-\displaystyle \frac{1}{2}{{\bf{r}}}^{T}{K}^{-1}{\bf{r}}-\displaystyle \frac{1}{2}\mathrm{lndet}(K)-N\mathrm{ln}(2\pi ),\end{array}\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&{\bf{r}}={\left({z}_{1}-M({t}_{1}^{{\rm{\Delta }}}),...,\,{z}_{2n}-M({t}_{2n}^{{\rm{\Delta }}})\right)}^{T},\end{eqnarray} \tag{ 9 }$$

and the elements of the covariance matrix K are ${[K]}_{{ij}}=k({t}_{i}^{{\rm{\Delta }}},{t}_{j}^{{\rm{\Delta }}})=R({t}_{i}^{{\rm{\Delta }}}-{t}_{j}^{{\rm{\Delta }}})$, with R(· ) defined in Equation (A5) of Appendix A.

The evaluation of a GP likelihood function given a realization of n data points requires the inversion of matrix K, which computationally scales as O(n³). This represents a prohibitive bottleneck that limits the computational scalability to large data sets. However, Foreman-Mackey et al. (2017) developed an algorithm exploiting the semiseparable structure of CARMA covariance matrices, enabling fast computation of the likelihood function with only linear complexity O(n). The likelihood function of a CARMA process may also be computed using the Kalman filter algorithm (Jones 1981; Jones & Ackerson 1990; Kelly et al. 2014), and although both methods have the same computational time complexity, the GP-based algorithm has been found to be an order of magnitude faster in practice (Foreman-Mackey et al. 2017). We use the EzTao package (Yu & Richards 2022; Yu et al. 2022) for the GP-based computation of our model likelihood. A python package, TD-CARMA, is publicly available to implement our proposed method.

The reader is referred to Williams & Rasmussen (2006), Foreman-Mackey et al. (2017), and Yu et al. (2022) for reviews on the GPs and GP-based computations for CARMA processes.

We adopt a Bayesian statistical approach to perform parameter inference, which allows us to quantify the uncertainties in the model parameters, including Δ, via their joint posterior distribution, given the observed data. From this we can compute point estimates and error bars, as required. Once the prior distributions are specified and the data are observed (and incorporated into the model via the likelihood function), Bayes' Theorem gives the expression for the posterior distribution of the parameters, i.e., the updated summary of information about the model parameters after observing the data

$$\begin{eqnarray}&&p({\rm{\Delta }},{\boldsymbol{\theta }},{\rm{\Omega }}| {\boldsymbol{D}})=\displaystyle \frac{L({\rm{\Delta }},{\boldsymbol{\theta }},{\rm{\Omega }}| {\boldsymbol{D}})p({\rm{\Delta }},{\boldsymbol{\theta }},{\rm{\Omega }})}{p({\boldsymbol{D}})},\end{eqnarray} \tag{ 10 }$$

where L(Δ, θ , Ω∣ D ) is the likelihood under the model defined in Equation (8) and is equal to p( D ∣Δ, θ , Ω), the sampling distribution of the data given the model parameters; p(Δ, θ , Ω) is the prior distribution of the parameters, and p( D ) is the marginal distribution of the data, otherwise known as the Bayesian evidence. The Bayesian evidence, hereafter denoted by ${ \mathcal Z }$, i.e., ${ \mathcal Z }=p({\boldsymbol{D}})$, can be used for model comparisons, and is defined as

$$\begin{eqnarray}&&{ \mathcal Z }={\int }_{{\rm{\Phi }}}p({\boldsymbol{D}}| \phi )p(\phi )d\phi ,\end{eqnarray} \tag{ 11 }$$

where ϕ is a generic (vector) parameter and Φ is the (multivariate) set of possible values for ϕ. The Bayesian evidence ${ \mathcal Z }$ is a measure of the "goodness of fit" of a model to the data D , so that models with higher evidence provide a better explanation of the data. The definition of ${ \mathcal Z }$ naturally incorporates Occam's razor, as more complex models, i.e., models defined on higher-dimensional parameter spaces, are penalized if they do not sufficiently improve the fit to the data (Jefferys & Berger 1991).

Suppose we fit each of several models, M₁,..., M_K and wish to identify the models that dominate in terms of their evidence. For example, we may wish to compare TD-DRW(m) and TD-CARMA(p, q, m) for each of several combinations of values of the hyperparameters, p, q, and m (auto-regressive, moving average, and microlensing orders).

For this, it is useful to compute the relative probability of each considered model

$$\begin{eqnarray}&&p({M}_{i})=\displaystyle \frac{{{ \mathcal Z }}_{i}}{{\sum }_{k=1}^{K}{{ \mathcal Z }}_{k}}.\end{eqnarray} \tag{ 12 }$$

In TD-CARMA we must specify the hyperparameter triplet (p, q, m) that best fits the data. Statistically, this is a model selection problem. Because CARMA(p, q) models are nonnested, likelihood ratio tests cannot be used for model selection (Protassov et al. 2002; Kelly et al. 2014). In the astrophysics literature, when fitting CARMA processes to light-curve data, (p, q) pairs are selected using information criteria such as the Akaike information criterion (AIC; Kelly et al. 2014) or deviance information criterion (DIC; Kasliwal et al. 2017). These may require additional computation (e.g., AIC requires the computation of the maximum likelihood estimator) and rely on asymptotic approximations. We instead use the principled Bayesian evidence, ${ \mathcal Z }$, which is a commonly used model selection approach in the Bayesian regime (mostly via Bayes' factors; e.g., Kass & Raftery 1995). The Bayesian evidence is conveniently computed directly by MultiNest ⁷ at the same time it generates a posterior sample, see Section 3.3.

We assume uniform prior distributions for each model parameter. Specifically, we set a uniform prior for the time delay Δ in the range [t₁ − t_n , t_n − t₁]. In instances where it is clear which light curve is delayed, the prior distribution for Δ can accordingly be reduced to [t₁ − t_n , 0] or [0, t_n − t₁]. The prior for the microlensing parameters, θ , is defined to be uniform on the range [−10⁵, 10⁵]. The average apparent magnitude μ is given a uniform prior on the range [−30, 30], which covers the magnitude range from that of the Sun to that of the faintest object observable by the Hubble Space Telescope (Tak et al. 2017). The other TD-CARMA parameters α , β , and σ are sampled on the logarithmic scale, with priors on their (natural) logs defined to be uniform on the range [−15, 15]. This largely covers the range of frequencies probed by the data. ⁸

We obtain a sample from the posterior distribution p(Δ, β , Ω∣ D ) in Equation (10) via the MultiNest implementation of (multimodal ellipsoidal) nested sampling (Feroz et al. 2009; Buchner et al. 2014), for both TD-DRW and TD-CARMA. In the case of time delay estimation, inference via nested sampling has a variety of advantages over traditional MCMC methods, including the computation of the Bayesian evidence to perform model selection (see Section 3.1), and the efficient sampling of multimodal posterior distributions.

The multimodality in this problem is complicated, because it manifests itself in two distinct ways.

First, the posterior distribution of CARMA(p, q) parameters is notoriously multimodal when p > 1 (Kelly et al. 2014). Moreover, even for fixed p and q there is potential multimodality in the time delay parameters Δ and θ (e.g., Tak et al.2017).

The sampling procedures in Tak et al. (2017) have difficulties navigating such posteriors even with the simpler TD-DRW model. Thus, we use MultiNest, which is specifically designed to sample efficiently from complicated and multimodal posterior distributions. MultiNest identifies the multiple modes in the posterior, partitions the parameter space into regions in which the separated modes are supported, and evaluates the local evidence of each identified mode.

Specifically, assuming MultiNest identifies mode i defined in region Φ_i of the parameter space Φ, the local evidence ${{ \mathcal Z }}_{i}$ of mode i is defined as

$$\begin{eqnarray}&&{{ \mathcal Z }}_{i}={\int }_{{{\rm{\Phi }}}_{i}}p({\boldsymbol{D}}| \phi )p(\phi )d\phi .\end{eqnarray} \tag{ 13 }$$

The local evidence ${{ \mathcal Z }}_{i}$ can be used to compute the relative probability of mode i, denoted as π_i and defined as

$$\begin{eqnarray*}&&{\pi }_{i}={\int }_{{{\rm{\Phi }}}_{i}}p(\phi | {\boldsymbol{D}})d\phi =\displaystyle \frac{1}{{ \mathcal Z }}{\int }_{{{\rm{\Phi }}}_{i}}p({\boldsymbol{D}}| \phi )p(\phi )d\phi =\displaystyle \frac{{{ \mathcal Z }}_{i}}{{ \mathcal Z }},\end{eqnarray*}$$

where ${ \mathcal Z }$ is the (total) Bayesian evidence.

In the event of a multimodal posterior distribution, computing the above quantities offers a principled way of quantifying the relative probability of the several modes.

As discussed in Section A.2, we use an alternative parameterization of the auto-regressive polynomial A(z) to enforce stationarity when fitting CARMA processes. Instead of sampling the auto-regressive parameters α , we sample the coefficients a = {a₁,..., a_p } (defined in Section A.2), on the natural log-scale to enforce nonnegativity (for stationarity of the resulting CARMA process).

The likelihood of a CARMA process is invariant to permutations in the indices of the roots of the auto-regressive polynomial. In order to avoid identifiability issues, Kelly et al. (2014) imposes ordering constraints on the auto-regressive roots. This, however, involves rejecting Monte Carlo draws that do not satisfy the ordering constraint, which renders sampling inefficient. Since swaps in the indices of the auto-regressive roots that give the same likelihood create mirror modes in the posterior distribution, we do not impose an ordering constraint but rather combine the modes and consider the full posterior in our analyses.

Finally, the polynomial regression coefficients θ = {θ₁,..., θ_m } are highly correlated in their posterior distributions, since the regression covariates $\{1,t-{\rm{\Delta }},{\left(t-{\rm{\Delta }}\right)}^{2},...,\,{(t-{\rm{\Delta }})}^{m}\}$ are themselves highly correlated. This yields a posterior distribution with complex geometry that is difficult to sample from efficiently. To circumvent this issue, we perform a QR decomposition ⁹ (e.g., Gander 1980) of the polynomial regression covariate matrix, a common computational trick in standard regression problems used to decorrelate covariates.

At each iteration, MultiNest proposes new values of the parameters by sampling from their joint prior distribution, constrained to regions of the parameter space corresponding to isocontours in the likelihood (Feroz et al. 2009). The resolution at which the parameter space is explored therefore depends (at least initially) on the number of points sampled from the prior at each iteration, denoted N_live (Feroz et al. 2009), and the range covered by the prior. We repeatedly fitted our models with different values of N_live, and found that N_live = 400 for our simulations and N_live = 1000 for our data set applications yield stable results.

In the event of a multimodal posterior distribution of the CARMA(p, q) model parameters, Ω = { α , β , σ}, one or more of the modes might correspond to a quasiperiodic frequency (QPO) in the posterior PSD (i.e., the PSD features a Lorentzian centered away from 0). For example, the left panel of Figure 2 shows a bimodal PSD, where mode 1 detects a frequency at ≈2 days. This posterior PSD was obtained by fitting TD-CARMA(3, 1, 3) to the doubly lensed quasar data set SDSS J1001+5027, the full analysis of which appears in Section 5.3.

Zoom In Zoom Out Reset image size

Following Kelly et al. (2014), we can attribute a detected frequency in the inferred PSD of the intrinsic light curve to either the signal or the noise by comparing the posterior distribution of the power for the detected frequency with the measurement noise level. ¹⁰ If the posterior distribution of power for the detected frequency is supported above the measurement noise level, we recognize the frequency as a feature of the light curve. Otherwise, this provides statistical evidence that the CARMA processes are essentially identifying a periodic artifact of the noise.

If a detected frequency were to be attributed to the intrinsic light curve, it would dramatically reduce the size of the error bars for the estimate of Δ, as shown in the right panel of Figure 2 (and Tables 12, 14, and 15 in Appendix C). Indeed, when a frequency is detected in the PSD for the intrinsic light curve, this naturally constrains Δ to multiples of the frequency. Thus, TD-CARMA can leverage a frequency that is not attributed to the noise to reduce the error bars on Δ dramatically. This would not be possible with TD-DRW, for instance.

We simulate doubly lensed synthetic data sets under nine different data generation models. The nine models were specified by three different continuous time series processes: DRW (i.e., CARMA(1, 0)), CARMA(2, 1), and CARMA(3, 1) crossed with three different orders of polynomial microlensing: m = 1, 2, 3. For each of the nine resulting data generating processes for the microlensed light curves, 300 replicate data sets are simulated, sharing all model parameters Ω related to the latent process. This yields 2700 data sets in total. In the next few paragraphs, we detail our procedure for the simulation of the synthetic data sets which, in an effort to mimic potential light curves from the Rubin Observatory Legacy Survey of Space and Time, ¹¹ roughly follows the simulation design presented in Liao et al. (2015).

All 2700 data sets share the same generative values for Δ, θ, and μ. We generate the shared value of Δ by sampling from a uniform distribution on [0, 200] days, and sample θ and μ from their priors, defined in Section 3.2. We then generate the time series parameters Ω for the three different intrinsic light-curve models DRW, CARMA(2, 1), and CARMA(3, 1). When an intrinsic light curve is generated under a DRW process, we set the timescale τ to 1000 days. For the CARMA(2, 1), we choose α such that the auto-correlation function (ACF) of the process features an exponentially damped sinusoid with a frequency equal to 1/100 days and an e-folding timescale of 1/500 days. For the CARMA(3, 1) process, we modify the α generated for the CARMA(2, 1) process to add an exponential decay with an e-folding timescale of 1/1000 days to the ACF of the CARMA(2, 1) process. ¹² The moving average parameter vector β is sampled uniformly from its prior, defined in Section 3.2. Finally, the σ² parameter is computed such that all three intrinsic light curves share the same marginal variance, R(0); expressions for R(0) are given in Equations (A1) and (A5) in Section A.2. We sample this shared value of R(0) uniformly from the [0, 10] range. Figure 7 shows the ACFs (defined in Sections A.1 and A.2) for the three intrinsic light-curve models resulting from this procedure.

Given the generative values of the parameters Ω, Δ, and θ , we create a doubly lensed synthetic data set as follows.

• 1.  
  Given Ω, generate an observation of the intrinsic light curve Z(t) in a high-resolution grid of measurement times, with 100,000 evenly spaced epochs over a 1000 day period.
• 2.  
  Given the values of Δ and θ , generate high-resolution observations of the light curves X(t) and Y(t), via Equation (3).
• 3.  
  Downsample the data to five observational seasons of length four months, and to a median observation cadence of three days. This yields approximately 200 data points for each light curve. The seasonal thinning is the same across all data set replicates, but the daily thinning is random.
• 4.  
  Generate measurement error standard deviations $\{{\delta }_{i}^{x}\}{}_{i=1}^{n}$ and $\{{\delta }_{i}^{y}\}{}_{i=1}^{n}$, sampled from a normal distribution with mean and standard deviation equal to $10 \% \sqrt{R(0)}$. The measurements ${\{{x}_{i}\}}_{i=1}^{n}$ and ${\{{y}_{i}\}}_{i=1}^{n}$ are then polluted with noise, and sampled at each observation time from the normal distributions $N(0,{{\delta }_{i}^{x}}^{2})$ and $N(0,{{\delta }_{i}^{y}}^{2})$.

For every generated data set, we fit nine different time delay models: TD-DRW, TD-CARMA(2, 1), and TD-CARMA(3, 1), each crossed with m = {1, 2, 3}.

We use two different classes of metrics to evaluate the performance of TD-CARMA in the context of this simulation study. First, we assess TD-CARMA purely in terms of its statistical estimator of Δ and use the classical measures of statistical fidelity of bias, variance, and root-mean-squared error (RMSE). Roughly speaking, these metrics describe the statistical properties of the error $\hat{{\rm{\Delta }}}-{\rm{\Delta }}$ of the estimator $\hat{{\rm{\Delta }}}$ (see Appendix B.1 for a rigorous definition).

To facilitate comparison with other time delay methods developed in the astrophysics literature, we also evaluate TD-CARMA as a time delay estimation technique using the accuracy, precision, and χ² metrics defined in Liao et al. (2015) for their TDC.

Specifically, suppose N simulated doubly lensed data sets with true time delays ${\{{{\rm{\Delta }}}_{i}\}}_{i=1}^{N}$ are fitted with TD-CARMA, which produces estimates ${\{{\hat{{\rm{\Delta }}}}_{i}\}}_{i=1}^{N}$ and posterior standard deviations ${\{{\sigma }_{i}\}}_{i=1}^{N}$. Then the accuracy A, precision P, and χ² are defined as

$$\begin{eqnarray}&&A=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{\hat{{\rm{\Delta }}}}_{i}-{{\rm{\Delta }}}_{i}}{{{\rm{\Delta }}}_{i}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&P=\displaystyle \frac{1}{N}\sum _{i=1}^{N}\displaystyle \frac{{\sigma }_{i}}{{{\rm{\Delta }}}_{i}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{\hat{{\rm{\Delta }}}}_{i}-{{\rm{\Delta }}}_{i}}{{\sigma }_{i}}\right)}^{2},\end{eqnarray} \tag{ 16 }$$

respectively.

We fit the nine different inference models (TD-DRW(m) and TD-CARMA(p, q, m) with m = {1, 2, 3}, p = {2, 3}, and q = 1) to each of the 2700 generated data sets, and for every data set we select the inference model with the highest Bayesian evidence for our estimate of the time delay. We then evaluate the performance of the 300 selected estimates for each of the nine data generation models using the metrics described in Section 4.2 and report them in Table 2. The results indicate that the highest evidence models (the overwhelming majority of which are TD-CARMA(3, 1, ·)) are able to estimate the true value of Δ to very high accuracy and precision. Irrespective of the generative model, the (absolute) accuracy of TD-CARMA is ∣A∣ = 0.00016 and precision is P = 0.00358, and the chi-squared metric is χ² = 0.63669. For comparison, entries in the 2015 TDC (Liao et al. 2015) were considered competitive if ∣A∣ < 0.03, P < 0.01, and χ² < 1.5. The bias, variance, and RMSE metric also illustrate the quality and robustness of the TD-CARMA estimator for time delay Δ. Out of the 2700 data sets, the TD-DRW model was only selected five times as the highest evidence model, which clearly shows the lack of flexibility of TD-DRW in modeling complex lens systems.

Table 3. Summary of our Analysis (Carried out in Section 5, of Six Doubly Lensed Quasars from SDSS, Observed and Previously Analyzed by the COSMOGRAIL Collaboration)

Data Set	COSMOGRAIL reference	COSMOGRAIL estimate	TD-CARMA
SDSS HS2209+1914	Eulaers et al. (2013)	−20 ± 5	−21.96 ± 1.448
SDSS J1001+5028	Kumar et al. (2013)	119.1 ± 3.3	120.93 ± 1.015
SDSS J1206+4332	Millon et al. (2020b)	113.0 ± 3	111.51 ± 1.452
SDSS J1515+1511	Millon et al. (2020b)	${210.5}_{-5.7}^{+5.5}$	210.80 ± 2.18
SDSS J1455+1447	Millon et al. (2020b)	${47.2}_{-7.8}^{+7.5}$	45.36 ± 1.93
SDSS J1349+1227	Millon et al. (2020b)	*	432.05 ± 1.950

Download table as:  ASCIITypeset image

To ensure that our error bars are well calibrated, we examine the frequency coverage of the posterior intervals. Specifically, we examine the proportion of simulated data sets for which the posterior intervals contain the true value of Δ, i.e., the coverage probability of the intervals. To be well calibrated from a frequency point of view, posterior intervals should also be confidence intervals (as defined in Equation (B4)), e.g., a (say) 95% confidence interval for a model parameter is defined to be an interval computed from the data in such a way that upon repeated sampling of the data set (and thus repeated computation of the interval), at least 95% of the intervals contain the true value of the parameter. We therefore aim to check that posterior intervals with posterior probability X% (i.e., the theoretical coverage) have (approximately) a coverage probability of X% (i.e., the actual coverage). The reader is referred to Appendix B for a review of frequency coverage analyses.

Figure 3 plots the theoretical coverage against the actual coverage probability observed in our simulation study, for p_n ∈ [0, 1]. The actual coverage is the proportion of p_n × 100% intervals (obtained from the highest evidence inference model for each of the 2700 simulated data sets) that actually contained the true value of Δ.

Zoom In Zoom Out Reset image size

The coverage of the 66% (1σ) and 95% (2σ) intervals, reported in Table 2 for each of the nine generative models, indicates overcoverage at the 66% level and good coverage at the 95% level. This is corroborated by Figure 3, where the estimated coverage line lies slightly above the theoretical coverage (45° line). This means that TD-CARMA produces conservative posterior intervals (i.e., slightly overestimated error bars), a desirable property in the context of time delay estimation. Nevertheless, our analyses of COSMOGRAIL quasar data sets in Section 5 show that TD-CARMA still produces much tighter error bars than other methods in the literature.

In conclusion, this simulation study demonstrates the flexibility of TD-CARMA in modeling lens systems of varying complexity, and the very high accuracy and precision at which TD-CARMA recovers the time delay parameter, Δ.

In this section, we apply TD-CARMA to a selection of doubly lensed quasar data sets, and compare the resulting estimates to those obtained in the literature. We have selected six quasars that appear to have particularly complex intrinsic light curves and on which TD-DRW failed to return reliable estimates of Δ (using the timedelay package of Tak et al. 2017). Thus, a DRW process may not be sufficiently robust to model the underlying light curves.

For each data set we fit a total of 28 models. Specifically, we fit TD-DRW with m = 3 and TD-CARMA with auto-regressive order p = {2, 3, 4}, moving average order q = {1, 2, 3}, and microlensing order m = {1, 2, 3}, such that p > q. In Appendix C, we report the full results of our analyses in tables that give summary statistics for the posterior distribution of the model parameters for each fitted model, including the posterior mean and variance for Δ, the Bayesian (log) evidence for the model, and the posterior mean of the frequency (computed from the posterior mean of the auto-regressive parameters α ) if the inferred PSD features a frequency (or frequencies). In the event of a multimodal posterior distribution, we report these statistics for each mode and give the value for their Bayesian local (log) evidence, i.e., the mode-specific evidence. ¹³

In all cases where the inferred PSD featured a frequency, the associated posterior distribution for the power of the frequency was below the measurement noise level. Therefore, following Section 3.5, we ignore all fits that included a frequency. We elaborate the numerical results for the frequencies analyses for all six examples in Section 5.6.

We select the model with highest Bayesian evidence (or relative probability), and compare its estimates for Δ with those derived by a baseline model, specifically our simplest model, TD-CARMA(2, 1, 3), ¹⁴ and those obtained by other methods in the literature. If the posterior distribution of the CARMA parameters under a particular model contains a mode with a frequency, we discard this mode and report the Bayesian evidence obtained by summing over the remaining modes (without frequencies). If the joint posterior is multimodal in the CARMA parameters, but none of the modes correspond to a frequency, we simply report the full Bayesian evidence given in Equation (11). The final estimates of Δ are reported in Table 3.

HS2209+1914 is a doubly lensed quasar observed as part of the COSMOGRAIL study (Eulaers et al. 2013) with four telescopes: Euler, Mercator, Maidanak, and HCT, during a total of 8.5 years (2004 to mid-2012). Eulaers et al. (2013) propose a time delay measurement of Δ = −20 ± 5 days (using a combination of four different techniques, referenced in Table 4), warning that the lack of obvious fast variability features in the observed light curves and the possible presence of microlensing made it a particularly difficult observation to study.

Our analysis also suggests the presence of a strong microlensing effect in this observation, as indicated by the dependence of the estimated time delay, $\hat{{\rm{\Delta }}}$, on the microlensing order, m. TD-CARMA models with m = {1, 2} microlensing produce estimates of Δ around −40 days, whereas those with m = 3 produce $\hat{{\rm{\Delta }}}\approx -22$ days (see Table 11); the latter is consistent with the literature (Eulaers et al. 2013; Kumar et al. 2015). (Our fitted $\hat{{\rm{\Delta }}}$ for all models that did not detect a frequency appear in Table 11 in Appendix C.) To conduct a more principled model comparison when models indicate different possibilities of Δ, we use the values of $\mathrm{ln}({ \mathcal Z })$ computed by MultiNest. This allows us to quantify the evidence in favor of models with m = 3, and the posterior probability of each model, via Equation (12). The cumulative probability of the m = 3 models is 0.999.

Our implementation of TD-DRW using MultiNest produced a (joint) posterior distribution with a total of 10 modes, and modes for Δ spanning from ≈−11 days to ≈43 days. The bulk of the probability mass of the posterior of Δ under TD-DRW is around ≈20 days, which indicates that TD-DRW fails even to get the sign of Δ correct. In addition, Tables 5, 11, and 12 show the Bayesian evidence $\mathrm{ln}({ \mathcal Z })$ in favor of TD-DRW is the lowest of all the considered models.

Table 4. Comparison of Our Highest Evidence Model, i.e., TD-CARMA(3, 2, 3), and our Baseline, TD-CARMA(2, 1, 3), to Estimates of Δ for HS2209+1914 in the Literature

Technique	$\hat{{\rm{\Delta }}}$	SD($\hat{{\rm{\Delta }}}$)
Dispersion a	−20.99	10.09
Difference a	−20.08	13.78
Spline a	−19.77	6.03
NMF a	−19.28	1.48
Combination of the four above a	−20	5
Difference smoothing (modified) b	−22.9	5.3
ROA c	−24.0	2.4
TD-CARMA(2, 1, 3)	−21.74	1.423
TD-CARMA(3, 2, 3)	−21.96	1.448

Notes.

^a Eulaers et al. (2013), where estimates for Δ are reported without the minus sign, as they swap labels for the two light curves compared to our notation. ^b Kumar et al. (2015). ^c Donnan et al. (2021).

Download table as:  ASCIITypeset image

The top left panel of Figure 4 compares the posterior distributions for Δ under the TD-CARMA model with the highest evidence, i.e., TD-CARMA(3, 2, 3), and the baseline model, i.e., TD-CARMA(2, 1, 3), with estimates found in the astrophysics literature. Figure 5 plots the data against the model predictions derived by TD-CARMA(3, 2, 3). Eulaers et al. (2013) propose an estimate of $\hat{{\rm{\Delta }}}=-20\pm 5$ days, and Kumar et al. (2015) give $\hat{{\rm{\Delta }}}=-22.9\pm 5.3$ days. TD-CARMA(3, 2, 3) and TD-CARMA(2, 1, 3) estimate $\hat{{\rm{\Delta }}}=-21.96\pm 1.448$ and $\hat{{\rm{\Delta }}}=-21.74\pm 1.423$, hence significantly reducing the error bars relative to existing methods (see Table 4). The simulation study in Section 4 demonstrates that our error bars are reliable and conservative.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. Posterior Mean and Standard Deviation for Δ under the Models with the Highest Bayesian Log-evidence, with Relative Probability π Summing to 0.999, for HS2209+1914

Model	$\hat{{\rm{\Delta }}}$	SD($\hat{{\rm{\Delta }}}$)	$\mathrm{ln}({ \mathcal Z })$	π
TD-DRW(3)	20.23	0.918	2536.00	0.000
TD-CARMA (3, 2, 3)	−21.96	1.448	2760.24	0.601
TD-CARMA (4, 2, 3)	−21.95	1.403	2759.83	0.399

Download table as:  ASCIITypeset image

SDSS J1001+5027 is a doubly lensed quasar, observed by the COSMOGRAIL collaboration (Kumar et al. 2013) with three telescopes: Mercator, Maidanak, and HCT for more than six years (from 2005 March to 2011 July), producing 443 observational epochs.

Our implementation of TD-DRW using MultiNest produced a (joint) posterior distribution with a total of 20 modes, and modes for Δ spanning values from ≈123 days and ≈135 days. The bulk of the probability mass of the posterior of Δ is around ≈132 days, almost two weeks away from the consensus value in the relevant literature, which is ≈120 days. The TD-DRW(3) model has $\mathrm{ln}({ \mathcal Z })=1803.30$, well below that of any of the TD-CARMA(p, q, m) (see Tables 6 and 13).

Table 6. Posterior Mean and Standard Deviation for Δ under the Five Models with the Highest Bayesian Log-evidence, for SDSS J1001+5027

Model	$\hat{{\rm{\Delta }}}$	SD($\hat{{\rm{\Delta }}}$)	$\mathrm{ln}({ \mathcal Z })$	π
TD-DRW(3)	132.71	0.770	1803.30	0.000
TD-CARMA(4, 3, 2)	120.93	1.015	2761.25	0.416
TD-CARMA(4, 3, 3)	120.21	1.018	2761.12	0.366
TD-CARMA(3, 2, 3)	120.16	1.003	2759.48	0.104
TD-CARMA(3, 2, 2)	120.88	1.013	2759.44	0.068
TD-CARMA(4, 2, 3)	120.21	0.970	2758.35	0.023

Download table as:  ASCIITypeset image

Table 6 shows the five TD-CARMA(p, q, m) with the highest Bayesian (log) evidences, as well as their relative probabilities, π_i , computed as detailed in Section 3.3. Together the five top models account for 97.7% of the total model probabilities. For given values of hyperparameters p and q, the increase in $\mathrm{ln}({ \mathcal Z })$ obtained from increasing m is small, indicating a weak effect of microlensing in this observation; see Tables 6 and 13 (the latter in Appendix C).

In Table 7, we compare the TD-CARMA(p, q, m) model with the highest Bayesian (log) evidence, i.e., TD-CARMA(4, 3, 2), with recently published analyses of SDSS J1001+5027. TD-CARMA(4, 3, 2) produces error bars that are between 30% and 83.63% lower than the existing methods; the baseline model, i.e., TD-CARMA(2, 1, 3), produces error bars that are between 48.34% and 87.92% lower. The larger error bars of TD-CARMA(4, 3, 2) relative to TD-CARMA(2, 1, 3) may be attributable to a bias–variance trade off: more highly parameterized models are more flexible and tend to provide better fits, but at the cost of higher levels of uncertainty.

Table 7. Comparison of Our Highest Evidence Model, i.e., TD-CARMA(4, 3, 2), and our Baseline, TD-CARMA (2, 1, 3), to Estimates of Δ for SDSS J1001+5027 in the Literature

Technique	$\hat{{\rm{\Delta }}}$	SD($\hat{{\rm{\Delta }}}$)
Dispersion-like a	120.5	6.2
Difference smoothing a	118.6	3.7
Regression difference a	121.1	3.8
Free-knot spline a	119.7	2.6
Combination of the four above a	119.1	3.3
GPs b	117.8	3.2
Difference smoothing (modified) c	119.7	1.8
ROA d	119.9	(−1.5, 1.4)
TD-CARMA(2, 1, 3)	120.18	0.749
TD-CARMA(4, 3, 2)	120.93	1.015

Note.

^a Eulaers et al. (2013). ^b Hojjati et al. (2013). ^c Kumar et al. (2015). ^d Donnan et al. (2021).

Download table as:  ASCIITypeset image

Kumar et al. (2013), the main reference in the estimation of the time delay of SDSS J1001+5027, argue that their four independent measurement techniques (listed in Table 7) each have their own limitations and thus combine the individual estimations from each method to produce their final estimate: $\hat{{\rm{\Delta }}}=119.3\pm 3.3$ days. Our final estimate is $\hat{{\rm{\Delta }}}\,=120.93\pm 1.015$ days, thus providing a consistent (since our posterior distribution covers their estimate) yet novel estimate of the time delay, that is three times more certain.

SDSS J1206+4332 is a doubly lensed quasar observed by the COSMOGRAIL collaboration (Eulaers et al. 2013) with three telescopes: Maidanak, Mercator, and HCT for a total of seven years (from 2005 to 2011).

For this data, TD-DRW has the lowest Bayesian evidence of the set of considered models, and has a relative probability of essentially 0 (up to numerical precision). TD-DRW produces a joint posterior distribution with six modes, with posterior modes for Δ ranging from ≈106 days to ≈117 days. Even though our computational procedure allows us to obtain a sample from the posterior (whereas the MCMC based implementation of TD-DRW in the timedelay package from Tak et al. 2017 fails to converge), TD-DRW does not appear to be sufficiently flexible for this data set.

Turning to fits obtained with TD-CARMA, Table 8 shows that the posterior distribution of Δ is sensitive to the value of m, which indicates the presence of a strong microlensing effect in this observation. Indeed, Figure 6 illustrates that third-order polynomial regression seems to pull the posterior distribution of Δ toward ≈112 days, whereas the posterior distributions of Δ obtained with m = 2 tend to center around ≈109 days. The posterior PSDs of the intrinsic light curve under both models are identical, which shows that this effect is indeed solely attributable to microlensing. In some cases, such as TD-CARMA(2, 1, 3), the posterior distribution covers both areas of likely values. This is also the case in Eulaers et al. (2013), where the numerical model fit (NMF) technique estimates $\hat{{\rm{\Delta }}}=113.80\pm 0.80$, whereas the spline technique finds $\hat{{\rm{\Delta }}}=111.31\pm 3.93$. This suggests that either the observational sampling does not allow a strong constraint for Δ, or that the effect of microlensing may be too complex for the currently employed modeling approaches.

Zoom In Zoom Out Reset image size

Table 8. Posterior Mean and Standard Deviation for Δ under the Five Models with the Highest Bayesian Log-evidence, for SDSS J1206+4332

Model	$\hat{{\rm{\Delta }}}$	SD($\hat{{\rm{\Delta }}}$)	$\mathrm{ln}({ \mathcal Z })$	π
TD-DRW(3)	114.68	1.413	1551.02	0.000
TD-CARMA(3, 2, 3)	111.51	1.452	1639.06	0.999
TD-CARMA(3, 2, 2)	109.27	1.351	1628.43	2.42e-5
TD-CARMA(2, 1, 3)	112.78	1.693	1624.65	5.52e-7
TD-CARMA(3, 1, 3)	112.77	1.670	1623.66	2.05e-7
TD-CARMA(2, 1, 2)	108.33	2.240	1616.64	1.83e-10

Download table as:  ASCIITypeset image

TD-CARMA(3, 2, 3) dominates the model space, with probability 0.999, and estimates $\hat{{\rm{\Delta }}}=111.51\pm 1.452$ days. This halves the uncertainty of Eulaers et al. (2013; 111.3 ± 3 days), and reduces the uncertainty derived in Kumar et al. (2015) by 23.58% (110.3 ± 1.9 days), as shown in Table 9. Although the estimate produced by Donnan et al. (2021; see "ROA Method" in Table 9) has less uncertainty than TD-CARMA for this particular data set, its similarity to the results obtained by the NMF technique (Eulaers et al. 2013) suggest that it may have left part of the parameter space unexplored.

Table 9. Comparison of our Highest Evidence Model, i.e., TD-CARMA(3, 2, 3), and our Baseline, TD-CARMA(2, 1, 3), to Estimates of Δ for SDSS J1206+4332 in the Literature

Technique	$\hat{{\rm{\Delta }}}$	SD($\hat{{\rm{\Delta }}}$)
Dispersion a	113.65	6.79
Difference a	109.73	8.28
Spline a	111.31	3.93
NMF a	113.80	0.90
Combination of the four above a	111.3	3
Splines+Regression difference b	111.8	(−2.7, 2.4)
Difference smoothing (modified) c	110.3	1.9
ROA d	113.0	(−1.2, 1.1)
TD-CARMA(2, 1, 3)	112.78	1.693
TD-CARMA(3, 2, 3)	111.51	1.452

Notes.

^a Eulaers et al. (2013). ^b Birrer et al. (2019). ^c Kumar et al. (2015). ^d Donnan et al. (2021).

Download table as:  ASCIITypeset image

In this section, we apply TD-CARMA to three doubly lensed quasars observed as part of the COSMOGRAIL study and recently analyzed in Millon et al. (2020b): SDSS J1515+1511, SDSS J1455+1447, and SDSS J1349+1227. For SDSS J1515+1511 and SDSS J1455+1447, we obtain estimates that are highly consistent with Millon et al. (2020b) and the relevant literature, and dramatically improve on the previously published uncertainties.

For SDSS J1515+1511, Shalyapin & Goicoechea (2017) propose a time delay measurement of Δ = 211 ± 5 days, and more recently Millon et al. (2020b) estimate ${\rm{\Delta }}={210.5}_{-5.7}^{+5.5}$ days. Our model with highest evidence, TD-CARMA(4, 2, 1), estimates Δ = 210.80 ± 2.18 days, dramatically reducing the uncertainty compared to previously published estimates. The full results are reported in Tables 17 and 18 in Appendix C.

For SDSS J1455+1447, Millon et al. (2020b) propose a time delay estimate of ${\rm{\Delta }}={47.2}_{-7.8}^{+7.5}$ days. Our highest evidence model, TD-CARMA(3, 2, 1), estimates Δ = 45.36 ± 1.93 days. Both our estimate and the one proposed on Millon et al. (2020b) are within one standard deviation of their respective error bars. Our estimate, however, is roughly four times less uncertain. The full results are reported in Tables 19 and 20 in Appendix C.

For SDSS J1349+1227, Kayo et al. (2010) predicted a time delay of the order of one year, and Millon et al. (2020b) could not determine the time delay of this lens system. We produced an estimate that is consistent with the one year order of magnitude value predicted in Kayo et al. (2010). For all values of the (p, q, m) hyperparameter triplet, TD-CARMA produces highly consistent estimates for Δ, all gravitating around 431–432 days. Our highest evidence model, TD-CARMA(3, 1, 1), estimates Δ = 432.05 ± 1.950 days. The full results are reported in Tables 21 and 22 in Appendix C.

As with Applications I–III, the TD-DRW(3) model produced highly multimodal posterior distributions for Δ. These estimates had much lower Bayesian evidence that our more complex TD-CARMA models, and were inconsistent with previously published analyses.

Tables 12, 14, 16, 18, 20, and 22 in Appendix C show that under some of the TD-CARMA(p, q, m), the posterior distribution of the CARMA parameters ( α , β , σ) correspond to PSDs with a frequency (i.e., a Lorentzian centered away from 0), either for specific modes or for the whole distribution. For instance, one mode of the posterior PSD under TD-CARMA(3, 1, 1) fitted to the SDSS J1206+4332 data set has a frequency of ≈1 day, and the posterior PSD under TD-CARMA(4, 1, 3) fitted to the SDSS J001+5027 data set finds a frequency of ≈2 days.

These frequencies might be an inherent property of the intrinsic light curves themselves, or simply a feature of the measurement noise.

For illustration, Table 10 shows the (posterior mean of the) detected frequencies, the associated 99% posterior intervals for the power, and the measurement noise level for the SDSS HS2209+1914, SDSS J1001+5027, and SDSS J1206+4332 data sets. In all three cases, the 99% posterior intervals for the power of the detected frequencies either cover or are well below the measurement noise level. ¹⁵ We therefore discard all models that detect a frequency in our analyses.

The ACF of a stochastic process $X(t),t\in {\mathbb{R}}$ is defined as E[(X(t_i ) − E(X(t_i ))(X(t_j ) − E(X(t_j ))], and gives the covariance of the process with itself at two different times. The ACF of a DRW process is given by

$$\begin{eqnarray}&&R({d}_{{ij}})=\displaystyle \frac{{\sigma }^{2}\tau }{2}{e}^{-\tfrac{{d}_{{ij}}}{\tau }},\end{eqnarray} \tag{ A1 }$$

where d_ij = t_i − t_j is the time difference between any two observation times t_i and t_j , with t_i > t_j (Kelly et al. 2009). The ACF in Equation (A1) features an exponentially decaying component with e-folding timescale τ which highlights the mean-reverting behavior of the DRW process.

The PSD, defined as the Fourier transform of the ACF, describes the distribution of power ¹⁶ over frequencies for continuous and stationary stochastic processes. The PSD of the DRW process, given by

$$\begin{eqnarray}&&P(f)={\sigma }^{2}\displaystyle \frac{1}{{\left(\tfrac{1}{\tau }\right)}^{2}+{\left(2\pi f\right)}^{2}},\end{eqnarray} \tag{ A2 }$$

is a Lorentzian centered at 0, with a break frequency at 1/2π τ² (Kelly et al. 2009). This illustrates the limited applicability of a DRW process in describing the stochastic variability in light curves. For example, quasiperiodic variability features in an AGN light curve can be represented as Lorentzians centered away from 0; hence a DRW process would be unable to characterize this type of variability. Moreover, the DRW process can only account for a single break frequency in the PSD, which renders it somewhat inflexible.

To ensure stationarity of the process, we must have q < p, and the roots of the auto-regressive polynomial A(z), defined as

$$\begin{eqnarray}&&A(z)=\displaystyle \sum _{k=0}^{p}{\alpha }_{k}{z}^{k},\end{eqnarray} \tag{ A3 }$$

must have negative real parts (Kelly et al. 2014; Tsai & Chan 2000; Tsai et al. 2011). In practice, this constraint is more easily enforced when considering an alternative parameterization of the auto-regressive polynomial (Jones 1981; Kelly et al. 2014), obtained by factorizing it into linear and quadratic terms as follows

$$\begin{eqnarray*}\begin{array}{ccc}A\left(z\right) & = & \left({a}_{1}+{a}_{2}z+{z}^{2}\right)\left({a}_{3}+{a}_{4}z+{z}^{2}\right)...\\ & & \times \,\left\{\begin{array}{ll}{a}_{p-1}+{a}_{p}z+{z}^{2} & {\rm{i}}{\rm{f}}\,p\,{\rm{i}}{\rm{s}}\,{\rm{e}}{\rm{v}}{\rm{e}}{\rm{n}},\\ {a}_{p}+z & {\rm{i}}{\rm{f}}\,p\,{\rm{i}}{\rm{s}}\,{\rm{o}}{\rm{d}}{\rm{d}}.\end{array}\right.\end{array}\end{eqnarray*}$$

The roots of the auto-regressive polynomial are negative if the coefficients a = {a₁..., a_p } are positive. We discuss computational strategies for stable implementation and inference in Section 3.4.

The PSD of a stationary CARMA(p, q) process is given by

$$\begin{eqnarray}&&P(f)={\sigma }^{2}\displaystyle \frac{{\left|{\sum }_{j=0}^{q}{\beta }_{j}{\left(2\pi {if}\right)}^{j}\right|}^{2}}{{\left|{\sum }_{k=0}^{p}{\alpha }_{k}{\left(2\pi {if}\right)}^{k}\right|}^{2}},\end{eqnarray} \tag{ A4 }$$

and its ACF is

$$\begin{eqnarray}&&R({d}_{{ij}})={\sigma }^{2}\displaystyle \sum _{k=1}^{p}\displaystyle \frac{\left[{\sum }_{l=0}^{q}{\beta }_{l}{r}_{k}^{l}\right]\left[{\sum }_{l=0}^{q}{\beta }_{l}{\left(-{r}_{k}\right)}^{l}\right]\exp ({r}_{k}{d}_{{ij}})}{-2\mathrm{Re}({r}_{k}){\prod }_{l=1,l\ne k}({r}_{l}-{r}_{k})({r}_{l}^{* }+{r}_{k})},\end{eqnarray} \tag{ A5 }$$

where r_k , k = 1,..., p are the roots of the auto-regressive polynomial A(z) (Tsai & Chan 2000; Tsai et al. 2011; Kelly et al. 2014).

The ACF of a CARMA process is therefore a weighted sum of exponentially damped sinusoidal functions (corresponding to complex roots) and exponential decays (corresponding to real roots; Kelly et al. 2014). A CARMA process has the flexibility to model exponentially decaying auto-correlations with multiple timescales (rather than a single one for the DRW process), and to incorporate periodic oscillations in the auto-correlations. Figure 7 illustrates the difference between the DRW and CARMA processes in terms of their auto-correlations and highlights the fact that even within a CARMA(p, q) family, CARMA processes with different orders of p and q may look quite different. This implies that the flexibility of CARMA processes makes them suitable to characterize the variability in a wide range of AGN light curves.

Zoom In Zoom Out Reset image size

Here we formally define the concepts of statistical bias, variance, and RMSE used in Section 4 in our numerical evaluation of the statistical veracity of TD-CARMA. These quantities are widely used to describe the frequency properties of statistical estimators. Denote by $\hat{{\rm{\Delta }}}| {\boldsymbol{D}}$ the distribution of the estimate for time delay parameter over replicates of data d ∼ D , obtained by repeatedly fitting a given model to randomly generated replicate data sets. The bias of an estimate $\hat{{\rm{\Delta }}}$ is defined as

$$\begin{eqnarray}&&\mathrm{Bias}(\hat{{\rm{\Delta }}})={{\mathbb{E}}}_{\hat{{\rm{\Delta }}}| {\boldsymbol{D}}}(\hat{{\rm{\Delta }}})-{\rm{\Delta }},\end{eqnarray} \tag{ B1 }$$

i.e., the average estimate of Δ across replicate data sets, from which the true value of Δ is subtracted. Similarly, the variance of the estimate $\hat{{\rm{\Delta }}}$ is defined as

$$\begin{eqnarray}&&\mathrm{Var}(\hat{{\rm{\Delta }}})={{\mathbb{E}}}_{\hat{{\rm{\Delta }}}| {\boldsymbol{D}}}\left[{\left({{\mathbb{E}}}_{\hat{{\rm{\Delta }}}| {\boldsymbol{D}}}(\hat{{\rm{\Delta }}})-\hat{{\rm{\Delta }}}\right)}^{2}\right].\end{eqnarray} \tag{ B2 }$$

Finally, the RMSE is defined as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{RMSE}(\hat{{\rm{\Delta }}}) & = & {{\mathbb{E}}}_{\hat{{\rm{\Delta }}}| {\boldsymbol{D}}}\left[{\left(\hat{{\rm{\Delta }}}-{\rm{\Delta }}\right)}^{2}\right]\\ & = & \sqrt{\mathrm{Bias}{\left(\hat{{\rm{\Delta }}}\right)}^{2}+\mathrm{Var}(\hat{{\rm{\Delta }}})}.\end{array}\end{eqnarray} \tag{ B3 }$$

Models with lower RMSEs are generally preferred as they return estimates that tend to be closer to the true value of the parameter.

A 100(1 − α)% confidence interval for parameter Δ is a random interval [a( D ), b( D )] such that

$$\begin{eqnarray}&&{\mathbb{P}}(a({\boldsymbol{D}})\lt {\rm{\Delta }}\lt b({\boldsymbol{D}}))=1-\alpha .\end{eqnarray} \tag{ B4 }$$

In other words, the interval [a( D ), b( D )] covers Δ with probability 1 − α. In theory, if we were to sample new data d ∼ D repeatedly and construct a 95% confidence interval [a(d), b(d)] every time, then approximately 95% of our intervals will contain the true value of Δ. The proportion of 100(1 − α)% confidence intervals that contain the true value of Δ is referred to as the empirical or actual coverage of the 100(1 − α)% confidence intervals, and 100(1 − α)% is referred to as the theoretical coverage. A model with perfectly calibrated confidence intervals (i.e., error bars) would have its actual coverage equal to the theoretical coverage for all α ∈ [0, 1]. Because of the approximate nature of the statistical model and the inferential procedure for the model parameters (nested sampling in our case), this may, however, not be the case in practice.

To test the "well calibratedness" of TD-CARMA's error bars, and to verify that the posterior distributions for Δ inferred by TD-CARMA on real data can be trusted, we compute the empirical coverage of 100(1 − α)% confidence intervals for all α ∈ [0, 1] as follows. We simulate N = 300 doubly lensed data sets under our model (as detailed in Section 4) with true time delay Δ_true. Then, we fit TD-CARMA to each of these data sets. For a given confidence level 100(1 − α)%, we compute the 100(1 − α)% confidence intervals derived by TD-CARMA on each simulated data set, and derive the empirical coverage as the proportion of those intervals that contain Δ_true. We repeat this procedure for every α ∈ [0, 1], and plot the empirical coverage against the corresponding α (see Figure 3). Finally, we visually inspect the resulting plot, and if the empirical coverage matches the theoretical coverage (i.e., if the plotted curve is "close enough" to the 45° line representing perfect calibration), we conclude that our model produces well-calibrated error bars.