In many cases, the detection of periodic signals in astronomical time series may be a difficult task due to the low number of available observations with an uneven sampling, to the existence of large gaps, or to the high level of noise. In such cases the amplitude spectrum analysis can be supplemented with Monte Carlo simulations in order to estimate the statistical significance of the highest peak [see e.g. Kuschnig et al. (1997); Pop et al. (2010) and references therein]. Assuming that the values of the parameters of the detected signal are those estimated through the nonlinear least-squares method [see e.g. Pop (1996, 1999)], an estimation methodology of the detectability of the detected periodic signal was also proposed (Pop & Vamoş, 2013). It relies on two different Monte Carlo simulations, each of them treated both through frequentist and Bayesian approaches.

A conceptually different analysis method which can be applied to such astronomical time series is the self-correlation (SC) method proposed by Percy et al. (1981) [see also Percy et al. (1993, 2003)]. Although it is an autocorrelation method, it must not be confused with the classical autocorrelation function. As Percy and his collaborators mentioned, the SC analysis ‘measures the cycle-to-cycle behaviour of the star averaged over all the data’ [e.g. Percy et al. (2006a); Percy & Palaniappan (2006)]. According to Percy et al. (2006a) the method “does not have the equivalent of the ‘aliasing’ problem found in the Fourier analysis” [see also Percy et al. (1993, 2002, 2003); Percy (2007) and Marinova & Percy (1999)]. Consequently, it was found to be useful in conjunction with the visual inspection of the light curves and the Fourier analysis in case of data displaying some degree of irregularity and having regular gaps of the same order of magnitude as the length of the involved (quasi)period in order to distinguish between true and alias periods [e.g. Percy et al. (1993, 2002, 2003, 2006a, 2006b, 2010); Percy & Mohammed (2004)].

The SC is defined as the average absolute value of the variations over a given time lag $h$ of a time series. If the time series is periodic with the period $P$, then its variation for lags equal to the multiples of $P$ vanishes and the SC reaches a minimum. If the periodic component is hidden in noise, then the shape of SC is distorted, but the positions of its minima can supply useful information concerning the respective periodicity. For example, Percy et al. (1993) noted that the level of the SC minima is the observational scatter $\sigma$ of the data, while that of the maxima is related to the amplitude of the periodic signal in the observational data. In addition, the SC method has the desirable property to be suitable for the analysis of time series arising from quasiperiodic and irregular phenomena.

Unlike the amplitude/power spectrum analysis, the SC method allows to perform the signal detection not in frequency domain, but in the time domain. The presence of a periodic component is signaled not by a single peak, but by a series of minima which are approximately equidistant. Depending on the length of the analysed time series and on the involved periodicity, it is possible to use more than one minimum in order to perform the periodicity detection.

The present paper has two main goals. First, we want to remove an important deficiency of the SC method already outlined by Percy et al. (2006a): the lack of information concerning the statistical significance of a periodicity or of a time scale detected using the SC function. The second goal is to prove the power of this method in the case of pathological data sets on variable stars featured by low number of unevenly sampled observations in the presence of a high noise level. We propose an improvement of the SC method by thickening the grid on which the SC is computed. Using it in Monte Carlo simulations we significantly improve the false alarm probability which allow us to claim the detection of the deterministic component in a noisy time series.

As an application we consider the orbital period variability of the eclipsing binary system ER Vulpeculae investigated in a recent paper by Pop & Vamoş (2013) for which contradictory results were obtained by different researchers. After rejecting some outlying points and averaging those corresponding to the same events, they obtained a data set consisting in 91 times of minimum light which cover a relatively short time base of 54.9 yr (Fig. 1). They found that the run of the $O-C$ residuals can be described by a statistically significant secular parabolic trend corresponding to a relatively slow increasing trend in the orbital period featured by a relative period change rate of $2.73\times {10}^{-8}$ yr^-1, and a 17.78 yr periodic component, with a low amplitude of about 3.5 min.

