\section{Introduction}\label{sec:intro}

The identification of periodicity in time series data is a fundamental task across scientific fields. 
Applications include the classical study of periodic mass extinctions~\citep{raup1986biological}, 
gene expression cycles~\citep{cho1998genome}, circadian rythms in biology~\citep{costa2013inference}, 
and seasonal oscillations in epidemiology~\citep{murray2011seasonal}.
Typically, observations times in these applications are unequally spaced, which 
introduces identification issues that standard statistical methods, be it frequentist or Bayesian, cannot handle.
When the hidden true periodicity cannot be point-identified, it is an open question how statistical inference should proceed, 
and whether it is possible to anticipate such identification issues through better observation designs.

Our work in this paper is motivated by the analysis of radial velocity data in exoplanet detection.
The idea behind this method is to look for periodic variations in velocity data from a star.
When these variations cannot be explained by known signals (e.g., rotation or magnetic activity of the star), they are attributed to an (unseen) exoplanet orbiting around the star. In recent years, the deployment of 
an extensive network of high-resolution, Earth-based observatories has led to astonishing 
exoplanet detections, including ``51 Pegasi b"~\citep{mayor1995jupiter} ---awarded the 2019 Nobel Prize in Physics--- and even a potential discovery in our immediate stellar neighborhood of  Proxima Centauri~\citep{anglada2016terrestrial}. 

The statistical methodology behind these exoplanet detections proceeds in two main steps. 
The first step is to test for any sign of periodicity. If this test succeeds, the second step is to estimate the hidden periodicity.

To test periodicity, statistical methods rely on the so-called {\em periodogram} of the data, which essentially quantifies the power of the Fourier spectrum of the data across different periods~\citep{schuster1898investigation}.
For instance, the classical method of~\citet{fisher1929tests} relies on a statistic 
defined as the ratio of the periodogram peak over its average value across all periods, which has a 
known distribution under a harmonic model with i.i.d.~normal errors.
This method has since been improved in terms of power~\citep{siegel1980testing, bolviken1983new, chiu1989detecting}, and  has been extended to more general hypotheses~\citep{juditsky2015detecting} 
and sparse alternatives~\citep{cai2016global}.
For adaptations of these tests to astronomical applications see also
\citep{baluev2008assessing, baluev2013detecting, delisle2020efficient, linnell1985test}, 
and the concept of ``false alarm probability".
In general, testing for periodicity presents no big challenges as current tests 
are robust and their properties are well-understood.

However, the problem of actually estimating the hidden periodicity is 
more subtle. For one, standard parametric results assume away identification issues, typically by assuming that the true periodicity lies in a constrained space allowing point-identification with an increasing number of observations~\citep{walker1971estimation, hannan1971non, hannan1973estimation, quinn1991fast, quinn2001estimation}; this is also related to the problem of aliasing and 
the Nyquist limit in statistical Fourier analysis~\citep[Section 2.5]{bloomfield2004fourier}.
Furthermore, standard methods that rely on Fourier analysis also require equally spaced observation times. In practice, however, observations are usually unequally spaced for practical reasons, 
such as the difficulty in gathering fossils to study biological extinctions or the logistics of scheduling 
astronomical observations with large telescopes.

To address the  problem of unequally spaced observation times,~\citet{lomb1976least} and \citet{scargle1982studies} independently developed a modification of Schuster's periodogram. 
Even though the Lomb--Scragle periodogram was motivated by Fourier analysis, 
it has an intuitive statistical interpretation as  it is  equivalent to  fitting the data with a harmonic model through least squares. As such, the Lomb--Scragle periodogram currently occupies a central place in exoplanet detection methodology~\citep[Section 1.1]{vanderplas2018understanding}, where the focus is on assessing whether the periodogram peak is statistically significant, 
typically via extreme value theory. 

A common flaw in current methodology, however, is that a significance result on the periodogram peak is interpreted to mean that the true unknown period is also ``somewhere near" the peak.
The implicit underlying assumption is, of course, that the finite-sample distribution of the periodogram peak is 
near the true period, and also normal in shape. This is unrealistic, however, 
because the likelihood in the harmonic models employed in practice is highly non-smooth and multimodal with respect to the underlying period. As such, without strong, unrealistic assumptions on the observation times, the periodogram peak is not ``well-behaved" neither in finite samples nor asymptotically, in the sense that it may deviate significantly from normality, and may even be inconsistent.
See also~\citep{hall2006using} for similar discussions on the challenges of estimating periodicity through periodograms.
Technically, one fundamental problem is that the separability conditions required for consistency of M-estimators~\citep[Theorem 5.7]{van2000asymptotic}, such as the peak of a periodogram, generally do not hold under arbitrary observation designs.  
Below, we present some examples to illustrate these issues and motivate our work.


\begin{example}[Synthetic data]
\label{example1}\em
Let  $t_i=i + 0.05 U_i$, $i=1, \ldots, 100$, be the observation times for data
$y_i = 1.5\cos(2\pi t_i/\sqrt{2}) + \varepsilon_i$, where $U_i\sim \mathrm{Unif}[-1, 1]$ and $\varepsilon_i\sim N(0, 1)$ i.i.d. The true periodicity is $\theta^\ast=\sqrt{2}\approx 1.414$.
The data and its generalized Lomb--Scragle periodogram~\citep{zechmeister2009generalised} are shown in Figure~\ref{fig0}~(left). In our sample, the peak is at $\hat\theta_n = 3.417$ and is  distinct 
from the peak around the true period. There are also multiple peaks 
associated with the true period harmonics. Using Fisher's test we can reject 
the null of no periodicity~($p$-value is 0).

\begin{figure}[t!]
\centering
\includegraphics[scale=0.36]{fig/fig0}
\caption{
\small
{\em Left:} Periodogram from one dataset from model of Example~\ref{example1}. The 
dataset was seeded to illustrate the identification issue.~{\em Right:} Sampling distribution of the periodogram peak from the same model over 1,000 replications. We see that 
the periodogram peak is not normal, which complicates inference.}
\label{fig0}
\end{figure}

Inference on $\theta^\ast$, however, is challenging as we need to know the sampling distribution of $\hat\theta_n$. 
As argued above, normality of $\hat\theta_n$ ---as it is commonly assumed in practice--- is not plausible without  strong assumptions on the observation times, especially in finite samples.
Indeed, in this example, we see that the simulated sampling distribution of $\hat\theta_n$ is far from normal, and instead looks discrete with three distinct modes, as shown in Figure~\ref{fig0}~(right). The empirical density 
over 1,000 simulations is shown below~(rounded to two significant digits):
\renewcommand{\arraystretch}{1.2}
\begin{table}[h!]
\begin{tabular}{c | lllllllll}
$\hat\theta_n$   & 0.37  & 0.59   & 0.77  & 1.41   & 1.42   & 3.40   & 3.41  & 3.42  & 3.43  \\
\hline
$\mathrm{pr}(\hat\theta_n | \theta_0)$ & 0.4\% & 15.3\% & 0.2\% & 55.6\% & 10.9\% & 1.6\% & 7.6\% & 7.8\% & 0.6\%
\end{tabular}
\end{table}
\end{example}

We see that the sampling distribution of our estimator is non-regular. 
This complicates inference as it is unclear, for example, how to construct confidence intervals for the true period. In this paper, we propose a method that can transform such sampling distributions into 
valid confidence sets that can still cover the true period with a target coverage probability.

\begin{example}[Extinction cycles]
\label{example2}
\em
In a series of influential articles,~\citet{raup1986periodic} and \citet{raup1986biological} analyzed 
data from biological extinction events, and attempted to estimate their periodicity.
Their approach did not use the periodogram, but it modeled the event times at the most prominent extinction events directly for a best-matching period. Significance was determined through a resampling scheme where all the  recorded extinction times (not just the prominent ones) were permuted.~Through this method,~\citet{raup1986periodic} determined that there exists a biological extinction cycle of roughly 26 million years.
However, this result may be an artifact of the small sample size (only 40 observations), and the unequal spacing of observation times.~Indeed, \citet{stigler1987substantial} showed that the same signal at 26 million years appears even in simulated data where 
the observation times are fixed but the actual data are random.
We analyze the extinction data in Section~\ref{sec:extinction} using the method proposed in this paper, and present a statistically principled complement to the heuristic analysis in~\citep{stigler1987substantial}.
\end{example}

\begin{example}[Planet around~$\alpha$ Centauri B?]
\label{example3}
\em
In a widely-publicized article,~\citet{dumusque2012earth} reported the detection of an Earth-mass planet around 
the star $a$ Centauri B at a period of 3.24 days. This detection was startling because $\alpha$ Centauri is the closest stellar system to our own system, implying 
that habitable exoplanets may be ubiquitous. However, the detection has not yet been confirmed, 
and there have been some early objections to it~\citep{hatzes2013radial}. Notably,~\citet{rajpaul2015ghost} followed a methodology similar to~\citep{stigler1987substantial} discussed in~Example~\ref{example2}, 
and showed that the 3.24-day signal also appears as a local peak in various simulations with quasi-random data and fixed observation times. In this paper, we analyze these data in Section~\ref{sec:exoplanet} and show that the statistical evidence is weak in the sense that the observed data could be explained by many other candidate signals. Moreover, in Section~\ref{sec:design} we 
discuss how the observation design could be improved to identify the 
3.24-day signal, assuming that the signal is true.
\end{example}

We emphasize that neither bootstrap- nor Bayesian-based methods have so far addressed this type of identification issues. 
On the one hand, bootstrap procedures require random observation times and convergence in distribution for $\hat\theta_n$, which requires strong, unrealistic assumptions on the design. Moreover, the small-sample performance of bootstrap procedures with time series is theoretically not well understood. This is particularly important in exoplanet detection as we usually have to work with only a few hundred observations.
Bayesian methods, on the other hand, have the problem of prior specification. 
For example, even through they are usually referred to as ``non-informative", uniform priors 
give preference to parameter regions that not only have high likelihood but  are also wide. 
In effect,  this sweeps the identification problem ``under the rug" as narrow parameter regions that cause identification problems~(such as the region around $\theta =0.59$ in Example~\ref{example1}) are washed away; see also~\citep[Section 1]{hall2003nonparametric} for similar arguments that invoke number theory, and~\citep[Section 3]{gelman2020holes} for a more general discussion on how 
uniform priors can unintentionally lead to sharp  inference.
In a Bayesian context, this issue also spills over to model selection through Bayes factors in the sense that 
 model decisions based on Bayes factors may strongly depend on  features that are esoteric to the selected models~(e.g., size or structure of the parameter space); see also~\citep[Section 6]{gelman2020holes} for details.
 
\subsection{Our contribution}
In this paper, we propose a set identification method to infer hidden periodicity. 
The main idea is to develop a sharp test for a null hypothesis on a given periodicity value, 
$H_0: \theta_0 = \theta$, and then invert the test to build a confidence set, $\Theta_{1-\alpha}$. Similar to a confidence interval, a confidence set 
is constructed so that it contains the true parameter with a pre-specified coverage probability:
\begin{align}\label{eq:thetan}
\mathrm{pr}(\theta^\ast \in \Theta_{1-\alpha} ) \ge 1-\alpha.
\end{align}
In contrast to confidence intervals, confidence sets are not  necessarily contiguous.
The use of confidence sets is therefore more suitable here since 
it does not require asymptotic normality, or even consistency, of periodogram peaks.
Furthermore, the ``topology" of $\Theta_{1-\alpha}$ can reveal identification issues based on, say, 
how wide or narrow the confidence set is.

We propose two constructions for $\Theta_{1-\alpha}$. 
The first construction is idealized and yields a finite-sample valid procedure, where~\eqref{eq:thetan} holds exactly in the sample~(Theorem~\ref{thm1}). 
This construction is computationally challenging, in practice, and so we also propose an alternative  construction, where~\eqref{eq:thetan} holds asymptotically~(Theorem~\ref{thm3}). This approximate construction may be conservative  but in empirical applications we show this not to be a major concern.
Our approach borrows ideas from randomization inference in statistics~\citep{lehmann2006testing, rosenbaum2010design}, where inference is  based on structural assumptions (e.g., symmetries) on the data, and is therefore robust~\citep{toulis2019life}. Our approach is also related to partial identification in econometrics~\citep{manski2003partial, romano2008inference, tamer2010partial}, where test inversion is used to build confidence sets for unknown parameters that cannot be point-identified.


While several heuristic techniques have been used to diagnose identification issues in periodicity estimation~(see Examples~\ref{example1} and~\ref{example2}), 
our paper is the first to develop a valid, principled statistical procedure for inference when such issues are present. 
The proposed procedure is described in detail in Section~\ref{sec:concrete}, and can in fact work with many different choices of the statistic (other than periodogram peaks), and can be parallelized for computational efficiency.
Moreover, our procedure does not require normality or i.i.d.~errors, 
and can seamlessly work with equally or unequally spaced data.
We note that while some important non-parametric methods with unequally spaced data have recently been proposed, they still require observation designs with equal spacings~\citep{sun2012nonparametric, vogt2014nonparametric} or
i.i.d.~spacings between observation times~\citep{hall2000nonparametric, hall2008nonparametric, hall2003nonparametric}, and emphasize 
asymptotic consistency over finite-sample performance.
That said, our procedure could leverage these methods in the choice of the test statistic.
This could combine the best of both worlds, namely the finite-sample validity and robustness from our randomization-based procedure, and the flexibility and power from nonparametric methods.
See~Remark~\ref{remark:statistic} and Section~\ref{sec:np} for more discussion.

Empirically, we illustrate our method in several examples from exoplanet detection using radial velocity data.
In this context, we show that our method gives sharp inference on the confirmed exoplanets in our sample, namely 51 Pegasi b~\citep[``51Pegb"]{mayor1995jupiter} and Gliese 436 b~\citep[``GJ436b"]{butler2004neptune}, showing that our inference is not conservative.
For other, yet unconfirmed, exoplanets in our sample our method gives more mixed results.
For example, we show that the statistical evidence is notably weak for a candidate discovery 
around $\alpha$ Centauri B~\citep{dumusque2012earth} in the sense that the confidence set, $\Theta_{1-\alpha}$, contains a large range of signals, some possibly non-planetary, including nuisance 1-day signals that are usually associated with the nightly cycle of astronomical observations. 
Testing another widely-publicized discovery in Proxima Centauri~\citep{anglada2016terrestrial}, 
we find that the discovery seems to be robust apart from a nuisance signal at 0.916 days that cannot be rejected at the 1\% level.
 Importantly, our method suggests ways to improve the observation designs for sharper inference. Specifically, we quantify how much additional randomness is needed in the observation times or, alternatively, how many more observations are needed to precisely identify a candidate signal that is suspected to be true.
 Following this approach, we propose designs that could address the aforementioned identification issues, and also help with future scientific discoveries.

\section{Problem Setup and Notation}\label{sec:setup}
We observe data $Y^{\nnn} = \{ y_1, \ldots, y_n\}\in\mathbb{R}^n$ at observation times  
$T^{\nnn} = \{t_1, t_2, \ldots, t_n\}$. For example, in exoplanet detection, the data
could represent measured radial velocity of a star at a time $t$ expressed as a barycentric Julian date.
We decompose the data model as $y_i = y^{\mathrm{p}}(t_i) + \varepsilon(t_i)$, where $y^{\mathrm{p}}(t)$ is an 
unobserved periodic, non-stochastic function with unknown period $\theta^\ast$, such that 
$$
y^{\mathrm{p}}(t+\theta^\ast) = y^{\mathrm{p}}(t),~\quad(t\in\mathbb{R}).
$$
The terms $\varepsilon(t_i)$ correspond to the unobserved errors.  We use 
$Y^{\nnn,\mathrm{p}}=(y^{\mathrm{p}}(t_1), \ldots, y^{\mathrm{p}}(t_n))$ and $\eps^{\nnn} = (\varepsilon(t_1), \ldots, \varepsilon(t_n))$ 
to denote the (unobserved) vector values of these two components.

Our main goal is to do inference on $\theta^\ast$, the true unknown period. 
We  will need the following additional assumptions.
First, we require a consistent estimator of $y^{\mathrm{p}}$, which could either be 
parametric  or nonparametric.
For simplicity, in the following section we start with the simple but popular trigonometric regression model for $y^{\mathrm{p}}$. 
Later, in Section~\ref{sec:np}, we also discuss a nonparametric approach.
Second, we assume that the errors, $\varepsilon(t)$, form a stochastic process satisfying the following 
invariance property. For any observation times $T^{\nnn} = \{t_1, \ldots, t_n\}$,  with $n$ finite, there exists 
a known algebraic group $\mathcal{G}^{\nnn}$ of $n\times n$ matrices such that 
\begin{align}\label{A1}
\gg\cdot \eps^{\nnn}\eqd \eps^{\nnn}~\mid T^{\nnn},~\quad(\gg\in\mathcal{G}^{\nnn}).\tag{A1}
\end{align}

This formulation follows the framework of randomization tests~\citep[Chapter 15]{lehmann2006testing}, 
where testing is based on structural assumptions rather than analytical.
To keep things simple, for the rest of this paper we assume that $\mathcal{G}^{\nnn} = \mathbb{D}^n$, the set of $n\times n$ diagonal matrices with $\pm 1$ in the diagonal. 
This definition implies that the regression errors, $\varepsilon(t)$, are assumed to have a sign symmetric distribution.
As such, our inference method not only covers Gaussian errors, which are frequently assumed in practice, 
but also extends to more general model classes where
the error distribution is unspecified, or even varies with $t$, as long as the symmetry in~\eqref{A1} holds.



\begin{remark}[Generality of~\eqref{A1}]
\label{rem21}
\em
Assumption~\eqref{A1} can also allow inference with complex error structures.
Suppose, for instance, that $\varepsilon(t)$ is a
 Gaussian process~(GP) with covariance matrix $C$. 
Then, we could define $\mathcal{G}^{\nnn} = \{ \gg' = C^{1/2} \gg C^{-1/2} : ~\gg\in\mathbb{D}^{n}\}$.
Of course, this requires $C$ to be known, or be consistently estimated.
To accomplish that, one idea would be to combine our approach with 
recent GP-based methods on radial velocity data~\citep{rajpaul2015gaussian, damasso2017proxima,faria2016uncovering,faria2020decoding}. We leave this for future work.~$\blacksquare$
\end{remark}

\section{Main Method}\label{sec:method}
Here, we follow a parametric approach to model the periodic component, $y^{\mathrm{p}}$. 
A nonparametric approach is discussed later, in Section~\ref{sec:np}.
All proofs can be found in the Supplement.

We start with the simple harmonic model 
\begin{align}\label{A2}
y^{\mathrm{p}}(t) \equiv y^{\mathrm{p}}(t; \theta^\ast, \psi^\ast) = \psi^\ast_{1} + \psi^\ast_{2} \cos(2\pi t/\theta^\ast) + \psi^\ast_{3} 
\sin(2\pi t/\theta^\ast),
\tag{A2}
\end{align}
where $\theta^\ast\in\Theta\subseteq\mathbb{R}^+$ is the true 
unknown period, and $\psi^\ast = (\psi^\ast_1, \psi^\ast_2, \psi^\ast_3) \in \Psi\subseteq\mathbb{R}^3$ are unknown nuisance parameters, 
where $\Psi$ is convex and compact. We note that we use this model for simplicity. 
For instance,  we could work with, say, a Keplerian model that includes additional parameters to describe various patterns of motion (e.g., eccentricity) of a body relative to another~\citep{wright2009efficient}.
All that is required for our theory to work is for the model of choice to be consistently estimated conditional on 
the true period, $\theta^\ast$.

To measure the fit of $y^{\mathrm{p}}$ we use the simple squared loss,
$$
L(\theta, \psi \mid Y^{\nnn}, T^{\nnn}) = \sum_{i=1}^n [y_i - y^{\mathrm{p}}(t_i | \theta, \psi)]^2 / \sigma_i^2,
$$
where the precisions $\sigma_i^2$ are known from instrument specifications, as is typical in exoplanet detection.
Note that the loss function~(i.e., negative log-likelihood under Gaussian errors) is non-convex and non-smooth with respect to $\theta$, and so standard statistical inference cannot be applied.

To proceed, let $L_{0n}(Y^{\nnn}, T^{\nnn}) = \sum_{i=1}^n (y_i - \bar y)^2/\sigma_i^2$ denote the baseline loss, where 
$\bar y=(1/n) \sum_{i=1}^n y_i$.
 Then, the {\em oracle} generalized periodogram, $A_n: \Theta\to\mathbb{R}^+$ can be defined as follows:
\begin{align}\label{eq:oracle_per}
A_n(\theta \mid Y^{\nnn}, T^{\nnn}, \psi) = \frac{L_{0n}(Y^{\nnn}, T^{\nnn})  - L(\theta, \psi \mid Y^{\nnn}, T^{\nnn}) } { L_{0n}(Y^{\nnn}, T^{\nnn})  }.
\end{align}
In practice, $\psi$ is not the target of inference, and so we use the profile likelihood
\begin{align}\label{eq:real_per}
\hat A_n(\theta \mid Y^{\nnn}, T^{\nnn}) \equiv A_n(\theta \mid Y^{\nnn}, T^{\nnn}, \hat\psi_\theta),
~\text{where}~\hat\psi_\theta = \arg\min_{\psi\in\Psi} L(\theta, \psi \mid Y^{\nnn}, T^{\nnn}).
\end{align}
$\hat A_n$ is known as the generalized periodogram, and is easy to calculate as $\hat\psi_\theta$ is known in closed form~\citep{zechmeister2009generalised}. Our first goal is to test a global null hypothesis of the form:
\begin{align}\label{eq:H01}
H_0^{\mathrm{full}}: \theta^\ast=\theta_0, \psi^\ast=\psi_0.
\end{align}
This hypothesis is not particularly interesting because, ultimately, we are only interested in the unknown period $\theta^\ast$. However, testing for $H_0^{\mathrm{full}}$ will serve as a good starting point for testing and inference on $\theta^\ast$.
Of course, to test any hypothesis we need a test statistic, $s_n(Y^{\nnn}, T^{\nnn})\in\mathbb{R}$. 
Many choices are possible here, and so we defer a detailed discussion about this choice to Section~\ref{sec:concrete}.
Here, we simply define
$$
 s_n(Y^{\nnn}, T^{\nnn}) = \max_{\theta \in \Theta} A_n(\theta \mid Y^{\nnn}, T^{\nnn}, \psi_0) - A_n(\theta_0 \mid Y^{\nnn}, T^{\nnn}, \psi_0).
$$
The idea behind this choice of the test statistic is that $s_n$ will tend to take small values when $\theta^\ast=\theta_0$, as we expect the periodogram to peak around $\theta_0$ under $H_0^{\mathrm{full}}$; 
when $H_0^{\mathrm{full}}$ is not true we expect $s_n$ to take large values as the periodogram will peak in 
periods different from $\theta_0$, so that the null can be rejected. 

Let $s_\obs = s_n(Y^{\nnn}, T^{\nnn})$ denote the observed value of the test statistic in the sample. 
To construct the null distribution of $s_n$, we simulate data as follows:
$$
Y^{\nnn,(i)} = Y^{\nnn,\mathrm{p}}_0+ G^{(i)} \cdot (Y^{\nnn} - Y^{\nnn,\mathrm{p}}_0),
$$
where $Y^{\nnn,\mathrm{p}}_0 = [y^{\mathrm{p}}(t_1; \theta_0, \psi_0), \ldots, y^{\mathrm{p}}(t_n; \theta_0, \psi_0)]$ are the periodic component values, which are known under the null hypothesis; $G^{(i)}\sim\mathrm{Unif}(\mathcal{G}^{\nnn})$ is a random error transformation from $\mathcal{G}^{\nnn}$. 

Then, a  $p$-value for $H_0^{\mathrm{full}}$ is
\begin{align}\label{eq:pval_full}
\mathrm{pval}(\theta_0, \psi_0) = E\{s_n(Y^{\nnn,(i)}, T^{\nnn}) \ge s_\obs \mid  Y^{\nnn}, T^{\nnn}\},
\end{align}
where the randomness is with respect to $G^{(i)}$ while $Y^{\nnn}, T^{\nnn}$ are fixed. 
The following result shows that this $p$-value is valid in finite samples for testing $H_0^{\mathrm{full}}$.

\begin{theorem}\label{thm1}
\ThmOne
\end{theorem}


In most settings, $\psi^\ast$ is a nuisance parameter, and so we only want to test for $\theta^\ast$:
\begin{align}\label{eq:H02}
H_0: \theta^\ast=\theta_0.
\end{align}
In this context, the result in Theorem~\ref{thm1} can still be used, in principle, to build a confidence set for 
$\theta^\ast$ using the following construction: 
\begin{align}\label{eq:Theta1}
\Theta_{1-\alpha} = \left\{\theta \in \Theta :   \max_{\psi\in\Psi}  \mathrm{pval}(\theta, \psi) > \alpha \right\}.
\end{align}

\begin{theorem}\label{thm2}
\ThmTwo
\end{theorem}
While $\Theta_{1-\alpha}$ is a valid confidence set, its construction is challenging due to the optimization problem 
in~\eqref{eq:Theta1} that requires maximizing the $p$-value over the entire nuisance parameter space, $\Psi$. This is especially challenging when $\Psi$ is multidimensional. 
To address this computational issue, a feasible procedure is  to estimate 
$\psi^\ast$ under the null hypothesis~$H_0$ ---for example, through the least-squares estimate $\hat\psi_{\theta_0}$--- and then plug  this estimate into~\eqref{eq:Theta1}. 
We note that this approach works because the estimator $\hat\psi_{\theta_0}$ is consistent for $\psi^\ast$  under the standard assumptions of least squares since, under fixed $\theta$, the loss function in the harmonic model is smooth and convex with respect to $\Psi$.
In contrast, estimating $\theta^\ast$ in a similar fashion (e.g., through least squares) does not work, 
precisely because of the identification issues discussed in Section~\ref{sec:intro}.



To test $H_0$ in a computationally manageable way, we therefore have to simply re-define the test components using the generalized periodogram~\eqref{eq:real_per} as follows:
\begin{align} 
 s_n(Y^{\nnn}, T^{\nnn}) &  = \max_{\theta \in \Theta} \hat A_n(\theta \mid Y^{\nnn}, T^{\nnn}) - \hat A_n(\theta_0 \mid Y^{\nnn}, T^{\nnn}),
 \label{eq:re1} \\
 \hat Y^{\mathrm{p}}_0 &  = [y^{\mathrm{p}}(t_1; \theta_0, \hat\psi_{\theta_0}), \ldots, y^{\mathrm{p}}(t_n; \theta_0,\hat\psi_{\theta_0})], \label{eq:re2}\\
 Y^{\nnn,(i)}  & = \hat Y^{\mathrm{p}}_0 + G^{(i)} \cdot [Y^{\nnn} - \hat Y^{\mathrm{p}}_0 ]. \label{eq:re3}.
\end{align}
The $p$-value for $H_0$ is then defined as
\begin{align}\label{eq:pval_theta}
\widehat\mathrm{pval}(\theta_0) = E\{s_n(Y^{\nnn,(i)}, T^{\nnn}) \ge s_\obs \mid Y^{\nnn}, T^{\nnn}\},
\end{align}
where, again, the randomness is with respect to $G^{(i)}$ while the data $(T^{\nnn}, Y^{\nnn})$ are 
fixed.
The following construction for the confidence set of $\theta^\ast$ is valid asymptotically:
\begin{align}\label{Theta2}
\hat\Theta_{1-\alpha} = \left\{\theta \in \Theta :  \widehat\mathrm{pval}(\theta) > \alpha \right\}.
\end{align}

\begin{theorem}\label{thm3}
\ThmThree
\end{theorem}



\subsection{Concrete procedure and practical considerations}\label{sec:concrete}
In summary, we propose the following procedure to estimate the true periodicity.

\vspace{-5px}
\MakeProcedure{1}{
\begin{enumerate}
\item Choose a grid of possible period values, $\Theta$ that contains $\theta^\ast$ w.p. 1.
Set $\hat\Theta_{1-\alpha} \gets \emptyset$.
\item Obtain data $(Y^{\nnn}, T^{\nnn})$, possibly after removing known 
stellar signals, e.g., rotational periods, magnetic cycles, etc.~\citep[Section 11.3.4]{feigelson2012modern}.
\item {\bf For} all $\theta_0\in\Theta$ {\bf do:}
\begin{enumerate}[(i)]
\item Estimate the nuisance parameters, $\hat\psi_{\theta_0} = \arg\min_{\psi\in\Psi} L(\theta_0, \psi \mid Y^{\nnn}, T^{\nnn})$, through weighted least squares.
\item Calculate the observed value, $s_\obs$, of the test statistic 
defined in~\eqref{eq:re1}.
\item With fixed $T^{\nnn}$, sample new data, $Y^{\nnn,(i)}$, where $i=1, \ldots, R$ for some fixed $R$, according 
to definitions~\eqref{eq:re2} and \eqref{eq:re3}.
\item Using the samples  from 3(iii), calculate the $p$-value~\eqref{eq:pval_theta}, and 
if it exceeds $\alpha$ then include $\theta_0$ in the confidence set;
i.e., set $\hat\Theta_{1-\alpha} \gets \hat\Theta_{1-\alpha} \cup \{\theta_0\}$ 
if $\widehat\mathrm{pval}(\theta_0) > \alpha$.
\end{enumerate}
\item {\bf Return} $\hat\Theta_{1-\alpha}$ as the $100(1-\alpha)\%$ confidence set of $\theta^\ast$.
\end{enumerate}
}

In words, Procedure~1 considers all possible candidate period signals in $\Theta$. 
For every candidate signal, $\theta_0$, the procedure, first, estimates the harmonic model assuming that the true period is $\theta_0$, and then 
it randomly flips the signs of the residuals to generate new synthetic data~(i.e., assumes symmetric distribution for the errors).
The observed value of the test statistic is then compared against the test statistic values in the synthetic data, and the resulting $p$-value determines whether $\theta_0$ should be in the confidence set or not.

We end this section with a few remarks on some practical aspects of implementing Procedure~1.
\begin{remark}[Choice of $\Theta$]
\label{remark:Theta}
\em
The choice of $\Theta$ can be crucial as the periodogram is generally sensitive to the average spacing in $\Theta$. In general, it is a good idea to choose spacings that are smaller than the expected width of the major periodogram peaks, and also induce a spacing in $\Theta$ that is uniform in the log scale~\citep[Section 6.1]{vanderplas2018understanding}.
Another useful technique is to keep refining $\Theta$ until the shape of the observed periodogram and its major peaks are stable.~$\blacksquare$
\end{remark}

\begin{remark}[Choice of test statistic, $s_n$]
\label{remark:statistic}
\em
Procedure 1 is valid for any choice of the test statistic. Our choice of test statistic in~\eqref{eq:re1} performs well under several controlled simulations~(see Section~\ref{sec:example}), but other options are certainly possible. For instance, Fisher's classical statistic to test periodicity 
can be defined as $s_n = \max_{\theta \in \Theta} \hat A_n(\theta) 
/ \bar A_n$, where $\bar A_n = |\Theta|^{-1}\sum_\theta A_n(\theta)$ is the average periodogram value; here, we suppressed the dependence on $Y^{\nnn}, T^{\nnn}$ to ease notation. Improvements using a trimmed mean in place of 
$\bar A_n$ have also been suggested~\citep{bolviken1983new, siegel1980testing, damsleth1982estimation}, including sparse alternatives~\citep{cai2016global}. See also~\citep{mcsweeney2006comparison} for some numerical comparisons.~$\blacksquare$
\end{remark}


\begin{remark}[Power]
\label{remark:power}
\em
The power of Procedure 1 depends on how sensitive the chosen test statistic, $s_n$, is in detecting violations of the null hypothesis. In general, statistics that lead to better estimators of $\theta^\ast$ will also lead to more powerful tests, and thus more narrow confidence sets. In this paper, we don't explore the power of our randomization tests theoretically to maintain the focus on our application---we refer the reader to~\citep{toulis2019life} and~\citep{dobriban2021consistency} for such theoretical investigations on the 
power of randomization tests, including those applied in a regression context like ours.
Some optimal results in testing periodicity are known but require strong assumptions on the errors. 
For instance, Fisher's classical statistic (described above) is optimal 
only under i.i.d. normal errors in the harmonic model. Empirically, our choice of test statistic in~\eqref{eq:re1} appears to be powerful compared to alternatives---see Section~\ref{sec:example} and Section~\ref{sec:exoplanet} for details.~$\blacksquare$
\end{remark}


\begin{remark}[Computation]
\label{remark:computation}
\em 
The complexity of Procedure~1 can be described, prima facie, as $O(|\Theta|^2 \cdot R \cdot C)$, 
where $C$ is the average computation time for a single weighted least-squares fit. 
For some ballpark numbers, we can take $|\Theta| = 10^4$, $R=10^3$, and $C=50\mu$s for an analysis on a conventional laptop of a time series with 200 observation times, using an implementation in {\tt R}.
The full computation requires  a total of 1,388 hrs.~of wall clock time~(approx. 58 days).
However, several reductions of computation time are possible. 
First, Procedure 1 can be fully parallelized since step 3 can be executed independently for each $\theta_0\in\Theta$. With a moderate-sized cluster of 100 computing nodes, the 
wall clock time thus drops to 14 hrs. Moreover, again in step 3, there is no need to consider all values in $\Theta$ but only a proportion of the periodogram peaks. In our applications, for example, 
we only consider those peaks of the periodogram that are at least 20\% as high as the highest peak.
This leads to a complexity $O(\gamma |\Theta|^2 \cdot R \cdot C)$ with $\gamma$ ranging between
$0.1\%$-$3\%$ in our applications.
As such, the required computation drops dramatically to roughly 30 mins.~of wall clock time.
Indeed, in our application, we are able to ramp up our samples to $R=100,000$ and still 
finish all exoplanet analyses within a couple of hours using a cluster of 400 computing nodes.
Further significant gains are possible if the method is implemented in a more efficient language than {\tt R}, such as 
{\tt C++} or {\tt Python}.~$\blacksquare$
\end{remark}

\section{Numerical Examples}\label{sec:example}

\subsection{Extinction data}\label{sec:extinction}
Here, we return to Example~\ref{example2} on estimating periodicity in biological extinction events.
The original data are shown in the left subplot of Figure~\ref{fig:extinction} below.
~\citet{raup1986periodic} argued that there is a 26 million year periodicity based on fitting 
a simple periodic model on the largest extinction events. Formally, for any given $\theta\in\Theta$, they fit a model $t_i = t_0 + k_\theta \theta$ for those timepoints that correspond to the eight largest peaks; here $t_0, k_\theta$ are free parameters and $k_\theta$ is an integer.
Using this model,~\citet{raup1986periodic} calculated which $\theta$ yields the smallest least squares fit, and then used a permutation-based test to assess significance.

The periodogram of the data is shown in the right subplot of Figure~\ref{fig:extinction}.
Indeed, using a set $\Theta$ ranging between 12.5-62.5 million years, as suggested by~\citet{raup1986periodic}, we find the largest periodogram peak to be around 26 million years.
However, we also see additional peaks around 30 million years and 35 million years.
Using Procedure~1, none of these values can be rejected as the true periodicity~($p$-values are 0.56 and 0.6, respectively).\footnote{All data and code for the numerical examples of this section, including the exoplanet detection applications of the next sections, are  provided with the Supplementary Code. 
An online repository with {\tt R} code to replicate all results and figures in the following sections can be found at \url{https://github.com/ptoulis/ri-exoplanet-detection}. 
A {\tt Python} implementation can be found at \url{https://github.com/afeinstein20/periodic-ri}.
} Moreover, a similar analysis with a larger range for $\Theta$~(not shown in the figure) reveals multiple peaks ---at periods less than 0.5 million years--- that are even larger than the peaks shown in Figure~\ref{fig:extinction}, and are accepted in the confidence set.
We conclude that there are significant identification issues with this dataset, and so the discovery cannot be confirmed. These issues were also highlighted by~\citet{stigler1987substantial}, albeit 
using a heuristic approach.

\begin{figure}[h!]
    \centering
    \begin{minipage}{0.50\textwidth}
        \centering
        \includegraphics[width=1.08\textwidth]{fig/extinction1}
    \end{minipage}\hfill
    \begin{minipage}{0.5\textwidth}
        \centering
        \includegraphics[width=1.08\textwidth]{fig/extinction2}
    \end{minipage}
          \caption{\small
          {\em Left:}~Original extinction data of~\citet{raup1986periodic}. 
          The data include 40 observations of extinction events measured as \% of estimated biomass lost during the event.~{\em Right:}~
The periodogram of the extinction data in the range 12.5-62.5 million years, as suggested by~\citet{raup1986periodic}.}
\label{fig:extinction}
\end{figure}

\subsection{Simulated data with realistic observation designs}
Here, we revisit Example~\ref{example1} but with a more realistic design for the observation times.
Specifically, we build upon  the  setup of~\citet[Section 3.2]{hall2006using} and define the data model:
\begin{align}\label{eq:sim1_Y}
y_i = 1-\cos(2\pi t_i/\theta^\ast) + N(0, \sigma^2),~\quad(i=1,\ldots, n).
\end{align}
Here, the design of observation times relies on a density $f_1$ on the interval $[0, 1]$ that represents a full 24-hour day, with $[0, 0.5]$ denoting day time and $[0.5, 1]$ denoting night time~($t=.75$ is midnight, and so on).  Adapting the setup of~\citet{hall2006using}, we consider three choices for $f_1$:

\begin{enumerate}[(i)]
\item  $f_1(t) = \mathbbm{1}(.5 \le t \le  1)$;
\item  $f_1(t) = \max\{0, -\pi \sin(2\pi t)\} \mathbbm{1}(0\le  t \le 1)$;
\item  $f_1(t) = \mathbbm{1}( | t- 0.75| \le 0.0208)$;
\end{enumerate}
The full observation design is then defined as:
\begin{align}\label{eq:sim1_T}
f(T^{\nnn}) \propto \prod_{i=1}^n \left(\sum_{j=0}^{n-1} f_1(t_i - j \theta^\ast_{\mathrm{obs}})\right).
\end{align}


We set $\theta^\ast=\sqrt{2}$ as the true signal period with $\sigma^2 =1.5^2$, and $\theta^\ast_{\mathrm{obs}}=1$ for
the observation period. Thus, the particular form of $f(T^{\nnn})$ in~\eqref{eq:sim1_T} represents a density that has daily periodicity in the interval $[0,n]$, which represents $n$ days of observations.
The definitions of $f_1$ are meant to represent certain recognizable classes of observation designs. 
Specifically, in design (i) the observation times $t_i$ are sampled uniformly during night time.
In design (ii), we define a more realistic scenario where observations happen only at night time, but are more probable during midnight. Finally, in design (iii) the observations are uniformly sampled within one hour window around midnight.

Figure~\ref{fig1} shows the sampling distribution of $\hat\theta_n = \arg\max_{\theta\in\Theta} \hat A_n(\theta \mid Y^{\nnn}, T^{\nnn}, \hat\psi_\theta)$, the periodogram peak, where $(T^{\nnn}, Y^{\nnn})$ are sampled, respectively, according to~\eqref{eq:sim1_T} and~\eqref{eq:sim1_Y}, 
and $\hat A_n$ is the generalized Lombe--Scragle periodogram in~\eqref{eq:real_per}.
As discussed in Section~\ref{sec:intro}, we see that the sampling distribution is far from normal across all designs, 
and is better described as an irregular, discrete distribution on $\Theta$.
The amount of irregularity depends on the observation design.
Under Design (i), for example, $\hat\theta_n$ has converged to a point distribution around $\theta^\ast$ at $n=500$
data points, but there is significant irregularity even at $n=100$.
Under Design (iii), the situation is much worse since the distribution of $\hat\theta_n$ is irregular 
even at $n=500$ with three distinct modes similar to the mode around $\theta^\ast$.
This example illustrates one of the key motivations of this paper, which is that standard asymptotic normal theory 
is inappropriate for inference of periodicity in time series due to the 
irregularity of periodic models.

\begin{figure}[t]
\centering
\includegraphics[scale=0.53]{fig/fig1_sampling2}
\caption{\small
Sampling distribution of periodogram peak~(approximated via 1,000 samples) in the model described by
~\eqref{eq:sim1_T} and~\eqref{eq:sim1_Y}.}
\label{fig1}
\end{figure}

\begin{figure}[h!]
\begin{floatrow}
\centering
\includegraphics[scale=0.315]{fig/fig_periodogram}

\capbtabbox{%
  \begin{tabular}{r | r cc}
  $\theta_0$ & $p$-value & $\hat\Theta_{0.95}$ & $\hat\Theta_{0.99}$ \\
   \hline
   0.1216  &  0.0645& yes & yes  \\
   0.1357 &  0.0465 &  no &  yes \\
   0.1385  & 0.0625 &   yes &   yes \\
   0.1998  & 0.0570 &   yes &   yes  \\
   0.2183  & 0.0680 &   yes &   yes  \\
   0.2792  & 0.0440 &  no &   yes  \\
   {\bf 0.5858}  & 0.6460 &   yes &   yes  \\
  0.7787  & 0.0760 &  yes &   yes  \\
 1.4136  & 0.0955 &   yes &   yes  \\
   3.1546  & 0.0675 &  yes &   yes  \\
 9.6474 &  0.0720 & yes & yes \\
 12.7177 & 0.0600 & yes & yes
   \end{tabular}
}{%
}
\end{floatrow} 
\caption{\small
{\em Left:} Periodogram for a problematic dataset sampled according to~\eqref{eq:sim1_T} and~\eqref{eq:sim1_Y}, where the peak $\hat\theta_n=0.586$ is not at the true period~$\theta^\ast=1.414$.
{\em Right:}~Inference of periodicity based on Procedure~1. The table 
shows the $p$-values for the hypothesis $H_0:\theta^\ast=\theta_0$ for values of $\theta_0$ that correspond 
to the highest peaks of the periodogram shown on the left.
}
\label{fig2}
\end{figure}


Figure~\ref{fig2} shows how to address this problem using  the identification method proposed in this paper. Specifically, following Example~\ref{example1}, we seed a dataset $(T^{\nnn}, Y^{\nnn})$ such that estimation is problematic 
in the sense that the observed periodogram peak is not at the true value, $\theta^\ast$~(see Figure~\ref{fig2}, left) ---in this dataset $\hat\theta_n=0.586$. 
The table on the right of Figure~\ref{fig2} shows the $p$-values for testing 
$H_0:\theta^\ast=\theta_0$ obtained from Procedure 1. For simplicity, we only show results 
for the 10 highest peaks in the observed periodogram.
We see that all these values cannot be rejected at the 1\% level and so they are included in the 
99\% confidence set, $\hat\Theta_{0.99}$~(see column ``$\hat\Theta_{0.99}$" in the table). 
By looking at the range of values included in this confidence set, we conclude that there are 
severe identification issues in the data. This was not immediately evident from the periodogram as the peak in Figure~\ref{fig2} is fairly distinct and far exceeds the 1\% false alarm probability threshold~(not shown above). We also note that, beyond the mode at $\theta_0=0.586$, our proposed method puts the highest confidence at the value $\theta_0=1.4136$~($p$-value = 0.096), which is close to the true signal.

Finally, we check whether the finite-sample coverage of $\hat\Theta_{1-\alpha}$ used in Procedure~1 
matches the asymptotic nominal coverage, as derived in Theorem~\ref{thm3}. 
For this simulation we consider all aforementioned designs ($f_1$), and also multiple sample sizes ($n$), and observation times restricted to a 6-month observation period. For each setting we calculate how often $\hat\Theta_{.95}$ covers the true period, $\theta^\ast$, over 4,000 replications for each $(f_1, n)$-pair. The results are shown in the table below. 
\renewcommand{\arraystretch}{1.1}
\begin{table}[h!]
\begin{tabular}{c | ccccc}
Design, $f_1$ & \multicolumn{5}{c}{Sample size, $n$} \\
& 25 & 50   & 100  & 250  & 500  \\
\hline
(i)  & 97.97\% &  97.32\%   & 96.05\% & 94.45\% &  95.27\% \\
(ii) & 97.82\% & 97.65\% & 95.42\%  & 94.50\%   & 97.50\%  \\
(iii) & 97.72\% & 98.85\% & 96.95\% & 93.40\% & 93.60\%
\end{tabular}
\label{tab:cover}
\end{table}


We note that coverage gets closer to the nominal value 
for larger sample sizes, which is consistent with Theorem~\ref{thm3}. Moreover, 
we see that better observation designs generally lead to better coverage.

\section{Application: Exoplanet Detection}\label{sec:exoplanet}
In this section, we use our method to analyze radial velocity data from four known star systems. 
Specifically, we first analyze data on two exoplanets, namely 51 Pegasi b~\citep[``51Pegb"]{mayor1995jupiter} and Gliese 436 b~\citep[``GJ436b"]{butler2004neptune}.
These exoplanets have been confirmed independently by other means, including
direct observations of transit.
We also analyze data from two candidate exoplanets potentially orbiting the stars $\alpha$ Centauri B~\citep{dumusque2012earth}, and Proxima Centauri~\citep{anglada2016terrestrial}, respectively.

We start our analysis with the confirmed exoplanets in our sample. The periodogram of the stellar radial velocity data for 51Pegb is depicted in Figure~\ref{fig3}. 
This shows a clear peak around $\hat\theta_n=4.23$ days with additional 
distinct peaks around 1-day and 0.5-day periods.
Our method correctly shows no identifiability issues here.
All periods corresponding to the highest peaks are strongly rejected ($p$-values numerically zero) except for 
$\theta_0=4.23$ (observed peak).
Next, we analyze data on GJ436b. Compared to 51Pegb, the periodogram for GJ436b shows  a peak around $\hat\theta_n=2.6441$, but this peak is less distinct than other peaks. Our method again correctly shows  that the true signal can be easily identified.
As in 51Pegb,  all periods corresponding to the highest peaks are strongly rejected ($p$-value = 0) except for the observed peak at $\theta_0=2.64$.

\begin{figure}[t!]
\begin{floatrow}
\centering
\includegraphics[scale=0.36]{fig/fig_51pegb}

\capbtabbox{%
  \begin{tabular}{r | r cc}
  $\theta_0$ & $p$-value & $\hat\Theta_{0.95}$ & $\hat\Theta_{0.99}$ \\
   \hline
0.3085 & 0.0000 & no & no \\ 
  0.5662 & 0.0000 & no & no \\ 
  0.8069 & 0.0000 & no & no \\ 
  0.8089 & 0.0000 & no & no \\ 
  0.8295 & 0.0000 & no & no \\ 
  1.3047 & 0.0000 & no & no \\ 
  1.3095 & 0.0000 & no & no \\ 
  3.7033 & 0.0000 & no & no \\ 
  4.1807 & 0.0000 & no & no \\ 
 {\bf  4.2311} & 1.0000 & yes & yes \\ 
  4.2821 & 0.0000 & no & no \\ 
  4.9331 & 0.0000 & no & no \\ 
   \end{tabular}
}{%
}
\end{floatrow} 
\caption{\small
{\em Left:}~Periodogram of radial velocity on exoplanet ``51Pegb". Here, $\Theta = \{0.1, \ldots, 1000\}$ is split uniformly in the log-space so that 
$|\Theta| = 25,000$. 
{\em Right:}~Inference of periodicity of 51Pegb~based on Procedure~1.
The table shows the $p$-values for the hypothesis $H_0:\theta^\ast=\theta_0$ for values of $\theta_0$ that correspond to high peaks of the periodogram shown on the left.
We see that there are no identification issues, as the 4.23-day signal is the only one accepted 
in the confidence sets.
\vspace{25px}
}
\label{fig3}
\end{figure}




\begin{figure}[h!]
\begin{floatrow}
\centering
\includegraphics[scale=0.36]{fig/fig_GJ436b}

\capbtabbox{%
  \begin{tabular}{r | r cc}
  $\theta_0$ & $p$-value & $\hat\Theta_{0.95}$ & $\hat\Theta_{0.99}$ \\
   \hline
0.4200 & 0.0000 & no & no \\ 
  0.6155 & 0.0000 & no & no \\ 
  0.7067 & 0.0000 & no & no \\ 
  0.7438 & 0.0000 & no & no \\ 
  1.3641 & 0.0000 & no & no \\ 
  1.5187 & 0.0000 & no & no \\ 
  1.6013 & 0.0000 & no & no \\ 
  1.6086 & 0.0000 & no & no \\ 
  1.7008 & 0.0000 & no & no \\ 
  2.4103 & 0.0000 & no & no \\ 
 {\bf  2.6441} & 1.0000 & yes & yes \\ 
  3.7092 & 0.0000 & no & no \\ 
   \end{tabular}
}{%
}
\end{floatrow} 
\caption{\small
{\em Left:}~Periodogram of radial velocity on exoplanet ``GJ436b". Here, $\Theta = \{0.1, \ldots, 1000\}$ is split uniformly in the log-space so that 
$|\Theta| = 30,000$. 
{\em Right:}~Inference of periodicity of GJ436b~based on Procedure~1.
The table shows the $p$-values for the hypothesis $H_0:\theta^\ast=\theta_0$ for values of $\theta_0$ that correspond to high peaks of the periodogram shown on the left.
We see that there are no identification issues, as the 2.64-day signal is the only one accepted 
in the confidence sets.
\vspace{20px}
}
\label{fig4}
\end{figure}


Our results are mixed for the other, yet unconfirmed exoplanets.
For the candidate exoplanet orbiting $\alpha$ Centauri B~\citep{dumusque2012earth}, in particular,
our analysis ---summarized in Figure~\ref{fig5}--- shows severe identification issues.
The periodogram shows a peak around $\hat\theta_n=3.23$ days but it is not as distinct 
from the other peaks. Moreover, there is significant variation in the periodogram. 
Indeed, our method reveals that many periods corresponding to the highest peaks cannot be rejected at the 1\% level, and there is even a period ($\theta_0=0.7622$) that cannot be rejected  at the 5\% level.
These results indicate that additional statistical evidence is necessary to claim that the particular exoplanet discovery is plausible.
As in the case of the extinction data, we note that  earlier work 
has raised concerns about this discovery as well~\citep[e.g.]{rajpaul2015ghost}, 
albeit using only  heuristic methods similar to those suggested by~\citet{stigler1987substantial}.



Finally, we discuss the data for a recent discovery of an exoplanet ostensibly orbiting 
Proxima Centauri~\citep{anglada2016terrestrial}.
We use a recent release of 63 observations of Proxima Centauri from the ESPRESSO spectograph~\citep{mascareno2020revisiting} that has increased precision.
The analysis and results are shown in Figure~\ref{fig6}. 
The periodogram shows  a peak around $\hat\theta_n=11.17$ days but also 
a substantial peak around the 1-day period signal, possibly  associated with the observation design.  Indeed, our method shows that out of all highest peaks in the periodogram, only 
$\theta_0=11.1739$ and $\theta_0=0.9164$ cannot be rejected at the 1\% level. 
Overall, while the finding of \citet{anglada2016terrestrial} seems to be robust, there is still a small chance of observing the observed periodogram pattern under the particular observation times, at least with respect to those aspects of the periodogram that are captured by our test statistic.
More statistical evidence is therefore needed to statistically validate this discovery, 
possibly through the use of a slightly more randomized observation design. 
We discuss such designs in the  following section.

\begin{figure}[t!]
\begin{floatrow}
\centering
\includegraphics[scale=0.35]{fig/fig_dumusque}

\capbtabbox{%
  \begin{tabular}{r | r cc}
  $\theta_0$ & $p$-value & $\hat\Theta_{0.95}$ & $\hat\Theta_{0.99}$ \\
   \hline
0.7622 & 0.0705 & yes & yes \\ 
  0.8882 & 0.0271 & no & yes \\ 
  1.0086 & 0.0174 & no & yes \\ 
  1.0678 & 0.0079 & no & no \\ 
  2.0292 & 0.0122 & no & yes \\ 
  3.2074 & 0.0163 & no & yes \\ 
  {\bf 3.2371} & 1.0000 & yes & yes \\ 
  3.2670 & 0.0178 & no & yes \\ 
  7.9394 & 0.0116 & no & yes \\ 
  8.1169 & 0.0175 & no & yes \\ 
  52.2242 & 0.0121 & no & yes \\ 
  61.1334 & 0.0226 & no & yes \\ 
   \end{tabular}
}{%
}
\end{floatrow} 
\caption{\small
{\em Left:}~Periodogram of radial velocity on candidate exoplanet orbiting $\alpha$ Centauri B~\citep{dumusque2012earth}.  Here, $\Theta = \{0.1, \ldots, 1000\}$ is split uniformly in the log-space, so that 
$|\Theta| = 10,000$. 
{\em Right:}~Inference of periodicity of the exoplanet~based on Procedure~1.
The table shows the $p$-values for the hypothesis $H_0:\theta^\ast=\theta_0$ for values of $\theta_0$ that correspond to high peaks of the periodogram shown on the left.
We see that there are severe identification issues as several signals other than the periodogram peak are accepted in the confidence sets.\vspace{20px}}
\label{fig5}
\end{figure}



\begin{figure}[h!]
\begin{floatrow}
\centering
\includegraphics[scale=0.35]{fig/fig_ESPRESSO}

\capbtabbox{%
  \begin{tabular}{r | r cc}
  $\theta_0$ & $p$-value & $\hat\Theta_{0.95}$ & $\hat\Theta_{0.99}$ \\
   \hline
0.1106 & 0.0007 & no & no \\ 
  0.3355 & 0.0022 & no & no \\ 
  0.3552 & 0.0055 & no & no \\ 
  0.4778 & 0.0025 & no & no \\ 
  0.5512 & 0.0047 & no & no \\ 
  0.7532 & 0.0052 & no & no \\ 
  0.8412 & 0.0059 & no & no \\ 
  0.9164 & 0.0173 & no & yes \\ 
  0.9266 & 0.0005 & no & no \\ 
  1.0957 & 0.0080 & no & no \\ 
  {\bf 11.1739} & 1.0000 & yes & yes \\ 
  12.8769 & 0.0006 & no & no \\ 
   \end{tabular}
}{%
}
\end{floatrow} 
\caption{
\small
{\em Left:}~Periodogram of radial velocity on candidate exoplanet Proxima Centauri b~\citep{anglada2016terrestrial}.  Here, $\Theta = \{0.1, \ldots, 1000\}$ split regularly in the log-space, so that 
$|\Theta| = 10,000$. 
{\em Right:}~Inference of periodicity of the exoplanet~based on Procedure~1.
The table shows the $p$-values for the hypothesis $H_0:\theta^\ast=\theta_0$ for values of $\theta_0$ that correspond to high peaks of the periodogram shown on the left.
We see that there are no severe identification issues, and the detection appears to be robust.
There is only only signal (at 0.9164 days) that cannot be rejected at the 1\% level, and thus causes concern.
An improved observation design could decisively rule out this nuisance signal; see Section~\ref{sec:design} for 
more details.\vspace{20px}
}
\label{fig6}
\end{figure}




\section{Discussion and Extensions}\label{sec:extensions}

\subsection{Observation design}\label{sec:design}

The importance of observation times in identifying a periodic signal is well understood~\citep{feigelson2012modern, vanderplas2018understanding, ivezic2014statistics}.
It may therefore come as a surprise that  there is little systematic work on optimizing observation designs.
One exception is perhaps~\citep{loredo2003bayesian}, which discusses a Bayesian approach to 
sequential experimental design of observations---see also~\citep{wynn1984jack,kiefer1975optimal, mitchell2000algorithm,smucker2018optimal} for a statistical background.
However, this approach also relies on regularity conditions for the underlying data model, and so it does not account for
the identification issues  discussed in Section~\ref{sec:intro}.

One way to address these issues is to use the identification method developed in this paper.  The idea is simply to synthesize data under alternative designs, and then pick the design 
that yields the most precise confidence set, $\hat\Theta_{1-\alpha}$. 
Conceptually, the main challenge with this approach is how to define an appropriate space for these alternative designs in a way that is both practical and computationally efficient.
Given a choice for the design space, we can then simulate synthetic data 
and choose from those designs that yield  narrow confidence sets.

In practice, we usually face a scenario where we suspect the peak $\hat\theta_n$ in the periodogram to be the true signal, but there is also another ``nuisance peak", say $\theta^{\mathrm{nuis.}}$. 
Typically, this nuisance peak appears near the 1-day signal as a by-product of a daily observation design.
For example, in the ESPRESSO data discussed in Section~\ref{sec:exoplanet} and Figure~\ref{fig6}, the problem is that 
even though there is very strong evidence for a signal at $\hat\theta_n=11.17$ days, 
there is also a signal at $\theta^{\mathrm{nuis.}} =0.9164$ that cannot be rejected at the 1\% level~($p$-value = 0.017). The observed peak could statistically be an artifact of the daily observation design. 
The key question is thus: 
\begin{quote}
``What kind of observation designs should we use to accurately 
identify the signal at $\hat\theta_n$, {\em assuming} this signal to be the true one?"
\end{quote}

To address this problem, we propose the following high-level procedure.

\begin{enumerate}
\item Select a space $\mathcal{D}$ for observation designs and a tolerance threshold $0 < \epsilon < |\hat\theta_n - \theta^{\mathrm{nuis.}}|$.
\item {\bf For} each design $D$ in the design space, {\bf do:}
\begin{enumerate}[(a)]
\item Assume $\theta^\ast= \hat\theta_n$ as the true signal, and then calculate $\hat\psi_{\theta^\ast}$ as defined in~\eqref{eq:real_per}.
\item {\bf For} $r=1, \ldots, R$, {\bf do:}
\begin{enumerate}[(i)]
\item Generate synthetic data $(\tilde T^{\nnn'}, \tilde Y^{\nnn'})$ according to the design and $(\theta^\ast, \hat\psi_{\theta^\ast})$ from step 2(a), where  we may 
introduce additional (synthetic) observations, such that $n' \ge n$.
\item  Calculate $\hat\Theta_{1-\alpha}$ based on synthetic data $(\tilde T^{\nnn'}, \tilde Y^{\nnn'})$.
\item Calculate how many periods in $\hat\Theta_{1-\alpha}$ are $\epsilon$-far from $\theta^\ast$:
$$
I_{D}^r = \mid \{ \theta\in\hat\Theta_{1-\alpha} : |\theta - \theta^\ast| > \epsilon \} \mid.
$$
\end{enumerate}
\end{enumerate}
\item Select the ``simplest" design possible that also yields sharp identification:
\begin{align}\label{eq:best_design}
\hat{D} = \arg\min_{D\in\mathcal{D}} 
\left\{ |D | : \sum_{r=1}^R I_{D}^r = 0\right\},
\end{align}
where $|D| \in\mathbb{R}$ denotes the complexity of the design~(explained below).
\end{enumerate}

In words, the design $\hat D$ in~\eqref{eq:best_design} is the simplest possible observation design 
for which the confidence sets include only values in the $\epsilon$-neighborhood of $\theta^\ast$.
In this case, we can accurately identify the underlying true signal. To implement the optimal design procedure outlined above, we only have to define the design space $\mathcal{D}$ along with an appropriate notion of design complexity, and also discuss choices for threshold $\epsilon$. Here, for simplicity, we consider two choices for the design space:
\begin{enumerate}[(A)]
\item {\bf Fixed observations, $n$}.  
Each design in this space is parameterized by $\delta>0$ such that 
$\tilde T^{\nnn} = T^{\nnn} + \delta U^{(n)}$, where $U^{(n)}$ is a vector of $n$ independent $\mathrm{Unif}([-1, 1])$, 
i.e., random variables uniformly distributed on $[-1,1]$.
Thus, the alternative observation time $\tilde t_i$ is uniformly distributed $\pm \delta$ days around 
the actual observation time $t_i$. 
Given $\tilde T^{\nnn}$, we may generate synthetic radial velocity data, $\tilde Y^{\nnn}$, through
$y(\tilde t_i) = y^{\mathrm{p}}(\tilde t_i; \theta^\ast, \hat\psi_{\theta^\ast}) + \varepsilon_i$, 
where $y^{\mathrm{p}}$ is the periodic component defined in~\eqref{A2}, and $\varepsilon_i \sim N(0, \sigma_i^2)$, i.i.d.
The complexity of a design $D$ with parameter $\delta$ can simply be defined as
$| D| = \delta$.

\item {\bf Additional observations, $n' - n$.}
Each design in this space is parameterized by $n'>0$ such that 
$\tilde T^{\nnn'} = (T^{\nnn}, \tilde T^{(n'-n)})$, where we sample $(n'-n)$ additional observation times in a hypothetical span of 3 months as follows:
$\tilde t_j = \max \{T^{\nnn}\} + \mathrm{Unif}(\{1,\ldots, 90\}) + 0.01 \cdot \mathrm{Unif}([-1,1])$, for $j=1, \ldots, (n'-n)$.
Given $\tilde T^{\nnn'}$, we can then generate synthetic radial velocity data $\tilde Y^{\nnn'}$ exactly as in Design (A).
The complexity of a design $ D$ can simply be defined as
$| D| = n'-n$.
\end{enumerate}

Intuitively, Design (A) adds some level of randomization (controlled by $\delta$) on top of the actual observation times $T^{\nnn}$. Roughly speaking, this determines the amount of randomization necessary for accurate identification in future observations. When there are identification issues, a large $\delta$ will generally be necessary for identification, and, conversely, when there are no identification issues a small $\delta$ (or zero) will be needed.
Design (B) adds new observation time points at a pre-specified level of noise. 
This level is set to $\pm 0.01$ on the day-scale which amounts to $2\times 0.01\times 24 \times 60 \approx$ 30 mins of randomized observation window around $n'-n$ additional daily observations.
Similar to Design (A), when there are identification issues, a large number of additional observations, $n'-n$, will be necessary for identification; and, conversely, when there are no identification issues, a small number of 
additional observations will be needed.

Regarding the hyper-parameters, we note that $R$ controls the number of synthetic datasets 
generated per design, and thus large $R$ leads to more accurate calculation of the optimal design in~\eqref{eq:best_design}. Practically, since we want to estimate small, near-zero proportion values, $R$ could be small, e.g., $R=1,000$ or  even $R=100$.
Moreover, a reasonable value for $\epsilon$ is in the range $[0.1, 0.5]$. Increasing that threshold to larger values could conflate real small-periodicity signals with the common 1-day period signal that is associated with the daily observation designs.
 
To illustrate these ideas, we apply the above optimal design procedure to the exoplanet detection problems of Section~\ref{sec:exoplanet}. The results are shown in Table~\ref{tab:design}. We see 
 that for 51 Pegasi b and Gliese 436 b, there is no need for improving the observation times. 
 Indeed, we need neither additional randomness in observation times ($\delta=0$), 
 nor any additional observations ($n'-n=0)$. 
 In contrast, for the candidate exoplanet around $\alpha$ Centauri B~\citep{dumusque2012earth}
 we need a substantial improvement to the observation design. 
Specifically, in order to identify precisely the observed signal at $\hat\theta_n=3.23$ at the 1\% level~(assuming it's true),  we need a large additional variation of $\pm 0.18$ days around the actual observation times 
(i.e., $\pm 4.32$ hrs./observation). Alternatively, we need 137 additional observations with a random variation of $\pm 15$ mins./observation, as specified by Design (B), almost equal to the size of the actual dataset.
Finally, the situation is much better for Proxima Centauri~\citep{anglada2016terrestrial}. 
Specifically, in order to identify the observed signal at $\hat\theta_n=11.17$ at the 1\% level ---and so eliminate the only nuisance signal at $\theta^{\mathrm{nuis.}}=0.9164$---  we need an additional variation of $\pm 0.06$ days~(i.e., $\pm 1.44$ hrs./observation) around the actual observation times. Alternatively, we only need 17 additional observations with a random variation of 
$\pm 15$ mins./observation, as specified by Design (B). 

In summary,  the results in Table~\ref{tab:design} illustrate the benefits of the design procedure outlined above in guiding the selection of observation times for future exoplanet detections. In practice, we can envision a two-stage approach where at the first stage the scientist analyzes the stellar signal through the periodogram and identifies the most plausible peak along with other nuisance peaks. In the second stage, the scientist performs an analysis as described here in order to quantify how much variation or how many more samples are needed for identification of the plausible signal at a desired significance level, and proceeds to make these observations.

Of course, one limitation of our discussion here is that we have only considered two simple classes of observations designs. In practice, we need to also take into account the particular logistics 
of Earth-based observations to determine more realistic observation designs. 
The optimal design procedure in~\eqref{eq:best_design} could then be used in conjunction with such more realistic designs. Overall, this is an intriguing statistical problem, which we suggest for future work.


\renewcommand{\arraystretch}{1.3}
\begin{table}[t!]
\small
\begin{tabular}{lccc}
& \multicolumn{2}{c}{Design (A)} & Design (B) \\
\hline
\multicolumn{1}{l|}{Candidate Exoplanet (or Star)}                                                  & \begin{tabular}[c]{@{}c@{}} \small{randomness needed} \\ \small{for identification} (best $\delta$)
\end{tabular} & \multicolumn{1}{c|}{ \small{$\pm$ hrs.}} & \begin{tabular}[c]{@{}c@{}} \small{\#additional obs.~needed} \\
\small{for identification~(best $n'-n$)}
\end{tabular} \\ \hline
\multicolumn{1}{l|}{51 Pegasi b~\footnotesize{\citep{mayor1995jupiter}  }   }      & 0                                                                         & \multicolumn{1}{c|}{0}                   & 0                                                                                  \\
\multicolumn{1}{l|}{Gliese 436 b~\footnotesize{\citep{butler2004neptune}}}           & 0                                                                         & \multicolumn{1}{c|}{0}                   & 0                                                                                  \\
\multicolumn{1}{l|}{$\alpha$ Centauri B~\footnotesize{\citep{dumusque2012earth}}}     & 0.18                                                                      & \multicolumn{1}{c|}{4.32}                & 137                                                                            \\
\multicolumn{1}{l|}{Proxima Centauri\footnotesize{~\citep{anglada2016terrestrial}}} & 0.06                                                                      & \multicolumn{1}{c|}{1.44}                & 17                                                                             
\end{tabular}
\caption{
\small Observation designs (A) and (B) to achieve identification 
in the exoplanet applications of Section~\ref{sec:exoplanet}.
Design (A) introduces additional randomness in the observation times, while design (B) introduces additional observations. 
We see that 51Pegb and GJ436b require no improvement in the observation times as the signal is well-identified.
To identify the candidate exoplanet around $\alpha$ Centauri B, however, 
  we need an additional variation of $\pm 0.18$ days around the actual observation times 
(i.e., $\pm 4.32$ hrs./observation). 
Alternatively, we need 137 additional observations with a random variation of $\pm 15$ mins./observation.
For the candidate exoplanet around Proxima Centauri~\citep{anglada2016terrestrial} we need an additional variation of $\pm 0.06$ days~(i.e., $\pm 1.44$ hrs./observation) on the actual observation times. Alternatively, we only need 17 additional observations with a random variation of 
$\pm 15$ mins./observation
}
\label{tab:design}
\end{table}



\subsection{Non-parametric approach}\label{sec:np}
The method developed in Section~\ref{sec:method} is parametric and employs the
harmonic model for the periodic component~\eqref{A2}. While this can be extended to include other periodic models (e.g., Keplerian model) it is interesting to consider non-parametric approaches as well. While nonparametric approaches have been considered in the literature before~\citep{feigelson2012modern, ivezic2014statistics, hall2000nonparametric, hall2006using, hall2008nonparametric}, their asymptotics require strong assumptions on the observation design (typically i.i.d observation times), which is unrealistic. As such, nonparametric methods face the  the same identification problems tied to irregular observation times as parametric methods.
For example, the identification of a 26-million year biological extinction cycle 
in Example~\ref{example2} used a nonparametric randomization test. However, as discused in Section~\ref{sec:extinction}, the period cannot be identified from this dataset.

Here, we discuss a nonparametric variant of our set identification method. As before, we assume 
that the signal is decomposed as $y_i = y^{\mathrm{p}}(t_i) + \varepsilon(t_i)$, but make no parametric assumption about $y^{\mathrm{p}}$. Before we define and test our main hypothesis, we need some notation.
For any two time points $t, t'$ and period $\theta\in\Theta$, let $t \equiv t'(\mod\theta)$ denote that 
$t-t' = k\theta$ for some integer $k$. Also, let $\mathsf{S}_n$ denote the symmetric group of permutations of $n$-element vectors.
For some fixed $\theta$, define:
$$
\Pi(T^{\nnn}; \theta) = \{\pi \in\mathsf{S}_n : t_{\pi(i)} \equiv t_i(\mod\theta),~i=1, \ldots, n\}.
$$
In words, $\Pi(T^{\nnn}; \theta)$ is the set of permutations of $\{1, \ldots, n\}$ such that  
any time $t_i$ is mapped only to an observation time that is equivalent to $t_i$ modulo $\theta$.
It is straightforward to see that $\Pi(T^{\nnn}; \theta)$ is a subgroup of $\mathsf{S}_n$ since the $\mod$operator defines an equivalence relationship. Let us assume that 
for any $T^{\nnn}$ and period $\theta\in\Theta$, the errors satisfy the invariance Assumption~\eqref{A1} with $\mathcal{G}^{\nnn} = \Pi(T^{\nnn}; \theta)$. This assumption is generally mild. It requires that the data are exchangeable modulo $\theta$, and actually allows certain correlated structures for the errors $\varepsilon(t)$.

We wish to test the following nonparametric null hypothesis of periodicity $\theta_0$:
\begin{align}
H_0^{\mathrm{np}}:~y^{\mathrm{p}}(t') = y^{\mathrm{p}}(t),~\text{for all}~t',t~\text{such that}~t'\equiv t(\mod\theta_0).\nonumber
\end{align}
To test $H_0^{\mathrm{np}}$ we can adapt Procedure 1 as follows.
\MakeProcedure{2}{
\begin{enumerate}
\item {\bf For} all $r = 1, \ldots, R$ {\bf do:}
\begin{enumerate}[(i)]
\item Sample $\pi\sim\mathrm{Unif}\big(\Pi(T^{\nnn}; \theta_0)\big)$.
\item Generate synthetic outcome data $Y^{\nnn,(r)} = \pi \cdot Y^{\nnn}$ obtained by permuting the data $Y^{\nnn}$ according to $\pi$ while the observation times, $T^{\nnn}$, are fixed.
\end{enumerate}
\item Using the samples  from 2(ii), calculate the $p$-value, say $\mathrm{pval}(\theta_0)$, as in~\eqref{eq:pval_theta}, and 
reject if the $p$-value is less than $\alpha$.
\end{enumerate}
}

In words, to test $H_0^{\mathrm{np}}$, Procedure~2, first, picks a permutation of the observation times in $T^{\nnn}$ such that 
the permuted times are equivalent to the corresponding original periods modulo $\theta_0$.
Then, it uses this permutation to generate a new outcome vector $Y^{\nnn,(r)}$. In the last step, the procedure recalculates the test statistic using data $(T^{\nnn}, Y^{\nnn,(r)})$, which can then be used 
in the calculation of the $p$-value.


\begin{theorem}\label{thm:np}
\ThmFour
\end{theorem}

According to Theorem~\ref{thm:np}, the test of Procedure~2 is exact without requiring a parametric model for the periodic component, $y^{\mathrm{p}}$. This is possible because under $H_0^{\mathrm{np}}$ and the invariance assumption on the errors, any two 
observations $y(t), y(t')$ are exchangeable when $t, t'$ are equivalent modulo $\theta_0$.

Despite the strength of this result, one practical difficulty with Procedure 2 is that it requires 
us to find several pairs of distinct time points in $T^{\nnn}$ that are periodically equivalent modulo $\theta$, for all $\theta\in\Theta$. This may be unrealistic in practical applications, including exoplanet detection, because the sample sizes can be small. Technically, the problem is that $\Pi(T^{\nnn}; \theta)$ may be a singleton group (including only the identity permutation) for many values of $\theta$, so that Procedure 2 leads to an exact by ``empty" test, and thus to confidence sets that are too conservative. 
To address this problem, it may be tempting to expand the definition of $\Pi(T^{\nnn}; \theta)$
to include permutations between $t,t'$ that are $\epsilon$-away from being periodically equivalent modulo $\theta$. Under this approach, however, $\Pi(T^{\nnn}; \theta)$ is not a group anymore.
An alternative approach would be to use the nonparametric estimators of $\theta^\ast$ 
developed by~\citet{hall2000nonparametric, hall2006using, hall2008nonparametric} together with a variation of Procedure~1 or Procedure~2. Both these procedures do not require regularity conditions on the observation times but only a consistent estimator for the periodic component, $y^{\mathrm{p}}$. We leave these directions for future work.

\section{Conclusion}
We developed a method of set identification for hidden periodicity in unequally spaced time series. 
Our approach is more appropriate than standard methods of statistical inference because common estimators, such as the periodogram peak, are not well-behaved and may even be inconsistent without 
strong (and practically unrealistic) assumptions on the observation design.
We illustrated empirically our method in examples from  exoplanet detection using radial velocity data.
In the confirmed exoplanets of our sample, our method can identify the hidden periodicity precisely. 
We also find that the statistical evidence in a recent candidate around  $\alpha$ Centauri B  is particularly weak, as several nuisance signals cannot be rejected even at the 5\% level. In another popular discovery around Proxima Centauri, the evidence appears to be robust but a nuisance signal due to the daily observational cycle can still not be rejected at the 1\% level.
Our method suggests ways to improve the observation designs, either by introducing 
some randomness in the observation times, or by just making more observations.
These improved designs could resolve the aforementioned identification issues, and also help in 
efficient scheduling of measurements for future discoveries.

\section{Acknowledgements}
We would like to thank Yang Chen, Eric Feigelson, Vinay Kashyap, Xiao-Li Meng, Aneta Siemiginowska, Katja Smetanina, Michael Vogt, and the participants at the Statistics Department Seminar at UC Berkeley and the ``Topics in Astrostatistics" Seminar at Harvard University for useful suggestions and feedback.

\small

\section*{A. Proof of Theorem~\ref{thm1}}

\begin{align}
\label{A1}
\gg\cdot \eps^{\nnn} & \eqd \eps^{\nnn}~\mid T^{\nnn}~\quad(\gg\in\mathcal{G}^{\nnn}).\tag{A1} \\
\label{A2}
y^{\mathrm{p}}(t; \theta^\ast, \psi^\ast) & = \psi^\ast_{1} + \psi^\ast_{2} \cos(2\pi t/\theta^\ast) + \psi^\ast_{3} 
\sin(2\pi t/\theta^\ast).
\tag{A2}
\end{align}

The definition of the $p$-value for $H_0^{\mathrm{full}}$ is as follows:
\begin{align}\label{eq:pval_full}
\mathrm{pval}(\theta_0, \psi_0) = E\{s_n(Y^{\nnn,(i)}, T^{\nnn}) \ge s_\obs \mid T^{\nnn}\},
\end{align}

\begin{theorem}\label{thm1}
\ThmOne
\end{theorem}
\begin{proof}
Recall that $Y^{\nnn,(i)} = Y^{\nnn,\mathrm{p}} + G^{(i)} \eps^{\nnn}$ where $G^{(i)}\sim \mathrm{Unif}(\mathcal{G}^{\nnn})$ and $\eps^{\nnn} = Y^{\nnn} - Y^{\nnn,\mathrm{p}}$.
Under $H_0^{\mathrm{full}}$ and Assumptions~\eqref{A1}-\eqref{A2} we obtain
$Y^{\nnn,(i)} \eqd Y^{\nnn} \mid T^{\nnn}$. It follows that
$$
s_n(Y^{\nnn,(i)}, T^{\nnn}) \eqd s_\obs \mid T^{\nnn}
$$
and so the $p$-value in~\eqref{eq:pval_full} is finite-sample valid, that is,
$\mathrm{E}\big[ \mathbbm{1}\{\mathrm{pval}(\theta_0, \psi_0) \le \alpha\} \mid H_0^{\mathrm{full}}, T^{\nnn}\big] = \alpha$.
It follows that the $p$-value is also valid {\em unconditionally}:
\begin{align}
\mathrm{E}\big[ \mathbbm{1}\{\mathrm{pval}(\theta_0, \psi_0) \le \alpha\} \mid H_0^{\mathrm{full}}\big]
= 
\mathrm{E}\big[~\mathrm{E}\big[ \mathbbm{1}\{\mathrm{pval}(\theta_0, \psi_0) \le \alpha\} \mid H_0^{\mathrm{full}}, T^{\nnn}\big] \mid T^{\nnn}\big]
= \alpha.\nonumber
\end{align}
\end{proof}


\section*{B. Proof of Theorem~\ref{thm2}}
\begin{align}
\label{eq:H02}
& H_0: \theta^\ast  =\theta_0. \\
\label{Theta1}
 & \Theta_{1-\alpha} = \{\theta \in \Theta :   \max_{\psi\in\Psi}  \mathrm{pval}(\theta, \psi) > \alpha \}.
\end{align}

\begin{theorem}\label{thm2}
\ThmTwo
\end{theorem}
\begin{proof}
Suppose that $H_0$ in~\eqref{eq:H02} holds., Then,
\begin{align}
\mathrm{pr}(\theta_0\notin\Theta_{1-\alpha} \mid H_0)
 & = \mathrm{E}[\mathbbm{1}\{\max_\psi \mathrm{pval}(\theta^\ast, \psi) \le \alpha\} ]~\nonumber\\
 & \le \mathrm{E}[\mathbbm{1}\{\mathrm{pval}(\theta^\ast, \psi^\ast) \le \alpha\} ]~\nonumber\\
 %
 & = \alpha.
 \end{align}
The first equality follows from the definition of $H_0$ in~\eqref{eq:H02}; the second line follows from the indicator identity
  $\mathbbm{1}(a \le b) \le \mathbbm{1}(c \le b)$ when $a\ge c$; 
  and the third line follows from Theorem~\ref{thm1}. Since $\theta^\ast\in\Theta$ w.p.~1 it follows that $\mathrm{pr}(\theta^\ast\in\Theta_{1-\alpha}) \ge 1-\alpha$.
\end{proof}

\newcommand{\mathrm{pval}(\theta^\ast, \psi^\ast)}{\mathrm{pval}(\theta^\ast, \psi^\ast)}
\newcommand{\numeq}[1]{\quad\quad~\mathrm{(#1)}}
\section*{Proof of Theorem~\ref{thm3}}
\begin{theorem}\label{thm3}
\ThmThree
\end{theorem}
\begin{proof}
Following a similar reasoning as above:
\begin{align}
\mathrm{pr}(\theta_0\notin\hat\Theta_{1-\alpha} \mid H_0)
 & = \mathrm{pr}\{ \mathrm{pval}(\theta^\ast, \hat\psi_{\theta^\ast}) \le \alpha \mid H_0 \} \numeq{i}
 \nonumber\\
 & = \mathrm{pr}\{\mathrm{pval}(\theta^\ast, \psi^\ast) +  \underbrace{\mathrm{pval}(\theta^\ast, \hat\psi_{\theta^\ast}) - \mathrm{pval}(\theta^\ast, \psi^\ast)}_{\Delta}  \le \alpha \mid H_0  \} \nonumber\\

& \le  \mathrm{pr}\{ \mathrm{pval}(\theta^\ast, \psi^\ast)  \le  \alpha+  |\Delta| \mid H_0 \} \numeq{ii}
\nonumber\\
& \le \mathrm{pr}(|\Delta | > \epsilon \mid H_0 ) + \mathrm{pr}\{ \mathrm{pval}(\theta^\ast, \psi^\ast) \le \alpha + \epsilon \mid H_0 \}
\numeq{iii}
\nonumber\\
& = \alpha + \epsilon + P(|\Delta | > \epsilon \mid H_0 ). 
\numeq{iv}
\nonumber\\
& \le \alpha + \min_\epsilon\{\epsilon + P(|\Delta | > \epsilon\mid H_0 )\}
\numeq{v}
\nonumber\\
& = \alpha + o_P(1).
\numeq{vi}
\end{align}
In these derivations, line (i) follows from definition of the confidence set and $H_0$; 
line (ii) follows from the indicator identity $\mathbbm{1}(A + B \le a) \le \mathbbm{1}(A\le a + |B|)$ for any scalar random variables $A, B$ and $a\in\mathbb{R}$; 
line (iii) follows from the 
indicator identity $\mathbbm{1}(A \le a + B) \le \mathbbm{1}(B > \epsilon) + \mathbbm{1}(A \le a + \epsilon)$ for any $\epsilon <0$, 
and $A, B, a$ as above; line (iv) follows from Theorem~\ref{thm1}; 
line (v) gives the tightest bound exploiting the fact that  $\epsilon>0$ is arbitrary.
Finally, line (vi) follows from consistency of $\hat\psi_{\theta^\ast}$.
\end{proof}


\subsection*{Proof of Theorem~\ref{thm:np}}
\begin{theorem}\label{thm:np}
\ThmFour
\end{theorem}
\begin{proof}
The proof follows the same steps as Theorem~\ref{thm1} 
as we only have to derive that $s_n \eqd s_\obs \mid T_n$.
The only difference is that this equivalence is implied 
immediately from the assumed invariance based on $\Pi(T^{\nnn}; \theta)$ rather than relying 
on the correct periodic component $y^{\mathrm{p}}$ as in Theorem~\ref{thm1}.
\end{proof}



