Improving the Lomb–Scargle Periodogram with the Thomson Multitaper

A common approach for characterizing the properties of time-series data that are evenly sampled in time is to estimate the power spectrum of the data using the periodogram. The periodogram as an estimator of the spectrum is (1) statistically inconsistent (i.e., its variance does not go to zero as infinite data are collected), (2) biased for finite samples, and (3) suffers from spectral leakage. In astronomy, time-series data are often unevenly sampled in time, and it is popular to use the Lomb–Scargle (LS) periodogram to estimate the spectrum. Unfortunately, from a statistical standpoint, the LS periodogram suffers from the same issues as the classical periodogram and has even worse spectral leakage. Here, we present an improvement on the LS periodogram by combining it with the Thomson multitaper approach. The multitaper spectral estimator is well established in the statistics and engineering literature for evenly sampled time series. It is attractive because it directly trades off bias and variance for frequency resolution, and is fast to compute: compared to an untapered spectral estimator, the multitaper adds no more than a couple of seconds for a time series with a million data points on a current desktop computer. Here, we describe an estimator that combines the multitaper with the LS periodogram. We show examples in which this new approach has improved properties compared to traditional approaches in the case of unevenly sampled time series. Finally, we demonstrate an application of the method to astronomy with an application to Kepler data.


Introduction
The analysis of time-series data is becoming increasingly common in astrophysics research, and will continue to grow in prevalence as data from telescopes such as the National Science Foundation Vera C.Rubin Telescope (LSST Science Collaboration et al. 2009) comes online in the 2020s. Data from spacecraft and ground-based telescopes such as Kepler (Koch et al. 2010), the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015), and the Zwicky Transient Facility (ZTF; Bellm et al. 2018;Graham et al. 2019) have and are improving our understanding of exoplanets and stars, but we are limited by issues such as unevenly sampled data, measurement uncertainties, and the statistical methods we employ. This paper is most concerned with the latter, and seeks to provide an improved statistical method for identifying and estimating periodic signals in time-series data.
Various methods exist for analyzing data sampled in time, but a common approach is called spectral analysis. Spectral analysis is used to estimate and characterize the process underlying timeseries data by expressing power as a function of frequency, and obtaining an estimate of the frequency spectrum or power spectrum. When performing spectral analysis, we obtain from the data a statistical estimate of the true power spectrum-we will never know to true power spectrum, but can hope to estimate it well using reliable methods.
Most spectrum estimators (e.g., the periodogram) depend heavily on the assumption that the data are regularly sampled in time. In astronomy, we do not always have the luxury of wellcontrolled observation times, and so time-series data are almost always irregularly sampled to some degree. Thus, the Lomb-Scargle (LS) periodogram (Lomb 1976;Scargle 1982), which accounts for unevenly sampled data, is a popular method for estimating the frequency spectra of astronomical time-series data.
The application of the LS estimator has become very popular in astronomy, despite its many drawbacks (see VanderPlas 2018 for a thorough discussion about LS). One method for spectral analysis that is widely used in other areas of science but less so in astronomy is the multitaper statistic invented by Thomson (1982). 5 One reason for this may be that the Thomson multitaper statistic, like many time-series analysis methods, relies on evenly sampled data. Closely related fields of helioseismology, solar physics, and geophysics have seen the use of the multitaper statistic (e.g., Thomson et al. 1996;Komm et al. 1999; Thomson & Vernon 2015, 2016. In this paper, we propose the multitaper-Lomb-Scargle (MTLS) statistic: a method for time-series analysis that combines the multitaper approach with the LS technique. The method is a rapid computation of the multitaper spectral statistic (Thomson 1982) under irregular sampling. We foresee this method as an easily implemented improvement to the most common approach to spectrum estimation when times are known but irregular.
We aim to make this paper accessible to those relatively new to time-series analysis and/or statistical notation for time-series analysis. If you are already well versed in both, then you may want to skip ahead to Section 4. The paper is organized as follows.
Section 2: a short introduction to mathematical notation for time-series analysis and a description of important assumptions and concepts such as stationarity, the autocovariance function, and the spectral density function.
Section 3: a review of the classical periodogram and its statistical issues, which necessitate the use of tapers.
Section 4: description of the discrete prolate spheroidal sequences (DPSS) as tapers, which are used in the multitaper spectral estimator.
Section 5: a well-known example that demonstrates the superior estimate of the frequency power spectrum provided by the multitaper statistic over the periodogram.
Section 6: issues with interpolating missing data in a time series.
Section 7: the MTLS statistic; the LS periodogram is briefly reviewed (Section 7.1) and then combined with the multitaper spectral statistic (Section 7.2).
Section 8: the MTLS is used in a preliminary analysis of time-series data from a red giant (RG) star as measured by Kepler.
Much of the work presented in this paper is part of Springford (2017), which also contains more details pertaining to time-series analysis with irregular and latent times (i.e., when there is a proxy for time).

Notation for Time-series Analysis
In time-series analysis, we consider a sequence of random variables X t i collected at the set of times , where we assume the set of times  are strictly increasing. Next, assume that the random variables have a joint distribution Here, t and r are two different times from , and E[ ] denotes the expected value of whatever is inside the square brackets. 6 A time series is a realization x drawn from X, at a set of times . For example, x might be the measurements of the flux of a star collected at a set of times .
A time series is regularly (or evenly) sampled if for all i, where c is a constant. Regularly sampled time series are the most common in other scientific disciplines, particularly when the observation times can be controlled by the observer.
Of primary interest in the analysis of time series is the covariance of the random variables {X t }. This differentiates time-series analysis from much of statistics in which random variables are often assumed to be independent or conditionally independent (spatial statistics has a similar interest in covariance structure). It also presents a problem for inference because typically we only obtain a single time series-that is, a single sample-making it impossible to characterize the covariance without simplifying assumptions about the covariance structure. Easily the most common assumption is that of stationarity, of which there are two flavors: 1. strong stationarity assumes that any finite collection of random variables ) if the times are continuous, and Î  s for integer-valued discrete times ( is the set of integers); and 2. weak or second-order stationarity assumes that m = E X t [ ] for every t and that the covariance ) does not depend on t, assuming that the variance of each X t is finite.
If {X t } follow a multivariate Gaussian distribution of finite dimension, weak stationarity implies strong stationarity because the multivariate Gaussian is uniquely determined by its mean and covariance. In practice, it is impossible to verify these stationarity assumptions, but some natural systems seem more likely than others to be non-stationary. For example, pulsars often exhibit obvious, predictable, and regular variation in flux over several years. On the other hand, the flux measured from a supernova changes quite a lot depending on when you start measuring! Assuming weak stationarity, we are interested in the autocovariance function m m We will assume for now that we have = ¼ - N 0, 1, 2, , 1, where N is the length of the time series. These time values will make the mathematics simpler compared to regular Julian dates, for example.
Under mild regularity conditions, the spectral density function of frequency f can be represented as which provides another approach to statistical inference via the spectral density or spectrum of the time series. The autocovariance and the spectrum are Fourier transform pairs, so that given the spectrum, the autocovariance can be represented as This is the Einstein-Wiener-Khintchine theorem (Einstein 1914;Wiener 1930;Khintchine 1934). The spectrum has the following important properties: The spectrum is non-negative The spectrum is periodic and has period equal to one = - The spectrum is an even function for real-valued x 6 Also known as the expectation in physics and astronomy literature.
The equivalence of the autocovariance R(s) and the spectrum S( f ) means that either can be used for inference. However, the spectrum has an arguably simpler interpretation and provides the opportunity to discover something unexpected (Tukey & Hamming 1949). The spectrum is precisely tuned to detect periodicities.
Another strong point in favor of the spectrum compared to the autocovariance is that at suitably spaced frequencies, the distribution of spectrum estimates are or can be made approximately independent. This independence simplifies interpretation and inference considerably.