The confidence level for the rejection of the null hypothesis that there is no periodic component in the ER Vul time series was a little over 96 per cent, smaller than the threshold usually taken into account (99.9 per cent). Using an improved SC algorithm, the enhanced self-correlation (ESC), we have obtained a confidence level of 99.98 per cent for the existence of a deterministic component of the variability phenomenon of ER Vul, much higher than that previously estimated for a periodic component by Pop & Vamoş (2013). The evaluation of the detectability of the above inferred periodicity supplied somewhat worse results with respect to those estimated by Pop & Vamoş (2013). However, the relatively good coincidence between the results of the analyses performed through two complementary methods is interpreted by us as a confirmation of the detection of the periodic modulation of the orbital period of the eclipsing binary ER Vul. It represents both an important step in understanding this system and an insight into the involved physical phenomena. Obviously, this result needs future confirmation relying on higher quality observational data.

In this section we analyse the properties of the SC function for some particular cases in which we gradually introduce the characteristics of an observed astronomical time series, like that of ER Vul. First of all, we analytically study the SC function for a continuous sinusoid. After that, we numerically study its discretization both on equidistant and nonequidistant grids. Finally, we analyse the influence of a superposition of a Gaussian white noise on this deterministic signal.

In Fig. 4 we plot $S_k$ of the $O-C$ residuals of ER Vul computed for three different values of the bins number $K$. The empirical criteria for the computations of $S_k$ listed in Section 2.2 are satisfied for $K\le 31$. But then the SC diagram is too smoothed even for $K=31$ (upper panel in Fig. 4) when the average number of $\mathrm{\Delta }{x}_{nm}$ per bin is 132. For $K=70$ this value reduces to 58, almost half of the bins contain less than 10 $\mathrm{\Delta }{x}_{nm}$ values and the SC diagram becomes irregular (lower panel in Fig. 4). Pop & Vamoş (2013) used the value $K=55$ for which the average number of $\mathrm{\Delta }{x}_{nm}$ per bin is 74 and only three bins at the end of the SC diagram contain less than 10 values.

The first minimum of $S_k$ located at $h_{\min }\approx 18$ yr is significantly higher than zero indicating a high noise level in the time series of ER Vul. We specify that thereafter, before the computation of the SC, every time series is always standardised, i.e. the average value is subtracted from it and then it is divided by its standard deviation.

From the SC graphs displayed in Fig. 4 one can see that although in the first half of the SC exists one maximum and one minimum, their shape varies quite a lot as a function of $K$. The minimum at $h_{\min }\approx 18$ yr is better defined for low values of $K$, while when $K$ increases, a secondary local maximum appears in the region of this minimum. In the same way, for low values of $K$ the maximum situated at $h_{\max }\approx 12$ yr is well defined. On the contrary, when $K$ increases, another secondary maximum appears, sometimes even higher than the primary one. In fact, the maximum must occur at approximately $h_{\min }/2\approx 9$ yr, but as we can see in Fig. 2b (continuous line with circle markers), even for the simple sine function the first maximum is strongly distorted and shifted to the right due to the uneven sampling. These fluctuations could explain the less defined second minimum which should be located at $2h_{\min }\approx 36$ yr.

In order to remove the shape fluctuations of $S_k$, we apply the stacking method used in seismology, which enhances the deterministic signal when the seismograms are very noisy (Shearer, 2009). Instead of a single value of $K$, we compute $S_k$ for several values, $K\in [K_{\min },K_{\max }]$, and by superposing the corresponding SC values we obtain the ‘enhanced self-correlation’ (ESC). For each $K$ we obtain another grid $h_k$, ($k=1,2,\mathrm{\dots },K+1$) on which we compute the corresponding SC values $S_k$. We construct a thickened grid by putting together the bins centers of the grids for all $K$. The new grid contains

values and we denote by $h_j^{\ast }$, $j=1,2,\mathrm{\dots },J$, its nodes arranged in ascending order and by $S_j^{\ast }$ the corresponding values of the SC. Figure 5a shows the ESC of the $O-C$ curve of ER Vul for $K_{\min }=30$ and $K_{\max }=75$ ($J=2415$). We chose the $K_{\min }$ value according to the criteria given in the previous section. The value $K_{\max }$ is chosen so that the ESC contains all the possible shapes of the extrema of individual SCs.

