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 Equation (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.  
  if f(t) ∈ V_m then for some integer n, f(t − 2^−mn) ∈ V_m,
• 2.  
  if f(t) ∈ V_m then f(2t) ∈ 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 _m(t) ∈ V_m is an approximation to the signal (t), then the next finer approximation _m+1(t) ∈ V_m+1 contains all the information encoded in _m(t) plus some additional detail d_m(t) defined as

$$\begin{equation} d_m(t) \equiv \epsilon _{m+1}(t) - \epsilon _{m}(t). \end{equation} \tag{ 15 }$$

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

$$\begin{eqnarray} \epsilon _M(t) & = & \epsilon _k(t)+\sum _{m=k}^{M} d_m(t). \end{eqnarray} \tag{ 16 }$$

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 ψ^m_n(t), obtained by translating and dilating some mother wavelet ψ(t). The orthogonal basis functions of V_m are denoted as ϕ^m_n(t), obtained by translating and dilating a so-called father wavelet ϕ(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.³

In Equation (16), the approximation _k(t) is a member of V_k, which is spanned by the functions ϕ^k_n(t), and d_m(t) is a member of W_m, which is spanned by the functions ψ^m_n(t). Thus we may rewrite Equation (16) as

$$\begin{eqnarray} \epsilon _M(t) & = & \sum _{n=-\infty }^{\infty } \bar{\epsilon }_n^k \phi _n^k (t) + \sum _{m=k}^{M} \sum _{n=-\infty }^{\infty } \epsilon _n^m \psi _n^m(t). \end{eqnarray} \tag{ 17 }$$

The wavelet coefficients ^m_n and the scaling coefficients $\bar{\epsilon }_n^m$ are given by

$$\begin{eqnarray} \epsilon _n^m &=& \int _{-\infty }^{\infty } \epsilon (t) \psi _n^m(t) dt, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} \bar{\epsilon }_n^m &=& \int _{-\infty }^{\infty } \epsilon (t) \phi _n^m(t) dt. \end{eqnarray} \tag{ 19 }$$

Equation (17) reduces to the wavelet-only Equation (13) for the case of a continuously sampled signal (t), when we have access to all resolutions m from −∞ to ∞.⁴

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

$$\begin{eqnarray} \phi (t) &=& \left\lbrace \begin{array}{@{}ll@{}} 1\;\;\; & \mbox{if $0 < t \le 1$} \\ 0\;\;\; & \mbox{otherwise} \end{array} \right., \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} \psi (t) &=& \left\lbrace \begin{array}{@{}ll@{}} 1\;\;\; & \mbox{if $-\frac{1}{2} < t \le 0$} \\ -1\;\;\; & \mbox{if $0 < t \le \frac{1}{2}$} \\ 0\;\;\; & \mbox{otherwise} \end{array} \right.. \end{eqnarray} \tag{ 21 }$$

The left panel of Figure 2 shows several elements of the approximation and detail bases for a Haar multiresolution analysis. The left panels of Figure 3 illustrate a Haar multiresolution analysis for an arbitrarily chosen signal (t), by plotting both the approximations _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 following section, for which the elements and the multiresolution analysis are illustrated in the right panels of Figures 2 and 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray} \epsilon _M(t) & = & \sum _{n=1}^{N_1} \bar{\epsilon }_n^1 \phi _n^1 (t) + \sum _{m=2}^{M} \sum _{n=1}^{N_m} \epsilon _n^m \psi _n^m(t), \end{eqnarray} \tag{ 22 }$$

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 _M(t_i) at a finite set of times t_i:

$$\begin{eqnarray} \epsilon (t_i) & = & \sum _{n=1}^{N_1} \bar{\epsilon }_n^1 \phi _n^1(t_i) + \sum _{m=2}^{M} \sum _{n=1}^{N_m} \epsilon _n^m \psi _n^m(t_i). \end{eqnarray} \tag{ 23 }$$