Statistical Properties of the Periodogram
A naïve estimator of the spectrum, known as the periodogram (Schuster 1898), is Recall that N is the length of the time series.) The estimator in Equation (4) assumes we have a time series x, a realization of X, observed at a set of known and discrete times. It also assumes that the observation times are regularly spaced with = ¼ - N 0, 1, 2, , 1. The periodogram estimator of the spectrum, S f P ( ), has some problems. It is inconsistent: the variance of the estimator does not tend to zero as N goes to infinity (first proved by Rayleigh 1903). In fact, the variance of the estimator does not decrease as the sample size increases.
The spectral window for the periodogram is Fejér's kernel = . This kernel has side lobes that contribute to the spectrum estimate at frequencies away from a given frequency, an effect known as spectral leakage. Spectral leakage can cause the periodogram to be badly biased (Thomson 1982).
Both the inconsistency and bias of the periodogram estimator can be displayed through an example first presented in Box & Jenkins (1970) and later in Percival & Walden (1993). The example begins by considering a time series generated by an autoregressive process, of order p=4, with coefficients j 1 =2.7607, j 2 =−3.8106, j 3 =2.6535, and j 4 =−0.9238. This autoregressive process (denoted AR(4)) is used to generate the time series shown in Figure 1. The evenly sampled time series is then used with the classical periodogram estimator to estimate the spectral density (or power spectrum). Figure 2 shows the power spectrum as estimated by the periodogram (blue curve) for the AR(4) time series of length N=1000 (top) and N=50,000 (bottom), and compares it to the true power spectrum that is known analytically (the red curve). Note that the classical periodogram estimator is biased at high frequency for the time series of length N=1000 (top). The spectrum estimate for the time series of length N=50,000 (bottom) is less biased, but the power spectrum remains statistically inconsistent-the variance does not go to zero as N increases.
These statistical issues with the periodogram are well known and so estimators that taper the data are preferred (Brillinger 1981). The tapered estimator of the spectrum, known as a direct estimate, is where h t is the data taper, which notably has the property that . Examining Equation (6), the taper reweights the time-series values, typically down-weighting values at the beginning and end of the time series. In the frequency domain, this results in a spectral window with lower side lobes than Figure 1. A single, evenly sampled time series generated from an AR(4) process with coefficients j 1 =2.7607, j 2 =−3.8106, j 3 =2.6535, and j 4 =−0.9238. Figure 2. Power spectrum estimate from the classical periodogram (blue curves) and the true power spectral density (red dashed curve) for an AR(4) process time series of length N = 1000 (top) shown in Figure 1. The bottom figure is the same but for N = 50,000. Note that the classical periodogram estimate of the power spectrum is biased at high frequency for finite N, and is inconsistent.
There are many kinds of data tapers h t (sometimes called windows) that can be used to mitigate the spectral leakage problem (Harris 1978). The family of cosine tapers are popular (of which the Hann window is a special case), as is the modified cosine taper or Hamming window. 7 In the next section, we describe the particular set of tapers used by the Thomson (1982) multitaper statistic and why they are an optimal choice.