Using the ESC we can determine the interval of $h_k$ values which contains the most significant information concerning the deterministic component of the ER Vul time series. In the classical approach of the SC method there are no rigorous criteria for the choice of the optimum number of extrema to be used in the analysis (Percy et al., 2004). In the case of ER Vul, the shape of $S_j^{\ast }$ displayed in Fig. 5a confirms the existence of the two extrema identified in Fig. 4. They appear together with weakly defined additional extrema and we want to establish if they are the signature of the existence of a deterministic component. We compare their shape with those of $S_j^{\ast }$ in the case of Gaussian white noise having the ER Vul sampling (Fig. 5b). We see that in the first half the variations of $S_j^{\ast }$ for the Gaussian white noise have lower amplitude compared to that of ER Vul time series and display a random shape. In the second half of the two ESCs, both curves have fluctuations with higher amplitude than those appearing in the first half. Therefore, we have to find a suitable criteria in order to choose the region of ESC over which we perform the comparison.

This separation of the properties of the ESC may be analysed through its variability on intervals $h_j^{\ast }\in [0,H]$, where $H$ is a given parameter. We denote by ${\sigma }_S(H)$ the standard deviation of the ESC over such an interval. We have computed $⟨{\sigma }_S(H)⟩$, the average of ${\sigma }_S(H)$ over a statistical ensemble of 1000 Gaussian white noise, having the time sampling of the ER Vul data. The results are displayed in Fig. 6 together with ${\sigma }_S(H)$ for ER Vul time series. We can make the following remarks:

• i.

  For the variability of the ESC of ER Vul and that of Gaussian white noises are comparable. This means that the initial region of the ESC of ER Vul does not contain the signature of a deterministic component.

• ii.

  The difference between the variability of the two ESCs becomes significant for $H\in (10,30)$, i.e. in the region where the two main extrema of the ER Vul ESC are situated.

• iii.

  For $H>30$, the variability of the ESC for Gaussian white noises increases much faster. In this region, the random fluctuations dominate both the ESC functions, and the difference between their variabilities diminishes.

Basically, the problem of the detection of a periodic signal in the presence of an additive noise is tackled through statistical tests [see e.g. Scargle (1982); Schwarzenberg-Czerny (1993, 1998, 2003); Kuschnig et al. (1997); Pop et al. (2010) (see also references therein), Appourchaux (2011)]. Usually, one claims the detection of a periodic signal if the power/amplitude of a suspected peak exceeds a level computed for a given threshold value $\alpha$ of the probability of false alarm $P_{FA}$, i.e. the probability of the type I error that the observed time series is purely random. Alternatively, one may estimate the false alarm probability $P_{FA}$ corresponding to the height of the considered spectral peak and compare it with a previously established threshold value $\alpha$, e.g. 0.1, 0.05, 0.03, 0.01, and 0.001 (Scargle, 1982; Black & Scargle, 1982; Horne & Baliunas, 1986; Kuschnig et al., 1997). Then, we may conclude that the observed time series consists not only of pure noise, but it contains a deterministic component. Moreover, due to the nature of the Fourier transform, from the amplitude/power spectrum analysis one may infer that the respective spectral peak is the signature of a periodic harmonic signal.

In the case of SC analysis we can design a statistical test for the existence of a deterministic component in the ER Vul time series, without any reference to its periodicity. The standard deviation ${\sigma }_S(H)$ of the ESC over the interval $[0,H]$ introduced in the previous section is an indicator for the existence of a deterministic component, not only for a periodic one. The difference of variability between the ESC of ER Vul and that of the Gaussian white noises occurs in the first half of ESC. Therefore, we choose as test statistic the standard deviation ${\sigma }_S(H_0)$ for $H_0=\mathrm{\Delta }{t}_{\max }/2=27.46$ yr.

In Fig. 7 we display the sample pdf of ${\sigma }_S(H_0)$ obtained from a statistical ensemble of 10000 Gaussian white noises. We also marked the value of ${\sigma }_S(H_0)$ corresponding to the $O-C$ curve for ER Vul (dashed vertical line). In very few cases (0.5 per cent from all) the variability of the ESC for Gaussian white noise is greater than that observed for ER Vul. Then $P_{FA}=0.005$ and we can reject the null hypothesis that the ER Vul time series is a Gaussian white noise at a confidence level of 99.5 per cent. This result does not significantly depend on the value of $H$. For $14.65\le H\le 31.12$, the above confidence level takes values between 99.18 per cent and 99.74 per cent.

