A method for the estimation of the significance of cross-correlations in unevenly sampled red-noise time series
W. Max-Moerbeck, J.L. Richards, T. Hovatta, V. Pavlidou, T.J. Pearson, A.C.S. Readhead
aa r X i v : . [ a s t r o - ph . I M ] A ug Mon. Not. R. Astron. Soc. , 1–24 (0000) Printed 28 August 2014 (MN L A TEX style file v2.2)
A method for the estimation of the significance of cross-correlationsin unevenly sampled red-noise time series
W. Max-Moerbeck , ⋆ , J. L. Richards , T. Hovatta , , V. Pavlidou , , , T. J. Pearson ,A. C. S. Readhead Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA National Radio Astronomy Observatory (NRAO), P.O. Box 0, Socorro, NM 87801, USA Department of Physics, Purdue University, West Lafayette, IN 47907, USA Aalto University Mets¨ahovi Radio Observatory, Mets¨ahovintie 114, 02540 Kylm¨al¨a, Finland Max-Planck-Institut f¨ur Radioastronomie, Auf dem H¨ugel 69, 53121 Bonn, Germany Department of Physics, University of Crete / Foundation for Research and Technology - Hellas, Heraklion 71003, Greece
Accepted 2014 August 19. Received 2014 August 18; in original form 2014 June 20
ABSTRACT
We present a practical implementation of a Monte Carlo method to estimate the significanceof cross-correlations in unevenly sampled time series of data, whose statistical propertiesare modeled with a simple power-law power spectral density. This implementation buildson published methods, we introduce a number of improvements in the normalization of thecross-correlation function estimate and a bootstrap method for estimating the significanceof the cross-correlations. A closely related matter is the estimation of a model for the lightcurves, which is critical for the significance estimates. We present a graphical and quan-titative demonstration that uses simulations to show how common it is to get high cross-correlations for unrelated light curves with steep power spectral densities. This demonstra-tion highlights the dangers of interpreting them as signs of a physical connection. We showthat by using interpolation and the Hanning sampling window function we are able to reducethe effects of red-noise leakage and to recover steep simple power-law power spectral den-sities. We also introduce the use of a Neyman construction for the estimation of the errors inthe power-law index of the power spectral density. This method provides a consistent way toestimate the significance of cross-correlations in unevenly sampled time series of data.
Key words: methods: data analysis — methods: statistical — techniques: miscellaneous
Studies of the variability in astronomical sources can reveal as-pects that are not accessible to imaging, which is limited by theangular resolution of current instruments. For example, variabil-ity can be used to set limits on the sizes of the emitting re-gions through causality arguments (e.g., Abdo et al. 2011), to de-termine the size of the broad line region in active galactic nu-clei (e.g., Peterson & Cota 1988), or to detect extrasolar planets(e.g., Charbonneau et al. 2000) among many other applications. Inthis paper we describe the practical implementation of a cross-correlation technique to determine the location of the gamma-rayemission site in blazars, by studying the relation between the vari-ability in the radio and gamma-ray bands. For this purpose we arecarrying out a blazar monitoring program with the Owens Val-ley Radio Observatory (OVRO) 40 meter telescope (Richards et al.2011) and the Large Area Telescope (LAT) on board of the
Fermi ⋆ E-mail: [email protected]
Gamma-ray Space Telescope ( Fermi , Atwood et al. 2009). Our ap-proach is to search for correlated variability between these two en-ergy bands, which would enable us to determine the location ofthe gamma-ray emission regions relative to the radio emission re-gions. The study of cross-correlations between two energy bandspresents a number of challenges from the data analysis and statis-tical point of view: among these are uneven sampling, non-equalerror bars, and short time duration of the light curves. The tech-niques we develop here should be useful for other applications.Related methods have been presented in the literature, for ex-ample the study of cross-correlations with unevenly sampled lightcurves has an extensive literature about its application to rever-beration mapping (e.g., Peterson 1993). These methods present adetailed treatment of the estimation of cross-correlations and timelags, but not of the estimation of significance of the observed cor-relations, a critical aspect for the interpretation of cross-correlationresults. The literature abounds with claims of statistically signifi-cant correlations that are not backed up by rigorous statistical anal-yses. c (cid:13) W. Max-Moerbeck et al.
This paper presents a detailed discussion of the meth-ods used for our investigation of time-correlation between ra-dio and gamma-ray activity in blazars, which is discussed inMax-Moerbeck et al. (2014). Here we present a description ofthe Monte Carlo method used to estimate the significance ofcross-correlations between unevenly sampled time series usingthe method of Edelson & Krolik (1988). In order to estimatethe distribution of cross-correlations in two uncorrelated datastreams we need a model for the light curves. A commonly usedmodel for time variability in blazars and other AGNs is a sim-ple power-law power spectral density (PSD ∝ /ν β ), as hasbeen measured for a small number of sources at various wave-lengths (e.g., Hufnagel & Bregman 1992; Edelson et al. 1995;Uttley et al. 2003; Ar´evalo et al. 2008; Chatterjee et al. 2008;Abdo et al. 2010). The results presented in Abdo et al. (2010) areof particular interest for the OVRO blazar monitoring program.In their paper, they find a value of β γ = 1 . ± . for brightBL Lacs and β γ = 1 . ± . for bright FSRQs in the gamma-ray band. In the radio band a number of publications have mea-sured β radio . It has been found that β radio = 2 . ± . for 3C279at 14.5 GHz (Chatterjee et al. 2008) using a fit to the PSD foran 11 year light curve. Additional indirect estimates for the PSDpower-law index are obtained by Hufnagel & Bregman (1992) us-ing structure function fits. For five sources, they obtain valuesof α = 0 . ± . to . ± . , where α is the exponent onthe structure function SF( τ ) ∝ τ α . The same method is usedfor 51 sources by Hughes et al. (1992) who found that most val-ues of α lie between 0.6 and 1.8, while a couple are closer to0. However, the often assumed relation between the exponents ofthe PSD and the structure function ( β = α + 1 ) is only validunder special conditions, not necessarily found in real data sets(Paltani 1999; Emmanoulopoulos et al. 2010). The structure func-tion has been widely used in blazar variability studies but its inter-pretation is not straightforward, as has been recently discussed byEmmanoulopoulos et al. (2010). These authors used simulations todemonstrate that many of the features in the structure function areassociated with the length and sampling patterns of the light curvesrather than anything of statistical significance. For these reasons,values obtained from the structure function can only be taken as arough measure of the properties of the time series, and thereforewe do not use them here. Instead we fit the PSDs directly.We start by giving a brief description of the data sets used(Section 2), and then provide detailed descriptions of the methodsin Sections 3 and 4. In Section 3 we describe our approach to thecritical problem of estimating a model for the light curves to usewith the Monte Carlo significance estimate. Here we describe animplementation of the method of Uttley et al. (2002) that containssome important modifications. In Section 4 we provide a descrip-tion of a number of modifications we propose to common methodsused to estimate the significance of cross-correlations. We give ajustification for the use of the local normalization (Welsh 1999)in the Edelson & Krolik (1988) method, demonstrate the strongdependence of the significance estimate on the model light curvesand introduce a bootstrap method to estimate the error in the cross-correlation significance estimates. We close the paper with a sum-mary of our main findings and recommendations for the use of thisand related techniques (Section 5).An important aspect of this work is the use of a statisticallywell-defined data set, where long light curves are used independentof the flaring state of the object. A fatal trap that many authors fallinto is that of ”cherry picking” the data by selecting small intervals of data. This approach can produce spurious levels of significancefor the cross-correlations, and hence cannot be used to draw con-clusions about blazar populations, or the long term behaviour ofindividual sources. The methods we discuss in this paper can be adapted to use withany data set, but since we are describing a particular implemen-tation, our simulated data sets are generated making some choicesrelated to the intended application. These choices are motivated bythe data sets associated with our blazar monitoring program in theradio and gamma-ray bands. These data sets are described in detailin Max-Moerbeck et al. (2014) and here we only summarize theirmain properties.A radio observation for each one of the monitored blazars isattempted twice per week, but because of the effects of weatherand other technical problems we obtain unevenly sampled lightcurves (Richards et al. 2011).
Fermi observes the whole sky onceevery three hours (Atwood et al. 2009), but because of the highlyvarying nature of blazars in gamma-rays, sources may sometimesfall below that detection threshold, resulting in upper limits fora given integration period. We consider gamma-ray light curveswith a time binning of one week, which allows us to detect about100 sources most of the time. We have upper limits for about30% of the data. At this level we find that treating the upperlimits as non-observations does not have important effects in themeasured time lags or significance of cross-correlations for thecases with interesting values of the cross-correlation significance(Max-Moerbeck et al. 2014). This behaviour could be a result ofthe particular properties of the light curves we considered in thatstudy, such as long time scales for the variability compared to thegaps created by ignoring the upper limits or the power spectraldensities, and it might not hold in other situations. We thus obtainunevenly sampled gamma-ray light curves as the ones we discusshere.
We begin with a brief summary of the standard methods usedfor the estimation of the PSD and then move to the uneven sam-pling and short time series cases. This discussion is based on themethod presented in Uttley et al. (2002) which is modified to suitour dataset and the range of PSDs we fit. Additional justification ofthe need for binning and interpolation of the light curves is givenin Section 3.2.2. We also present an example of the applicationof our method to a simulated light curve, and a number of testsusing real data sampling for simulated light curves that demon-strate the accuracy of the fitting procedure under different condi-tions. The real data sampling is based on the data set presented inMax-Moerbeck et al. (2014), which have 4 year 15 GHz radio lightcurves from the OVRO 40 meter telescope blazar monitoring pro-gram, and 3 year gamma-ray light curves from the LAT on boardof
Fermi . A study of the effect of increasing the number of simu-lations in the fitting procedure is performed to guide our choice ofparameters for the data analysis. A summary of the method, withemphasis on the improvements we add to the original formulationis given in Section 5. c (cid:13)000
Fermi . A study of the effect of increasing the number of simu-lations in the fitting procedure is performed to guide our choice ofparameters for the data analysis. A summary of the method, withemphasis on the improvements we add to the original formulationis given in Section 5. c (cid:13)000 , 1–24 ignificance of cross-correlations for uneven sampling We define a time series as a time ordered sequence of triplets ( t i , f i , e i ) , where t i is the observation time, f i is the measuredvalue of the quantity of interest (e.g., flux density, photon flux,etc.), and e i is an estimate of the observational error associatedwith the measurement. We assume that the time series is sorted intime and i = 1 , ..., N . An estimation of the PSD can be obtained through the peri-odogram, which is conventionally defined as the squared modulusof the discrete Fourier transform: P ( ν k ) = " N X i =1 f i cos(2 πν k t i ) + " N X i =1 f i sin(2 πν k t i ) (1)where the periodogram is evaluated at the discrete set of fre-quencies ν k = k/T for k = 1 , ..., N/ for N even, or k =1 , ..., ( N − / for N odd, ν Nyq = N/ T is the Nyquist fre-quency and T = N ( t N − t ) / ( N − , see footnote .Estimating the PSD in this way requires sampling a continu-ous time series at discrete times for a finite amount of time. Thesampling operation is equivalent to multiplication of the time se-ries by a Dirac comb, while sampling for a finite time correspondsto a multiplication by a rectangular observing window. These twomultiplications appear as convolutions in frequency space: theoriginal spectrum is convolved with the Fourier transform of theDirac comb and of the rectangular window. As a final step we onlylook at a discrete set of frequencies which is equivalent to multi-plication by a Dirac comb in frequency space. Ignoring the effect of sampling with a Dirac comb in the fre-quency domain, and omitting normalization factors, we find thatthe periodogram is given by P ( ν ) = | W ( ν ) ∗ III t ( ν ) ∗ F ( ν ) | , (2)where F ( ν ) is the Fourier transform of the time series ( t i , f i ) , III t ( ν ) is the Fourier transform of the Dirac comb with samplinginterval ∆ t , and W ( ν ) is the Fourier transform of the samplingwindow function, which is by default a rectangular window, and ∗ denotes convolution.As a result of the convolution with the Dirac comb, we do nothave access to the original spectrum but a modified version that re-peats periodically. Another distortion comes from the convolutionwith the sampling window function, which modifies the shape ofthe original spectrum, and finally we only look at discrete set offrequencies. All these factors have to be taken into account whenanalysing data and interpreting the results. The periodic repetitionof the spectrum gives rise to aliasing, in which high frequencycomponents are mistaken as low frequency components. Convo-lution with a window function can be a serious problem when thesidelobes of the frequency window function lie on regions of thespectrum where the power is much higher than at the frequency ofinterest – this is the origin of the red-noise leakage problem. Hav-ing the spectrum sampled at a number of discrete frequencies canbe problematic if we are searching for narrow spectral componentswhich can be smeared or missed. In what follows we use ν for the frequency and f i for time series data,e.g., flux density, photon flux, etc. This choice of T is consistent with the definition of the discrete Fouriertransform (Brigham 1988) and allows us to make use of the Fast FourierTransform algorithm to increase the speed of the computations. A graphical representation of these operations can help the reader un-derstand their effect. See Figure 6.1 in Brigham (1988) or elsewhere.
For the case of evenly sampled time series, PSD estimationamounts to using the discrete Fourier transform (DFT) along withperiodogram or frequency averaging to decrease the noise whichis distributed as a χ for a single frequency component. Each ofthese averaging processes can reduce the variance at the price ofreduced spectral resolution. For example, in the case of frequencyor periodogram averaging of M components the resulting distribu-tion is χ M , which reduces the variance by a factor of /M withrespect to the non-averaging case.The application of these methods is straightforward in thecase of long time series, where a good estimate of the PSD can beobtained at the expense of reduced frequency resolution. Nonethe-less, problems of aliasing and red-noise leakage can still compli-cate the analysis of broadband signals like the simple power-lawPSDs we fit to our data ( P ( ν ) ∝ /ν β ), for the reasons out-lined below. For relatively flat spectra ( β from 0 to 2) aliasingcan be a problem as high frequency power above the Nyquist fre-quency contaminates low frequencies. This problem is less seriousfor steep spectra ( β > ), that have relatively small amounts ofpower at high frequencies. But in this case red-noise leakage canflatten the high frequency part of the spectrum: power from lowfrequencies contaminates the low amplitude high frequency partsof the spectrum through sidelobes on the sampling window func-tions. To reduce the effects of these problems a combination offilters and sampling window functions can be used (e.g., Brigham1988; Press et al. 1992; Shumway and Stoffer 2011). When working with time series data, problems often arise be-cause the time series is unevenly sampled and relatively short.Uneven sampling requires the use of a different estimate of theperiodogram: the best known alternatives are the Deeming pe-riodogram (Deeming 1975) and the Lomb-Scargle periodogram(Scargle 1982). The Lomb-Scargle periodogram is well suited tothe detection of periodic signals in white noise, because its statis-tical properties are well understood. For the analysis of broadbandsignals, the Deeming periodogram is often used for reasons that aremainly historical as it does not present any real advantages. Thesetwo methods allow us to obtain an estimate of the periodogram forunevenly sampled time series directly, but do not provide a wayto correct for distortions produced by the sampling window func-tions, which can modify the shape of the periodogram significantlyas explained below.
The method we present here was originally developed and de-scribed in detail in Uttley et al. (2002). We describe the main stepshere to highlight the differences between theirs and our implemen-tation.(i) Obtain the periodogram for the light curve and bin it in fre-quency to reduce scatter. The periodogram is given by a frequencybinned version of the following expression P ( ν k ) = 2 TN " N X i =1 f i cos(2 πν k t i ) + " N X i =1 f i sin(2 πν k t i ) (3)where the frequencies are ν k = k/T for k = 1 , ..., N/ for N c (cid:13) , 1–24 W. Max-Moerbeck et al. even, or k = 1 , ..., ( N − / for N uneven. The minimum fre-quency is ν min = 1 /T , the maximum frequency is the Nyquistfrequency ν Nyq = ( N/ /T ) , and T = N ( t N − t ) / ( N − .The multiplicative factor is a normalization, such that the integralfrom ν i to ν f is equal to the variance contributed to the light curvein this frequency range. The evenly sampled time series ( t i , f i ) is obtained from the original one by interpolation onto a regulargrid. This interpolated time series is first multiplied by an appro-priate sampling window in order to reduce red-noise leakage. Ajustification of these steps is given in Sections 3.2.2 and 3.2.3.(ii) Choose a PSD model to test against the data. In this casewe are fitting power-laws of the form P ( ν ) ∝ /ν β but this canbe generalized to any functional dependence. For the given modelsimulate M time series, where M is a large number that allows usto represent a variety of possible realizations of this PSD model.These and all the simulated light curves used in this work are gen-erated using the method described in Timmer & Koenig (1995),which randomizes both the amplitude and phase of the Fouriertransform coefficients consistent with the statistical properties ofthe periodogram.(iii) For each simulated light curve apply the same sampling,add observational noise, and interpolate into the same even grid.Calculate the periodogram for each one. From these M peri-odograms determine the mean periodogram and its associated er-ror as the scatter at each frequency bin.(iv) Using the mean periodogram and errors obtained in step(iii) construct a χ -like test, defined by χ = ν max X ν = ν min [ P sim ( ν ) − P obs ( ν )] ∆ P sim ( ν ) (4)where P obs ( ν ) is the periodogram of the observed light curve, P sim ( ν ) and ∆ P sim ( ν ) are the mean and scatter in the peri-odogram obtained from the simulated light curves, and χ isthe χ of the observed light curve when compared to the simu-lations for a given PSD model. This χ is then compared to thesimulated distribution of χ for which we can obtain M samples, χ , i , by replacing the P obs ( ν ) term by the periodogram of thesimulated light curves, P sim , i ( ν ) , in Equation 4. The fraction ofthe distribution for which χ , i > χ is the significance levelat which the tested PSD model can be rejected, also known as the p -value. Thus a high value of this percentage represents a good fit,while a low one corresponds to a poor fit.The process described above can be repeated for a number ofmodels with different parameters. The final step consists of select-ing the best model fit as the one with the highest value of p . As withany statistical procedure, a measurement of the uncertainty in theparameters of the model needs to be given. In this point we departfrom the original formulation and provide uncertainties based onMonte Carlo simulations of the model fitting process (as describedin Section 3.2.6).The most significant differences with the original implemen-tation are the use of sampling window functions to reduce red-noise leakage and the Monte Carlo estimation of fitting uncertain-ties. Another difference is that we simulate the effects of aliasingby simulating light curves with high frequency components with asampling period of 1 day, instead of adding a constant noise termto the PSD of the simulated light curves as in the original formu-lation. The high frequency cut at 1 day − is justified in our imple-mentation by the small amount of power seen at higher frequen-cies specially in the radio band. Intraday variable sources showvariability in short time scales (Wagner & Witzel 1995), but even in these cases the amplitude of the variability is only a few percentfor most sources (Quirrenbach et al. 1992). At gamma-rays this isnot necessarily true as fast variability has been observed, but giventhat gamma-ray photon fluxes correspond to mean values of longintegrations of at least a week for most blazars, the effects of fastvariability are less important as they are averaged out. Other ap-plications where fast variability is expected might require a higherfrequency cut, making this approach impractical. Another impor-tant difference, although less important conceptually, is the use ofthe Fast Fourier Transform to perform the computations, whichsubstantially decreases computing time. Further discussions of themost important elements of the method are given below. This step is very important when estimating steep PSDs. It is easyto get mislead by intuition developed from the behavior of win-dow functions for evenly sampled time series, but it turns out thatwindow functions for unevenly sampled data do not behave in thesame way. An example is presented in Figure 1, where we show thefrequency response of a uneven sampling pattern with rectangularand Hanning windows, for the periodogram of power-law PSDswith different values of β from 0 to 5. These window functions are w rect ( t ) = ( , t T , otherwise (5) w Hanning ( t ) = ( cos( π ( t − T/ T ) , t T , otherwise (6)From Figure 1 it can be seen that even though we can calcu-late the periodogram directly for an unevenly sampled time seriesthe results we obtain are very noisy and do not vary much amongdifferent values of β . The main problem is that all the PSDs with β > look very similar, showing almost the same slope whenfitted with a linear function after a log-log transformation. This isproblematic as the fitting procedure relies on the differences be-tween different PSD power-law indices to choose the best model.Doing the same exercise for a time series with the same timelength and number of data points but with even sampling we obtainthe results shown in Figure 2. In this case the results are much lessnoisy and the estimated PSDs look different from each other evenfor very steep PSDs. This allows for better discrimination and isrequired to find an upper limit to the source power-law exponentof the PSD.The problems associated with the window functions becomeevident when trying to apply the fitting method using unevenlysampled data, and show up as an inability to find an upper limit tothe power-law exponent β due to the lack of difference between theestimated PSDs for the simulated data. This problem can be solvedby the use of interpolation and an appropriate window function, asubject that is discussed in Section 3.2.3.Figures 1 and 2 illustrate the limited use we can make of di-rect PSD fitting, even for the case of long time series. In this case,red-noise leakage makes it impossible to recover the right power-law index for steep PSDs.The subject of windowing of unevenly sampled data is brieflydiscussed in Scargle (1982). In particular figure 3 in Scargle (1982)shows a few example window functions for the cases of even anduneven sampling using the classic periodogram. That figure illus-trate the very different sidelobe structure that is obtained for the c (cid:13) , 1–24 ignificance of cross-correlations for uneven sampling Figure 1.
Effect of the use of window functions for uneven sampling cases using the rectangular (blue) and Hanning window (green). Each figure shows theresult of simulating 1000 light curves with a given simple power-law PSD ∝ /ν β , with β given in each figure title. The data points are the mean PSD andthe error is the standard deviation in the simulation, while the units of power (vertical axis) and frequency (horizontal axis) are arbitrary. Also included aredirect fits of the slopes of the mean PSDs for the simulated data for each window (individual panels legend). Notice how the linear fits can hardly discriminatebetween different slopes and how all the estimated PSDs look very similar. Figure 2.
Effect of the use of window functions for even sampling cases using the rectangular (blue) and Hanning window (green). Each figure shows theresult of simulating 1000 light curves with a given simple power-law PSD ∝ /ν β , with β given in each figure title. The data points are the mean PSD andthe error bar is the standard deviation in the simulation, while the units of power (vertical axis) and frequency (horizontal axis) are arbitrary. Also includedare direct fits of the slopes of the mean PSDs for the simulated data for each window (individual panels legend). In this case, the shape of the PSDs is lessnoisy and the estimated PSDs for steep cases look different from each other. Even in this case, direct linear fitting of the PSD produces biased results.c (cid:13) , 1–24 W. Max-Moerbeck et al. uneven sampling case, which is at the root of the problem de-scribed here.To clarify this point, we also include the window functions forour test data along with the results of applying the Hanning win-dow. An examination of Figure 3 helps us understand the resultsdescribed below. In conventional Fourier analysis, window func-tions change the frequency response of the sampling, changing thesidelobe structure and thus helping mitigate the effects of red-noiseleakage and aliasing. This behavior can be seen when using evenlysampled data sets, where the sidelobe structure is regular and de-cays as frequency increases. The case for uneven sampling is verydifferent: the shapes of the window functions explains the strongred-noise leakage seen in the simulations and the increased noise.In the case of even sampling we recover the results of conven-tional Fourier analysis, with all the known properties of windowfunctions.For the reasons described above we use linear interpolationand rebinning to interpolate the unevenly sampled light curves toa regular grid, thus allowing for the PSD fitting.
One fundamental difference between the implementation of themethod of Uttley et al. (2002) and ours is that we use windowfunctions to reduce the effects of red-noise leakage. We found thatthis is necessary when dealing with steep power spectral densities,like those found in blazar studies. In our first attempts to fit thePSDs we found that with a rectangular window we were not ableto set an upper limit to the value of β and were only able to set alower limit. The upper limit on β is necessary to constrain the sig-nificance of cross-correlations, as will be described in Section 4. Inthis section we explain the origin of that problem and the solutionwe implemented.For broadband time series a big problem is the leakage ofpower through far sidelobes of the spectral window response. Thisproblem is evident when dealing with high dynamic range PSDs,such as steep power-laws. For these power-law PSDs, it is seenas a flattening of the high frequency part of the periodogram dueto power leaking from low frequency part which has much higherpower. In practical terms, it means that after some critical value ofthe power-law index all the periodograms have a flat slope whichdoes not depend strongly on the PSD (Figures 1 and 2). Most ofthis high frequency power is actually coming from low frequen-cies through sidelobes of the window function. One way to dealwith this problem is by using window functions with low levelsidelobes; some details about their application to our data set arepresented below. Spectral window functions for our data sets
There is a greatvariety of window functions, which differ mainly in the width oftheir main lobe, the maximum level, and the fall-off rate of thesidelobes. The ideal window function will depend on the appli-cation and some experimentation might be necessary. Propertiesof various window functions can be found elsewhere (e.g., Harris1978)We tried a number of them and compared their performancein recovering steep PSDs. We found that among the ones testedthe most suitable one was the Hanning window, which is able torecover a steep spectrum in a range that allows us to fit our lightcurves. Among the special characteristics of this window are itslow sidelobe level, more than dB below main lobe, and the fast Table 1.
Properties of selected window functionsWindow Sidelobe Level Sidelobe Fall-Off 3-dB BW(dB) (dB/oct) (bins)Rectangular − − . Triangle or Bartlett − −
12 1 . ( x ) or Hanning − −
18 1 . fall-off at − dB/decade. As a downside the Hanning windowhas a broader main lobe at 3 dB ( . · /T ) when compared tothe rectangular window ( . · /T ), where T is the length of thetime series.The effect of different windows is illustrated in Figure 4,which shows the periodogram for a series of steep PSDs. Fromthe figure it is also clear why other window functions fail to dis-tinguish between steep PSDs, and thus are not suitable to use withthis method.The results of Figure 4 can be understood by comparing theproperties of the window functions shown in Table 1 (Harris 1978).The reduction of the red-noise leakage when using the Hanningwindow allows us to discriminate between different steep power-law indices of the PSD, and is due to the low level and fast fall-offits sidelobes.Windowing is good for fitting a featureless PSD, but it canbe a source of problems if the goal is to find narrow spectral com-ponents. This is because of the well-known trade-off between res-olution and sidelobe level: tapering the window function in orderto decrease the sidelobe level must reduce the resolution. This hasto be considered when searching for periodic components, a casewhich is outside of the scope of the current analysis. The windowing technique is able to solve the problem with red-noise leakage, but another method that can be used is filtering inthe time domain followed by a correction to the frequency domainresult. The goal of filtering is to eliminate the low-frequency com-ponents that produce the red-noise leakage before computing theperiodogram. Since this changes the spectrum of the time series, ithas to be compensated in the final periodogram by the applicationof a frequency filter.One of these techniques is called pre-whitening and post-darkening by first differencing. In this case, the original timeseries ( t i , f i ) with even sampling is transformed to ( t i , g i ≡ f i − f i − ) . In the frequency domain this is equivalent to filter-ing with | H ( ν ) | = 2[1 − cos(2 πν )] . Higher order filtering ispossible, for example by the application of first order differencingmultiple times (Shumway and Stoffer 2011).Figure 5 shows the result of applying this procedure to simu-lated data with even sampling and a range of power-law slopes ofthe PSD. It can be seen that this method has problems recoveringflat PSDs with β and very steep PSDs with β > . We alsotested it with the sampling of the OVRO data set and found thatin a large number of cases it was not able to provide good upperlimits for β and was outperformed by windowing with the Han-ning window. We therefore use Hanning windowing for the dataanalysis. c (cid:13) , 1–24 ignificance of cross-correlations for uneven sampling Figure 3.
Spectral window functions in the cases of even and uneven sampling. In the uneven sampling case the rectangular (blue curve) and Hanning (redcurve) windows have a response with a relatively high sidelobe level, that does not decay as the frequency increases. For the even sampling case with thesame time length and number of data points we see that the rectangular (green curve) and Hanning window (cyan curve) behave as expected in the usualcase, with a regular sidelobe structure whose amplitude decreases as the frequency increases.
Figure 4.
Comparison of windowed periodogram for power-law PSDs for evenly sampled data. Each figure shows the result of simulating 1000 light curveswith a given simple power-law PSD ∝ /ν β , with β given in each figure title. The data points are the mean PSD and the error is the standard deviation inthe simulation, while the units of power (vertical axis) and frequency (horizontal axis) are arbitrary. Also included are direct fits of the slopes of the meanPSDs for the simulated data in each case using a rectangular (blue), triangular (red) and Hanning (green) windows.c (cid:13) , 1–24 W. Max-Moerbeck et al.
Figure 5.
Effect of the pre-whitening and post-darkening filters in the PSD of evenly sampled time series. Each figure shows the result of simulating 1000light curves with a given simple power-law PSD ∝ /ν β , with β given in each figure title. The data points are the mean PSD and the error is the standarddeviation in the simulation, while the units of power (vertical axis) and frequency (horizontal axis) are arbitrary. Also included are direct fits of the slopes ofthe mean PSDs for the simulated data in each case using first difference (blue curve) and second difference (green curves). A final issue is the addition of noise to simulated light curves, anecessary step to consider the effect of observational uncertaintiesin our ability to measure the PSD. This is not a serious problem forthe radio light curves, which in most cases have very high signal-to-noise ratio. But it is important for most gamma-ray light curves,which have moderate signal-to-noise ratios.In order to add the observational noise to the light curves wefirst need to normalize the simulated data to match the observa-tions. One way to obtain an approximate normalization is by usingParseval’s theorem, which with the normalization we use impliesthat σ = ν max X ν min P ( ν )∆ ν (7)We can estimate the variance for the observations and thesimulations and use a constant factor to make them equal, thusgetting an approximate normalization of the PSD. One problem isthat the data already contain observational noise added to the sig-nal, so for each data point we have d i = s i + n i , where d is thedata, s the signal and n the noise. We estimate the variance to ob-tain σ = σ + σ , under the assumption that the noise and signalare uncorrelated.The variance of the noise can be obtained from the observa-tional uncertainty by σ ≈ ¯ e i , where e i is the σ observationaluncertainty associated with the i -th measurement. The final nor-malization equation is σ = A ( σ − ¯ e i ) (8)We multiply the originally arbitrarily normalized simulateddata by A − , to get a normalization equivalent to the one in the observations. In practice, we use A to transfer the observationalerror bars to the simulations, to which we add Gaussian observa-tional noise to the time domain signal such that e sim , i = A e i . Inthe original formulation the noise is applied to the periodogram,but we choose to apply it directly to the time series in order to ac-count for the different magnitudes of the observational uncertain-ties. The assumption of Gaussian error bars is only approximatefor the gamma-ray data, which have a Poisson distribution. Sincein this analysis we are only considering highly significant gamma-ray detections, we usually have at least 5 photons in each integra-tion and in most cases many more. In this regime, the differencebetween Poisson and Gaussian distributed errors is negligible. Uttley et al. (2002) defined the region of confidence for the fittedmodel parameters as the region for which p (ˆ θ ) > p conf , where p (ˆ θ ) is the p -value for a given set of parameters ˆ θ . For examplea 68.3% confidence interval has p conf = 0 . , while a 95.5%confidence interval has p conf = 0 . . One problem with this ruleis that it is not possible to get 68.3% confidence intervals for fitsin which p < . . This contrast with the usual approach tomeasure uncertainties from χ fits that defines a 68.3% (or anyother level) confidence interval by the region of parameter spacefor which χ ( θ ) − χ ∆ χ , where ∆ χ depends on thenumber of interesting parameters being fit and the confidence level(Avni 1976; Press et al. 1992; Wall et al. 2003). In this widely usedmethod a confidence interval can be obtained independently of thevalue of χ for the fit, thus effectively decoupling the goodnessof fit estimate from the estimation of confidence intervals.For these reasons we decided to estimate frequentist best fit c (cid:13) , 1–24 ignificance of cross-correlations for uneven sampling confidence intervals by using the method of the Neyman construc-tion. These intervals are constructed to include the true value ofthe parameter with a probability greater than a specified level, asdemonstrated in Beringer et al. (2012) and James (2006). Here weonly describe the mechanics of obtaining confidence intervals andrefer the reader to the references for a formal demonstration. Theprocedure requires that we know the probability of a given ex-perimental result β fit , as a function of the value of the unknownparameter β . The distribution for the fitted value of the power-lawindex for a given value of β is estimated with a Monte Carlo sim-ulation resulting in the distribution of β fit as a function of β . Ateach value of β we can construct a confidence interval at a desiredlevel, the results of these confidence intervals are summarized ina figure with β in the vertical axis and β fit in the horizontal axis(see the lower panel of Figure 7 for an example). For each value of β we can draw a horizontal line going from the lower to the upperlimit of the confidence interval. We join the confidence intervalsfor a range of values of β to get a confidence band. We then fitthe PSD to the data, draw a vertical line at β fit , and determine theconfidence interval as the intersection of the vertical line and con-fidence band. This procedure requires that for each PSD fit we runa large number of fits to simulated data, which increases the com-putational time. This is feasible when fitting a single power-lawindex but it can be prohibitive when fitting a larger number of pa-rameters. In principle we are required to have a confidence intervalfor each possible value of β , but to reduce the computational timewe discretize a reasonable parameter range and use linear interpo-lation to fill the gaps in the confidence band.An example of the application of this method is presented inSection 3.3.1. This section starts with an example of the application of themethod to a simulated light curve of known PSD. We present fourtests intended to validate the procedure by fitting a large numberof simulated data with known PSD, using sampling patterns takenfrom the OVRO program, and with observational errors consis-tent with our data. This section ends with a study of the effects ofchanging the number of simulated light curves ( M as defined inSection 3.2.1), when fitting simulated data in one of them, and anexample light curve from the OVRO program in the other. The goalis to get an indication of the associated uncertainties by changing M , as it can have a large impact on the computational time. A simulated light curve with a power-law PSD with β = 2 . , noobservational noise, and sampled in the same way as the sourceJ1653+3945 is shown in Figure 6; along with the periodogram andbest fit.The results of the fitting procedure are summarized by a plotof p vs β (Figure 7). The best fit corresponds to β = 1 . ± . , where the errors were obtained with a Neyman construction,whose resulting 68.3% confidence band is also shown in the figure.In what follows all the errors are obtained in this way. This errorcan be compared with the original error prescription, which in thiscase produces a value of ± . , more than twice the value estimatedfrom the simulations. Figure 6.
Example of the PSD fit method applied to simulated data. Upperpanel is the simulated light curve with a PSD ∝ /ν and no noise. Lowerpanel is the data periodogram binned in frequency (black line) and themean PSD and scatter for the best fit with β = 1 . ± . (black dots anderror bars). In order to validate the implementation we tested it with simulateddata sets of known PSD. Typical sampling patterns and variousrelative amounts of noise are considered to investigate the behav-ior of the method under different conditions. In each of the testswe use M = 1000 to get the mean PSD and scatter at each trialvalue of β . We use trial values of β from 0.0 to 3.5 in steps of 0.05.Besides validating the method, these tests illustrate how our abil-ity to measure the PSD and the associated error in the power-lawexponent depend on the sampling and noise for the particular lightcurves. OVRO sampling pattern 1 and no noise
In this test we use thesampling pattern for the source J1653+3945, the OVRO data areshown in Figure 8 as reference. Note that for all the tests in thissection, the fitted data are simulations, with only the samplingtaken from the OVRO observations. The results of the fit for sim-ulated data as a distribution of best fit values are shown in Figure8. We find that in all cases we are able to recover the true β with atypical uncertainty of 0.2. OVRO sampling pattern 1 and noise
In this test we use the sam-pling pattern for the source J1653+3945 and error bars consistentwith the noise in this source. The results of the fit for simulateddata as a distribution of best fit values are shown in Figure 9. In c (cid:13) , 1–24 W. Max-Moerbeck et al.
Figure 7.
Example of the fitting method applied to simulated data of knownPSD. Upper panel is p versus β for the different model power-law ex-ponents we tested. The peak at 1.85 indicates the best fit. The error onthe fit is obtained from the confidence band which is shown in the lowerpanel. The intersection of the vertical line with the confidence band giveus β = 1 . ± . . this case, the large measurement errors make recovering the PSDexponent very hard and the fitting procedure fails to yield a mean-ingful constraint. OVRO sampling pattern 2 and noise
In this test we use the sam-pling pattern for the source J0423 − β except for the case of β = 3 . . If necessary this could be handledby the use of a different window function. OVRO sampling pattern 3 and noise
In this test we use the sam-pling pattern for the source J2253+1608 and error bars consistentwith the noise in this source. The OVRO data are shown in Figure11. The results of the fit for simulated data as a distribution of bestfit values are shown in Figure 11. In this last case we are also ableto constrain β with an uncertainty of about 0.2. Here we study the effect of varying the number of simulated lightcurves ( M ) on the repeatability of the results. We then establish Figure 8.
Upper panel shows the OVRO data used to get the time sampling.The flux density uncertainties are not used in this test and we assume a per-fect measurement. Four lower panels are the distribution of best fit valuesfor 1000 simulated light curves in each case. In each case this distributiongives an estimation of the error on the fit and is used to construct the con-fidence band. Top left is for β sim = 0 . and β fit = 0 . +0 . − . , top right isfor β sim = 1 . and β fit = 1 . ± . , lower left is for β sim = 2 . and β fit = 2 . +0 . − . , and lower right is for β sim = 3 . and β fit = 3 . +0 . − . .In the case of β sim = 0 . we report the mode and dispersion about thatvalue. All the other cases use the median and dispersion. a criterion to select M for the data analysis and to get an idea ofpossible errors associated with that choice. We test this by fittingthe same simulated data set used in Section 3.3.1, 100 times using M =100, M =1,000 and M =10,000 simulated light curves at eachtrial power-law exponent of the PSD. The distribution of best fitvalues is used to estimate the repeatability of the fitting process.The second test does the same but in this case it fits the OVRO datafor J0423 − M increases, as we would expect. The scatter is reduced byhalf when going from M =100 to M =10,000. We also note that inthe case of the OVRO data we get a big increase (a factor of 1.9)in accuracy when going from M =100 to M =1,000, but a much c (cid:13) , 1–24 ignificance of cross-correlations for uneven sampling Figure 9.
Upper panel shows the OVRO data used to get the time samplingand the flux density uncertainties. Four lower panels are the distribution ofbest fit values for 1000 simulated light curves in each case. In each casethis distribution gives an estimation of the error on the fit and is used toconstruct the confidence band. Top left is for β sim = 0 . and β fit =0 . +0 . − . , top right is for β sim = 1 . and β fit = 1 . +1 . − . , lower left isfor β sim = 2 . and β fit = 1 . +0 . − . , and lower right is for β sim = 3 . and β fit = 3 . +0 . − . . In the cases of β sim = 0 . and . we report themode and dispersion about that value. All the other cases use the medianand dispersion. smaller one (a factor of 1.2) by going to M =10,000. This showsthat exceeding M =1,000 is not necessary for these data and cansave a significant amount of time when studying large samplesof sources. Similar tests can be performed for other data sets todetermine the required number of simulations. Here we deal with the statistical problem of quantifying the signif-icance of the cross-correlation between two time series, in the caseof uneven sampling and non-uniform measurement errors. The twotime series are assumed to contain no upper or lower limits.
Figure 10.
Upper panel shows the OVRO data used to get the time sam-pling and the flux density uncertainties. Four lower panels are the distri-bution of best fit values for 1000 simulated light curves in each case. Ineach case this distribution gives an estimation of the error on the fit andis used to construct the confidence band. Top left is for β sim = 0 . and β fit = 0 . +0 . − . , top right is for β sim = 1 . and β fit = 1 . ± . ,lower left is for β sim = 2 . and β fit = 2 . ± . , and lower right isfor β sim = 3 . and β fit = 3 . ± . . In the case of β sim = 0 . wereport the mode and dispersion about that value. All the other cases use themedian and dispersion. Table 2.
Repeatability of fitted parameters as a function of number of sim-ulated light curvesTest β β βM =100 M =1,000 M =10,000Simulated, known PSD . ± .
08 1 . ± .
05 1 . ± . OVRO data with noise . ± .
13 2 . ± .
07 2 . ± . c (cid:13) , 1–24 W. Max-Moerbeck et al.
Figure 11.
Upper panel shows the OVRO data used to get the time sam-pling and the flux density uncertainties. Four lower panels are the distri-bution of best fit values for 1000 simulated light curves in each case. Ineach case this distribution gives an estimation of the error on the fit andis used to construct the confidence band. Top left is for β sim = 0 . and β fit = 0 . +0 . − . , top right is for β sim = 1 . and β fit = 1 . +0 . − . ,lower left is for β sim = 2 . and β fit = 2 . ± . , and lower right is for β sim = 3 . and β fit = 3 . +0 . − . . In the case of β sim = 0 . we report themode and dispersion about that value. All the other cases use the medianand the 15.86 and 84.15 percentiles to measure the σ dispersion. Our basic data sets are two time series we call A and B. These timeseries are time ordered sequences of triplets ( t ai , a i , σ ai ) with i =1 , ..., N and ( t bj , b j , σ bj ) with j = 1 , ..., P . In both cases t xi isthe observation time, x i is the measured value of a quantity ofinterest (e.g. flux density, photon flux, etc.) and σ xi an estimate ofthe observational error associated with the measurement.Since the time interval between successive samples is not uni-form and the A and B time series are not sampled simultaneously,we need to resort to some kind of time binning in order to measurethe cross-correlation. The cross-correlation between two unevenlysampled time series can be measured using a number of differ-ent approaches. The usual approach is to generalize a standardmethod and use time binning to deal with the uneven sampling. Here we consider two methods that are commonly encountered inthe literature: the discrete correlation function (Edelson & Krolik1988) and the local cross-correlation function (e.g., Welsh 1999).A number of other alternatives have been used to handle the prob-lem of measuring the correlation between unevenly sampled timeseries. Among them are the interpolated cross-correlation func-tion (ICCF; Gaskell & Peterson 1987), inverse Fourier transformof the cross-spectrum (Scargle 1989) and the z -transformed cross-correlation function (Alexander 1997). We do not explore thesealternative methods in this work. These methods provide a wayof estimating the cross-correlation coefficients, but do not providean estimate of the associated statistical significance, which is dis-cussed in Section 4.2.The two most commonly found alternatives are presented be-low. The discrete correlation function was proposed byEdelson & Krolik (1988) and developed in the context of re-verberation mapping studies. For two time series a i and b j , wefirst calculate the unbinned discrete correlation for each of thepairs formed by taking one data point from each time series UDCF ij = ( a i − ¯ a )( b j − ¯ b ) σ a σ b (9)where ¯ a and ¯ b are the mean values for the time series, and σ a and σ b are the corresponding standard deviations. This particularvalue, UDCF ij , is associated with a time lag of ∆ t ij = t bj − t ai . The discrete cross-correlation is estimated within independenttime bins of width ∆ t , by averaging the unbinned values withineach bin, DCF( τ ) = 1 M X UDCF ij (10)The uncertainty in the binned discrete cross-correlation isgiven by the scatter in the unbinned values for each time bin, andis given by σ DCF ( τ ) = 1 M − (cid:16)X [UDCF ij − DCF( τ )] (cid:17) / (11)In the expressions above the sum is over the M pairs for which τ ∆ t ij < τ + ∆ t , where τ is the time lag, and all the binshave at least two data points in order to get a well-defined error.In practice it is recommended to choose M much larger than 2 toreduce the effect of statistical fluctuations.In this case, the mean and standard deviation use all the datapoints in a given time series, but the DCF for a given time lag onlyincludes overlapping samples. This particular choice for normal-ization produces values of the DCF which are not restricted to theusual [ − , interval of standard correlation statistics. This imme-diately challenges the interpretation of the amplitude of the DCF asa valid measure of the cross-correlation and invalidates the use ofstandard statistical tests developed for other correlation statistics,forcing us to find alternative ways to estimate the significance ofcorrelations. A modification that corrects this normalization prob-lem but not the significance evaluation issue is described below. Motivated by the normalization problems presented by the DCF,some authors have proposed a different prescription (e.g., Welsh c (cid:13) , 1–24 ignificance of cross-correlations for uneven sampling ∆ t . Hence, we have LCCF( τ ) = 1 M P ( a i − ¯ a τ )( b j − ¯ b τ ) σ aτ σ bτ (12)where the sum is over the M pairs of indices ( i, j ) such that τ ∆ t ij < τ + ∆ t . The averages ( a τ and b τ ) and standard deviations( σ aτ and σ bτ ) are also over the M overlapping samples only.The main motivation for using this expression instead ofthe DCF is that we recover cross-correlation coefficients that arebound to the [ − , interval. This latter property is a result of us-ing only the overlapping samples to compute the means and stan-dard deviations, which in effect reduces the problem to a stan-dard cross-correlation bounded to [ − , , as a consequence ofthe Cauchy-Schwarz inequality. Additionally, Welsh (1999) showsthat the LCCF can determine time lags more accurately than theDCF in simulated data sets. These are certainly desirable proper-ties, but as explained in Section 4.2, they do not solve the estima-tion of significance problem. In Section 4.3 we perform a series of tests designed to help uscompare the detection efficiency of the DCF and LCCF. In lookingat these results, it is useful to consider the relation between thosetwo correlation measures.From our previous discussions, we can see that the only dif-ference between the DCF and LCCF is in the values used for themeans and standard deviations. In the case of the DCF, the meanand standard deviation are calculated from the complete time se-ries ( ¯ a, ¯ b for the means and σ a , σ b for the standard deviations),while for the LCCF only the overlapping samples at each time lagare used ( ¯ a τ , ¯ b τ for the means and σ aτ , σ bτ for the standard devi-ations). It can be shown that the two are related at a given time lagby DCF( τ ) = LCCF( τ ) σ aτ σ bτ σ a σ b + ( ¯ a τ − ¯ a )( ¯ b τ − ¯ b ) σ a σ b . (13)This linear relation has coefficients that depend on the samplingpattern and the overlap between the two time series at differenttime lags. For long stationary time series, the means and variancesof the overlapping and complete time series will be identical andthe DCF will equal the LCCF. For short or non-stationary timeseries, the coefficients will make the DCF different from the LCCF.Deviations of the multiplicative coefficient ( σ aτ σ bτ ) / ( σ a σ b ) from 1 change the amplitude of the DCF, while deviations of theadditive coefficient (( ¯ a τ − ¯ a )( ¯ b τ − ¯ b )) / ( σ a σ b ) from 0 changethe zero-point of the DCF. The combination of these variationsexplains why the DCF is not bounded to the [ − , interval as isthe LCCF, and can also explain why they have different detectionefficiencies. The standard method used by the reverberation mapping com-munity (Peterson et al. 1998), uses bootstrapping and randomiza-tion to generate slightly modified versions of the original data set,in order to quantify the uncertainty in the location of the cross-correlation peak. A modified data set is constructed by the appli-cation of two procedures. The first is “random subset selection”, in which a bootstrapped light curve is constructed by randomlyselecting with replacement samples from the original time series.In the second, we perturb the selected flux measurements by “fluxrandomization”, in which normally distributed noise with a vari-ance equal to the measured variance is added to the measuredfluxes. Each of these modified data sets is cross-correlated usingthe method of choice and a value for the cross-correlation peak ofinterest is measured. By repeating this for many randomized datasets, a distribution of measured time lags for the cross-correlationpeaks is obtained. This distribution is used to construct a confi-dence interval for the position of the peak.
There has been some discussion in the literature about the effectsof detrending the light curves in order to improve the accuracy ofthe time lag estimates. Welsh (1999) strongly recommended re-moving at least a linear trend from the light curves. His results arebased on simulations with even sampling and do not directly applyto uneven sampling as shown by Peterson et al. (2004). They findthat detrending does not improve accuracy in unevenly sampleddatasets, and produces large errors in some cases. Based on thatfinding, we have decided not to detrend our light curves.We emphasize that care must be taken when correlating timeseries where long term trends are present, as these are guaranteedto produce large values of the cross-correlation coefficient. Ourstudies are mostly concerned with the correlation between periodsof high activity in different energy bands for light curves that ap-pear to have a detectable “quiescent” level. This is generally truefor gamma-ray light curves, but is not always true for radio lightcurves. Radio light curves showing a single dominant increasingor decreasing linear trend should be analyzed with care, as theycan produce spurious correlations. In our opinion, the best remedyfor those cases is to collect longer light curves.
A complete quantification of the cross-correlation needs an esti-mate of its statistical significance. In our case, we need to con-sider the intrinsic correlation between adjacent samples of a giventime series, which are produced by the presence of flare-like fea-tures; a distinctive characteristic of blazar light curves. This be-havior can be modeled statistically by red-noise stochastic pro-cesses (e.g., Hufnagel & Bregman (1992) in the radio and opti-cal, Lawrence & Papadakis (1993) in the X-rays, and Abdo et al.(2010) in gamma-rays). Red-noise processes are characterized bytheir PSD, show variability at all time scales, and appear as timeseries in which flare-like features are a common phenomenon. Thefrequent appearance of flares means that high correlation coeffi-cients between any two energy bands are to be expected, even inthe absence of any relation between the processes responsible fortheir production. To illustrate this point Figure 12 shows simu-lated light curves with power-law power spectral densities (PSD ∝ /ν β ).In fact, every time we cross-correlate two time series, eachof which has a flare, we will get a peak in the cross-correlation atsome time lag. Then quantifying the chances of such peak beingjust a random occurrence is of critical importance. The problemis further complicated by the uneven sampling and non-uniformerrors, so the only feasible method is to use Monte Carlo simula-tions. Standard methods are not suitable for this analysis, as they c (cid:13) , 1–24 W. Max-Moerbeck et al.
Figure 12.
Illustration of the time domain characteristic of simulated light curves with different power-law PSD. In all panels the horizontal axis is time andthe vertical one is amplitude, both in arbitrary units. Top panels 1, 2 and 3 for PSD ∝ /ν , central panels 4, 5 and 6 for ∝ /ν and lower panels 7, 8 and9 for ∝ /ν . The light curves with steeper PSD show more flare-like features that can induce high values of the cross-correlation coefficient as shown inFigure 13. assume that the individual data points are uncorrelated. The ef-fect of ignoring the correlations will lead to an overestimate of thesignificance of the cross-correlations and to an erroneous physicalinterpretation.In Figure 13, we show the results of cross-correlating the in-dependently simulated light curves from Figure 12, which havedifferent values of the power-law exponent for the PSD. It can beseen that correlating light curves with steep PSD, which show fre-quent flare-like features, can result in high cross-correlation coef-ficients that have nothing to do with a physical relation betweenthe light curve pairs. The results illustrate how common it is to gethigh cross-correlations for unrelated light curves with steep PSDsand the dangers of interpreting them as signs of a physical connec-tion. Standard statistical tests that assume uncorrelated data areequivalent to the case of white noise time series (PSD ∝ /ν ),which is illustrated in the upper panels of Figure 13. Since blazarlight curves are more similar to simulated light curves with steepPSDs, it is easy to see how misleading it is to use statistical teststhat ignore the long term correlations in the individual time series. To estimate the significance of the cross-correlation coefficients,we use a Monte Carlo method to estimate the distribution of ran-dom cross-correlations, that uses simulated time series with statis-tical properties similar to the observations. These and related ideashave been applied by several authors (e.g., Edelson et al. 1995;Uttley et al. 2003; Ar´evalo et al. 2008; Chatterjee et al. 2008). Thedetails of the procedure vary from author to author, so we provide a detailed description of our implementation to enable others toevaluate and reproduce our analysis.The algorithmic description of the method we use to measurethe significance of the time lags is as follows:(i) We calculate the cross-correlation coefficients between theunevenly sampled time series using one of the methods describedin Section 4.1.(ii) Using an appropriate model for the PSDs at each energyband, we simulate time series with the given noise properties andsampled exactly as the data. The resulting flux densities are per-turbed by adding noise according to the observational errors. Wecalculate the cross-correlation coefficients of the simulated lightcurve pairs using the same method as for the real data.(iii) We repeat the previous step for a large number ofradio/gamma-ray simulated light curve pairs and accumulate theresulting cross-correlation coefficients for each time lag.(iv) For each time lag bin, the distribution of the simulatedcross-correlation coefficients is used to estimate the significancelevels of the data cross-correlation coefficients.An additional detail is that the gamma-ray time series are theresult of long integrations, so each simulated data point is gen-erated by averaging the required number of samples to replicatethe time binning. For the radio light curves, the integrations areso short that the closest sample can be chosen. Figure 14 showsthe application of the method for an example using simulated datawith the sampling pattern from our monitoring program. We use β = 2 in both bands for the DCF and the LCCF. In both cases,the cross-correlation coefficient at each time lag is represented bythe black dots and the distribution of random cross-correlations by c (cid:13) , 1–24 ignificance of cross-correlations for uneven sampling Figure 13.
Examples of the cross-correlation of simulated light curves shown in Figure 12 using the DCF (upper figure) and LCCF (lower figure). In allpanels the horizontal axis is time lag in arbitrary units and the vertical one is the amplitude of the cross-correlation. Upper panels, cross-correlation ofindependent β = 0 . light curves. Central panels, cross-correlation of independent β = 1 . light curves. Lower panels, cross-correlation of independent β = 2 . light curves. The pair of numbers on the upper left corner of each panel are the light curve numbers from Figure 12 which are correlated in eachcase. The light curve pairs have been simulated independently and yet show large peaks in the discrete cross-correlation function for the cases of β = 1 . and . . The existence and amplitude of peaks in the cross-correlation appears to increase for steeper power spectral densities, independently of the methodused.c (cid:13) , 1–24 W. Max-Moerbeck et al. the colored dotted lines. A time lag τ > means the gamma-rayemission lags the radio, while τ < represents the opposite. Thered lines contain 68.27% of the random cross-correlations, so werefer to them as the σ lines, the orange lines contains 95.45%( σ ), and the green lines contains 99.73% ( σ ) . The colored con-tours provide a quick way to evaluate the significance of the cross-correlation and are used for this purpose throughout this paper. Inthis case, although the amplitudes are relatively high for both theDCF and LCCF, the significance is not even 2 σ indicating onlymarginal evidence of a correlation. We use both DCF and LCCF in our tests, to determine quantita-tively which is the best for the problem of detecting significantcorrelations between two time series. The comparison is made interms of detection efficiency of correlations, at a given significancelevel, and a maximum time lag error. For the tests, we simulate atime series with a very fine time resolution and make two copies,one for each band, in which the only difference is a known timelag and the different sampling pattern, which is taken from exam-ple light curves from our monitoring program.In all the cases, we bin the cross-correlation with ∆ t = 10 days and model the time series with a PSD ∝ /ν , which is alsoused for the Monte Carlo evaluation of the significance. We use M = 1000 uncorrelated time series to estimate the distribution ofrandom cross-correlations and significance, and use these resultsto estimate the significance of cross-correlation for 1000 correlatedtime series. This enables us to determine the significance of thecorrelations and the error in the recovered time lag.This corresponds to the ideal case of a perfect intrinsic corre-lation, which is only distorted by the time lag and different sam-pling of the two time series. The case is also ideal with respectto the significance evaluation, as we perfectly know the model forthe light curves. It is important to keep these points in mind and torealize that the actual detection efficiencies could be much lowerthan what we find through these tests. As a check of the method and to help the reader understand theresults, we first test our ability to detect correlations in a very sim-ple case. In this case a time series with a uniform sampling periodof 3 days is correlated with a copy of itself without any delay ornoise. An example of the simulated data set along with the resultsfor the DCF and LCCF is shown in Figure 15. The same proce-dure is repeated for all simulated time series with known time lagand correlation properties, and the fraction of detected lags at theknown lag ( ± ∆ t ) with a given significance level is reported as anefficiency in Figure 16.In this case, we recover most of the time lags at the right valueand the behaviors of the DCF and LCCF are very similar. Thevalues of the coefficients of the linear relation for τ = 0 (Equation13), are very close to the case when the DCF and LCCF are equal(Figure 17). In what follows we refer to them as the 1, 2 and 3 σ lines or significancelevels Figure 15.
Example of simulated data with PSD ∝ /ν , uniform andidentical sampling, zero time lag and no noise. The upper panel shows thetwo time series which overlap perfectly in this case. The lower panel hasthe results of the DCF and LCCF for this case. The vertical lines showthe position of the most significant peak with color corresponding to themethod used. Horizontal color lines mark the amplitude of the most signif-icant peak for each method. The most striking difference between the twomethods is the normalization which is not restricted to the [ − , intervalin the case of the DCF.
65 70 75 80 85 90 95 100significance0.00.20.40.60.81.0 e ff i c i e n c y DCFlocal CCF
Figure 16.
Detection efficiency versus significance for both methods forthe case of uniform and identical sampling for both time series, zero lagand no noise . In this case close to 95% of the lags are recovered at the rightvalue and σ significance. Fermi -LAT.
We now study a case with sampling taken from the OVRO 40 me-ter blazar monitoring program and
Fermi -LAT data set. Again weadd no noise to the simulations and have zero lag between thetwo light curves, so the only difference is in the sampling pat-tern. In this case, a source was observed for two years with theOVRO 40 meter telescope at 15 GHz with a nearly twice per week c (cid:13) , 1–24 ignificance of cross-correlations for uneven sampling −1000 −500 0 500 1000 1500τ [days]−1.0−0.50.00.51.0 D C F simsource −1000 −500 0 500 1000 1500τ [days]−1.0−0.50.00.51.0 L CC F simsource Figure 14.
Example of cross-correlation significance results. Upper panel shows simulated radio data with arbitrary flux units (top) and simulated gamma-ray data with arbitrary units (bottom). Both light curves use a typical sampling pattern from our monitoring program. Lower left panel is for the DCF andlower right panel for the LCCF. The black dots represent the cross-correlation for the data, while the color contours show the distribution of random cross-correlations obtained by the Monte Carlo simulation with σ (red), σ (orange) and for σ significance (green). A time lag τ > indicates the gamma-rayemission lags the radio and τ < the opposite. sampling (Richards et al. 2011). The gamma-ray data for the samesource has one observation per week and a one year time duration(Abdo et al. 2010).An example of simulated data with this sampling is shownin Figure 18 (upper panel), along with the results for the cross-correlation (lower panel). In this case the radio sampling (bluedots) covers a longer time span than the gamma-ray one (red dots).Figure 19 shows that in this case we only recover a fraction ofthe time lags at a σ significance. This is because the DCF oftenfinds the most significant peak at a lag different from zero (Fig.20). Moreover some of those spurious lags are of high statisticalsignificance. We still get some significant peaks at lags differentfrom zero for the LCCF, but at a much smaller rate.To understand how we can get small values of the DCF atzero lag while still having large values of the LCCF, we can takea look at the distributions of the coefficients of the linear relation(Equation 13), shown in Figure 21. The multiplicative coefficientshould be one in the ideal case, but instead it has a broad distri-bution (upper panel). The additive coefficient should be zero inthe ideal case, but it also has a broad distribution (lower panel).This can effectively reduce the value of the correlation coefficientor make its distribution broader, either way reducing its discrim- inating power. This effect is seen in Figure 22, which shows thedistribution of cross-correlation coefficients at τ = 0 days. In thefigure, the distribution of random cross-correlations is representedwith a dotted line and the one for correlated data with a solid line.The upper panel is for the DCF and the lower panel for the LCCF.The vertical green line represents the σ significance thresholdamplitude for cross-correlation coefficients. The fraction of cross-correlations for correlated data (solid line) that is to the right ofthe green line is approximately equal to the detection efficiency . It can be seen that this fraction is much larger for the LCCF,as a result of increased scatter in the distribution of the DCF whencompared to the LCCF of correlated data, for the reasons presentedearlier. The equality is only approximate because a peak with larger significancemight have appeared in a lag different than τ = 0 . These cases are notexcluded from the histogram.c (cid:13) , 1–24 W. Max-Moerbeck et al.
Figure 17.
Distribution of coefficients of the linear relation between DCFand LCCF for τ = 0 day, for the case of uniform and identical samplingfor both time series, zero lag and no noise . Upper panel is the multiplicativefactor, which is very close to 1 in most cases. Lower panel is the additiveconstant which is very close to 0. These values make DCF ≈ LCCF whichmakes the results of both methods very similar as can be seen in Figure 16.
Fermi -LAT.
We make the same comparison using a data set with radio lightcurves of 4 years time duration sampled about twice a week, andgamma-ray light curves with a 3 year time duration and weeklysampling. We again consider the case with no noise and zero lagbetween the two light curves, so the only difference is in the sam-pling pattern. An example of a simulated data set with this sam-pling is shown in Figure 23 (upper panel), along with the resultsfor the cross-correlation (lower panel). Comparison of the resultsof this section with the shorter dataset test (Section 4.3.2), can giveus an idea of the variation of the relative power to detect correla-tions in different data sets.As in the case of the “short data set”, we find that the effi-ciency of detection strongly depends on the method used. Figure24 shows that the LCCF recovers the right time lag at high signif-icance for all the cases, while the DCF does so in only about 15%
Figure 18.
Example of simulated data with PSD ∝ /ν , for the case“short data set”. Upper panel shows the two time series, which have somesmall differences produced by the different sampling at each waveband.Lower panel has the results of the DCF and LCCF for this case. The ver-tical lines show the position of the most significant peak with color cor-responding to the method used. Horizontal color lines mark the amplitudeof the most significant peak for each method. In this example the LCCFrecovers the right time lag, but the DCF finds a spurious time lag.
65 70 75 80 85 90 95 100significance0.00.20.40.60.81.0 e ff i c i e n c y DCF local CCF
Figure 19.
Detection efficiency versus significance for both methods forthe case of the “short data set”. In this case the efficiencies differ signifi-cantly between both methods, with the LCCF being the more efficient. of the cases. An examination of Figure 25 shows that the DCFproduces spurious correlation peaks with a wide distribution. Asin the case of the “short data set”, some of those spurious peakshave high statistical significance.A comparison of Figures 19 and 24 shows that the perfor-mance of both methods improves as expected when using longertime series. However, as can be seen from Figure 25, the DCFproduces a large fraction of spurious statistically significant corre-lation peaks, while the LCCF recovers a significant correlation at τ = 0 in all cases.Figure 26 shows the distribution of the coefficients for the lin- c (cid:13) , 1–24 ignificance of cross-correlations for uneven sampling (cid:0) (cid:0) (cid:0)
200 0 200 400 600 (cid:1) N DCFLocal CCF
Figure 20.
Distribution of most significant peaks in the correlation for bothmethods for the case of the “short data set”. Upper panels show the lag andsignificance of the most significant peak for both methods. The lower panelis a histogram for the distribution of lags for the most significant peak. ear relation between the DCF and LCCF (Equation 13). We againsee that they significantly differ from the ideal case of a station-ary time series. This provides an explanation for the differencebetween these two estimators of the correlation. As for the caseof the “short data set”, we also look at the distribution of cross-correlation coefficients for the uncorrelated and correlated datasets at τ = 0 (Figure 27). We again see the broad distribution ofcorrelation coefficients for the DCF of correlated data sets, while amuch narrower distribution for the LCCF, demonstrating the betterdiscriminating power of the LCCF. Additional tests were performed introducing various time lags forthe time series and measuring the efficiency of detection for theDCF and LCCF. They all show the same qualitative informationand are thus not included here. In all cases the LCCF outperformsthe DCF and the efficiency of detection improves when using alonger time duration dataset. These results demonstrate that the
Figure 21.
Distribution of coefficients of the linear relation between DCFand LCCF for τ = 0 day, for the case of the “short data set”. Upper panelis the multiplicative factor, which has a very broad distribution, differentfrom 1 in most cases. Lower panel is the additive constant which also hasa very broad distribution, different from the ideal case of 0. These valuesshow the DCF to be different from the LCCF and have a role in producingspurious highly significant peaks in the correlation. LCCF is the more efficient method for recovering time lags withhigh significance.
In this section, we describe some additional issues that should beconsidered when estimating the significance of cross-correlationsusing the Monte Carlo test we have devised, or similar methods.The error on the significance estimate has been mostly ignored inthe literature, while the dependence of the significance estimate onthe model light curves - when not fully appreciated - can lead tosignificance tests that are not consistent with the basic statisticalproperties of blazar light curves.Another effect not considered here has recently been raisedby Emmanoulopoulos et al. (2013). In their paper, they proposea method to simulate light curves that reproduces not onlythe power spectral density, but also the power density func- c (cid:13) , 1–24 W. Max-Moerbeck et al.
Figure 22.
Distribution of the cross-correlation coefficient for both meth-ods at τ = 0 day, for the case of the “short data set”. Upper panel is forthe DCF and lower panel for the LCCF. Both panels show the distribu-tion of random cross-correlations with a dotted line and for correlated datawith a solid line. Points with cross-correlation coefficient to the right of thevertical green line have a significance of at least σ . Figure 23.
Example of simulated data for the “long data set”. Upper panelshows the two time series, which have some small differences producedby the different sampling at each waveband. Lower panel has the resultsof the DCF and LCCF for this case. The vertical lines show the positionof the most significant peak. In this example the LCCF recovers the righttime lag, but the DCF finds an spurious time lag. tion of the flux measurements. This method is an improvementon Timmer & Koenig (1995), that produces Gaussian distributedfluxes and can provide a better approximation to light curves thathave non-Gaussian probability density functions. For our data, thiscan be the case for gamma-ray light curves but it is not of muchconcern for the radio light curves.
65 70 75 80 85 90 95 100significance0.00.20.40.60.81.0 e ff i c i e n c y DCFlocal CCF
Figure 24.
Detection efficiency versus significance for both methods, forthe case of the “long data set”. In this case the efficiencies differ signifi-cantly between both methods, with the LCCF being the more efficient.
As illustrated in Figure 13, the distribution of random cross-correlation coefficients will depend on the model used for the sim-ulated light curves. In order to better appreciate that dependence,we have estimated the significance of the cross-correlation for anexample using simulated data with the sampling pattern from ourmonitoring program (same as Figure 14). We have used 10,000simulated light curves with PSD ∝ /ν β for β = 0 , and . Fig-ure 28 presents the results in the form introduced in Figure 14. Asin Figure 13, we observe an increase in the amplitude of the ran-dom cross-correlation when steeper power spectral densities areused in the simulations. This manifests as increased scatter in thedistribution of random cross-correlations and a lower significanceestimate for the cross-correlations. The dependence of the resultson the particular model of the light curves illustrates the impor-tance of a proper characterization of the light curves variability, asubject we discussed in section 3. It is expected that the precision of the significance estimates willincrease as the number of simulated light curve pairs increases.In order to get an estimate on the expected error in our signifi-cance estimate, due to the finite number of simulations, we havedivided a full simulation with 100,000 simulated light curve pairsinto independent subsets, and provide independent estimates foreach of them. The idea is to observe the scatter when a small num-ber of simulations is used and compare its variation as more sim-ulations are included. The original simulation is divided in twohalves which are subsequently divided into two. The process isrepeated until the number of simulations in each subset is smallenough that results have a very large scatter, and do not give us re-liable significance estimates. For all sources we find that the resultsof a test with smaller number of simulations is less precise than theone using all the simulations. In all cases, the average gives the re-sult of the complete simulation, an expected result since together c (cid:13) , 1–24 ignificance of cross-correlations for uneven sampling (cid:2)
200 0 200 400 600 800 1000 1200 1400 1600 (cid:3) N DCFLocal CCF
Figure 25.
Distribution of most significant peaks in the correlation for bothmethods, for the case of the “long data set”. Upper panels show the lag andsignificance of the most significant peak for both methods. The lower panelis a histogram for the distribution of lags for the most significant peak. they encode the same information. As expected, the scatter is muchsmaller when a large number of simulations is used. An exampleis presented in Figure 29, which clearly shows the reduction in thescatter as the number of simulated light curve pairs is increased.With less than 1,000 simulations the scatter is of a few percentagepoints, and gets to about 0.4% for more than 10,000 simulations.The process described above could in principle be used to ob-tain an error estimate, but instead we compute a more conventionalbootstrap estimate of the standard error, following the proceduredescribed below (this is applied in Max-Moerbeck et al. 2014). Forthe time lag of interests, we have N values of the random cross-correlations obtained from the N simulated light curves. Fromthese N random cross-correlations 1,000 bootstrap samples areobtained, each one giving a different significance estimate. Thesample standard deviation of these bootstrap replications is usedas the error in the significance estimate. An example of the distri-bution of bootstrapped estimates is shown in Figure 30. We thinkthis error estimate is a required step of any Monte Carlo estimate ofthe significance, and we recommend the adoption of this or equiv- Figure 26.
Distribution of coefficients of the linear relation between DCFand LCCF for τ = 0 day, for the case of the “long data set”. Upper panelis the multiplicative factor, which has a very broad distribution, differentfrom 1 in most cases. Lower panel is the additive constant which also hasa very broad distribution, different from the ideal case of 0. These valuesshow the DCF to be different from the LCCF and have a role in producingspurious highly significant peaks in the correlation. alent procedures - an issue that has surprisingly been up to nowignored by all authors. We presented a description of a Monte Carlo method to estimatethe significance of cross-correlations between two unevenly sam-pled time series. We demonstrated the dependence of the signif-icance estimates on the model of the light curves, and presenteda method based on Uttley et al. (2002), that allow us to determinethe best fit for a simple power-law power spectral density modelfor a light curve. An improved way of dealing with the effects ofred-noise leakage is implemented. This method uses interpolationand windowing with a Hanning window, and provides the ability tofit steep PSDs like those found in our data sets. We demonstratedthat windowing is essential to obtain an upper limit on the value of c (cid:13) , 1–24 W. Max-Moerbeck et al.
Figure 27.
Distribution of the cross-correlation coefficient for both meth-ods at τ = 0 day, for the case of the “long data set”. Both panels show thedistribution of random cross-correlations with dotted line and the one forcorrelated data with solid line. Points with cross-correlation coefficient tothe right of the vertical green line have a significance of at least σ . Upperpanel is for the DCF and lower panel for the LCCF. the PSD power-law index. An upper limit is required for meaning-ful cross-correlation significance estimates, which depend on themodel used for the light curves. The method used for error estima-tion of the best fit was modified for one which decouples the good-ness of fit estimate from the estimation of confidence intervals, andthat can indicate the presence of biases in the fitting procedure.The method was evaluated using simulated data sets and found tobe accurate with a typical error in β of less than ± . , for cases inwhich the signal power is large compared to observational noise.The performance of the method is degraded when fitting time se-ries in which the signal power is comparable to the observationalnoise. In these cases, the procedure fails to provide a reliable con-straint on the shape of the PSD, a situation we can consider whenanalysing our data set by using the Neyman construction to obtainconfidence intervals. We also checked the repeatability of the bestfit value when running the procedure multiple times, and find thatit improves when using a large number of simulated light curves( M ). For an example using the OVRO data set, we find that big im-provements are expected when going from M =100 to M =1,000,but any further increase provides a small improvement, and mightnot be worth the increased computational time.Finally, we described the problem of estimating the cross-correlation for unevenly sampled time series. We have shown thathigh values of the cross-correlation coefficients for red-noise timeseries are ubiquitous, and that any method that aims at quanti-fying the significance of correlation coefficients for light curveshaving flare-like features needs to take this into account. We havedescribed a general Monte Carlo method to estimate the signifi-cance of cross-correlation coefficients between two wavebands. Anumber of tests aimed at measuring the effectiveness of a particu-lar cross-correlation method have been performed to compare theLCCF and the DCF. Given the absence of a physical model for theexpected correlations, the method cannot be used to give a defini-tive value of the detection efficiency, but it can be used to comparedifferent alternatives. The main result is that the LCCF has a muchlarger detection efficiency than the DCF when trying to recover a −1000 −500 0 500 1000 1500τ [days]−1.0−0.50.00.51.0 L CC F simsource −1000 −500 0 500 1000 1500τ [days]−1.0−0.50.00.51.0 L CC F simsource −1000 −500 0 500 1000 1500τ [days]−1.0−0.50.00.51.0 L CC F simsource Figure 28.
Example of cross-correlation significance results for simulateddata using the typical sampling from our monitoring program (same lightcurves as Figure 14). We use β = 0 (upper panel), β = 1 (central panel)and β = 2 (lower panel). The black dots represent the LCCF for the data,while the color contours the distribution of random cross-correlations ob-tained by the Monte Carlo simulation with red for σ , orange for σ andgreen for σ . The increased amplitude of random cross-correlations is ev-ident for steeper PSDs. c (cid:13) , 1–24 ignificance of cross-correlations for uneven sampling Figure 29.
Example of scatter in the significance estimate for independentsubsets of the full simulation using different numbers of light curve pairs.The horizontal axis shows the number of simulations used to get each esti-mate and the vertical the significance. Black dots represent each of the in-dependent subsets of the full simulation. The empty circles and error barsrepresent the mean and standard deviation for subsets of a given number ofsimulations. The horizontal segmented line corresponds to the results us-ing the whole simulation. As expected, the scatter of the estimates obtainedusing smaller number of simulations is larger.
Figure 30.
The distribution of the significance estimates for the bootstrapsamples is represented as a histogram. The solid line represents the valueobtained using the whole simulation. The segmented line is the mean ofthe distribution, with the dotted lines the one standard deviation upper andlower limits. linear correlation. The DCF has the additional problem of produc-ing a large fraction of spurious high significance time correlations,which could be mistaken as real correlations. This problem is lessimportant for the LCCF especially when long time series are used.The origin of the difference, and the lack of discriminatingpower for the DCF, seems to originate in the short duration or non-stationarity of the time series involved. In conclusion, we recom-mend the use of the LCCF as a tool to search for correlations. We also show that the significance of the cross-correlationcoefficients is strongly dependent on the power-law slope of thePSD, which makes characterization of the light curves critical.We investigate the error on the estimated significance by repeat-ing the analysis using different numbers of simulations. Especiallyin cases where high significances are claimed, we suggest using abootstrap estimate of the error on the significance and reporting itsvalue as part of the analysis results. The results of the applicationof this method to a data set combining data from the OVRO mon-itoring program and
Fermi
Large Area Telescope are presented inMax-Moerbeck et al. (2014).
ACKNOWLEDGMENTS
The OVRO program is supported in part by NASA grantsNNX08AW31G and NNX11A043G and NSF grants AST-0808050 and AST-1109911. Support from MPIfR for upgradingthe OVRO 40-m telescope receiver is acknowledged. W.M. thanksJeffrey Scargle, James Chiang, Iossif Papadakis and Glenn Jonesfor discussions. The National Radio Astronomy Observatory isa facility of the National Science Foundation operated under co-operative agreement by Associated Universities Inc. TH was sup-ported in part by the Jenny and Antti Wihuri foundation and by theAcademy of Finland project number 267324. We thank the anony-mous referee for constructive comments that greatly improved thepresentation of some sections of this paper.
REFERENCES
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 722,520Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJL, 733,L26Alexander, T. 1997, Astronomical Time Series, 218, 163Aller, M. F., Aller, H. D., & Hughes, P. A. 2009, arXiv:0912.3176Angelakis, E., Fuhrmann, L., Nestoras, I., et al. 2012, Journal ofPhysics Conference Series, 372, 012007Ar´evalo, P., Uttley, P., Kaspi, S., et al. 2008, MNRAS, 389, 1479Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ,697, 1071Avni, Y. 1976, ApJ, 210, 642Beringer, J., et al.(PDG), PR D86, 010001 (2012)Brigham, E. O. 1988,
The fast Fourier transform and its applica-tions . Englewood Cliffs, N.J.: Prentice-Hall, 1988Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M.2000, ApJL, 529, L45Chatterjee, R., Jorstad, S. G., Marscher, A. P., et al. 2008, ApJ,689, 79Coles, W., Hobbs, G., Champion, D. J., Manchester, R. N., &Verbiest, J. P. W. 2011, MNRAS, 418, 561Deeming, T. J. 1975, Ap&SS, 36, 137Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646Edelson, R., Krolik, J., Madejski, G., et al. 1995, ApJ, 438, 120Emmanoulopoulos, D., McHardy, I. M., & Uttley, P. 2010, MN-RAS, 404, 931Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2013,MNRAS, 433, 907Gaskell, C. M., & Peterson, B. M. 1987, ApJS, 65, 1Harris, F. J., Proc. IEEE, 66, 51Hufnagel, B. R., & Bregman, J. N. 1992, ApJ, 386, 473 c (cid:13) , 1–24 W. Max-Moerbeck et al.
Hughes, P. A., Aller, H. D., & Aller, M. F. 1992, ApJ, 396, 469James, F. 2006, Statistical Methods in Experimental Physics: 2ndEdition. Published by World Scientific Publishing Co., Pte. Ltd.,Singapore, 2006.Lawrence, A., & Papadakis, I. 1993, ApJL, 414, L85Max-Moerbeck, W., Hovatta, T., Richards, J. L., et al. 2014, sub-mittedNieppola, E., Tornikoski, M., Valtaoja, E., et al. 2011, A&A, 535,A69Paltani, S. 1999, BL Lac Phenomenon, ASP Conference Series,159, 293Peterson, B. M., & Cota, S. A. 1988, ApJ, 330, 111Peterson, B. M. 1993, PASP, 105, 247Peterson, B. M., Wanders, I., Horne, K., et al. 1998, PASP, 110,660Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ,613, 682Press, W. H. 1978, Comments on Astrophysics, 7, 103Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P.1992, Cambridge: University Press, —c1992, 2nd ed.Quirrenbach, A., Witzel, A., Kirchbaum, T. P., et al. 1992, A&A,258, 279Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011,ApJS, 194, 29Scargle, J. D. 1982, ApJ, 263, 835Scargle, J. D. 1989, ApJ, 343, 874Shumway, R. H. and Stoffer, D. S. (2011).
Time Series Analysisand Its Applications .Springer, New York.Timmer, J., & Koenig, M. 1995, A&A, 300, 707Uttley, P., McHardy, I. M., & Papadakis, I. E. 2002, MNRAS,332, 231Uttley, P., Edelson, R., McHardy, I. M., Peterson, B. M., &Markowitz, A. 2003, ApJL, 584, L53van der Klis, M. 1989, Fourier techniques in X-ray timing. InTiming Neutron Stars (eds. H. ¨Ogelman and E. P. J. van denHeuvel) pp. 27–69. New York: Kluwer Academic / Plenum Pub-lishersWagner, S. J., & Witzel, A. 1995, ARA&A, 33, 163Wall, J. V., Jenkins, C. R., Ellis, R., et al. 2003, Practical statisticsfor astronomers, by J.V. Wall and C.R. Jenkins. Cambridge ob-serving handbooks for research astronomers, vol. 3. Cambridge,UK: Cambridge University Press, 2003Welsh, W. F. 1999, PASP, 111, 1347White, R. J., & Peterson, B. M. 1994, PASP, 106, 879This paper has been typeset from a TEX/ L A TEX file prepared by theauthor. c (cid:13)000