Parameter Estimation from Time-Series Data with Correlated Errors:
A Wavelet-Based Method and its Application to Transit Light Curves

We will see shortly that some extra terms are required in Eqn. (14) for real signals with some minimum and maximum resolution. To explain those terms it is useful to describe the wavelet transform as a multiresolution analysis, in which we consider successively higher-resolution approximations of a signal. An approximation with a resolution of $2^{m}$ samples per unit time is a member of a resolution space $V_{m}$. Following Wornell (1996) we impose the following conditions:

1. 1.

   if $f(t)\in V_{m}$ then for some integer $n$, $f(t-2^{-m}n)\in V_{m}$

2. 2.

   if $f(t)\in V_{m}$ then $f(2t)\in V_{m+1}$.

The first condition requires that $V_{m}$ contain all translations (at the resolution scale) of any of its members, and the second condition ensures that the sequence of resolutions is nested: $V_{m}$ is a subset of the next finer resolution $V_{m+1}$. In this way, if $\epsilon_{m}(t)\in V_{m}$ is an approximation to the signal $\epsilon(t)$, then the next finer approxmation $\epsilon_{m+1}(t)\in V_{m+1}$ contains all the information encoded in $\epsilon_{m}(t)$ plus some additional detail $d_{m}(t)$ defined as

We may therefore build an approximation at resolution $M$ by starting from some coarser resolution $k$ and adding successive detail functions:

The detail functions $d_{m}(t)$ belong to a function space $W_{m}(t)$, the orthogonal complement of the resolution $V_{m}$.

With these conditions and definitions, the orthogonal basis functions of $W_{m}$ are the wavelet functions $\psi_{n}^{m}(t)$, obtained by translating and dilating some mother wavelet $\psi(t)$. The orthogonal basis functions of $V_{m}$ are denoted $\phi_{n}^{m}(t)$, obtained by translating and dilating a so-called “father” wavelet $\phi(t)$. Thus, the mother wavelet spawns the basis of the detail spaces, and the father wavelet spawns the basis of the resolution spaces. They have complementary characteristics, with the mother acting as a high-pass filter and the father acting as a low-pass filter.³³3More precisely, the wavelet and scaling functions considered here are “quadrature mirror filters” (Mallat 1999).

In Eqn. (16), the approximation $\epsilon_{k}(t)$ is a member of $V_{k}$, which is spanned by the functions $\phi_{n}^{k}(t)$, and $d_{m}(t)$ is a member of $W_{m}$, which is spanned by the functions $\psi_{n}^{m}(t)$. Thus we may rewrite Eqn. (16) as

The wavelet coefficients $\epsilon_{n}^{m}$ and the scaling coefficients $\bar{\epsilon}_{n}^{m}$ are given by

Eqn. (17) reduces to the wavelet-only equation (13) for the case of a continuously sampled signal $\epsilon(t)$, when we have access to all resolutions $m$ from $-\infty$ to $\infty$.⁴⁴4The signal must also be bounded in order for the approximation to the signal at infinitely coarse resolution to vanish, i.e., $\lim_{k\rightarrow-\infty}\epsilon_{k}(t)=0$.

There are many suitable choices for $\phi$ and $\psi$, differing in the tradeoff that must be made between smoothness and localization. The simplest choice is due to Haar (1910):

The left panel of Fig. 2 shows several elements of the approximation and detail bases for a Haar multiresolution analysis. The left panels of Fig. 3 illustrate a Haar multiresolution analysis for an arbitrarily chosen signal $\epsilon(t)$, by plotting both the approximations $\epsilon_{m}(t)$ and details $d_{m}(t)$ at several resolutions $m$. The Haar analysis is shown for pedagogic purposes only. In practice we found it advantageous to use the more complicated fourth-order Daubechies wavelet basis, described in the next section, for which the elements and the multiresolution analysis are illustrated in the right panels of Fig. 2-3.

Real signals are limited in resolution, leading to finite $M$ and $k$ in Eqn. (17). They are also limited in time, allowing only a finite number of translations $N_{m}$ at a given resolution $m$. Starting from Eqn. (17), we truncate the sum over $n$ and reindex the resolution sum such that the coarsest resolution is $k=1$, giving

where we have taken $t=0$ to be the start of the signal. Since there is no information on timescales smaller than $2^{-M}$, we need only consider $\epsilon_{M}(t_{i})$ at a finite set of times $t_{i}$:

Eqn. (28) is the inverse of the Discrete Wavelet Transform (DWT). Unlike the continuous transform of Eqn. (13), the DWT must include the coarsest level approximation (the first term in the preceding equation) in order to preserve all the information in $\epsilon(t_{i})$. For the Haar wavelet, the coarsest approximation is the mean value. For data sets with $N=n_{0}2^{M}$ uniformly spaced samples in time, we will have access to a maximal scale $M$, as in Eqn. (28), with $N_{m}=n_{0}2^{m-1}$.