Equation (23) is the inverse of the DWT. Unlike the continuous transform of Equation (13), the DWT must include the coarsest level approximation (the first term in the preceding equation) in order to preserve all the information in (t_i). For the Haar wavelet, the coarsest approximation is the mean value. For data sets with N = n₀2^M uniformly spaced samples in time, we will have access to a maximal scale M, as in Equation (23), with N_m = n₀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₀2^M returning n₀(2^M − 1) wavelet coefficients and n₀ scaling coefficients for some n₀ > 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^γ noise for 0 < γ < 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 Figure 2, illustrating the rather unusual functional form of $_4 \!{\cal D}$. The right panel of Figure 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 Section 3, the wavelet transform acts as a nearly diagonalizing operator for the covariance matrix in the presence of 1/f^γ noise. The wavelet coefficients ^m_n of such a noise process are zero mean, nearly uncorrelated random variables. Specifically, the covariance between scales m, m', and translations n, n' is (Wornell 1996, p. 65)

$$\begin{eqnarray} \big\langle \epsilon _n^{m} \epsilon _{n^{\prime }}^{m^{\prime }} \big\rangle & \approx & \left(\sigma _r^22^{-\gamma m} \right) \delta _{m,m^{\prime }} \delta _{n,n^{\prime }}. \end{eqnarray} \tag{ 24 }$$

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:

$$\begin{eqnarray} \epsilon (t) & =& \epsilon _0(t) + \epsilon _{\gamma }(t), \end{eqnarray} \tag{ 25 }$$

where ₀(t) is a Gaussian white-noise process (γ = 0) with a single noise parameter σ_w, and _γ(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, σ_w is a known constant), while the γ ≠ 0 term represents the "rumble" on many timescales due to instrumental, atmospheric, or astrophysical sources. For the noise process of Equation (25), the covariance between wavelet coefficients is

$$\begin{eqnarray} \big\langle \epsilon _n^{m} \epsilon _{n^{\prime }}^{m^{\prime }} \big\rangle & \approx & \left(\sigma _r^2 2^{-\gamma m} + \sigma _w^2 \right) \delta _{m,m^{\prime }} \delta _{n,n^{\prime }}, \end{eqnarray} \tag{ 26 }$$

and the covariance between the scaling coefficients $\bar{\epsilon }_n^m$ is

$$\begin{eqnarray} \big\langle \bar{\epsilon }_n^m \bar{\epsilon } _n^m\big\rangle & \approx & \sigma _r^22^{-\gamma m} g(\gamma) + \sigma _w^2, \end{eqnarray} \tag{ 27 }$$

where g(γ) is a constant of order unity; for the purposes of this work, g(1) = (2ln 2)⁻¹ ≈ 0.72 (Fadili & Bullmore 2002). Equations (26) and (27) 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 a wavelet basis and on the spectral index γ. 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 fourth-order basis instead of the Haar basis. In the numerical experiments described in Section 4, we found the covariances to be negligible for the purposes of parameter estimation. In addition, the compactness of the Daubechies fourth-order basis reduces artifacts arising from the assumption of a periodic signal that is implicit in the FWT.

Given an observation of noise (t) that is modeled as in Equation (25), we may estimate the γ ≠ 0 component by rescaling the wavelet and scaling coefficients, and filtering out the white component:

$$\begin{equation} \epsilon _{\gamma }(t) = \sum _{n=1}^{N_1} \left(\frac{\sigma _r^2 2^{-\gamma } g(\gamma)}{\sigma _r^2 2^{-\gamma }g(\gamma)+\sigma _w^2}\right)\bar{\epsilon }_n^1 \phi _n^1(t) \end{equation} \tag{ 28 }$$

$$\begin{equation} +\sum _{m=2}^{M} \sum _{n=1}^{N_m} \left(\frac{\sigma _r^2 2^{-\gamma m}}{\sigma _r^2 2^{-\gamma m}+\sigma _w^2}\right) \epsilon _n^m\psi _n^m(t). \end{equation} \tag{ 29 }$$

We may then proceed to subtract the estimate of the correlated component from the observed noise, ₀(t) = (t) − _γ(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 Equation (6) in the wavelet domain. First we transform the residuals $r_i \equiv y_i-f(t_i; \vec{p})$, giving

$$\begin{eqnarray} r_n^m & = & y_n^m-f_n^m(\vec{p}) = \epsilon _{\gamma,n}^m+\epsilon _{0,n,}^m \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} \bar{r}_n^1 & = & \bar{y}_n^1-\bar{f}_n^1(\vec{p}) = \bar{\epsilon }_{\gamma,n}^1+\bar{\epsilon }_{0,n,}^1 \end{eqnarray} \tag{ 31 }$$

where y^m_n 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₀ scaling coefficients of the data and the model. Given the diagonal covariance matrix shown in Equations (26) and (27), the likelihood ${\cal L}$ is a product of Gaussian functions at each scale m and translation n:

$$\begin{eqnarray} {\cal L} &=& \left\lbrace \prod _{m=2}^{M} \prod _{n=1}^{n_0 2^{m-1}} \frac{1}{\sqrt{2 \pi \sigma _W^2} }\exp \left[ -\frac{\left(r_n^m \right)^2}{2 \sigma _W^2}\right] \right\rbrace \nonumber \\ &&\times \left\lbrace \prod _{n=1}^{n_0} \frac{1}{\sqrt{2 \pi \sigma _S^2} }\exp \left[ -\frac{\left(\bar{r}_n^1\right)^2}{2 \sigma _S^2}\right] \right\rbrace, \end{eqnarray} \tag{ 32 }$$

where

$$\begin{eqnarray} \sigma _W^2 & = & \sigma _r^22^{-\gamma m}+\sigma _w^2, \end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray} \sigma _S^2 & = & \sigma _r^22^{-\gamma } g(\gamma)+\sigma _w^2, \end{eqnarray} \tag{ 34 }$$

are the variances of the wavelet and scaling coefficients, respectively. For a data set with N points, calculating the likelihood function of Equation (32) 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 γ, σ_r, and σ_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 Equation (32) to estimate the model parameters; or, in a single step we could maximize Equation (32) with respect to all of γ, σ_r, σ_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 Equation (32). 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 2. 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 simplification is that the beginning and end of a time series are artificially associated with each other. This is one reason why we chose the fourth-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 γ, σ_r, and σ_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^γ) 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 Equation (6) for the white analysis, or Equation (32) 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:

$$\begin{eqnarray} {\cal N} & \equiv (\hat{p}_k-p_k)/\hat{\sigma }_{p_k}. \end{eqnarray} \tag{ 35 }$$

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 1, 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 1, 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 mid-transit. 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 mid-transit 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 Equations (26) and (27) for 1023 wavelet coefficients over nine 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 γ, σ_r, and σ_w. We also needed to find the single parameter σ for the white-noise analysis; it is not simply related to the parameters γ, σ_r, and σ_w. In practice, we found σ by calculating the median sample variance among 10⁴ unique realizations of a noise source with fixed parameters γ, σ_r, and σ_w.

For the transit model, we used the analytic formulae of Mandel & Agol (2002), with a planet-to-star ratio of R_p/R_⋆ = 0.15, a normalized orbital distance of a/R_⋆ = 10, and an orbital inclination of i = 90°, as appropriate for a gas giant planet in a close-in orbit around a K star. These correspond to a fractional loss of light δ = 0.0225, duration T = 1.68 hr, and partial duration τ = 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. 2008). Each simulated light curve spanned 3 hr centered on the mid-transit time, with a time sampling of 11 s, giving 1024 uniformly spaced samples. A noise-free light curve is shown in Figure 4.

Zoom In Zoom Out Reset image size

