FAST AND RELIABLE TIME DELAY ESTIMATION OF STRONG LENS SYSTEMS USING THE SMOOTHING AND CROSS-CORRELATION METHODS

Over the last few decades, cosmology has emerged as a precision science with access to a variety of data from different resources. Multiple images of strong lensing systems have proven to be a rich source of information for extracting some precise knowledge about the expansion history of the universe and for breaking the degeneracies between different cosmological models (Linder 2004; Suyu et al. 2013).

In strong gravitational lensing, not only do we have multiple images of the same source, but since the path of light is different for different images there will be a time delay between these images. When the source is variable, we might be able to estimate the time delays between different images. If we can also obtain proper information about the strong gravitational lens system and model it appropriately (to estimate the gravitational potentials different images have passed through), then we can use all of this information to measure the cosmic distances between the source, the strong gravitational lens system, and our position; this information can then be used to reconstruct the expansion history of the universe (Refsdal 1964; Kochanek 2002). There have been recent advances on both fronts: in modeling of strong lens systems (Oguri 2007; Suyu et al. 2009; Fadely et al. 2010; Suyu et al. 2010), and there have also been programs for long-term observations of multiple images of some strong lensing systems in order to estimate the time delays between the multiple images (Courbin et al. 2011; Eulaers et al. 2013; Rathna Kumar et al. 2013; Tewes et al. 2013b). The importance of analyzing strong lens systems in order to reconstruct the expansion history of the universe and their broad cosmological implications are becoming increasingly evident to the cosmological community. This is due to the fact that one can use these systems to measure the expansion history of the universe in a completely different way from standardized candles and rulers such as Type Ia supernovae or baryon acoustic oscillations. Hence, strong lensing systems can be used as a complementary approach for studying the expansion of the universe and even to calibrate the other probes of the expansion history (Linder 2004, 2011; Suyu et al. 2012).