A crucial point is the availability of the Fast Wavelet Transform (FWT) to perform the DWT (Mallat 1989). The FWT is a pyramidal algorithm operating on data sets of size $N=n_{0}2^{M}$ returning $n_{0}(2^{M}-1)$ wavelet coefficients and $n_{0}$ scaling coefficients for some $n_{0}>0$, $M>0$. The FWT is a computationally efficient algorithm that is easily implemented (Press et al. 2007) and has $O(N)$ time-complexity (Teolis 1998).

Daubechies (1988) generalized the Haar wavelet into a larger family of wavelets, categorized according to the number of vanishing moments of the mother wavelet. The Haar wavelet has a single vanishing moment and is the first member of the family. In this work we used the most compact member (in time and frequency), $\psi=_{4}\!{\cal D}$ and $\phi=_{4}\!{\cal A}$, which is well suited to the analysis of $1/f^{\gamma}$ noise for $0<\gamma<4$ (Wornell 1996). We plot ${}_{4}\!{\cal D}_{n}^{m}$ and ${}_{4}\!{\cal A}_{n}^{m}$ in the time-domain for several $n$, $m$ in Fig. 2, illustrating the rather unusual functional form of ${}_{4}\!{\cal D}$. The right panel of Fig. 3 demonstrates a multiresolution analysis using this basis. Press et al. (2007) provide code to implement the wavelet transform in this basis.

As alluded in § 3, the wavelet transform acts as a nearly diagonalizing operator for the covariance matrix in the presence of $1/f^{\gamma}$ noise. The wavelet coefficients $\epsilon_{n}^{m}$ of such a noise process are zero-mean, nearly uncorrelated random variables. Specifically, the covariance between scales $m$, $m^{\prime}$ and translations $n$, $n^{\prime}$ is (Wornell 1996, p. 65)

The wavelet basis is also convenient for the case in which the noise is modeled as the sum of an uncorrelated component and a correlated component,

where $\epsilon_{0}(t)$ is a Gaussian white noise process ($\gamma=0$) with a single noise parameter $\sigma_{w}$, and $\epsilon_{\gamma}(t)$ has ${\cal S}(f)=A/f^{\gamma}$. In the context of transit photometry, white noise might arise from photon-counting statistics (and in cases where the detector is well-calibrated, $\sigma_{w}$ is a known constant), while the $\gamma\neq 0$ term represents the “rumble” on many time scales due to instrumental, atmospheric, or astrophysical sources. For the noise process of Eqn. (30) the covariance between wavelet coefficients is

and the covariance betwen the scaling coefficients $\bar{\epsilon}_{n}^{m}$ is

where $g(\gamma)$ is a constant of order unity; for the purposes of this work $g(1)=(2\ln 2)^{-1}\approx 0.72$ (Fadili & Bullmore 2002). Eqns. (31) and (32) are the key mathematical results that form the foundation of our algorithm. For proofs and further details, see Wornell (1996).

It should be noted that the correlations between the wavelet and scaling coefficients are small but not exactly zero. The decay rate of the correlations with the resolution index depends on the choice of wavelet basis and on the spectral index $\gamma$. By picking a wavelet basis with a higher number of vanishing moments, we hasten the decay of correlations. This is why we chose the Daubechies 4th-order basis instead of the Haar basis. In the numerical experiments decribed in § 4, we found the covariances to be negligible for the purposes of parameter estimation. In addition, the compactness of the Daubechies 4th-order basis reduces artifacts arising from the assumption of a periodic signal that is implicit in the FWT.

Given an observation of noise $\epsilon(t)$ that is modeled as in Eqn. (30), we may estimate the $\gamma\neq 0$ component by rescaling the wavelet and scaling coefficients and filtering out the white component:

We may then proceed to subtract the estimate of the correlated component from the observed noise, $\epsilon_{0}(t)=\epsilon(t)-\epsilon_{\gamma}(t)$ (Wornell 1996, p. 76). In this way the FWT can be used to “whiten” the noise.

Armed with the preceding theory, we rewrite the likelihood function of Eqn. (6) in the wavelet domain. First we transform the residuals $r_{i}\equiv y_{i}-f(t_{i};\vec{p})$, giving

where $y_{n}^{m}$ and $f_{n}^{m}(\vec{p})$ are the discrete wavelet coefficients of the data and the model. Likewise, $\bar{y}_{n}^{1}$ and $\bar{f}_{n}^{1}(\vec{p})$ are the $n_{0}$ scaling coefficients of the data and the model. Given the diagonal covariance matrix shown in Eqns. (31) and (32), the likelihood ${\cal L}$ is a product of Gaussian functions at each scale $m$ and translation $n$:

where

are the variances of the wavelet and scaling coefficients respectively. For a data set with $N$ points, calculating the likelihood function of Eqn. (37) requires multiplying $N$ Gaussian functions. The additional step of computing the FWT of the residuals prior to computing ${\cal L}$ adds $O(N)$ operations. Thus, the entire calculation has a time-complexity $O(N)$.

For this calculation we must have estimators of the three noise parameters $\gamma$, $\sigma_{r}$ and $\sigma_{w}$. These may be estimated separately from the model parameters $\vec{p}$, or simultaneously with the model parameters. For example, in transit photometry, the data obtained outside of the transit may be used to estimate the noise parameters, which are then used in Eqn. (37) to estimate the model parameters. Or, in a single step we could maximize Eqn. (37) with respect to all of $\gamma$, $\sigma_{r}$, $\sigma_{w}$ and $\vec{p}$. Fitting for both noise and transit parameters simultaneously is potentially problematic, because some of the correlated noise may be “absorbed” into the choices of the transit parameters, i.e., the errors in the noise parameters and transit parameters are themselves correlated. This may cause the noise level and the parameter uncertainties to be underestimated. Unfortunately, there are many instances when one does not have enough out-of-transit data for the strict separation of transit and noise parameters to be feasible.

In practice the optimization can be accomplished with an iterative routine [such as AMOEBA, Powell’s method, or a conjugate-gradient method; see Press et al. (2007)]. Confidence intervals can then be defined by the contours of constant likelihood. Alternatively one can use a Monte Carlo Markov Chain [MCMC; see, e.g., Gregory (2005)], in which case the jump-transition likelihood would be given by Eqn. (37). The advantages of the MCMC method have led to its adoption by many investigators (see, e.g., Holman et al. 2006, Burke et al. 2007, Collier Cameron et al. 2007). For that method, computational speed is often a limiting factor, as a typical MCMC analysis involves several million calculations of the likelihood function.

Some aspects of real data do not fit perfectly into the requirements of the DWT. The time sampling of the data should be approximately uniform, so that the resolution scales of the multiresolution analysis accurately reflect physical timescales. This is usually the case for time-series photometric data. Gaps in a time series can be fixed by applying the DWT to each uninterrupted data segment, or by filling in the missing elements of the residual series with zeros.

The FWT expects the number of data points to be an integral multiple of some integral power of two. When this is not the case, the time series may be truncated to the nearest such boundary; or it may be extended using a periodic boundary condition, mirror reflection, or zero-padding. In the numerical experiments described below, we found that zero-padding has negligible effects on the calculation of likelihood ratios and parameter estimation.

The FWT generally assumes a periodic boundary condition for simplicity of computation. A side effect of this simplication is that information at the beginning and end of a time series are artificially associated in the wavelet transform. This is one reason why we chose the 4th-order Daubechies-class wavelet basis, which is well localized in time, and does not significantly couple the beginning and the end of the time series except on the coarsest scales.

In this section we consider the case in which the noise parameters $\gamma$, $\sigma_{r}$, and $\sigma_{w}$ are known with negligible error. We have in mind a situation in which a long series of out-of-transit data are available, with stationary noise properties.

We generated transit light curves with known transit parameters $\vec{p}$, contaminated by an additive combination of a white and a correlated ($1/f^{\gamma}$) noise source. Then we used an MCMC method to estimate the transit parameters and their 68.3% confidence limits. (The technique for generating noise and the MCMC method are described in detail below.) For each realization of a simulated light curve, we estimated transit parameters using the likelihood defined either by Eqn. (6) for the white analysis, or Eqn. (37) for the wavelet analysis.

For a given parameter $p_{k}$, the estimator $\hat{p}_{k}$ was taken to be the median of the values in the Markov chain and $\hat{\sigma}_{p_{k}}$ was taken to be the standard deviation of those values. To assess the results, we considered the “number-of-sigma” statistic

In words, ${\cal N}$ is the number of standard deviations separating the parameter estimate $\hat{p}_{k}$ from the true value $p_{k}$. If the error in $p_{k}$ is Gaussian, then a perfect analysis method should yield results for ${\cal N}$ with an expectation value of $0$ and variance of $1$. If we find that the variance of ${\cal N}$ is greater than one, then we have underestimated the error in $\hat{p}_{k}$ and we may attribute too much significance to the result. On the other hand, if the variance of ${\cal N}$ is smaller than one, then we have overestimated $\sigma_{p_{k}}$ and we may miss a significant discovery. If we find that the mean of ${\cal N}$ is nonzero then the method is biased.