For the noise model, we chose σ_w = 1.35 × 10⁻³ and γ = 1, and tried different choices for σ_r. We denote by α the ratio of the rms values of the correlated noise component and the white-noise component.⁵ The example in Figure 4 has α = 1/3. As α is increased from zero, the correlated component becomes more important, as is evident in the simulated data plotted in Figure 5. Our choice of σ_w corresponds to a precision of 5.8 × 10⁻⁴ 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 × 10⁻⁴ and 4.0 × 10⁻⁴ 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 Figure 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.

Zoom In Zoom Out Reset image size

For a given choice of α, we made 10,000 realizations of the simulated transit light curve with 1/f noise. We then constructed two MCMCs 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 Equation (6). The other chain was for the wavelet analysis, using Equation (32) instead. Both chains used the Metropolis–Hastings jump condition, and employed perturbation sizes such that ≈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.

Figure 6 shows the resulting distributions of ${\cal N}$, for the particular case α = 1/3. Table 1 gives a collection of results for the choices α = 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 (α = 1), the white analysis gave an estimate of t_c that differs from the true value by more than 1σ nearly 80% of the time. The factor by which the white analysis underestimates the error in t_c appears to increase linearly with α. In contrast, for all values of α, the wavelet analysis maintains a unit variance in ${\cal N}$, as desired.

Zoom In Zoom Out Reset image size

Table 1. Estimates of Mid-transit Time, t_c, from Data with Known Noise Properties

Method	α	$\langle \hat{\sigma }_{t_c} \rangle$ (s)	$\langle {\cal N} \rangle$	$\sigma _{{\cal N}}$	Prob$({\cal N} > 1)\; (\%)$	Prob(best)a (%)
White	0	4.1	+0.004	0.95	29	50
	1/3	4.3	−0.005	1.93	61	39
	2/3	5.0	+0.005	3.04	75	35
	1	5.9	−0.036	3.82	79	34
Wavelet	0	4.0	+0.005	0.95	29	50
	1/3	7.2	−0.004	0.93	28	61
	2/3	11.5	−0.004	0.94	28	65
	1	16.0	−0.001	0.95	29	66

Note. ^aThe probability that the analysis method (white or wavelet) returns an estimate of t_c that is closer to the true value than the other method.

Download table as:  ASCIITypeset image

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 α = 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 γ, σ_r, and σ_w, but for most of the experiments presented here we restricted the problem to the case γ = 1. This may be a reasonable simplification, given the preponderance of natural noise sources with γ = 1 (Press 1978). Some experiments involving noise with γ ≠ 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 PSD [${\cal S}(f) \propto 1/f$], and drew the phases from a uniform distribution between −π and π. 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 α times the specified σ_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 α, 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 σ were estimated using the likelihood defined via Equation (6). For the wavelet analysis, we estimated t_c and the two noise parameters σ_r and σ_w using the likelihood defined in Equation (32).

Table 2 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 Section 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 α, 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 σ 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.

Table 2. Estimates of t_c from Data with Unknown Noise Properties

Method	α	$\langle \hat{\sigma }_{t_c} \rangle$ (s)	$\langle {\cal N} \rangle$	$\sigma _{{\cal N}}$	Prob$({\cal N} > 1)\; (\%)$	Prob(better)a (%)
White	0	4.0	−0.011	0.97	31	⋅⋅⋅
	1/3	4.2	+0.010	1.70	57	⋅⋅⋅
	2/3	4.9	+0.012	2.69	73	⋅⋅⋅
	1	5.8	+0.023	3.28	78	⋅⋅⋅
Wavelet	0	4.5	−0.009	0.90	26	50
	1/3	6.9	−0.003	1.03	33	56
	2/3	11.2	−0.005	1.07	35	57
	1	15.7	−0.007	1.09	36	57
Time-averaging	0	4.4	−0.006	0.88	26	50
	1/3	6.8	+0.009	1.15	36	50
	2/3	11.6	−0.012	1.24	40	50
	1	17.6	+0.007	1.21	38	50
Residual-permutation	0	3.5	−0.012	1.16	37	50
	1/3	6.6	+0.013	1.24	37	50
	2/3	11.8	−0.014	1.28	38	49
	1	17.3	+0.008	1.30	38	48