From an observational point of view, the COSMOGRAIL project (http://cosmograil.org) is will yield tens of time delays for strong lens systems, and in the future the Large Synoptic Survey Telescope (LSST; LSST Science Collaboration et al. 2009; Ivezic et al. 2008) will monitor several thousand time delay lens systems (Oguri & Marshall 2010; LSST Dark Energy Science Collaboration 2012). This ensures that there will be sufficient information in the near future to perform precision cosmology using strong lensing systems if we can appropriately model the lens systems and precisely estimate the time delays. From a theoretical point of view, while there have been important advances in the modeling of strong lens systems, there have also been various attempts to provide efficient and robust algorithms in order to estimate the time delays between strong lens system images with high precision and accuracy (Hojjati et al. 2013; Suyu et al. 2013; Tewes et al. 2013a, 2013b). These proposed algorithms are all very different in nature and have both advantages and weaknesses. In an effort to move toward more accurate and precise algorithms of time delay estimation, a group of scientists organized the Time Delay Challenge (TDC) to encourage different scientific groups to modify or propose new algorithms for time delay estimation (Dobler et al. 2015; Liao et al. 2015). In this challenge, some simulated data were prepared in two stages, TDC0 and TDC1, and participants were asked to perform their own analysis and report the estimated time delays and corresponding uncertainties. As one of the participating groups, we joined the challenge and developed our own fully automated and reliable algorithm to estimate time delays (for more details about TDC1 and the proposed algorithms by participating groups, see Liao et al. 2015). Our algorithm is also fast when analyzing light curves. The entire processing time for a pair of light curves using a typical PC is of the order of a few minutes. We use the R programing language (R Core Team 2013) in both computational and plotting tasks. In this paper, we introduce and explain the algorithm that we used in the TDC1 challenge.

In the following, Section 2 describes the TDC1 simulated data. In Section 3, we explain the methodology which contains the smoothing method and cross-correlation and the strategy we follow to estimate the time delays. We discuss our fast and reliable approach for estimating the uncertainties of the recovered time delays, and finally in Section 4 we present our results and end with a conclusion.

The TDC1 simulated data is provided in five different categories (rungs). Each rung contains the light curves of 720 Double and 152 Quad image systems. Tables 1 lists the properties of the data in each rung. For each rung and within all of the data samples, the number of data segments (segments are separated where there is a gap of more than 100 days between two neighboring data points.), the minimum and maximum number of data points, the minimum and maximum of the length of data, and the minimum and maximum of the shortest and longest gaps between two data points are tabulated. We see that the Doubles and Quads in each rung have nearly the same properties. Figure 1 depicts the typical light curves for each rung. From Figure 1 and the tables, we note that the data are not equispaced in time and also suffer from gaps in some intervals, and in some cases there are evident microlensing effects. These specifications ensure that there is a need for sophisticated and perhaps complicated statistical approaches in order to correctly estimate the time delays. In Section 3, we elaborate on the methodology that we employ to deal with this data.

Zoom In Zoom Out Reset image size

To find the time delay between two light curves, we somehow need to make a comparison between the curves. We apply cross-correlation to find the lag time between the two light curves. We argue that we should obtain the highest cross-correlation between the two light curves if one of them is shifted in time coordinate by approximately the actual time delay between the curves. This seems trivial since all of the light curves come from the same source and should have similar features. However, as we described in the previous section, the light curve data are not equispaced in time, with several seasonal gaps. This sparseness of data can in fact cause inaccurate or even incorrect estimates of the time delay. Therefore, in practical cases, one cannot apply cross-correlation directly to the data to estimate the time delays in strong lens systems. To tackle this issue, before using cross-correlation to compare the light curves in each system, we obtain the proper smoothed forms of the light curves. Then, we derive the cross-correlations between the various light curves and the smoothed forms. In this section, we first introduce the smoothing method and the cross-correlation which we use to estimate the time delays and we then propose a fast and efficient algorithm for error estimation.

3.1. Smoothing

We adopt a model-independent smoothing method introduced by Shafieloo et al. (2006), Shafieloo (2007), Shafieloo & Clarkson (2010), and Shafieloo (2012). In this method, the smooth light curve, A^s(t), is obtained at any required time t by an iterative procedure:

$$\begin{eqnarray}&&{{A}^{s}}(t)={{A}^{g}}(t)+N(t)\mathop{\sum }\limits_{i}\frac{{{A}^{d}}\left( {{t}_{i}} \right)-{{A}^{g}}\left( {{t}_{i}} \right)}{\sigma _{d}^{2}\left( {{t}_{i}} \right)}{\rm exp} \left[ -\frac{{{\left( {{t}_{i}}-t \right)}^{2}}}{2{{{\Delta }}^{2}}} \right],\ \;\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&N{{(t)}^{-1}}=\mathop{\sum }\limits_{i}{\rm exp} \left[ -\frac{{{\left( {{t}_{i}}-t \right)}^{2}}}{2{{{\Delta }}^{2}}} \right]\frac{1}{\sigma _{d}^{2}\left( {{t}_{i}} \right)},\end{eqnarray} \tag{ 2 }$$

where ${{A}^{d}}\left( {{t}_{i}} \right)$ and σ_d(t_i) represent the ith data point and its uncertainty, A^g(t) is the guess model, N(t) is the normalization factor, and Δ is the smoothing width. In this procedure, the first iteration begins with an initial guess model, A^g0(t), which results in the smooth fit, A^s1(t). In the next iteration, the obtained smooth fit is used as the next guess model ${{A}^{g1}}(t)={{A}^{s1}}(t)$, which yields the new smooth fit A^s2(t). The procedure can be continued through future iterations. It has been shown before that results are independent of the assumption of the initial guess model (Shafieloo 2007). We should also note that the size of Δ represents the minimum size of the fine structure in the data that will finally be presented in the smoothed reconstruction. In this procedure, the smoothing width Δ also has a direct connection to the speed of convergence to the final required smooth fit. In fact, in order to reduce the time of computation, it is important to find a desirable combination of Δ and the number of required iterations in the procedure. We select the optimum value of Δ and the number of iterations by simulating data based on some of the TDC0 data. In our procedure, we set Δ = 8 days and 3 iterations. By choosing a small value of Δ we reduce computation time by having a very low number of iterations in the procedure. Based on our simulations, we noted that choosing a smaller value of Δ will result in noise fitting rather than finding actual features in the light curves, while on the other hand choosing a much larger value of Δ could smooth out the light curves significantly and wash out the important features in the data. Our simulation results indicate that while Δ = 8 seems to be a reasonable value, in the case of the TDC0 and TDC1 data choosing 3 < Δ < 12 would be acceptable.

3.2. Cross-correlation

The cross-correlation between the light curves (two time series of data X, Y) in a system measures the correlation coefficient ρ_XY in different time lags. This results in a set of correlation coefficients corresponding to different lags. The lag corresponding to the maximum coefficient gives the estimated time delay (Peterson 2001). To derive these correlation coefficients at any time lag, we need to have the data on both light curves corresponding to the assumed time lag, which cannot be the case since the data is not equispaced and is sparse in some cases. However, as discussed in Section 3.1, we tackle this issue by reconstructing the smooth forms of the light curves so that we can measure the correlation coefficients at any given time lag.

The Pearson correlation coefficient has a value between -1 and 1, and shows the linear dependency between the two random variables (Wasserman 2004):

$$\begin{eqnarray}&&{{\rho }_{XY}}=\frac{\operatorname{cov}(X,Y)}{{{\sigma }_{X}}{{\sigma }_{Y}}},\end{eqnarray} \tag{ 3 }$$

where X and Y indicate two time series with standard deviations σ_X and σ_Y, respectively, and cov(X, Y) is their covariance. The Pearson correlation coefficient can be estimated through the sample correlation coefficient, r_xy,

$$\begin{eqnarray}&&{{r}_{xy}}=\frac{\mathop{\sum }\limits_{i=1}^{n}\left( {{x}_{i}}-\bar{x} \right)\left( {{y}_{i}}-\bar{y} \right)}{\sqrt{\mathop{\sum }\limits_{i=1}^{n}{{\left( {{x}_{i}}-\bar{x} \right)}^{2}}\mathop{\sum }\limits_{i=1}^{n}{{\left( {{y}_{i}}-\bar{y} \right)}^{2}}}},\end{eqnarray} \tag{ 4 }$$

where ${{x}_{i}},{{y}_{i}},(i=1,2,..,n)$ are data samples with sample means $\bar{x}$ and $\bar{y}$ of two time series X and Y, respectively (Peterson 2001). It can be also shown that the linear transformation of random variables X and Y does not affect the correlation coefficient (Wasserman 2004). Hence, ${{\rho }_{XY}}={{\rho }_{X^{\prime} Y^{\prime} }}$ if $X^{\prime} ={{a}_{x}}X+{{b}_{x}},$ $Y^{\prime} ={{a}_{y}}Y+{{b}_{y}}$, where a_x, b_x, a_y, and b_y are constant. This property plays an important role in accurately estimating the time delay when data suffers from microlensing, which we discuss in Section 3.3.

3.3. Time Delay Estimation

Given a pair of light curves ${{A}_{1}},{{A}_{2}}$, we first apply the smoothing method (Section 3.1) to obtain the smooth fits $A_{1}^{s},A_{2}^{s}$. Figure 2 depicts a pair of typical light curves and their smooth fits. We can see that the data are in segments and that there are gaps between them. We should note that we only consider the data segments and their corresponding smooth fit. To calculate the cross-correlation, we compare individual segments of the two light curves together when there are at least three data points in the overlapped region. Next, we add together the correlation coefficients corresponding to each time lag and divide by the number of segments to obtain the average correlation coefficient.

Zoom In Zoom Out Reset image size

To maintain good control of the errors, we calculate the cross-correlation between A₁ and A^s₂ as well as between A₂ and A^s₁. The best time delay in each case corresponds to the maximum correlation. This leads to two close values (absolute) for the time delay, $|{{\tilde{{\Delta }t}}_{{{A}_{1}};A_{2}^{s}}}|$ and $|{{\tilde{{\Delta }t}}_{{{A}_{2}};A_{1}^{s}}}|$. If the maximum average correlation coefficient in each case becomes more than 0.6, then we take their mean value, $|{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}|$, and the standard deviation,³ $\sigma _{{{A}_{1}}{{A}_{2}}}^{{\rm ini}}$, as the estimated time delay and its initial error:

$$\begin{eqnarray}&&|{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}|=\frac{|{{\tilde{{\Delta }t}}_{{{A}_{1}};A_{2}^{s}}}|+|{{\tilde{{\Delta }t}}_{{{A}_{2}};A_{1}^{s}}}|}{2}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \sigma _{{{A}_{1}}{{A}_{2}}}^{{\rm ini}} \\ =\sqrt{{{\left( |{{\tilde{{\Delta }t}}_{{{A}_{1}};A_{2}^{s}}}|-|{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}| \right)}^{2}}+{{\left( |{{\tilde{{\Delta }t}}_{{{A}_{2}};A_{1}^{s}}}|-|{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}| \right)}^{2}}} \\ =\sqrt{2}\times \frac{\left| {{\tilde{{\Delta }t}}_{{{A}_{1}};A_{2}^{s}}}|-|{{\tilde{{\Delta }t}}_{{{A}_{2}};A_{1}^{s}}} \right|}{2}. \\ \end{array}\end{eqnarray} \tag{ 6 }$$

As an example, Figure 3 (the upper left and upper right panels) illustrates the cross-correlation of segments of the light curve data and smooth fits (illustrated in Figure 2) obtained by comparing A₁ with A^s₂ and A₂ with A^s₁, respectively. The corresponding best time delays are obtained ${{\tilde{{\Delta }t}}_{{{A}_{1}};A_{2}^{s}}}=+103.15$ days and ${{\tilde{{\Delta }t}}_{{{A}_{2}};A_{1}^{s}}}=-103.59$ days at maximum average correlation coefficients of 0.87 and 0.88, respectively (Figure 3, lower panels). Using Equations (5) and (6), the time delay and its initial error are estimated as $|{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}|=103.37$ days and $\sigma _{{{A}_{1}}{{A}_{2}}}^{{\rm ini}}=0.31$ days respectively.⁴

Zoom In Zoom Out Reset image size

Our rejection criterion ρ < 0.6 is inferred from the feedback of our different submissions for simulated light curves in the TDC0 challenge where we tried to achieve zero outliers and to avoid any single incorrect estimation of a time delay. We will elaborate on the error estimation in Section 3.4.

Micolensing can generally cause distortions in light curves which cause difficulties in finding the time delay. Since we compare the light curves in segments, the effect of microlensing can be considered as a linear distortion in different segments of the light curves. As we discussed in Section 3.2, the correlation coefficient is not changed under the linear transformation of variables. Therefore, there is no concern for microlensing in our methodology and this is one of the important advantages of our proposed method. The results of the TDC0 and TDC1 challenge have shown that we were correct in our argument as we could submit time delay estimates consistently for more than 30% of the cases for all data rungs.

3.4. Error Estimation

Apart from the derived initial uncertainty (Section 3.3) of the estimated time delay for a pair of light curves (standard deviation of the two derived values of time delays comparing smooth fits and actual light curves), we need to consider other uncertainties. We should be very careful to avoid the effects of random fluctuations as they may mimic features or signals. Performing some extensive simulations, we found that even when $|{{\tilde{{\Delta }t}}_{{{A}_{1}},A_{2}^{s}}}|$ is very close to $|{{\tilde{{\Delta }t}}_{{{A}_{2}},A_{1}^{s}}}|$, we might still be considerably off from the actual time delay in some cases. This happens when the data is sparse with large gaps, and can result in a huge source of error when estimating time delays. For each specific case, one can perform a number of simulations to understand the effect of the data distribution on the uncertainties of the time delay, and this seems to be a proper approach. However, since we had limited time for our analysis (we joined the challenge two months before its deadline), we chose a novel approach for our error estimation. We estimated the errors for different light curve pairs using an estimation of the errors based on the rich information in the Quad systems. By using a consistency relation between the various light curves of the Quad systems, we calibrated our error estimation for different rungs with various time delays. In other words, we can extract more information from Quad systems which can indirectly be used in our error estimation.

Suppose the four light curves of a Quad image are labeled A₁, A₂, B₁, and B₂. For every Quad system, we should have

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} {{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}-\left( {{\tilde{{\Delta }t}}_{{{A}_{1}}{{B}_{1}}}}+{{\tilde{{\Delta }t}}_{{{B}_{1}}{{A}_{2}}}} \right) \\ \;\pm \sqrt{{{\left( \sigma _{{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}}^{{\rm ini}} \right)}^{2}}+{{\left( \sigma _{{{\tilde{{\Delta }t}}_{{{A}_{1}}{{B}_{1}}}}}^{{\rm ini}} \right)}^{2}}+{{\left( \sigma _{{{\tilde{{\Delta }t}}_{{{B}_{1}}{{A}_{2}}}}}^{{\rm ini}} \right)}^{2}}}\equiv 0 \\ \end{array},\end{eqnarray} \tag{ 7 }$$

where ${{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}$ is the estimated time delay between the light curves A₁ and A₂, and $\sigma _{{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}}^{{\rm ini}}$ is the corresponding initial error (and the same notation is used for other pairs). In short,

$$\begin{eqnarray}&&{{T}_{{\rm dif}}}\pm {{\sigma }_{{{T}_{{\rm dif}}}}}\equiv 0,\end{eqnarray} \tag{ 8 }$$

where ${{T}_{{\rm dif}}}={{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}-\left( {{\tilde{{\Delta }t}}_{{{A}_{1}}{{B}_{1}}}}+{{\tilde{{\Delta }t}}_{{{B}_{1}}{{A}_{2}}}} \right)$ and ${{\sigma }_{{{T}_{{\rm dif}}}}}=\sqrt{{{\left( \sigma _{{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}}^{{\rm ini}} \right)}^{2}}+{{\left( \sigma _{{{\tilde{{\Delta }t}}_{{{A}_{1}}{{B}_{1}}}}}^{{\rm ini}} \right)}^{2}}+{{\left( \sigma _{{{\tilde{{\Delta }t}}_{{{B}_{1}}{{A}_{2}}}}}^{{\rm ini}} \right)}^{2}}}$. Now, if $\left| {{T}_{{\rm dif}}} \right|\leqslant {{\sigma }_{{{T}_{{\rm dif}}}}}$, then Equation (8) is valid and we can therefore assume that all time delays and their corresponding errors are estimated consistently. If not, then the estimated errors are smaller than they should be. Hence, we add an error correction σ_ec such that Equation (8) will be satisfied. The error correction would be calculated as follows:

$$\begin{eqnarray}&&\sigma _{{\rm ec}}^{2}={{\left| {{T}_{{\rm dif}}} \right|}^{2}}-\sigma _{{{T}_{{\rm dif}}}}^{2}.\end{eqnarray} \tag{ 9 }$$

We add this error correction quadratically to the initial error. Therefore, the new uncertainty for ${{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}$ would be

$$\begin{eqnarray}&&\sigma _{{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}}^{{\rm new}}=\sqrt{{{\left( \sigma _{{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}}^{{\rm ini}} \right)}^{2}}+\frac{\alpha }{3}\sigma _{{\rm ec}}^{2}}.\end{eqnarray} \tag{ 10 }$$

To be more conservative, we choose α = 2 instead of α = 1 (which is an optimum consideration) in the above equation. The above argument and derivation is also valid for the error estimation of ${{\tilde{{\Delta }t}}_{{{B}_{1}}{{B}_{2}}}}$.

Based on the above error estimation for Quad pairs, we establish a profile to estimate the error for Double systems. We plot the estimated errors and relative errors versus estimated time delays (see Figure 4) from the Quad systems for all rungs. This information suggests how large the uncertainties should be for different systems in different rungs depending on their estimated time delays. According to this information, we suggest that the uncertainties σ_R for Double systems should be as follows.

Zoom In Zoom Out Reset image size

For Rung 0,

$$\begin{eqnarray*}&&{\rm if}\;|\tilde{{\Delta }t}|\leqslant 20\Rightarrow {{\sigma }_{R}}=0.06\times |\tilde{{\Delta }t}|,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm if}\;|\tilde{{\Delta }t}|\gt 20\Rightarrow {{\sigma }_{R}}=1.2.\end{eqnarray*}$$

For Rung 1,

$$\begin{eqnarray*}&&{\rm if}\;|\tilde{{\Delta }t}|\leqslant 20\Rightarrow {{\sigma }_{R}}=0.06\times |\tilde{{\Delta }t}|,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm if}\;|\tilde{{\Delta }t}|\gt 20\Rightarrow {{\sigma }_{R}}=1.3.\end{eqnarray*}$$

For Rung 2,

$$\begin{eqnarray*}&&{\rm if}\;|\tilde{{\Delta }t}|\leqslant 30\Rightarrow {{\sigma }_{R}}=0.07\times |\tilde{{\Delta }t}|,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm if}\;|\tilde{{\Delta }t}|\gt 30\Rightarrow {{\sigma }_{R}}=1.3.\end{eqnarray*}$$

For Rung 3,

$$\begin{eqnarray*}&&{\rm if}\;|\tilde{{\Delta }t}|\leqslant 30\Rightarrow {{\sigma }_{R}}=0.08\times |\tilde{{\Delta }t}|,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm if}\;|\tilde{{\Delta }t}|\gt 30\Rightarrow {{\sigma }_{R}}=1.5.\end{eqnarray*}$$

For Rung 4,

$$\begin{eqnarray*}&&{\rm if}\;|\tilde{{\Delta }t}|\leqslant 25\Rightarrow {{\sigma }_{R}}=0.08\times |\tilde{{\Delta }t}|,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm if}\;|\tilde{{\Delta }t}|\gt 25\Rightarrow {{\sigma }_{R}}=1.5.\end{eqnarray*}$$

One should note that we used these simple boundaries for simplification and in principle one can look directly at the profile of the results from Quad systems to estimate the expected uncertainties for a Double system with a particular time delay in a specific rung. To be conservative, these estimated errors have been added to the initial errors quadratically. So finally, the errors we considered for our submissions in the TDC1 challenge were (for any pair A₁A₂)

$$\begin{eqnarray}&&{{\sigma }_{{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}}}=\sqrt{{{\left( \sigma _{{{\tilde{{\Delta }t}}_{{{A}_{1}}{{A}_{2}}}}}^{{\rm ini}} \right)}^{2}}+\sigma _{R}^{2}}\end{eqnarray} \tag{ 11 }$$

To evaluate the results of the challenge, Evil team used the following metrics (Liao et al. 2015):

f is the fraction of the submitted time delays,

$$\begin{eqnarray}&&f\equiv \frac{{{N}_{{\rm submitted}}}}{N};\end{eqnarray} \tag{ 12 }$$

χ² shows the average weighted distance between the estimated and true time delays,

$$\begin{eqnarray}&&{{\chi }^{2}}=\frac{1}{fN}\mathop{\sum }\limits_{i}{{\left( \frac{\tilde{{\Delta }{{t}_{i}}}-{\Delta }{{t}_{i}}}{{{\sigma }_{i}}} \right)}^{2}};\end{eqnarray} \tag{ 13 }$$

P is the average precision metric,

$$\begin{eqnarray}&&P=\frac{1}{fN}\mathop{\sum }\limits_{i}\left( \frac{{{\sigma }_{i}}}{|{\Delta }{{t}_{i}}|} \right);\end{eqnarray} \tag{ 14 }$$

and A reflects the average bias of estimation,

$$\begin{eqnarray}&&A=\frac{1}{fN}\mathop{\sum }\limits_{i}\left( \frac{\tilde{|{\Delta }{{t}_{i}}|}-|{\Delta }{{t}_{i}}|}{|{\Delta }{{t}_{i}}|} \right),\end{eqnarray} \tag{ 15 }$$

where N is the total number of light curves for the analysis, N_submitted is the number of submitted time delays, Δt_i is the true time delay, and $\tilde{{\Delta }{{t}_{i}}}$ is the estimated time delay with uncertainty σ_i. The criteria for passing the TDC0 stage were f > 0.3, P < 0.15, A < 0.15, and χ² < 2. We passed the TDC0 stage in our first attempt, and working with the TDC1 data, since there was no fixed criteria, our rational expectation was that we should set the TDC0 limits of the metrics as the lower boundaries of the TDC1 criteria. Hence, we focused on improving our results in all four criteria.

Table 2 shows the feedback of our results given by the TDC1 challenge Evil team. For all rungs we achieved f > 0.3, which was one the main criteria for passing the TDC0 stage of the challenge. This table clearly shows that despite having a considerably large f factor, our results do not include numerous outliers and our algorithm estimates the time delays reliably and precisely. It is worth emphasizing that all of the steps of our algorithm are fully automated and there is no need for direct human evaluation at any stage. The results of this table show that we could achieve a very good f factor with f > 0.3 for all rungs and with precise measurements of time delays and their uncertainties with χ² < 1 and P < 0.06 for all rungs. However, we can see that there is a slight bias in our results, as reflected in the A metric.

Figures 5–7 depict the χ², P, and A individually for all of the entries that have χ² < 10 (the entries of rung 0, 1, 2, 3, and 4 are colored black, red, blue, green, and gold, respectively). We choose the χ² < 10 cut-off to be consistent with the Evil team settings in the TDC1 paper (Liao et al. 2015). Our results show that there are only five items with ${{\chi }^{2}}\geqslant 10$ out of a few thousand entries (within all rungs), which is indeed a very small fraction and statistically insignificant. We have plotted the histograms of the metrics in Figure 8. In Figure 9, the estimated time delays versus the true time delays are plotted. For comparison, the 45° line and its 2% deviations are also plotted. This plot clearly shows the consistency and power of our method where there is no outlier in any of the rungs and in almost all of the cases we had a reliable and precise estimate of the time delays.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As mentioned earlier, our results had a small bias (A was derived to be approximately 0.018 at all rungs) which indicates a need for further calibration. However, after the true time delay values of the TDC0 and TDC1 challenges were revealed, we noticed that our time delay estimations in all rungs had a constant 0.5 day bias. This constant bias appeared to be due to the simple fact that the TDC0 rung 0 simulated data was equispaced with one day gaps between neighboring data points and we calibrated the method to find time delays on these integer dates. So, in fact, the bias in our results is due to a calibration issue rather than a bias in the method. Hence, by simply adding 0.5 days to all estimated time delays, A was improved to 0.1%–0.6%, which is quite acceptable.

On the other hand, we noticed that our average χ² is very small in all of the rungs, which indicates that we could have a smaller P metric if we had smaller errors, while χ² could be enlarged slightly but still within an acceptable range. Hence, by dividing all of our estimated errors by a factor of $\sqrt{2}$, the average P value improved to 0.027 to ∼0.041 while the average χ² reached around 1. The calibrated results are tabulated in Table 3. These results show that with a minor modification, the algorithm could yield a reliable and precise time delay estimate for various forms of the data and in a fully automated manner.

In this paper, we present the algorithm that we used to participate in the TDC0 and TDC1 Strong Lens Time Delay Challenges (Dobler et al. 2015; Liao et al. 2015). In our algorithm, we combined various statistical methods of data analysis in order to estimate the time delay between the different light curves of the strong lens systems. At different stages of our analysis, we used iterative smoothing, cross-correlation, simulations and error estimation, bias control, and significance testing to prepare our results. Given the limited timeframe, we had to make some approximations in our error analysis. However, the novel approach we used to estimate the uncertainties has been shown to perform satisfactorily. In our approach to estimating the time delay between a pair of light curves A₁ and A₂, we first smoothed over both light curves using an iterative smoothing method (Shafieloo et al. 2006; Shafieloo 2007; Shafieloo & Clarkson 2010; Shafieloo 2012), producing the smoothed light curves A^s₁ and A^s₂. During smoothing, we recorded the ranges with no data available (which would have resulted in unreliable smoothing). The algorithm was set to detect such ranges automatically. Then, we calculated the cross-correlation between A₁ and A^s₂ and also between A₂ and A^s₁ for different time delays, and found the maximum correlations. These two maximum correlations should be for the same time delays (that is, the absolute values of the time delays should be consistent with each other). The difference between these two estimated time delays (with maximum correlations) was part of the total uncertainty considered for each pair (in the estimated time delay). To estimate the error on each derived time delay, we simulated many realizations of the data for each rung and for various time delays. Knowing the fiducial values, we derived the expected uncertainties in the estimated values of the time delays. In the case of the Quad sample, we used different combinations of smoothed and raw light curves to test the internal consistency of the results and the relative errors. These internal consistency relations were used to adjust the estimated error bars for each pair. In our final submission, we used these internal consistencies between the estimated time delays within Quad systems for our error estimates. For the bias control, since long time delays have a limited data overlap between the two light curves, we only passed systems with cross-correlations with a correlation coefficient getter than 0.6. Throughout this challenge, one of our main concerns was to achieve a high f value without any outliers. This was achieved with f > 0.3 for all five rungs. Our conservative approach yielded average χ² values of around 0.5–0.9 for different rungs with P of approximately 0.03–0.06. Since χ² and P are correlated, by simply dividing all estimated errors by a factor of $\sqrt{2}$, χ² of ∼1 and P of 0.02–0.04 could be achieved trivially. Since our method was initially calibrated only on TDC0 rung 0 data, due to a lack of time, a calibration bias of 0.5 days for all of the submissions was discovered, resulting in A ≈ 1.8%–2.5%. By adding this calibration correction of 0.5 days to all of our submissions' time delay estimates, the bias was removed, improving A substantially to only 0.1%–0.6%. It has recently been noted (Hojjati & Linder 2014) that it is very important to have a good bias control in time delay estimates in order to have a reasonable cosmology determination. In summary, our proposed method seems promising in both reliability and precision, and is automated in all steps. It is also fast with a processing time on the order of a few minutes (using a typical PC) for a pair of light curves. The method presented here is exactly the same algorithm we used to participate in the TDC challenges, and there are ways to improve our error estimation by performing appropriate simulations for each set of light curves separately. Significant modifications of this algorithm will be reported in future publications.

The authors thank Eric Linder for various discussions throughout this work. The authors thank Shi Won Ahn for computational support and Alireza Hojjati, Stephen Appleby, Hyungsuk Tak, and Tommaso Treu for their comments on our manuscript. We also thank the organizers of the Strong Lens Time Delay Challenge "Evil Team" (K. Liao, G. Dobler, C. D. Fassnacht, P. Marshall, N. Rumbaugh, T. Treu) for providing the simulated light curves (available at http://timedelaychallenge.org) and feedback during the challenges. We also thank both the Evil and Good teams of the Time Delay Challenge for useful discussions. A.A. and A.S. acknowledge support from the Korea Ministry of Education, Science and Technology, Gyeongsangbuk-Do and Pohang City for Independent Junior Research Groups at the Asia Pacific Center for Theoretical Physics. A.S. acknowledges the support of the National Research Foundation of Korea (NRF-2013R1A1A2013795).