2 Three stage period analysis (TSPA)

A sample of n data points y(t_i) is fitted to a model $g(t,\bar{\beta})$ in our period search. The Q free parameters are $\beta_1$, $\beta_2$, ..., $\beta_Q$. The most probable periodicity is searched for by determining the absolute minimum for the sum of the squared residuals

$${\mathrm{RSS}} = \sum_{i=1}^{n} \left[ y(t_{\mathrm{i}})-g(t_i,\beta_1,...,\beta_Q) \right]^2.$$

Let us assume that the model can be expressed as a linear sum

$$g(t,\beta_1,...,\beta_Q)= \sum_{r=1}^{\mathrm{R}} \beta_r g_r(t,\beta_{R+1},...,\beta_Q)$$

of separate model components $g_r(t,\beta_{R+1},...,\beta_Q)$, which depend nonlinearly on the periods (frequencies) or other parameters (e.g. decay rates, frequency shifts). In this case the RSS can be evaluated as a function of the nonlinear parameters only, because the linear parameters $\beta_1,...,\beta_R$ can be estimated by ordinary regression methods for any fixed set of nonlinear parameters $\beta_{R+1},\!...,\beta_Q$. This gives the Least Squares Spectrum

$${\mathrm{LSS}} (\beta_{R+1},\!...,\beta_Q) \! \! = \! \! \! \sum_... ...a_1, \!...,\hat \beta_R, \beta_{R+1},\!...,\beta_Q) \! \right]^{\! 2} \!\!\!\!,$$

where the linear parameter estimates $\hat \beta_1,...,\hat \beta_R$ are also functions of the nonlinear parameters. The LSS minima reveal the most probable nonlinear parameter combinations. If the model is correctly specified, the global LSS minimum should give just the set of parameter estimates searched for. In this way the problem of model parameter estimation can be converted to the problem of solving the global minimum or minima of LSS. This approach has been applied to astronomical data already by Vanicek (1969, 1971) and Taylor & Hamilton (1972), or more recently by Martinez & Koen (1994), Wilcox & Wilcox (1995) and Bossi & La Franceschina (1995). The LSS analysis is not trivial for semirandomly distributed observations, like the astronomical data that often contain periodic gaps. We formulate a flexible and computationally efficient approach to this problem by dividing the analysis into several stages, where the resolution of the spectra increases at each subsequent stage.

   
2.1 Multistage search

The first stage of the analysis, the Pilot Search, provides crude estimates for the nonlinear parameters. A standard Grid Search on a fine mesh centered at these crude estimates is performed during the second stage. Finally, traditional parameter refinement techniques are applied to obtain the highest possible precision (Refined Search). While the grid search and the refinement of the model parameters represent commonly used astronomical data processing techniques, the pilot search method is not so well known.

As a particular case of the general problem, we consider the LSS of the simple model

$$g(f)\!=\! \hat B_0(f) \! + \! \sum_{r=1}^R \! \left[ \hat B_r (f) \cos(r 2 \pi f t_i) \! + \! \hat C_r(f) \sin(r 2 \pi f t_i) \! \right]\! ,$$

where the frequency f is the only nonlinear free parameter. The linear free parameters for the mean ($\hat B_0$) and the amplitudes ($\hat B_r$, $\hat C_r$) are functions of f. The grid search over the full nonlinear parameter space tests a discrete frequency set $f_l\!=\!f_{\mathrm{min}} \!+\! (l\!-\!1)\Delta f$, where $l\!=\!1,2,...,L$ and $f_L\!=\!f_{\mathrm{max}}$. The chosen f_min and f_max limits may represent the physically possible frequency range or some other interval of interest. The fixed frequency step $\Delta f$ must be chosen carefully. If the data cover an interval of $\Delta T \!= \!t_n\!-\!t_1$, it is reasonable constrain the maximum phase shift $\Delta\phi$ that can occur within $\Delta T$ when any tested f_l changes by $\Delta f$. This gives $\Delta \phi \!=\! \Delta f \Delta T$. Sufficient values are $\Delta \phi \! \leq \! 0.1$, like the overfilling factor of $G\!=\![\Delta \phi]^{-1}\!\geq\!10$in Sect. 3.2. Hence the corresponding number of tested frequencies becomes $L\!=\!\Delta T (f_{\mathrm{max}}\!-\!f_{\mathrm{min}}) [\Delta \phi]^{-1}$. If the tested periods are much shorter than $\Delta T$, then L becomes large. Having fixed $\Delta\phi$, one can reduce L by shortening $\Delta T$ or $f_{\mathrm{max}}\!-\!f_{\mathrm{min}}$. The following general scheme of a simple multistage search reduces $\Delta T$ on the first, and $f_{\mathrm{max}}\!-\!f_{\mathrm{min}}$ on the second stage, respectively.

1.

Divide the data into k separate parts with a length of $\Delta T_k$. Compute the LSS spectrum within each part using a frequency step $\Delta f_{\mathrm{crude}}\!=\! \Delta\phi[\Delta T_k]^{-1}$, i.e. the number of tested frequencies is $L_k\!=\!\Delta T_k (f_{\mathrm{max}}\!-\!f_{\mathrm{min}}) [\Delta \phi]^{-1}$. Use the average of all spectra for obtaining the crude periodicity estimates.

2.

Compute the fine grid search spectra for all data within narrow frequency intervals centered at these crude estimates.

Apply a frequency step of $\Delta f_{\mathrm{fine}} \!=\! \Delta \phi [\Delta T]^{-1}$. A suitable width for the tested frequency interval is about $10 \langle \Delta f_{\mathrm{crude}} \rangle$, where $\langle \Delta f_{\mathrm{crude}} \rangle$ is the mean frequency step that was used during the 1st stage. With this particular width, the number of tested frequencies is $L \!=\! 10 \Delta T \Delta f_{\mathrm{crude}} [\Delta \phi]^{-1}$.

3.

Refine the precision of the frequency estimates with the standard Marquardt iteration (see e.g. Press et al. 1988).

The simplest multistage search techniques were first described by Evans (1961). Since then the computation of long spectra with powerful computers has become more popular than the application of complex multistage algorithms.

   
2.2 Nonparametric dispersion estimation

The average of the local LSS spectra gave the crude frequency estimates during the first stage of the simple multistage search described above. Pelt (1983, 1997) has introduced an alternative to compute these estimates. The idea is that if the data form a continuous periodic curve, the dispersion around this curve can be estimated without solving the curve itself. This dispersion is measured with the nonparametric test statistic

$$\Theta (f)={{\sum_{i=1,j=i+1}^{n-1,n} Z(t_i,t_j,f) \left[ y(t_i)-y(t_j) \right]^2} \over {\sum_{i=1,j=i+1}^{n-1,n} Z(t_i,t_j,f),}}$$

where $Z(t_i,t_j,f)\!\geq\!0$ deviates from zero only if

$$t_i-t_j \approx k f^{-1}, ~k=0,\pm 1,\pm 2,...$$

The $\Theta (f)$ test statistic is a good approximation of the LSS(f) spectrum. A spectral window can be introduced like in the standard Fourier estimation. This is achieved with the more general formulation

$$\Theta(f)={{\sum_{i=1,j=i+1}^{n-1,n} Z(t_i,t_j,f) W(t_i,t_j) \lef... ..._i)-y(t_j) \right]^2} \over {\sum_{i=1,j=i+1}^{n-1,n} Z(t_i,t_j,f) W(t_i,t_j)}}$$

that allows additional smoothing with the data window

The value of $D_{\mathrm{max}}\!=\!\Delta T$gives the maximum resolution for the $\Theta (f)$ spectrum, because $W(t_i, t_j)\!=\!1$ for all pairs of observations. Adjusting this algorithm parameter to $D_{\mathrm{max}} \!<\! \Delta T$enables a faster computation of a smoother spectrum. This reduces the number of tested frequencies during the first stage of the analysis, because a smaller D_max allows a longer step $\Delta f$ in the tested frequencies.

Although $\Theta (f)$ involves double summation for every tested frequency, this apparently complex test statistic can be computed effectively by rearranging these sums and using trigonometric approximations. Such computational techniques for different realizations of the algorithm can be found in (Pelt 1980, 1992). The application of this method speeds up the complete multistage analysis, because the $\Theta (f)$ spectrum for a full grid of tested frequencies can usually be computed as efficiently as the fast Fourier transform.

One particular version of a multistage analysis that can also utilize the additional information of the measurement errors is presented below.