Note. ^aThe probability that the analysis method returns an estimate of t_c that is closer to the true value than the white analysis.

Download table as:  ASCIITypeset image

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 σ_w, but it does depend on the time sampling. In particular, the width grows larger as the time sampling becomes finer (see Table 3). 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 3. Effect of Time Sampling on the White Analysis

Na	Cadence (s)	$\sigma _{{\cal N}}$
256	42.2	1.72
512	21.1	2.04
1024	10.5	2.69
2048	5.27	3.49
4096	2.63	4.39

Note. ^aThe number of samples in a 3 hr interval.

Download table as:  ASCIITypeset image

Table 4 gives the results of additional experiments with γ ≠ 1. In those cases, we created simulated noise with γ ≠ 1 but in the course of the analysis we assumed γ = 1. The correlated noise fraction was set to α = 1/2 for these tests. The results show that even when γ 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.

Table 4. Estimates of t_c from Data with Unknown Noise Properties

Method	γa	$\langle \hat{\sigma }_{t_c} \rangle$ (s)	$\langle {\cal N} \rangle$	$\sigma _{{\cal N}}$	Prob$({\cal N} > 1)\; (\%)$	Prob(best)b (%)
White	0.5	4.5	−0.025	1.34	47	50
	1.5	4.6	+0.020	3.10	77	32
Wavelet	0.5	6.7	−0.021	0.97	30	50
	1.5	6.9	+0.002	1.17	39	68

Notes. ^aThe spectral exponent of the PSD, S(f) ∝ 1/f^γ. ^bThe probability that the analysis method (white or wavelet) returns an estimate of t_c that is closer to the true value than the other method.

Download table as:  ASCIITypeset image

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 Equation (8).