Slepian Tapers
Thomson's multitaper statistic makes use of not one taper, but multiple tapers. Specifically, it uses a special set of sequences called the DPSS (colloquially referred to as the Slepians in reference to Slepian 1978).
The Slepians are usually denoted as v N W , , where k is the order of the sequence. They are sequences of length N (i.e., the length of the time series) that have the smallest possible spectral leakage outside of a defined frequency band of width 2W. As an example, the Slepian sequences of order k=0 through k=3 for a time series of length N=60 and frequency band of width 2W=1.2 are shown in Figure 3.
The higher-order DPSS can also have nearly minimal spectral leakage. In fact, the first K NW 2 ⪅ have nearly minimal spectral leakage, so that for any k will provide an estimate that is characteristically similar (note here that we have replaced Slepian's notation with an equivalent shorthand: v t k ( ) ). Moreover, the estimates S f k mt ( ) ( ) are approximately uncorrelated because the DPSS tapers are orthogonal (Thomson 1982). This suggests averaging the to not only control spectral leakage, but also to reduce the variance of the combined estimate.
The result is the multitaper spectrum estimate where K is the maximum order of the DPSS taper chosen in forming the estimate. An adaptively weighted average can be used to further improve the in-band concentration (Thomson 1982). It is also common to play it safe with respect to spectral leakage and use = - (Thomson 1982).
By using multiple tapers, the Thomson multitaper spectrum estimate reduces bias, reduces variance, and reduces spectral leakage. In the next section, we continue the example from Percival & Walden (1993) that helped bolster the application of the multitaper estimator in other disciplines.

A Motivating Example for the Multitaper Spectral Estimator
Continuing with the example in Section 3, we now apply the multitaper spectrum estimator to the same time series (AR(4) process) presented in Figure 1. Figure 4 shows the true spectrum for the AR(4) process (red curve), the spectral estimate from the periodogram (blue curve), and the spectral estimate from the Thomson multitaper approach (black curve; this example appears in Percival & Walden 1993). The Thomson multitaper estimator clearly reduces bias and variance, and accurately predicts the true power spectrum. Note that one way to reduce the variance of the periodogram is to apply a smoothing filter after computing the estimate. We hope that this example makes obvious the shortcoming of this approach: although the variance of the estimate would be reduced, the bias would remain.
The Thomson multitaper estimator is meant for regularly sampled time-series data, and Stoica & Sundin (1999) proved that when the spectrum is relatively smooth, the multitaper estimate is equivalent to the maximum likelihood. In the rest of this paper, we introduce a way to combine the multitaper estimator with the LS periodogram in order to make this approach feasible and useful for irregularly sampled timeseries data.  . Comparison of the power spectrum estimate made by the classical periodogram (blue curve) to that made by the multitaper spectral estimator (black curve). The true power spectrum for the AR(4) process is shown as a red dashed curve.

Irregular Times
There exists a version of the multitaper for irregular time series (Bronez 1988). The approach involves computing generalized DPSS for each frequency band of interest by solving a generalized matrix eigenvalue problem. By definition, the Bronez solution achieves optimal in-band concentration, just like Thomson (1982) in the case of regularly sampled times. However, the estimator is extremely computationally expensive and numerically difficult to compute. These qualities make it ill-suited to problems in which a quick and reliable solution is necessary for Markov Chain Monte Carlo (MCMC) sampling or in which there are thousands of time series to be analyzed (e.g., forthcoming data from the Legacy Survey of Space and Time; LSST Science Collaboration et al. 2009).
Still others have developed solutions for otherwise regular series with gaps (e.g., Fodor & Stark 2000;Smith-Boughner & Constable 2012). Most recently, Chave (2019) introduces a generalization of the multitaper for gappy time-series data that are otherwise regularly sampled. However, this method is not intended for data with arbitrary sampling times, and when applied to time series with many small gaps does not work well (A. Chave 2020, personal communication).
Approaches that interpolate the autocovariance matrix using slotted estimators based on discretization of time lags are also possible (Stoica & Sandgren 2006), but again would seem to be more suited to series that are otherwise regular or comprised of a small set of discrete time increments. A review of various parametric, non-parametric, and heuristic approaches is available in Babu & Stoica (2010).
Given the computational advantage afforded by the spectral estimators under regular sampling, a natural question is whether we can use methods of spectrum estimation that are designed for regularly sampled time series in cases where the times are irregularly sampled.
A seemingly popular approach is to interpolate the irregularly sampled time series to a mesh of regular times, then proceed with spectrum estimation as though the interpolated time series were what was originally observed. However, the consequences of the interpolation on the spectrum estimate and the inferred characteristics of the underlying spectrum should be considered. For example, Lepage (2009) demonstrated that spurious peaks in the spectrum can occur as a result of interpolation of an irregularly sampled time series to regular spacing.
Consider the following highly simplified example. Suppose we have a white noise process that follows for any t, where  represents a normal distribution. We draw a realization from this process x at observation times that are irregular, and in addition at observation times on a regular mesh for comparison. Next, we use linear interpolation to interpolate the times for the irregularly sampled time series. Figure 5 compares the multitaper spectrum estimate obtained by linear interpolation of the irregularly sampled series to the spectrum estimate that would have been obtained had the series been regularly sampled.
The linear interpolation has a noticeable effect on the overall shape of the spectrum. Although the underlying process is white, and should have resulted in a flat spectrum estimate, interpolation has reddened the estimate: high frequencies have a lower associated spectrum estimate, and low frequencies have a higher associated spectrum estimate.
Here we might like to know whether the estimate of the power spectrum in the bottom panel of Figure 5 would lead to a false detection of signals (as was seen in Lepage 2009). Mann & Lees (1996) describe a popular method for estimation of background noise and signal detection in climatic time series. The method consists of partitioning the spectral representation of a series into a background autoregressive spectrum and a set of signals. The signals are then tested for significance by comparison to a null distribution. Figure 6 shows the application of the Mann & Lees (1996) method to the interpolated spectrum estimate. Although it appears that the background spectrum does not do a good job of describing the   Mann & Lees (1996) method to the interpolated white-noise series spectrum seen in the bottom panel of Figure 5. The light gray is the locally smoothed spectrum estimate used to fit the autoregressive spectrum in solid black. The dotted lines are 90th, 95th, and 99th percentiles of the background spectrum. Estimated spectrum values exceeding these percentiles would be interpreted as significant signals. spectrum estimate at higher frequencies, several low to moderate frequency signals are detected by the procedure. In cases where the method is applied in "black box" fashion, false detection could easily become a common problem.
Given that interpolation can have misleading effects on estimates of the spectrum, it seems that computational efficiency alone does not justify interpolation to regular times. Thus, there is motivation to develop spectral estimators that use the unevenly sampled data directly without interpolation.
The LS periodogram is widely used and considered to be a solution to the spectrum estimation problem for irregular sampling. However, like the periodogram, it has statistical issues, namely: 1. it is inconsistent, 2. it has spectral leakage that is often worse than the classical periodogram, and 3. it has spectral leakage that depends on the pattern of observation times and that is frequency specific.
In Scargle (1982), the author is concerned primarily with the detection and identification of a single periodic component from a background noise process that is of no interest. Because of this narrow focus, Scargle argues that inconsistency and spectral leakage are not of particular concern, but that the estimator could be improved by applying a window function to reduce spectral leakage if desired.
Here, we propose using interpolated DPSS as windows for the LS periodogram, and averaging the resulting spectral statistics in the same way as for the regular sampling multitaper. Although not theoretically optimal in the sense of Bronez (1988), the statistic can be made fast to compute using a nonequispaced fast Fourier transform (nfft; Press & Rybicki 1989).
The properties of the multitaper spectral likelihood for the evenly sampled case do not necessarily transfer to the irregular sampling case. In particular, the spectral window under irregular sampling varies as a function of the actual sample times and differs as a function of frequency (Scargle 1982). In the next section, we describe an approximate multitaper spectral statistic based on the nfft, the LS periodogram, and the DPSS tapers.

The Multitaper LS Statistic
In this section, we review the LS periodogram and describe how it can be computed using a fast algorithm. Next, we define a multitaper version of the LS periodogram and demonstrate that it has spectral windows that exhibit reduced spectral leakage.

Fast Computation of the LS Periodogram
The LS periodogram estimator (Lomb 1976;Scargle 1982 N is the length of the time series, ω=2πf is the angular frequency, x i is the ith observation, and t i is the sampling time of the ith observation.
The LS estimator has the following properties.
1. The estimator is invariant to time translation. That is, if the t i are replaced witht t i * where t * is a constant, the estimator is unchanged. 2. The estimator has the same distribution for a Gaussian white-noise process as the regularly sampled periodogram. 3. The estimator is equivalent to a least-squares fit of a sinusoid at a given frequency f to the data x in the time domain.
Leroy (2012) points out that the nfft can be used to compute the LS periodogram using the decomposition of Press & Rybicki (1989). Following Press & Rybicki, the LS periodogram can be represented using the following sums: (11) and (12)  where wt cos 2 LS and sin 2ω τ LS can be computed from Equation (12), and cosω τ LS and wt sin LS can be computed using half-angle formulae from wt cos 2 LS and wt sin 2 LS . The LS periodogram is thus expressed as (11) and (12), we use the nfft library (Keiner et al. 2009), which provides a fast Fourier transform for irregularly sampled data. The nfft library defines the nonequispaced discrete Fourier transform (NDFT) as

To compute the sums in Equations
where the x are input complex Fourier coefficients, and x t j * are the time series as observed at times t j *, standardized so that the times are in the interval [−1/2, 1/2). Writing the NDFT in matrix notation, we have (Note that k here represents the index over which the summation occurs in Equation (18), not the kth taper.) To make the nfft library useful for time-series analysis, we need to consider the x as an input and the Fourier representation as an output. The adjoint (conjugate transpose) of the matrix A multiplied by the time-series vector x is Thus, the adjoint NDFT performs the necessary Fourier transform of the input time series x, except that it seems to return the complex conjugate of the usual transform. In fact, because of the centering of the times to the interval [−1/2, 1/2], there is a further rotation of kπ due to the −1/2 offset. This additional rotation requires a further transformation of z k , This subtlety is not mentioned by Leroy (2012).

MTLS: Multitapering the LS Periodogram
Scargle (1982) suggests that the spectral leakage in the irregular sampling case can be improved through data tapering. The Thomson multitaper limits the spectral leakage to a negligible amount in the case of regularly sampled times, so that multitaper spectral statistics spaced far enough apart in frequency are nearly uncorrelated. Here we generalize this approach to unevenly sampled times.
Our approach for multitapering using the DPSS for irregular time series is as follows.
1. We compute the Slepians at a mesh corresponding to an average sampling interval and N is the number of times in the series. The DPSS on the regular mesh are denoted as before: v N W , , where k is the order of the sequence.

The taper weights at intermediate points corresponding to
the irregular times are obtained by interpolation using a cubic spline. We denote these interpolated weights

The interpolated weights for each DPSS order are renormalized so that
The approximation in step 2 is best when the set of times is close to regular-for example, there are no long gaps-which we assume is true for many astronomical data sets (e.g., from the Kepler or TESS spacecraft). The resulting multitaper statistics are then used as in the regular sampling case described previously.
The MTLS is a compromise between the optimal quadratic estimator of Bronez (1988), which is difficult and costly to compute, and the untapered LS periodogram, which can be made fast to compute but does not achieve minimum spectral leakage out of band. Because the spectral leakage performance of the MTLS depends on the observation times  as well as the frequency band, it is crucial that the spectral pseudowindows (see Section 7.3) be examined across a range of frequencies for a given set of irregular times for which the statistic is to be computed.
We now revisit the earlier example of an irregularly sampled white-noise process, interpolated to regular sampling (bottom of Figures 5 and 6 in Section 6). We saw that interpolation to regular sampling induced an overall shape to the spectrum estimate and could lead to false signal detection. As shown in Figure 7, the MTLS spectrum estimate is able to overcome the irregular sampling without introducing any shape or other spurious features. (1982) defines and examines pseudowindows and shows that in some cases, spectral leakage can be significant even at remote frequencies. This suggests that checking pseudowindows is good practice in all cases of spectrum estimation under irregular sampling, whether using the MTLS, LS, the method of Bronez, or any other.

Scargle
A pseudowindow is the frequency-domain response of the estimator to a pure sinusoidal input at a particular frequency. In order to check pseudowindows, it suffices to create an artificial signal at a given frequency, sample it at the same times as the original irregular series, and examine the resulting spectral statistics. This reveals the pattern of spectral leakage out of the band defined by the Slepians.
When using the MTLS for irregular times, the DPSS are interpolated, and the resulting spectral leakage depends on both the particular set of irregular times and the frequency band. Although the spectral leakage properties of the MTLS should be close to optimal for nearly regular sampling, it is a good idea to check the pseudowindows for a range of bands to avoid errors in the interpretation of the MTLS statistics.

Application to Kepler Data
It is hoped that the MTLS spectral statistic can be used as an approximate direct estimator of the spectrum in cases where the observation times are known but irregular. In particular, we believe this estimator may be usefully applied to study timevarying processes in astronomical data.
Here, we perform a preliminary analysis of Kepler timeseries data from an RG star, Kepler-91. Figure 8 shows the time Figure 7. MTLS estimate for a white noise process sampled irregularly. The spectrum estimate is comparable in shape-that is, flat-to the estimate that would have been obtained under regular sampling ( Figure 5, top). Moreover, the average spectrum value and variation in the estimate are similar. series of the relative flux measurements measured by the Kepler spacecraft (preprocessing done by T. Bedding 2020, private communication), and Figure 9 shows detail by zooming in on the first 3000 observations. The star has a transiting planet, Kepler-91b, with an orbital period of approximately 6.25 days (Batalha et al. 2013;, in addition to other periodic variations in flux indicative of stellar oscillations. This system has been extensively studied previously (e.g., Lillo- Barclay et al. 2015;Placek et al. 2015;Budding et al. 2016).
By using measurements of total flux over time from the Kepler spacecraft and estimating the power spectrum, asteroseismologists can unravel stellar properties (e.g., helium core burning in RG stars; Bedding et al. 2011). This inference is based on the peak mode and the spacing of adjacent line components in the spectrum of the time series of the flux, which relies on accurate identification of signals from noise in the spectrum estimate.
The data in Figure 8 are sampled irregularly due to the orbit of the spacecraft around the Sun and correction to the heliocentric barycenter. There are also some longer time gaps in the series (Figure 8). For this preliminary analysis, the series was analyzed in its entirety, across these longer gaps.
We obtain estimates of the power spectrum of the relative flux using MTLS and three different frequency bandwidth values (corresponding to NW=4, 10, or 20), and also using the LS periodogram for comparison. Figure 10 shows the LS periodogram and the MTLS estimate for increasing numbers of   tapers k. The increasing number of tapers with increasing bandwidth has a noticeable effect on the variance of the estimate. In the extreme case of the LS periodogram (top), signals appear to be less obvious than for even the narrowbandwidth MTLS (second panel).
Comparing the three MTLS spectrum estimates, the one with a time-bandwidth parameter of 10 appears to retain the relative amplitude of any signals present while providing a good separation relative to the noise in the estimate. Figure 11 shows a zoom-in of the top and third panel in Figure 10 for comparison.
Although some work has been done on the choice of bandwidth parameter under some assumptions about the underlying true spectrum (Haley & Anitescu 2017), there is no universally applicable method for determining the optimal choice for this parameter.
The MTLS estimate of the spectrum appears to confirm previous estimates of the orbital period of the transiting planet Kepler-91b (Figure 12(b)). In particular, there is a peak in the spectrum estimate at a period that matches previous estimates (i.e., the planet orbital period of 6.25 days; Batalha et al. 2013;, as well as a series of harmonics, which is to be expected due to the shape of the light curve attenuation as the planet passes in front of the star. There are also a number of apparent signals or line components in the frequency band from 8 to 12 cycles per day, approximately 92-140 μHz (Figure 12(c)). The high-frequency portion of the MTLS estimate (seen in Figures 12(a) and (d)) should not be trusted, for reasons explained next.
The apparently large spectrum values at high frequency are potentially due to leakage from elsewhere in the spectrum. To investigate this possibility, we generated synthetic signals across the frequency range and calculated the response of the MTLS statistic for each.
The frequency response shows a pattern of aliasing beyond about 10 cycles per day (Figure 13). This pattern was verified by computing the LS periodogram of the same synthetic signals using an independent computer code (Ruf 1999). The aliasing effect is a result of the particular pattern of times.

Discussion
Unlike many branches of statistics in which observations are assumed to be realizations of independent or conditionally independent random variables, time-series analysis is preoccupied with the dependence of sequences of random variables. Unfortunately for time-series analysts, the random variables underlying the observations are fleeting, and it is impossible to obtain repeated measures of a random variable corresponding to a moment in time. Instead, methods of time-series analysis make simplifying assumptions about the nature of the underlying random variables and their covariance structure. The assumption of stationarity, for example, allows measurement to begin at any point in time, because if the assumption holds then the joint distribution of random variables depends only on their relative position in time. The stationarity assumption allows the covariance to be summarized as a function of time separation.
The autocovariance can be obtained from the spectrum if desired, meaning that spectral statistics or spectrum estimates can be used instead of estimates of the autocovariance. Inference in the spectral domain is convenient because many natural processes are approximately periodic, and show up as signals in the spectrum of a time series. This makes nonparametric spectrum estimation of the kind described here an important contributor to the exploratory data analysis toolbox. In cases where little is known beforehand, a good estimate of the spectrum can provide some insight into the underlying processes, while suggesting parametric approaches that could be successful in describing the data. Like other exploratory approaches, analysis might stop at the spectrum estimate, or might proceed to additional modeling.
Moreover, spectral statistics can be made to be approximately independent using data tapers, so that the total likelihood of a set of spectral statistics can be approximated using a product likelihood. This could greatly simplify inference for parametric models of stellar processes that imply a theoretical spectrum-given a model with unknown parameters, the theoretical spectrum could be used to calculate the likelihood for each of the spectral statistics. Because of their approximate independence, the total likelihood is the product of each, and maximum-likelihood estimates or Bayesian posterior distributions for model parameters are easily obtained. A more general example of this approach is described in Springford (2017).
For irregularly sampled time series, the LS spectrum estimate is widely used. Scargle (1982) was motivated by the case of a single periodic component embedded in a white noise process, and for this type of process the LS periodogram has its place. However, the spectral leakage properties of the estimator Figure 11. Zoom-in of the top and third panel of Figure 10 to illustrate the reduced variance in the MTLS-estimated power spectrum. can be extremely poor, and if the underlying spectrum is not of the type envisioned then it is possible to be misled by the result. Scargle (1982) suggests tapering the time series to control spectral leakage before computing his periodogram, a suggestion that we take to the same conclusion as Thomson (1982) by multitapering using interpolated DPSS tapers.
The generalized LS periodogram (GLS; Zechmeister & Kürster 2009) and the Bayesian GLS (Mortier et al. 2015) are increasingly popular in astronomy. While these methods offer improvements to the standard LS, these improvements are complementary to those offered by tapering, primarily: (1) measurement of error-derived weights and (2) separate mean values by frequency. Analysts who are interested in employing the GLS or BGLS for these reasons can also benefit by tapering or multitapering.
The MTLS statistic we introduce here matches the Thomson multitaper under regular sampling. Deviations from regular sampling introduce frequency-specific differences that depend on the particular set of observation times. For a given set of times, careful examination of the pseudowindows is recommended. If there appear to be strong signals or a large dynamic range in the spectrum, frequencies with high associated spectral statistics should be examined for their leakage properties because they are liable to contribute to spectral statistics at other frequencies in the spectrum.
Unlike in the case of regular sampling, leakage (or aliasing) to distant frequencies is possible for irregular sampling. In particular, high-frequency parts of the spectrum may not be trustworthy and, if judged to be suspect, should be ignored.
In the example of the RG star data, there appear to be many signals present, and the overall shape of the spectrum at high frequencies seems to be affected by a kind of aliasing. This limits the effective Nyquist frequency to something much smaller than D = t N T 1 2 2 . One possibility is that the aliasing is due to regular gaps in the series. If this is the case, an MTLS spectrum estimate could be computed for each contiguous portion of the series, and averaged as in a non-overlapped Welch estimate (Welch 1967). Although this does not resolve the higher-frequency components of the spectrum estimate, if the aliasing is removed it makes interpretation of the estimate more straightforward. Regardless, the technique of separating time series according to the location of long gaps could be helpful in certain cases where the long gaps make the high-frequency portion of the spectrum estimate unreliable.
There are obvious extensions to the current MTLS as an estimator of the spectrum, based on analogies to the use of the multitaper estimator in the regular sampling case.
The eigencoefficients (Thomson 1982) are the complexvalued result of taking the discrete Fourier transform of the DPSS-windowed time series. The eigencoefficients are the basis of several extensions to the multitaper method in the regular sampling case. The Thomson F-test (Thomson 1982) is a test for the presence of periodic components that is based on a regression estimate of the eigencoefficients. The variance explained by the regression is compared to the residual variance of the estimate to form the F-test. In principle a similar test could be constructed for the MTLS, but this might require an adjustment to the method of computation described in Section 7.2. A targeted approach in which frequencies of interest are tested for periodicity might be the most efficient if the fast computing advantage is lost.
Coherence, the frequency-domain analog of correlation, would also be a useful tool in cases where an estimate of the relationship between two or more time series is of interest. Jackknife error estimates for the spectrum estimate, periodic F-test, and coherence are also possible, and would provide insight into the uncertainties and biases associated with each (Thomson 2007).

Conclusion
We have described a multitapered version of the popular LS estimator of the spectral density of an irregularly sampled time series, which we call the MTLS. The MTLS was shown to have limited spectral leakage, but because the leakage depends on the particular pattern of observation times, examination of the pseudowindows is necessary in any particular application.
As an estimator, the MTLS was applied to a time series of photometry data and displayed the reduced variance and spectral leakage characteristic of the multitaper estimator for regularly spaced data. However, examination of the frequency response (pseudowindows) revealed aliasing and a much smaller effective Nyquist frequency than might have been expected.
The subfield of time-series analysis in astronomy is growing with data sets related to exoplanets, asteroseismology, and transients through Kepler, TESS, and ZTF. The subfield is destined to become even more prevalent as survey projects such as the LSST on the Vera C. Rubin Observatory begin in the 2020s. When scientific questions require the identification of more than a single mode in the power spectrum, the LS estimator is a weak tool because it suffers from spectral leakage, inconsistency, and bias. The MTLS estimator, however, shows promise for applications in astronomy where multiple modes are of interest.
In a future work, we will explore other applications of this method to astronomical data and compare it to other recent advances for estimating the power spectrum. Equally important is the development of an F-test for the MTLS statistic, which determines the probability that a frequency signal in the spectrum estimate originates from a genuine periodic process. The F-test provides a way to test for the presence of real periodic signals in time-series data that may contain evidence of, e.g., exoplanets, asteroseismic modes, and stellar rotation.
T. Bedding for supplying the preprocessed Kepler data and for providing comments that helped improve this manuscript. We also thank thank D. Riegert for helpful discussions, and an anonymous reviewer whose comments and suggestions significantly improved the manuscript.
Software: R Statistical Software Environment (R Development Core Team 2012), including the following packages: