# Detecting periodicity in experimental data using linear modeling techniques

## I Introduction

Periodic, and nearly periodic, behavior is a common feature of many biological and physical systems and there exist several widely-used techniques to estimate the period of a behavior, for example, spectral estimation , autocorrelation , spectrographs, band pass (comb) filters  and wavelet transforms . All of these standard techniques either employ, or are related to, or are a generalization of, Fourier series.

In this paper we propose an alternative method of detecting periodicity that is not so closely related to Fourier series. This new technique applies ideas from information theory to linear autoregressive models of time series to extract evidence of periods.

The basic principle is the following. Given a time series $`\{y_t\}_{t=1}^N`$ one can propose a linear autoregressive model $`AR(n)`$ by

$`y_t`$ $`=`$ $`a_1y_{t1}+a_2y_{t2}+a_3y_{t3}+\mathrm{}`$ (2)   
$`+a_ny_{tn}+e_tt=n+1,n+2,\mathrm{},N.`$   
where $`e_t`$ are assumed to be independent and identically distributed random variables, which are interpreted as the modeling errors . Under these assumptions the maximum likelihood estimate of the parameters $`a_1,a_2,\mathrm{},a_n`$ can be written in terms of a covariance function, and are therefore related to the autocorrelation function and Fourier spectrum. It is common practice to determine the optimal size $`n`$ of the model by using either the Akaike  or the Schwarz  information criteria; this is done to avoid over-fitting of the time series . It has recently been observed that a further optimization of an $`AR(n)`$ model may be possible by deleting some of the terms to obtain a model

$`y_t`$ $`=`$ $`a_0+a_1y_{t\mathrm{}_1}+a_2y_{t\mathrm{}_2}+a_3y_{t\mathrm{}_3}+\mathrm{}`$ (4)   
$`+a_ky_{t\mathrm{}_k}+e_t,`$   
where,

$`1\mathrm{}_1<\mathrm{}_2<\mathrm{}_3<\mathrm{}<\mathrm{}_kn`$

for $`\mathrm{}_i𝐙^+i=1,2,3,\mathrm{},k`$. The hope is to obtain a model that fits the time series equally well, but has far fewer parameters. Profound theoretical arguments, which are a codification of Occam’s razor, imply that if a reduced autoregressive model (RARM) is suitably optimized, then it is superior to an equivalent autoregressive model $`AR(n)`$. The key principle of this paper is that if one has an optimized RARM, that is the RARM that has been reduced to only the essential terms, then the parameters $`\mathrm{}_1,\mathrm{}_2,\mathrm{}_3,\mathrm{},\mathrm{}_k`$, often called lags, provide information about the periodicity of the time series.

A practical procedure for obtaining an optimal reduced autoregressive model (RARM) has been described by Judd and Mees . This procedure was introduced in the more general context of nonlinear modeling, but in the following section we describe briefly the underlying theory in the context of RARM.

The major part of this paper is aimed at presenting evidence that examining the lags of an optimal RARM provides a more robust and accurate means of detecting periods in time series than traditional spectral techniques. That is, the proposed technique unambiguously identifies periodicities even when spectral methods fail to do so, and furthermore, does not falsely suggest the presence of periods when none are present. The evidence presented is a combination of theoretical argument and numerical procedures, which are illustrated with both artificial and experimental data.

An important numerical procedure that will be used to establish that the proposed technique does not falsely identify periods is  surrogate data analysis. The principle of surrogate data analysis is the following. From experimental data one generates artificial data that are “similar” to the experimental data and satisfy a given hypothesis. One then calculates a test statistic for each surrogate data set, and hence obtains an ensemble of statistic values that estimate the distribution of the test statistic under the assumption that the original data is consistent with the given hypothesis. One then compares the statistic value of the original data with the estimated distribution of the surrogates. If the data has an atypical statistic value then the hypothesis will be rejected, otherwise it should be accepted. In this paper we employ this technique to ensure that RARM procedures do not spuriously identify periodicities in temporally uncorrelated surrogate data.

### A Minimum description length

The criteria we use for determining the optimal RARM is the minimum description length. Occam’s razor recommends that the best description of a phenomenon is the shortest description. This principle can be made rigorous using information theory, and the principle was independently developed by Wallace  and Rissanen .

Operationally the principle is applied as follows. Suppose you have a time series $`\{y_t\}_{t=1}^N`$ given to a certain fixed accuracy and that you wish to communicate the data to a colleague. To send the raw data would require a certain number of bits. Alternatively, one could build a predictive model, of the form (4) for example, and then send the model parameters (to some precision), the initial $`\mathrm{}_k`$ observations, and the differences between the model’s predictions and actual observations. Given this information your colleague can reconstruct the original data. If the model of the time series is good, then the total number of bits required for parameters, initial conditions and prediction errors is less than the number of bits of raw data, because the differences between the predicted and actual observations are smaller than the observations. The total number of bits sent in the second case is called the description length, and the model that achieves the minimum description length is the one recommended by the application of Occam’s razor. The dogma is that this model achieves the best prediction of the data without over-fitting.

In practice it is usually sufficient to estimate the description length of a model, rather than calculate it in detail. An estimate will usually have the form

$`\left(\text{description length}\right)(\text{number of data})`$   
$`\times \mathrm{log}\left(\text{sum of squares of prediction errors}\right)`$   
$`+\left(\text{penalty for number and accuracy of parameters}\right).`$   
Following Judd and Mees  the description length of a RARM can be estimated as follows. Given a time series $`\{y_t\}_{t=1}^N`$ define a set of vectors $`\{V_i\}_{i=1}^n`$ by

$`V_0=(1,1,\mathrm{},1)^T,`$

$`V_1=(y_n,\mathrm{},y_{N1})^T,`$

$`V_2=(y_{n1},\mathrm{},y_{N2})^T,`$

$`\mathrm{}`$

$`V_j=(y_{nj+1},\mathrm{},y_{Nj})^T,`$

$`\mathrm{}`$

$`V_n=(y_1,\mathrm{},y_{Nn})^T,`$

and define

$`y=(y_{n+1},\mathrm{},y_N)^T.`$

Observe that if the model (4) is appropriate for the time series one can write

$$y=\underset{i=1}{\overset{k}{}}a_iV_\mathrm{}_i+e$$
(5)   
$`=V_Ba_B+e,`$

where $`B=(\mathrm{}_1,\mathrm{}_2,\mathrm{},\mathrm{}_k)`$, $`V_B=[V_\mathrm{}_1|V_\mathrm{}_2|\mathrm{}|V_\mathrm{}_k]`$ is a matrix, and $`a_B=(a_1,a_2,\mathrm{},a_k)^T`$. The maximum likelihood estimates of $`a_B`$, that is, the values that minimize $`e^Te`$, are given by

$`a_B=(V_B^TV_B)^1V_B^Ty.`$

Now each parameter $`a_j`$ must be sent to some precision $`\delta _j`$, that is, the maximum likelihood estimate of $`a_j`$ is “rounded-off” by an amount $`\delta _j`$. It can be shown  that the optimal precisions $`\delta =(\delta _1,\delta _2,\mathrm{},\delta _k)`$, that is, the round-off for each $`a_j`$ that gives the minimum description length, satisfy

$`(Q\delta )_j=1/\delta _j`$

where

$`Q={\displaystyle \frac{NV_B^TV_B}{(a_BV_By)^T(a_BV_By)}}.`$

Consequently, it can be shown  that the approximate description length of the RARM (4) is

$$\frac{N}{2}\left(1+\mathrm{ln}\frac{2\pi e^Te}{N}\right)+\left(\frac{1}{2}+\mathrm{ln}\gamma \right)k\underset{j=1}{\overset{k}{}}\mathrm{ln}\delta _j,$$
(6)   
where $`\gamma `$ is a constant depending on the overall scale of the data.

Armed with this estimate of the description length of a RARM one can search over all combinations of lags $`B=(\mathrm{}_1,\mathrm{}_2,\mathrm{},\mathrm{}_k)`$ to obtain the optimal RARM, however, Judd and Mees  describe a fast and efficient method of doing this optimization.

## II Detecting periodicity using optimal RARM

A function $`f`$ is periodic with period $`\tau `$ if $`f(t)=f(t+\tau )`$ for all $`t`$. A time series (assumed stationary) has an (approximate) periodicity of period $`\tau `$ if $`y_ty_{t+\tau }`$ for all $`t`$, or, equivalently, the autocorrelation $`\rho `$ has a local maximum at $`\tau `$. The reduced autoregressive model (4) predicts the current value of a time series $`y_t`$ as a weighted average of the previous values, that is, at the time steps $`\mathrm{}_1`$, $`\mathrm{}_2`$, $`\mathrm{}`$, and $`\mathrm{}_k`$ previous to $`t`$. If a time series has periodic behavior, then the lags $`\mathrm{}_1,\mathrm{}_2,\mathrm{},\mathrm{}_k`$ should be (multiples of) the periods.

We claim that one can detect in time series a periodicity of period $`n_{MAX}`$ by the following procedure, called the RARM procedure. For $`n=0,1,2,3,\mathrm{},n_{MAX}`$ build optimal reduced autoregressive models of the form (4) using the algorithm described by Judd and Mees . For each model in this sequence calculate its description length (6) and take as the overall optimal model that model with the smallest description length. We claim that if the overall optimal RARM is non-trivial, then the lags $`\mathrm{}_1,\mathrm{}_2,\mathrm{},\mathrm{}_k`$ should be (multiples of) the periods $`n_{MAX}`$ in the original time series if the time series is sufficient long.

In order to establish our claim we must demonstrate that

1. if the time series contains a period then the RARM procedure detects this periodic behavior, and
2. if the RARM procedure detects a period then there is periodic behavior in the time series.

In section II A we provide a theoretical argument to establish the forward implication (i). In section II B we discuss an essential procedure for ensuring (ii).

### A Forward implication (i)

The argument to establish the forward implication proceeds as follows. First, we observe that a period in a time series will (regardless of whether it is linear or nonlinear) produce a local maximum in the autocorrelation function $`\rho (\tau )`$. Next it is shown below that, in the optimization of a RARM of given maximum size $`n`$, the criterion for inclusion of a particular term $`a_jy_{t\mathrm{}_j}`$ in (4) is closely related to the magnitude of the autocorrelation at $`\mathrm{}_j`$, $`\rho (\mathrm{}_j)`$. Hence, if $`n`$ is large enough, the optimal RARM will include a term corresponding to this periodicity. Rissanen’s minimum description length criterion guarantees that provided the time series is sufficiently long this will always be the case and so the RARM procedure will always detect periods that are present in a time series, provided the time series is sufficiently long.

The remainder of this section elaborates on the detail of this argument. A period $`\tau `$ in a time series $`\{y_t\}_{t=1}^N`$ of $`N`$ scalar measurements is a strong positive correlation between values separated by $`\tau `$ time steps, i.e. the autocorrelation

$$\rho (\tau )=\frac{(y\overline{y})^T(V_\tau \overline{y})}{_{n+1}^N(y\overline{y})^2}$$
(7)   
has a local maximum at $`\tau `$. Without loss of generality we may assume that $`\overline{y}=0`$, and therefore (7) reduces to

$$\rho (\tau )=\frac{V_\tau ^Ty}{y^Ty}.$$
(8)   
Let the set of lags for the optimal RARM of size $`k`$ be denoted by $`B_k=(\mathrm{}_1^{(k)},\mathrm{}_2^{(k)},\mathrm{},\mathrm{}_k^{(k)})`$. The vector $`B_k`$ uniquely determines the least squares model

$`y={\displaystyle \underset{i=1}{\overset{k}{}}}a_i^{(k)}V_{\mathrm{}_i^{(k)}}+e.`$

Define

$`L(\tau )`$ $`=`$ $`\left|V_\tau ^Ty{\displaystyle \underset{i=1}{\overset{k}{}}}a_i^{(k)}V_\tau ^TV_{\mathrm{}_i^{(k)}}\right|`$ (9)   
$`=`$ $`y^Ty\left|\rho (\tau ){\displaystyle \underset{i=1}{\overset{k}{}}}a_i^{(k)}\rho (\tau \mathrm{}_i^{(k)})\right|`$ (10)   
According to the algorithm of Judd and Mees , given $`B_k`$ and $`a_B^{(k)}`$, the next best term to add to the model has the lag $`\tau `$ that maximizes $`L(\tau )`$. However, identity (10) implies that such a $`\tau `$ is a local maximum of $`\rho (\tau )`$.

Rissanen’s minimum description length ensures that, for sufficiently large $`N`$, “if there is any machinery behind the data, which restricts the future observation in a similar manner as the past and which can be captured by the selected class of parametric functions, then we will find that machinery” . The argument in the preceding paragraphs demonstrates that RARM are a sufficiently broad class of parametric functions to capture “machinery” behind the data, including observed periodicities. Thus, if periodicity is present in the data then RARM techniques will detect it — provided $`N`$ is sufficiently large. This ensures the forward implication (i).

### B Reverse implication (ii): Surrogate data techniques

In order to establish that the RARM techniques does not falsely identify a period when none is present, the numerical procedure of surrogate data analysis can be used. The technique of surrogate data was originally introduced by Theiler and colleagues . They suggest three surrogate generation techniques to address three different hypotheses about a time series, but for our purposes we only use Theiler’s algorithm $`0`$ surrogates.

In the present case we are interested in whether a time series contains periodicities, or said in another way, we wish to test the null hypothesis that the time series contains no periodicities, that is, has no temporal correlation. Theiler’s algorithm $`0`$ generates surrogate time series having no temporal correlation by simply shuffling the original time series, or put another way, the surrogates are i.i.d. noise having the same same rank distribution as the original time series .

Our proposal is to use optimal RARM as the test for periodicity, that is, if the optimal RARM is non-trivial in that $`k>0`$ in (4), then periods are present in the time series. To believe the validity of this test one must require that if the optimal RARM detects a period in a time series, then it must not detect any period in algorithm $`0`$ surrogates . This surrogate test must be applied to each data set for which an optimal RARM has been constructed to ensure that the structure detected in each data set is genuine. That is, we propose that an algorithm $`0`$ surrogate test is a necessary part of the procedure of detecting periodicity using an optimal RARM. If RARM methods identify periodicity in the surrogates then this is clear evidence of false identification of periodicity in the data. However, if the RARM algorithm detects no periodicity in the surrogates then periodicity identified in the original data is genuine. To ensure the reverse implication (ii) holds one need only apply an algorithm $`0`$ surrogate calculation.

## III Calculations

In this section we demonstrate with artificial and experimental data that RARM detects periodic behavior (i) if and (ii) only if it is present in the original time series. To demonstrate that RARM detects periodic behavior if it is present in the data we construct artificial data contaminated with noise and demonstrate the effectiveness of the RARM algorithm. We compare the RARM results to traditional Fourier spectral and autocorrelation techniques. We repeat these calculations for some experimental data comparing the RARM algorithm and traditional techniques. To demonstrate that our RARM algorithm detects periodic behavior only if it is present in the data we apply the method of surrogate data.

In section III A we describe the application of these techniques to detect periodicities in recordings of infant respiratory patterns during natural sleep. Section III B applies these methods to artificial data sets to demonstrate the effectiveness of these techniques compared to traditional methods. Section III C describes the application of these same methods to global climatic data.

### A Infant respiratory data

Using inductance plethysmography we have collected measurements of cross-sectional area of the abdomen of infants during natural sleep. From these measurements we extract a measure that can be related to the breath volume . Figure 1 gives an example of data collected in this way.

We applied our RARM procedure to the data illustrated in figure 1 and obtained a model of the form

$$y_t=a_0+a_1y_{t1}+a_2y_{t6}+e_t$$
(11)   
where $`a_02.945206`$, $`a_10.300739`$ and $`a_20.202056`$. Figure 2 shows the result of analysis of this data set with a fast Fourier transform algorithm (MATLAB’s spectrum command.) and an estimate of the autocorrelation function. Both these techniques yield small peaks at the same value (that is, $`6`$) and are consistent with the results of our RARM algorithm. However, the results are not as unambiguous as the results of the RARM algorithm. That is, the RARM detects a periodicity that is not strong enough to be unambiguously identified by spectral methods.

For many time series of breath size  we have computed autocorrelation and Fourier spectral estimates. We have applied our RARM algorithm to each data set and compared this to the result of applying traditional techniques. For these data the period of periodic behavior detected by the RARM algorithm is consistent with the periods detected by autocorrelation. That is, if RARM detects periodic behavior, then it is of the same period as that detected by the autocorrelation estimate (if the autocorrelation detects periodic behavior). Furthermore, if RARM does not detect periodic behavior, then neither does the autocorrelation estimate. The traditional techniques will often fail to detect periodic behavior when the RARM algorithm does detect it.

We have provided experimental evidence that the RARM technique detects periodic behavior when it does occur. Now we will demonstrate that the RARM technique does not lead to spurious identification of periodic behavior. That is, we will show that if the RARM algorithm detects periodic behavior, then there is periodic behavior in the data. To do this we apply a surrogate data algorithm which will ensure that false indications of periodicities can always be identified.

For the data illustrated in figure 1, none of $`100`$ surrogates generated by shuffling the data exhibited periodic behavior of any period. This calculation was repeated with another $`48`$ data sets . In all $`49`$ cases the RARM failed to detect periodic behavior in the surrogate data in at least $`99`$ (of $`100`$) surrogates of each data set. This indicates that the RARM algorithm does not identify periodicities not present in the data.

### B Artificial data

In this section we use the optimal RARM from section III A as a basis for generating noisy artificial data with a known periodicity. From (11) we use the model

$$y_t=a_0+a_1y_{t1}+a_2y_{t6}+e_t$$
(12)   
(where $`a_02.945206`$, $`a_10.300739`$ and $`a_20.202056`$, as above) to generate an artificial data set $`y`$. To this data we add observational noise $`ϵ_t`$ and apply the above analysis to the series $`z`$, $`z_t=y_t+ϵ_t`$. Figure 3 demonstrates the result of this technique for an artificial data set of the same length as the data and normal observational noise with standard deviation $`1`$ ($`e_t,ϵ_tN(0,1)`$). Figure 4 is the result of the same technique for a longer data set ($`5000`$ data points) and more observational noise ($`e_tN(0,1)`$ and $`ϵ_tN(0,2)`$). In both cases RARM clearly identified periodic behavior with period $`6`$. For the time series in shown in figures 3 and 4 we constructed $`100`$ algorithm $`0`$ surrogates. none of them exhibited periodicity detected by RARM.

The traditional Fourier spectral and autocorrelation techniques identify the same period as the RARM technique for the shorter, but less noisy data illustrated in figure 3. However, for the data shown in figure 4 the RARM technique has identified periodicities that are not obvious from traditional techniques. Furthermore, it should be noted that in all cases that the results of the autocorrelation and spectral methods are not clear cut. For reasonably long, but extremely noisy data sets the RARM algorithm still provides a decisive and accurate estimate of the period of periodic behavior present in data.

### C Global climatic data

In this section we describe the application of these techniques with noisy physical data. The time series we use here is monthly deviations from monthly mean global air temperatures over the period 1856–1997 . These global air temperature measurements are obtained by averaging observations at many spatially separated sites on the globe. Figure 5 shows the complete data set. A more detailed discussion of this data may be found in . Analysis using the methods described in this paper demonstrates the presence of periodic fluctuation over periods of $`7`$ months, $`2`$ years and $`45`$ months. Fourier spectral and autocorrelation estimates were also applied (after de-trending this time series) and the results are illustrated in figure 5. From $`100`$ algorithm $`0`$ surrogates RARM did not detect periodicity in $`99`$ of them. These results demonstrate the presence of genuine periodic fluctuation in this time series and that the fluctuation is difficult to detect with traditional techniques. An advantage of the RARM technique is that no de-trending is required. The results of the RARM algorithm are not effected by trends or non-stationarity.

## IV Conclusions

We have provided theoretical and experimental evidence to support the use of RARM techniques to detect periodic behavior in noisy experimental time series. The concept of minimum description length ensures that a RARM built with an MDL modeling criterion will detect any periodicities present in the data. We provided numerical evidence using experimental and artificial data to support this. Moreover these calculations have demonstrated that the RARM algorithm provides an accurate and decisive method of detecting periodicities that is more sensitive than Fourier spectrum or autocorrelation methods.

By applying surrogate data techniques we have demonstrated that the RARM algorithm did not identify periodicities in temporally uncorrelated surrogates. This is strong experimental evidence that the RARM algorithm is robust against identification of false periodicities. It does not identify behavior not present in the original system. However this result has only been supported by numerical evidence and does not imply that true identification with arbitrary data. To guard against false positives we recommend application of surrogate data tests, as discussed in this paper. Periodicity detected using RARM are genuine provided RARM detects no periodicity in i.i.d. surrogates.

## Acknowledgments

We wish to thank Madeleine Lowe and Stephen Stick of Princess Margaret Hospital for Children for supplying the infant respiratory data, and for physiological guidance. We also thank Tiempo Climatic Cyberlibrary for making the global climatic data used in this article easily available.