For now, we consider only the single parameter $t_{c}$, the time of midtransit. The $t_{c}$ parameter is convenient for this analysis as it is nearly decoupled from the other transit parameters (Carter et al. 2008). Furthermore, as mentioned in the introduction, the measurement of the midtransit time cannot be improved by observing other transit events, and variations in the transit interval are possible signs of additional gravitating bodies in a planetary system.

The noise was synthesized as follows. First, we generated a sequence of $N=1024$ independent random variables obeying the variance conditions from Eqns. (31) and (32) for $1023$ wavelet coefficients over $9$ scales and a single scaling coefficient at the coarsest resolution scale. We then performed the inverse FWT of this sequence to generate our noise signal. In this way, we could select exact values for $\gamma$, $\sigma_{r}$, and $\sigma_{w}$. We also needed to find the single parameter $\sigma$ for the white-noise analysis; it is not simply related to the parameters $\gamma$, $\sigma_{r}$, and $\sigma_{w}$. In practice, we found $\sigma$ by calculating the median sample variance among $10^{4}$ unique realizations of a noise source with fixed parameters $\gamma$, $\sigma_{r}$, and $\sigma_{w}$.

For the transit model, we used the analytic formulas of Mandel & Agol (2002), with a planet-to-star ratio of $R_{p}/R_{\star}=0.15$, a normalized orbital distance of $a/R_{\star}=10$, and an orbital inclination of $i=90\arcdeg$, as appropriate for a gas giant planet in a close-in orbit around a K star. These correspond to a fractional loss of light $\delta=0.0225$, duration $T=1.68$ hr, and partial duration $\tau=0.152$ hr. We did not include the effect of limb darkening, as it would increase the computation time and has little influence on the determination of $t_{c}$ (Carter et al. 2009). Each simulated light curve spanned 3 hr centered on the midtransit time, with a time sampling of 11 s, giving 1024 uniformly spaced samples. A noise-free light curve is shown in Fig. 4.

For the noise model, we chose $\sigma_{w}=1.35\times 10^{-3}$ and $\gamma=1$, and tried different choices for $\sigma_{r}$. We denote by $\alpha$ the ratio of the rms values of the correlated noise component and the white noise component.⁵⁵5We note that although $\sigma_{w}$ is the rms of the white noise component, $\sigma_{r}$ is generally not the rms of the correlated component. The notation is unfortunate, but follows that of Wornell (1996). The example in Fig. 4 has $\alpha=1/3$. As $\alpha$ is increased from zero, the correlated component becomes more important, as is evident in the simulated data plotted in Fig. 5. Our choice of $\sigma_{w}$ corresponds to a precision of $5.8\times 10^{-4}$ per minute-equivalent sample, and was inspired by the recent work by Johnson et al. (2009) and Winn et al. (2009), which achieved precisions of $5.4\times 10^{-4}$ and $4.0\times 10^{-4}$ per minute-equivalent sample, respectively. Based on our survey of the literature and our experience with the Transit Light Curve project (Holman et al. 2006, Winn et al. 2007), we submit that all of the examples shown in Fig. 5 are “realistic” in the sense that the bumps, wiggles, and ramps resemble features in actual light curves, depending on the instrument, observing site, weather conditions, and target star.

For a given choice of $\alpha$, we made 10,000 realizations of the simulated transit light curve with $1/f$ noise. We then constructed two Monte Carlo Markov Chains for $t_{c}$ starting at the true value of $t_{c}=0$. One chain was for the white analysis, with a jump-transition likelihood given by Eqn. (6). The other chain was for the wavelet analysis, using Eqn. (37) instead. Both chains used the Metropolis-Hastings jump condition, and employed perturbation sizes such that $\approx$40% of jumps were accepted. Initial numerical experiments showed that the autocorrelation of a given Markov chain for $t_{c}$ is sharply peaked at zero lag, with the autocorrelation dropping below $0.2$ at lag-one. This ensured good convergence with chain lengths of $500$ (Tegmark et al. 2004). Chain histograms were also inspected visually to verify that the distribution was smooth. We recorded the median $\hat{t}_{c}$ and standard deviation $\hat{\sigma}_{t_{c}}$ for each chain and constructed the statistic ${\cal N}$ for each separate analysis (white or wavelet). Finally, we found the median and standard deviation of ${\cal N}$ over all 10,000 noise realizations.