The false alarm probability value $P_{FA}=0.005$ is one order of magnitude higher than that resulted from the Monte Carlo simulations relying on the amplitude spectrum performed by Pop & Vamoş (2013). When instead of the ESC method we have used the simple SC computed for a single $K$ value, we have obtained a $P_{FA}$ value a few times larger. Although we have improved the result of the statistical test, however we have not reached the significance level of 0.001 which allows us to claim the detection of a deterministic signal in the $O-C$ residuals of ER Vul. The statistical test described above may be used indifferently of the shape of the deterministic component possibly existing in the analysed time series. If there is no indication that the deterministic component is periodic or at least quasi-periodic, then it has to be determined by a detrending algorithm. The resulting trend strongly depends on the noise level in the analysed time series (Vamoş & Crăciun, 2012).

We can further improve the detection of the deterministic component found in the $O-C$ residuals of ER Vul by introducing additional information concerning the periodicity of the detected signal. Pop & Vamoş (2013) proposed the following model for the cyclic variability of the $O-C$ curve of ER Vul

where $j=1,2,\mathrm{\dots },N_d$, $N_d$ is the number of the analysed data, $n_j$ is the number of orbital cycles, and $\{z_j\}$ is a Gaussian white noise with zero mean and unit variance. The values of the model parameters were estimated via nonlinear least-squares method: ${\tau }_{11}=0.00246\pm 0.00048$ d, ${\nu }_{01}=0.0001075\pm 0.0000032$, and ${\varphi }_{11}=6.17\pm 0.19$ rad. The final residuals obtained after removing the estimated periodic component have a standard deviation ${\sigma }_z=0.00257$ d.

This model was used by Pop & Vamoş (2013) to evaluate the detectability of the estimated periodic signal through Monte Carlo simulations. In order to approach this problem relying on the self-correlation method we have to take into account the effect emphasized by Welsh (1997). The presence of a Gaussian additive noise in the analysed time series determines an overestimation of the amplitude of a periodic harmonic signal existing in the data, either in the case of least-squares method, or in that of the amplitude spectrum computation. In the previous study of ER Vul the amplitude ${\tau }_{11}$ was computed by imposing the condition that the value of the highest peak in the spectrum of the observed data and in the averaged spectrum of a statistical ensemble of synthetic time series of the type (5) are equal (Pop & Vamoş, 2013). They obtained a lower value for the amplitude of the periodic component ${\tau }_{11}=0.00142$ d.

This time we impose the condition that the shapes of the ESC of the observational data and of the artificial time series are similar. We generate 100 artificial time series with a fixed value of the ratio $p={\sigma }_z/{\tau }_{11}$. The numerical generation is performed in two steps. First, we use Eq. (5) with ${\tau }_{11}=1$ and ${\sigma }_z=p$ and then we standardise it. For each artificial time series we compute the ESC $S_{num}^{\ast }$ for $h\in (10,H_0)$, the interval where we have established that the ESC of the observed data differ most from a Gaussian white noise (see Section 3 and Fig. 5). After that we compute the average over the statistical ensemble $⟨S_{num}^{\ast }⟩$. We perform the same computation for several values of $p$ separated by a step $\delta p=0.1$.

In Figure 8 we display the averaged ESC for four values of the ratio $p$. One notices that the difference between the extrema of ESCs decreases as the noise amplitude increases. The square norm of the difference between the ESC of the observed data $S_{obs}^{\ast }$ and $⟨S_{num}^{\ast }⟩$ is a measure of the similarity of the two graphs. The minimum of this norm is obtained for $p=0.8$. Taking into account that the ER Vul $O-C$ residuals are standardised using the standard deviation ${\sigma }_{obs}$, it follows that the coefficients appearing in Eq. (5) are ${\tau }_{11}={\sigma }_{obs}/⟨S_{num}^{\ast }⟩=0.00287$ d and ${\sigma }_z=p{\sigma }_{obs}/⟨S_{num}^{\ast }⟩=0.00229$ d.

In Section 4 we have established only the existence of a deterministic component inside the ER Vul time series. Now we may use the numerical model presented in Section 5 for the estimation of the existence of a periodic signal, not only a deterministic one. We use the information related to the positions of the first two ESC extrema of the artificial time series with a periodic component to design a new statistical test. The position of these two extrema can be established from Fig. 8 in the regions where the differences between the ESC graphs for different values of ${\sigma }_z$ are the largest. These regions have unequal widths due to the uneven sampling of the ER Vul data (in Fig. 2b the width of the first maximum is half of that of the first minimum).

Using this criterion we choose the interval $I_{\max }=[11.44,13.35]$ yr which contains the values belonging to the first maximum of $⟨S^{\ast }⟩$ (the first two dashed vertical lines in Fig. 8). For the first minimum of $⟨S^{\ast }⟩$ we chose the interval $I_{\min }=[15.40,20.01]$ yr (the last two dashed vertical lines in Fig. 8). We define the average peak-to-peak amplitude of ESC as the difference between the average of the ESC values belonging to the corresponding two intervals

where $N_{\max }$ and $N_{\min }$ are the number of values within the intervals $I_{\max }$ and $I_{\min }$, respectively.

In the following we use the average amplitude $A$ of the ESC as a test statistic in order to evaluate the detectability of a periodic signal hidden in an additive noise. The statistical ensemble is composed of 10000 Gaussian white noises identical with those used in Section 4. The probability of false alarm is $P_{FA}=0.0002$, i.e. one order of magnitude smaller than that obtained in Section 4 (Fig. 9). Hence the $O-C$ residuals of ER Vul analysed using the ESC method contain, with a very high confidence level (99.98 per cent), a deterministic component having a characteristic time scale of about 18 yr.

We have also generated statistical ensembles composed by 10000 artificial data sets relying on the model described by Eq. (5). The sample pdf of the average amplitude $A$ of the ESC for $p={\sigma }_z/{\tau }_{11}=0.8$ is plotted in Fig. 9. The value of the ‘observed’ mean amplitude of the ESC computed for the $O-C$ residuals of ER Vul allows the estimation of detectability of the inferred periodic signal [Scargle (1982); Schwarzenberg-Czerny (1998, 2003); Levy (2008); Sturrock & Scargle (2009); Appourchaux (2004, 2011); Johnson (2013); Pop & Vamoş (2013), see also references therein]. Assuming well-defined values of the parameters of the periodic model and for those describing the noise distribution, for a given ‘observed’ value of the considered statistic one can estimate the corresponding values of the probability of false alarm $P_{FA}$ and of the probability of detection (or detection efficiency) $P_D$. As it is well-known, these two quantities allow us to evaluate the detectability of the assumed periodic signal. The probability of detection represented as a function of the probability of false alarm is known as the ROC (Receiver Operating Characteristic) curve.

The ROC curves for the noise levels corresponding to those used in Fig. 8 are plotted in Fig. 10. The detectability of a periodic signal in ER Vul data depends strongly on the noise level $p$ and the chosen significance level. For the optimal model presented in Section 5 with $p=0.8$, one remarks the following representative points on the run of the ROC curve ($P_{FA}$,$P_D$): (0.001, 0.681), (0.01, 0.884), and (0.05, 0.960). For the ‘observed’ value of the average amplitude estimated from the ESC diagram one has $P_{FA}=0.0002$ and $P_D=0.509$. Hence the high noise level in the data led to a low value of the detection probability, i.e. there is a high chance to obtain a non-detection result in the presence of the inferred signal. For lower noise levels the detection probability increases significantly, becoming almost ideal for $p=0.4$ (Fig. 10).

In the Bayesian approach we have the positive predicted value (the Bayesian detection rate)

and the negative predictive value

where $\alpha =P_{FA}$ and $\beta =1-P_D$ are the type I and the type II errors, respectively [see e.g. Levy (2008)]. According to Bayesian approach, the high $PPV$ indicates that the alternative hypothesis is very likely to be true, i.e. we have a very high chance to detect the inferred signal when it is really present in the analysed data. The relatively low $NPV$ indicates that, because of the high noise level in the data, there is a significantly high chance to obtain a detection result even when the inferred signal lacks from the data.

We have investigated the properties of the self-correlation function when it is used to search for periodicities in astronomical time series. We have considered the case of a harmonic periodic component existing in different types of time series with increasing complexity with respect to the sampling and to the presence of additive Gaussian white noise. We have begun with the ideal case of a continuous sinusoid defined on an infinite domain. In this case the minima of the SC function vanish and are situated at integer multiples of the period. The maxima occur at odd multiples of the half-period and are equal with $A=4a/\pi$, $a$ being the sinusoid amplitude. The second analysed case is that of a continuous sinusoid defined on a finite domain. The minima of the SC function are the same as in the previous case but the maxima are distorted, the distortion increasing with the time lag. The relation between the peak-to-peak amplitude $A$ of the SC function and that of the sinusoid $a$ remains approximately the same as above.

The following complication which we have introduced is the discretization of the finite continuous sinusoid. An equidistant grid does not distort the shape of the SC diagram, however, the level of minima of the SC does not vanish if the sampling interval of the time series is not a divisor of the period. For a non-equidistant grid (as the sampling of the ER Vul data) the positions of minima remain approximatively the same, but their levels become nonzero. The non-equidistant grid also induces large fluctuations in the region of maxima, so that their shapes can be determined only numerically. The proportionality between $A$ and $a$ is preserved, but its explicit form has to be also numerically studied.

Finally, we considered the case of realistic data by including in the last type of artificial time series the Gaussian white noise. We analytically derived the mathematical expression which allows us to estimate the standard deviation of the noise $\sigma =\sqrt{\pi }/2{⟨S_k⟩}_{\min }$. The amplitude $A$ of the SC decreases when the noise level $\sigma$ increases following a power law.

We have also investigated the statistical properties of the SC method especially concerning its ability to detect harmonic periodic signals hidden in noise and in unevenly sampled data. We propose an improvement of the SC method by applying ‘the stacking method’ used in seismology, which we call enhanced self-correlation (ESC) method. It achieves an enrichment of the statistics which allows us to obtain a better resolution of the SC diagram. Using ESC we have obtained through Monte Carlo simulations a value of the false alarm probability low enough to reject the null hypothesis that the $O-C$ residuals of ER Vul are only noise. The very low value of $P_{FA}$ indicates that almost surely the observed time series plotted in Fig. 1 is not a Gaussian white noise. But the high level of the noise prevents the precise identification of the characteristics of its deterministic component. It could be a simple harmonic periodic signal or a multiple periodic one dominated by a single period or a quasiperiodic phenomenon.

For the evaluation of the detectability of the inferred deterministic component, we take into account the previous result of Pop & Vamoş (2013) which describes it by a periodic harmonic model, the involved periodicity being of about 18 yr. For this approach we considered a test statistic with a higher specificity with respect to the assumed model, i.e. the average peak-to-peak amplitude in the ESC estimated from the first two extrema. Taking into account the average peak-to-peak amplitude of the ESC estimated for the standardised ER Vul time series, and using the binary decision making formalism within both the frequentist and Bayesian approaches, we performed Monte Carlo simulations which led us to the quantities featuring the detectability of the assumed periodic modulation.

Our conclusion is that the presence of a periodic signal was correctly detected at very low $P_{FA}$ values () and at very high values of $PPV$ ($PPV>0.999$), but due to the high noise level in the observational data, there is also a high risk to get false detections (relatively low values of $P_D$ and $NPV$) when the periodic component lacks and the time series contains only noise. These results obtained through ESC method are in good agreement with those obtained by Pop & Vamoş (2013). The present methodology is complementary to that relying on the amplitude spectrum analysis, and supplies a further substantiation of the orbital period modulation of about 18 yr occurring in ER Vul. However, because of the unsatisfactory detectability of this periodicity, the new observational data on this eclipsing binary system have to be collected with great care in order to reduce as much as possible the level of the observational noise.

We thank the referee, Professor J.R. Percy, for useful comments.