The time-domain calculation and the use of the covariance matrix raised two questions. First, how well can we estimate the autocovariance R(τ) 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(τ) 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²) 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 χ² of Equation (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 ≲ 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 a significant advantage.

We addressed these questions by repeating the experiments of the previous sections using a likelihood function based on the truncated χ² statistic. We assumed that the parameters of the noise model were known, as in Section 4.1. The noise was synthesized in the wavelet domain, with γ = 1, σ_w = 0.00135, and α set equal to 1/3 or 2/3. The parameters of the transit model and the time series were the same as in Section 4.1. We calculated the "exact" autocovariance function R(l) at integer lag l for a given α by averaging sample autocovariances over 50,000 noise realizations. Figure 7 plots the autocorrelation [R(l)/R(0)] as a function of lag for α = 1/3, 2/3. We constructed the stationary covariance Σ_ij = R(|i − j|) and computed its inverse (Σ⁻¹)_ij for use in Equation (9).

Zoom In Zoom Out Reset image size

Then we used the MCMC method to find estimates and errors for the time of mid-transit, and calculated the number-of-sigma statistic ${\cal N}$ as defined in Equation (35). In particular, for each simulated transit light curve, we created a Markov chain of 1000 links for t_c, using χ²(L) in the jump-transition likelihood. We estimated t_c and $\sigma _{t_c}$, and calculated ${\cal N}$. We did this for 5000 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. Figure 8 shows the dependence of $\sigma _{{\cal N}}$ upon the maximum lag L.

Zoom In Zoom Out Reset image size

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^−β with β = 0.15, 0.25 for α = 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 ≈ 50 for α = 1/3 or L ≈ 75 for α = 2/3. In our implementation, the calculation based on the truncated covariance matrix (Equation (9)) took 30–40 times longer than the calculation based on the wavelet likelihood (Equation (32)).

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 Equation (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

$$\begin{eqnarray} \hat{\sigma }_{n}^2 &=&\frac{ \hat{\sigma }_{1}^2}{n} \left(\frac{m}{m-1}\right). \end{eqnarray} \tag{ 36 }$$

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 Δn corresponding to timescales that are judged to be most important. In the case of transit photometry, the duration of ingress or egress is the most relevant timescale (corresponding to averaging timescales on the order of tens of minutes, in our example light curve). A white analysis is then performed, using the noise parameter σ = βσ₁ instead of σ₁. This causes the parameter errors $\hat{\sigma }_{p_k}$ to increase by β but does not change the parameter estimates $\hat{p}_k$ themselves.⁶

A second method is the "residual-permutation" method that has been used by Jenkins et al. (2002), 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 χ² (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 5000 realizations of a noise source with γ = 1 and a given value of α (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 5000 noise realizations. Then we found the median and standard deviation of $\hat{\beta }/\beta$ over all 5000 realizations. Table 2 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 β. 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^γ 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 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^γ noise, and we also make contact between our method and the large body of statistical literature on ARMA models.

For 1/f^γ 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-χ² likelihood (Equation (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; ψ)) model. An AR(1; ψ) process (t_i) is defined by the recursive relation

$$\begin{eqnarray} &&\epsilon (t_i) = \eta (t_i)+\psi \epsilon (t_{i-1}), \end{eqnarray} \tag{ 37 }$$

where η(t_i) is an uncorrelated Gaussian process with the width parameter σ and ψ is the sole autoregressive parameter. The autocorrelation γ(l) for an AR(1; ψ) process is

$$\begin{eqnarray} &&\gamma (l) = \frac{\sigma ^2}{1-\psi ^2} \psi ^l. \end{eqnarray} \tag{ 38 }$$

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

Zoom In Zoom Out Reset image size

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^γ noise, it is still useful in the presence of noise with shorter-range correlations.

Table 5. Estimates of t_c from Data with Autoregressive Correlated Noise

Method	$\langle \hat{\sigma }_{t_c} \rangle$ (s)	$\langle {\cal N} \rangle$	$\sigma _{{\cal N}}$	Prob$({\cal N} > 1)\; (\%)$	Prob(better)a (%)
White	4.5	−0.010	2.50	70	...
Wavelet	8.7	−0.016	1.33	44	51
Time-averaging	9.9	−0.010	1.25	40	49
Residual-permutation	10.2	−0.010	1.23	38	51

Note. ^aThe probability that the analysis method returns an estimate of t_c that is closer to the true value than the white analysis.

Download table as:  ASCIITypeset image

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 PSD, 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 Section 2. Short-memory, ARMA-like processes can be identified by large autocorrelations at small lags or by finite PSD at zero frequency. Long-memory processes (1/f^γ) 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^γ processes. In particular, fractional autoregressive integrated moving-average models (ARFIMA) describe "nearly" 1/f^γ 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^γ 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 mid-transit times of a known transiting planet. Typically, an observer would fit the mid-transit times t_c,i, to a model in which the transits are strictly periodic:

$$\begin{eqnarray} t_{c,i} & = & t_{c, 0}+E_i P, \end{eqnarray} \tag{ 39 }$$

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 χ²-distribution, prob(χ², N_dof), where

$$\begin{eqnarray} &&\chi ^2 = \sum _i \left[ \frac{ t_{c,i} - (t_{c,0} + E_i P)}{\sigma _{t_{c,i}}} \right]^2 \end{eqnarray} \tag{ 40 }$$

and N_dof = N − 2 is the number of degrees of freedom. Values of χ² with a low probability of occurrence indicate that 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_⋆ = 0.084, a/R_⋆ = 12.25, i = 8594, and P = 2.644 d. This gives δ = 0.007, T = 1 hr, and τ = 0.24 hr. We chose limb-darkening parameters as appropriate for the Sloan Digital Sky Survey (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 hr centered on the transit time. Noise was synthesized in the Fourier domain (as in Section 4.2), with a white component σ_w = 0.001 and a 1/f component with an rms of 0.0005 (α = 1/2). These 10 simulated light curves are plotted in Figure 10. Visually, they resemble the best light curves that have been obtained for this system.

Zoom In Zoom Out Reset image size

To estimate the mid-transit time of each simulated light curve, we performed a wavelet analysis and a white analysis, allowing only the mid-transit time and the noise parameters to vary while fixing the other parameter values at their true values. We used the same MCMC technique that was described in Section 4.2. Each analysis method produced a collection of 10 mid-transit times and error bars. These 10 data points were then fitted to the linear model of Equation (39). Figure 11 shows the residuals of the linear fit (observed − calculated). Table 6 gives the best-fit period for each analysis (wavelet or white), along with the associated values of χ².

Zoom In Zoom Out Reset image size

Table 6. Linear Fits to Estimated Mid-transit Times

Method	Fitted period/True period	$\hat{\chi }^2/N_{\rm dof}$	Prob$(\chi ^2 < \hat{\chi }^2)\; (\%)$
White	1.00000071 ± 0.00000043	2.25	98
Wavelet	1.00000048 ± 0.00000077	0.93	51

Download table as:  ASCIITypeset image

As was expected from the results of Section 4.2, the white analysis gave error bars that are too small, particularly for epochs 4 and 7. As a result, the practitioner of the white analysis would have rejected the hypothesis of a constant orbital period with 98% confidence. In addition, the white analysis gave an estimate for the orbital period that is more than 1σ away from the true value, which might have complicated the planning and execution of future observations. The wavelet method, in contrast, neither underestimated nor overestimated the errors, giving χ² ≈ N_dof in excellent agreement with the hypothesis of a constant orbital period. The wavelet method also gave an estimate for the orbital period within 1σ of the true value.

Thus far we have focused exclusively on the determination of the mid-transit time, in the interest of simplicity. However, there is no obstacle to using the wavelet method to estimate multiple parameters, even when there are strong degeneracies among them. In this section, we test and illustrate the ability of the wavelet method to solve for all the parameters of a transit light curve, along with the noise parameters.

We modeled the transit as in Sections 4.1 and 4.2. The noise was synthesized in the frequency domain (as in Section 4.2), using σ_w = 0.0045, γ = 1, and α = 1/2. The resulting simulated light curve is the upper time series in Figure 12. We used the MCMC method to estimate the transit parameters {R_p/R_⋆, a/R_⋆, i, t_c} and the noise parameters {σ_r, σ_w} (again fixing γ = 1 for simplicity). The likelihood was evaluated with either the wavelet method (Equation (32)) or the white method (Equation (6)).

Zoom In Zoom Out Reset image size

Figure 13 displays the results of this analysis in the form of the posterior distribution for the case of t_c, and the joint posterior confidence regions for the other cases. The wavelet method gives larger (and more appropriate) confidence regions than the white analysis. In accordance with our previous findings, the white analysis underestimates the error in t_c and gives an estimate of t_c that differs from the true value by more than 1σ. The wavelet method gives better agreement. Both analyses give an estimate for R_p/R_⋆ that is smaller than the true value of 0.15, but in the case of the white analysis, this shift is deemed significant, thereby ruling out the correct answer with more than 95% confidence. In the wavelet analysis, the true value of R_p/R_⋆ is well within the 68% confidence region. Both the wavelet and white analyses give accurate values of a/R_⋆ and the inclination, and the wavelet method reports larger errors. As shown in Figure 13, the wavelet method was successful at identifying the parameters (α and σ_w) of the underlying 1/f noise process.

Zoom In Zoom Out Reset image size

Figure 12 shows the best-fitting transit model, and also illustrates the action of the "whitening" filter that was described in Section 3.4. The jagged line plotted over the upper time series is the best estimate of the 1/f contribution to the noise, found by applying the whitening filter (Equation (29)) to the data using the estimated noise parameters. The lower time series is the whitened data, in which the 1/f component has been subtracted. Finally, in Figure 14 we compare the estimated 1/f noise component with the actual 1/f component used to generate the data. Possibly, by isolating the correlated component in this way, and investigating its relation to other observable parameters, the physical origin of the noise could be identified and understood.

Zoom In Zoom Out Reset image size