Fig. 6 shows the resulting distributions of ${\cal N}$, for the particular case $\alpha=1/3$. Table 1 gives a collection of results for the choices $\alpha=0$, $1/3$, $2/3$, and $1$. The mean of ${\cal N}$ is zero for both the white and wavelet analyses: neither method is biased. This is expected, because all noise sources were described by zero-mean Gaussian distributions. However, the widths of the distributions of ${\cal N}$ show that the white analysis underestimates the error in $t_{c}$. For a transit light curve constructed with equal parts white and $1/f$ noise ($\alpha=1$), the white analysis gave an estimate of $t_{c}$ that differs from the true value by more than $1~\sigma$ nearly $80\%$ of the time. The factor by which the white analysis underestimates the error in $t_{c}$ appears to increase linearly with $\alpha$. In contrast, for all values of $\alpha$, the wavelet analysis maintains a unit variance in ${\cal N}$, as desired.

The success of the wavelet method is partially attributed to the larger (and more appropriate) error intervals that it returns for $\hat{t}_{c}$. It is also partly attributable to an improvement in the accuracy of $\hat{t}_{c}$ itself: the wavelet method tends to produce $\hat{t}_{c}$ values that are closer to the true $t_{c}$. This is shown in the final column in Table (1), where we report the percentage of cases in which the analysis method (white or wavelet) produces an estimate of $t_{c}$ that is closer to the truth. For $\alpha=1$ the wavelet analysis gives more accurate results $66\%$ of the time.

In this section we consider the case in which the noise parameters are not known in advance. Instead the noise parameters must be estimated based on the data. We did this by including the noise parameters as adjustable parameters in the Markov chains. In principle this could be done for all three noise parameters $\gamma$, $\sigma_{r}$, and $\sigma_{w}$, but for most of the experiments presented here we restricted the problem to the case $\gamma=1$. This may be a reasonable simplification, given the preponderance of natural noise sources with $\gamma=1$ (Press 1978). Some experiments involving noise with $\gamma\neq 1$ are described at the end of this section.

We also synthesized the noise with a non-wavelet technique, to avoid “stacking the deck” in favor of the wavelet method. We generated the noise in the frequency domain, as follows. We specified the amplitudes of the Fourier coefficients using the assumed functional form of the power spectral density [${\cal S}(f)\propto 1/f$], and drew the phases from a uniform distribution between $-\pi$ and $\pi$. The correlated noise in the time domain was found by performing an inverse Fast Fourier Transform. We rescaled the noise such that the rms was $\alpha$ times the specified $\sigma_{w}$. The normally-distributed white noise was then added to the correlated noise to create the total noise. This in turn was added to the idealized transit model.

For each choice of $\alpha$, we made 10,000 simulated transit light curves and analyzed them with the MCMC method described previously. For the white analysis, the mid-transit time $t_{c}$ and the single noise parameter $\sigma$ were estimated using the likelihood defined via Eqn. (6). For the wavelet analysis we estimated $t_{c}$ and the two noise parameters $\sigma_{r}$ and $\sigma_{w}$ using the likelihood defined in Eqn. (37).

Table 3 gives the resulting statistics from this experiment, in the same form as were given in Table 1 for the case of known noise parameters. (This table also includes some results from § 4.4, which examines two other methods for coping with correlated noise.) Again we find that the wavelet method produces a distribution of ${\cal N}$ with unit variance, regardless of $\alpha$; and again, we find that the white analysis underestimates the error in $t_{c}$. In this case the degree of error underestimation is less severe, a consequence of the additional freedom in the noise model to estimate $\sigma$ from the data. The wavelet method also gives more accurate estimates of $t_{c}$ than the white method, although the contrast between the two methods is smaller than it was with for the case of known noise parameters.

Our numerical results must be understood to be illustrative, and not universal. They are specific to our choices for the noise parameters and transit parameters. Via further numerical experiments, we found that the width of ${\cal N}$ in the white analysis is independent of $\sigma_{w}$, but it does depend on the time sampling. In particular, the width grows larger as the time sampling becomes finer (see Table 2). This can be understood as a consequence of the long-range correlations. The white analysis assumes that the increased number of data points will lead to enhanced precision, whereas in reality, the correlations negate the benefit of finer time sampling.

Table 4 gives the results of additional experiments with $\gamma\neq 1$. In those cases we created simulated noise with $\gamma\neq 1$ but in the course of the analysis we assumed $\gamma=1$. The correlated noise fraction was set to $\alpha=1/2$ for these tests. The results show that even when $\gamma$ is falsely assumed to be unity, the wavelet analysis still produces better estimates of $t_{c}$ and more reliable error bars than the white analysis.

Having established the superiority of the wavelet method over the white method, we wish to show that the wavelet method is also preferable to the more straightforward approach of computing the likelihood function in the time domain with a non-diagonal covariance matrix. The likelihood in this case is given by Eqn. (8).

The time-domain calculation and the use of the covariance matrix raised two questions. First, how well can we estimate the autocovariance $R(\tau)$ from a single time series? Second, how much content of the resulting covariance matrix needs to be retained in the likelihood calculation for reliable parameter estimation? The answer to the first question depends on whether we wish to utilize the sample autocorrelation as the estimator of $R(\tau)$ or instead use a parametric model (such as an ARMA model) for the autocorrelation. In either case, our ability to estimate the autocorrelation improves with number of data samples contributing to its calculation. The second question is important because retaining the full covariance matrix would cause the computation time to scale as $O(N^{2})$ and in many cases the analysis would be prohibitively slow. The second question may be reframed as: what is the minimum number of lags $L$ that needs to be considered in computing the truncated $\chi^{2}$ of Eqn. (9), in order to give unit variance in the number-of-sigma statistic for each model parameter? The time-complexity of the truncated likelihood calculation is $O(NL)$. If $L\lesssim 5$ then the time-domain method and the wavelet method may have comparable computational time-complexity, while for larger $L$ the wavelet method would offer significant advantage.

We addressed these questions by repeating the experiments of the previous sections using a likelihood function based on the truncated $\chi^{2}$ statistic. We assumed that the parameters of the noise model were known, as in § 4.1. The noise was synthesized in the wavelet domain, with $\gamma=1$, $\sigma_{w}=0.00135$, and $\alpha$ set equal to $1/3$ or $2/3$. The parameters of the transit model and the time series were the same as in § 4.1. We calculated the “exact” autocovariance function $R(l)$ at integer lag $l$ for a given $\alpha$ by averaging sample autocovariances over 50,000 noise realizations. Fig. 7 plots the autocorrelation [$R(l)/R(0)$] as a function of lag for $\alpha=1/3$, $2/3$. We constructed the stationary covariance $\Sigma_{ij}=R(|i-j|)$ and computed its inverse $(\Sigma^{-1})_{ij}$ for use in Eqn. (9).

Then we used the MCMC method to find estimates and errors for the time of midtransit, and calculated the number-of-sigma statistic ${\cal N}$ as defined in Eqn. (40). In particular, for each simulated transit light curve, we created a Markov chain of $1,000$ links for $t_{c}$, using $\chi^{2}(L)$ in the jump-transition likelihood. We estimated $t_{c}$ and $\sigma_{t_{c}}$, and calculated ${\cal N}$. We did this for $5,000$ realizations and determined $\sigma_{{\cal N}}$, the variance in ${\cal N}$, across this sample. We repeated this process for different choices of the maximum lag $L$. Fig. (8) shows the dependence of $\sigma_{{\cal N}}$ upon the maximum lag $L$.

The time-domain method works fine, in the sense that when enough non-diagonal elements in the covariance matrix are retained, the parameter estimation is successful. We find that $\sigma_{{\cal N}}$ approaches unity as $L^{-\beta}$ with $\beta=0.15$, $0.25$ for $\alpha=1/3$, $2/3$, respectively. However, to match the reliability of the wavelet method, a large number of lags must be retained. To reach $\sigma_{{\cal N}}$ = 1.05, we need $L\approx 50$ for $\alpha=1/3$ or $L\approx 75$ for $\alpha=2/3$. In our implementation, the calculation based on the truncated covariance matrix [Eqn. (9)] took $30$–$40$ times longer than the calculation based on the wavelet likelihood [Eqn. (37)].

This order-of-magnitude penalty in runtime is bad enough, but the real situation may be even worse, because one usually has access to a single noisy estimate of the autocovariance matrix. Or, if one is using an ARMA model, the estimated parameters of the model might be subject to considerable uncertainty as compared to the “exact” autocovariance employed in our numerical experiments. If it is desired to determine the noise parameters simultaneously with the other model parameters, then there is a further penalty associated with inverting the covariance matrix at each step of the calculation for use in Eqn. (9), although it may be possible to circumvent that particular problem by modeling the inverse-covariance matrix directly.

In this section we compare the results of the wavelet method to two methods for coping with correlated noise that are drawn from the recent literature on transit photometry. The first of these two methods is the “time averaging” method that was propounded by Pont et al. (2006) and used in various forms by Bakos et al. (2006), Gillon et al. (2006), Winn et al. (2007, 2008, 2009), Gibson et al. (2008), and others. In one implementation, the basic idea is to calculate the sample variance of unbinned residuals, $\hat{\sigma}_{1}^{2}$, and also the sample variance of the time-averaged residuals, $\hat{\sigma}_{n}^{2}$, where every $n$ points have been averaged (creating $m$ time bins). In the absence of correlated noise, we expect

In the presence of correlated noise, $\hat{\sigma}_{n}^{2}$ differs from this expectation by a factor $\hat{\beta}_{n}^{2}$. The estimator $\hat{\beta}$ is then found by averaging $\hat{\beta}_{n}$ over a range $\Delta n$ corresponding to time scales that are judged to be most important. In the case of transit photometry, the duration of ingress or egress is the most relevant time scale (corresponding to averaging time scales on the order of tens of minutes, in our example light curve). A white analysis is then performed, using the noise parameter $\sigma=\beta\sigma_{1}$ instead of $\sigma_{1}$. This causes the parameter errors $\hat{\sigma}_{p_{k}}$ to increase by $\beta$ but does not change the parameter estimates $\hat{p}_{k}$ themselves.⁶⁶6Alternatively one may assign an error to each data point equal to the quadrature sum of the measurement error and an extra term $\sigma_{r}$ (Pont et al. 2006). For cases in which the errors in the data points are all equal or nearly equal, these methods are equivalent. When the errors are not all the same, it is more appropriate to use the quadrature-sum approach of Pont et al. (2006). In this paper all our examples involve homogeneous errors.

A second method is the “residual permutation” method that has been used by Jenkins et al. (2002), Moutou et al. (2004), Southworth (2008), Bean et al. (2008), Winn et al. (2008), and others. This method is a variant of a bootstrap analysis, in which the posterior probability distribution for the parameters is based on the collection of results of minimizing $\chi^{2}$ (assuming white noise) for a large number of synthetic data sets. In the traditional bootstrap analysis the synthetic data sets are produced by scrambling the residuals and adding them to a model light curve, or by drawing data points at random (with replacement) to make a simulated data set with the same number of points as the actual data set. In the residual permutation method, the synthetic data sets are built by performing a cyclic permutation of the time indices of the residuals, and then adding them to the model light curve. In this way, the synthetic data sets have the same bumps, wiggles, and ramps as the actual data, but they are translated in time. The parameter errors are given by the widths of the distributions in the parameters that are estimated from all the different realizations of the synthetic data, and they are usually larger than the parameter errors returned by a purely white analysis.

As before, we limited the scope of the comparison to the estimation of $t_{c}$ and its uncertainty. We created 5,000 realizations of a noise source with $\gamma=1$ and a given value of $\alpha$ (either $0$, $1/3$, $2/3$, or $1$). We used each of the two approximate methods (time-averaging and residual-permutation) to calculate $\hat{\beta}$ and its uncertainty based on each of the 5,000 noise realizations. Then we found the median and standard deviation of $\hat{\beta}/\beta$ over all 5,000 realizations. Table (3) presents the results of this experiment.

Both methods, time-averaging and residual-permutation, gave more reliable uncertainties than the white method. However they both underestimated the true uncertainties by approximately 15-30%. Furthermore, neither method provided more accurate estimates of $t_{c}$ than did the white method. For the time-averaging method as we have implemented it, this result is not surprising, for that method differs from the white method only in the inflation of the error bars by some factor $\beta$. The parameter values that maximize the likelihood function were unchanged.

We have shown the wavelet method to work well in the presence of $1/f^{\gamma}$ noise. Although this family of noise processes encompasses a wide range of possibilities, the universe of possible correlated noise processes is much larger. In this section we test the wavelet method using simulated data that has correlated noise of a completely different character. In particular, we focus on a process with exclusively short-term correlations, described by one of the aforementioned autoregressive moving-average (ARMA) class of parametric noise models. In this way we test our method on a noise process that is complementary to the longer-range correlations present in $1/f^{\gamma}$ noise, and we also make contact between our method and the large body of statistical literature on ARMA models.

For $1/f^{\gamma}$ noise we have shown that time-domain parameter estimation techniques are slow. However, if the noise has exclusively short-range correlations, the autocorrelation function will decay with lag more rapidly than a power law, and the truncated-$\chi^{2}$ likelihood [Eqn. (9)] may become computationally efficient. ARMA models provide a convenient analytic framework for parameterizing such processes. For a detailed review of ARMA models and their use in statistical inference, see Box & Jenkins (1976). Applications of ARMA models to astrophysical problems have been described by in Koen & Lombard (1993), Konig & Timmer (1997) and Timmer et al. (2000).

To see how the wavelet method performs on data with short-range correlations we constructed synthetic transit data in which the noise is described by a single-parameter autoregressive [AR(1; $\psi$)] model. An AR(1; $\psi$) process $\epsilon(t_{i})$ is defined by the recursive relation

where $\eta(t_{i})$ is an uncorrelated Gaussian process with width parameter $\sigma$ and $\psi$ is the sole autoregressive parameter. The autocorrelation $\gamma(l)$ for an AR(1; $\psi$) process is

An AR(1;$\psi$) process is stationary so long as $0<\psi<1$ (Box & Jenkins 1976). The decay length of the autocorrelation function grows as $\psi$ is increased from zero to one. Figure (9) plots the autocorrelation function of a process that is an additive combination of an AR(1; $\psi=0.95$) process and a white noise process. The noise in our synthetic transit light curves was the sum of this AR(1; $\psi=0.95$) process, and white noise, with $\alpha=1/2$ (see Fig 9). With these choices, the white method underestimates the error in $t_{c}$, while at the same time the synthetic data look realistic.

We proceeded with the MCMC method as described previously to estimate the time of mid-transit. All four methods assessed in the previous section were included in this analysis, for comparison. Table 5 gives the results. The wavelet method produces more reliable error estimates than the white method. However, the wavelet method no longer stands out as superior to the time-averaging method or the residual-permutation method; all three of these methods give similar results. This illustrates the broader point that using any of these methods is much better than ignoring the noise correlations. The results also show that although the wavelet method is specifically tuned to deal with $1/f^{\gamma}$ noise, it is still useful in the presence of noise with shorter-range correlations.

It is beyond the scope of this paper to test the applicability of the wavelet method on more general ARMA processes. Instead we suggest the following approach when confronted with real data [see also Beran (1994)]. Calculate the sample autocorrelation, and power spectral density, based on the out-of-transit data or the residuals to an optimized transit model. For stationary processes these two indicators are related as described in § 2. Short-memory, ARMA-like processes can be identified by large autocorrelations at small lags or by finite power spectral density at zero frequency. Long-memory processes ($1/f^{\gamma}$) can be identified by possibly small but non-vanishing autocorrelation at longer lags. Processes with short-range correlations could be analyzed with an ARMA model of the covariance matrix [see Box & Jenkins (1976)], or the truncated-lag covariance matrix, although a wavelet-based analysis may be sufficient as well. Long-memory processes are best analyzed with the wavelet method as described in this paper.

It should also be noted that extensions of ARMA models have been developed to mimic long-memory, $1/f^{\gamma}$ processes. In particular, fractional autoregressive integrated moving-average models (ARFIMA) describe “nearly” $1/f^{\gamma}$ stationary processes, according to the criterion described by Beran (1994). As is the case with ARMA models, ARFIMA models enjoy analytic forms for the likelihood in the time-domain. Alas, as noted by Wornell (1996) and Beran (1994), the straightforward calculation of this likelihood is computationally expensive and potentially unstable. For $1/f^{\gamma}$ processes, the wavelet method is probably a better choice than any time-domain method for calculating the likelihood.

We present here an illustrative calculation that is relevant to the goal of detecting planets or satellites through the perturbations they produce on the sequence of midtransit times of a known transiting planet. Typically an observer would fit the midtransit times $t_{c,i}$, to a model in which the transits are strictly periodic:

for some integers $E_{i}$ and constants $t_{c,0}$ and $P$. Then, the residuals would be computed by subtracting the best-fit model from the data, and a test for anomalies would be performed by assessing the likelihood of obtaining those residuals if the linear model were correct. Assuming there are $N$ data points with normally-distributed, independent errors, the likelihood is given by a $\chi^{2}$-distribution, prob$(\chi^{2}$, $N_{\rm dof})$, where

and $N_{\rm dof}=N-2$ is the number of degrees of freedom. Values of $\chi^{2}$ with a low probability of occurrence indicate the linear model is deficient, that there are significant anomalies in the timing data, and that further observations are warranted.

We produced 10 simulated light curves of transits of the particular planet GJ 436b, a Neptune-sized planet transiting an M dwarf (Butler et al. 2004, Gillon et al. 2007) which has been the subject of several transit-timing studies (see, e.g., Ribas et al. 2008, Alonso et al. 2008, Coughlin et al. 2008). Our chosen parameters were $R_{p}/R_{\star}=0.084$, $a/R_{\star}=12.25$, $i=85.94$ deg, and $P=2.644$ d. This gives $\delta=0.007$, $T=1$ hr, and $\tau=0.24$ hr. We chose limb-darkening parameters as appropriate for the SDSS $r$ band (Claret 2004). We assumed that 10 consecutive transits were observed, in each case giving $512$ uniformly-sampled flux measurements over $2.5$ hours centered on the transit time. Noise was synthesized in the Fourier domain (as in § 4.2), with a white component $\sigma_{w}=0.001$ and a $1/f$ component with rms $0.0005$ ($\alpha=1/2$). The 10 simulated light curves are plotted in Fig. (10). Visually, they resemble the best light curves that have been obtained for this system.