Periodic Variable Stars Modulated by Time-varying Parameters

RR Lyrae stars are important astrophysical tools for the measurement of distances and studies of the astrophysical properties of old stellar populations. They are moderately bright, evolved low-mass stars, currently in the core helium-burning phase, also known as the horizontal branch. Their periods are typically in the range of between about 0.2 and 1.0 day, which together with their characteristic light-curve shapes, allow them to be relatively easily identified in time series photometric surveys. An overview of their properties can be found in the monographs by Smith (1995) and Catelan & Smith (2015).

In spite of their astrophysical importance, RR Lyrae stars are still not fully understood. Indeed, one of the longest-standing problems in stellar astrophysics is also one that specifically affects RR Lyrae stars: the so-called "Blazhko effect" (Blažko 1907). It is a long-term modulation of an RR Lyrae's light curve, over timescales ranging from a few to hundreds of days (for recent reviews, see Catelan & Smith 2015; Gillet et al. 2019). The Blazhko effect is particularly common among fundamental-mode (ab-type) pulsators (e.g., Plachy et al. 2019), but is also present, to a lesser extent, in first-overtone (c-type) RR Lyrae stars (e.g., Netzel et al. 2018).

Over the decades since it was first described, the Blazhko effect has persistently defied theoretical explanations as to its cause (e.g., Gillet et al. 2019). Gradual strengthening and weakening of turbulent convection in the stellar envelope (Stothers 2006), a 9:2 resonance between the fundamental and ninth-overtone radial modes (Buchler & Kolláth 2011), and interaction between fundamental and first-overtone modes in the "either-or" region of the instability strip (Gillet 2013) are the most recent candidates, but no consensus has yet been reached as to the root cause of the Blazhko effect, due in large part to the difficulties involved in the nonlinear hydrodynamical modeling of the phenomenon.

In this article, we introduce a model for time series observations of variable stars having a smoothly time-varying trend and amplitudes. More precisely, we develop a semi-parametric method for unequally spaced time series measuring the brightness of a modulated variable star. Our approach is flexible because it avoids assumptions about the functional form of the trend and amplitudes. The estimation of our time-varying curves translates into the estimation of time-invariant parameters that can be performed by ordinary least squares, with the following two advantages: modeling and forecasting can be implemented in a parametric fashion, and we are able to cope with missing observations. We also study the spectral density of the residuals obtained from the fit of our novel model.

In order to detect serial correlation in the residuals, in this paper we derive the definition of the spectral density for unequally spaced time series. There are many reasons why astronomical time series are not sampled equidistantly, and the gaps can be either regular or random. From the Earth, stars cannot be observed during the day, which introduces regular gaps in the time series. Also, for about half a year, most objects become unobservable, as they are up on the sky at the same time as the Sun, which introduces yearly gaps. There could be clouds or high winds, forcing the closure of telescopes, producing random gaps. There could be high-priority alerts overriding the observations, or the telescope could be available only on certain nights.

In some cases, observations are unevenly spaced due to missing values. Astronomical data sets often contain missing values, and this limitation is sometimes due to incomplete observations or varying survey depths. Even if telescopes are recording and storing information systematically, that is, at a regular cadence, there are a few things that can alter the regular sampling. For example, an astronomer might decide to increase the exposure time if there are clouds obscuring the target, to try to increase the signal-to-noise ratio. Conversely, if the observing conditions are excellent, the astronomer might decide to decrease the exposure times (and hence the cadence) to avoid saturating the detector, for example. Missing values are usually handled via imputation, that is, the gap generated by the missing value is "filled in" by an estimated value. If observations are missing because of the survey, imputation can be performed using statistical models. For example, Feigelson et al. (2018) apply ARIMA models to fill in missing data in astronomical time series.

However, in astrostatistics, missing value problems are sometimes inherently brought about by the manner in which physical processes are recorded. In particular, telescopes are not located in the center of the solar system. Since the speed of light is finite, this results in a time delay between the arrival times of signals at our position and at the center of the solar system. This is typically corrected for by referring the times of observations to either Heliocentric Julian Dates (HJD) or Barycentric Julian Dates (BJD), which refer to the center of the Sun or the entire solar system, respectively. Thus, even if telescopes record data strictly evenly according to the local time at the observatory (e.g., one observation performed every night at local midnight), this correction will slowly change between observations, modifying what was initially a regular grid to an irregular one. Also, this correction is different for every source on the sky, even though sources close to each other may have very similar corrections. Therefore, for some astronomical data sets where missing values may arise from the manner in which observations of a physical process are collected, or even the nature of the physical process itself (e.g., sudden, extreme dimming events that may occasionally render an object impossible to detect for a certain amount of time), the imputation method may not be applicable (see Chattopadhyay 2017). Our novel approach, which involves the classical periodogram, has the advantage of not relying on any imputation method.

We shall divide the present study into six main sections. In Section 2 we introduce our novel model and clarify analogies and differences as compared with previous approaches. In Section 3, the present status of important ingredients of amplitude and frequency modulations is critically discussed. In Section 4 we present the method we adopt to estimate the time-varying parameters. In Section 5 we present a new method to estimate the spectral density of unequally spaced times series, which is needed for the analysis of the residuals. Section 6 provides simulation results, whereas Section 7 illustrates the advantages of using our novel method by means of an application to an RR Lyrae variable star. Finally, our main conclusions are summarized in Section 8.

Through the paper we use bold uppercase letters to denote matrices, and bold slanted to denote vectors. We denote by I_m the identity matrix of size m, by 0_n a column-vector of zeros of length n, by $\mathrm{tr}\{{\bf{A}}\}$ the trace of A, by A^⊤ the transpose of A, by ∥A∥ the Frobenius norm $\parallel {\bf{A}}\parallel ={[\mathrm{tr}\{{{\bf{A}}}^{\top }{\bf{A}}\}]}^{1/2}$, and by A⁻¹ the inverse of the square matrix A, that is, the square matrix A⁻¹ such that A⁻¹ A = AA⁻¹ = I.

Light curves of variable stars are typically fitted using harmonic models with a linear (or constant) trend and time-invariant amplitudes (see Equations (1) and (5) in Richards et al. 2011). This type of model would be inappropriate when the underlying trend and amplitudes change over time in a more complex way. Eilers et al. (2008) proposed a model with one harmonic component (K = 1) where the trend and amplitudes vary smoothly over time. In this paper, we extend the model by Eilers et al. (2008) to the case of K ≥ 1 harmonic components, where the amplitudes associated with each harmonic component vary smoothly over time. We estimate our model by means of P-splines (Eilers & Marx 1996), which are a combination of B-splines and penalties. The estimation of the time-varying curves translates into the estimation of time-invariant parameters that can be performed by the least-squares method, with the following three advantages: it is computationally fast, forecasting can be implemented in a parametric fashion, and we can cope with missing observations.

Compared to local smoothers (such as kernel smoothers), the main advantage of regression spline in the context of time series is that the unknown parameters are time-invariant and thus they can be estimated globally rather than locally. As a consequence, forecasting only requires good estimates of the global unknown parameters. We can think of regression splines with B-splines as a semi-parametric model in the sense that it contains parametric as well as nonparametric components. The parametric component is given by a finite number of parameters, whereas the nonparametric component is given by the basis functions. Another advantage of parametric and semi-parametric models over nonparametric models is the computing speed, as many nonparametric models are computationally intensive. Finally, the use of B-splines in regression allows us to rewrite the estimation problem as a least-squares fit, avoiding the use of numerical methods—such as Newton−Raphson—which can be time consuming.

Let $\{{Y}_{i}\equiv {Y}_{{t}_{i}},i=1,\,\ldots ,\,N\}$ be a set of observations occurring at certain discrete times t₁, ..., t_N . In the case of equally spaced observations, t_i = t₀ + iΔ where i is an integer, and Δ > 0 is the constant data spacing. Then ∣t_i − t_k ∣ = Δ∣i − k∣, and typically Δ = 1. Astronomical light curves are often observed unequally in time, that is, the data spacing of observation times is not constant.

We decompose the observed light curve into the sum of a deterministic signal μ and a random noise z. The deterministic part $\mu \left(t\right)$ consists of a trend m(t) and a modulated periodic signal. The modulated periodic signal is a linear combination of K cosines and sines, with weights given by the modulating functions g(t):

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{Y}_{i} & = & \mu ({t}_{i})+{z}_{i},\\ i & = & 1,\,\ldots ,\,N,\,\{{z}_{i}\}\sim {\rm{W}}{\rm{N}}(0,{\sigma }_{z}^{2}),\\ \mu ({t}_{i}) & = & m({t}_{i})+\displaystyle \sum _{k=1}^{K}\{{g}_{1,k}({t}_{i})\cos ({w}_{k}{t}_{i})\\ & & +\,{g}_{2,k}({t}_{i})\sin ({w}_{k}{t}_{i})\},\end{array}\end{array}\end{eqnarray} \tag{ 1 }$$

or in matrix notation Y = μ + z , where ${\boldsymbol{Y}}={\left({Y}_{1},\,\ldots ,\,{Y}_{N}\right)}^{\top }$ is the vector of observations at time ${\boldsymbol{t}}={\left({t}_{1},\,\ldots ,\,{t}_{N}\right)}^{\top }$, ${\boldsymbol{\mu }}\,={[\mu ({t}_{1}),\,\ldots ,\,\mu ({t}_{N})]}^{\top }$ is the expectation of Y , m(t_i ) is the smooth time-varying trend at time t_i , the g(t_i )'s are smooth time-varying amplitudes of the cosine and sine waves at time t_i , respectively, w_k = 2π f_k is the angular frequency, and f_k is the ordinary frequency. Since the errors are zero mean, the expectation of the observed brightness at time t_i is equal to the deterministic part of the signal at time t_i , that is, ${\rm{E}}\left[{Y}_{i}\right]=\mu ({t}_{i})$.

We refer to m( · ) as the "trend," that is, the (typically) aperiodic change in the mean of the light curve. By contrast, we call "amplitudes" the functions g( · )'s that weigh the periodic variation (of this average brightness) of cosine and sine waves. In Appendix A we summarize the standard mathematical definitions of amplitude modulation and frequency modulation. Both trend and amplitudes are typically restricted to be sinusoidal, whereas in this paper our trend m(t) and our amplitude functions g(t)'s are general smooth functions and not necessarily sinusoidal. In Section 3 we clarify the mathematical connection between the standard modulation models and our novel modulation model in Equation (1).

The error vector ${\boldsymbol{z}}={\left({z}_{1},\,\ldots ,\,{z}_{N}\right)}^{\top }$ is a white noise (WN) process with mean zero and variance ${\sigma }_{z}^{2}$. That is, each error z_i , i = 1, ..., N, follows a zero-mean WN process with variance ${\sigma }_{z}^{2}$:

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{E}}\left[{z}_{i}\right] & = & 0\quad {\rm{and}}\quad {\rm{E}}\left[{z}_{i}\,{z}_{j}\right]\\ & = & {\delta }_{\{i=j\}}{\sigma }_{z}^{2},\quad {\rm{for}}\,{\rm{all}}\ i,j=1\ldots ,N,\end{array}\end{eqnarray*}$$

where δ_{i=j} = 1 if i = j and zero otherwise.

Our model in Equation (1) is defined in discrete time, and it focuses on the time domain. Kelly et al. (2014) adopt the continuous-time autoregressive moving average (CARMA) models to estimate the variability features of a light curve in the frequency domain. More specifically, Kelly et al. (2014) use the power spectral density (PSD) of CARMA models to account for irregular sampling and measurement errors. A stationary CARMA(p, q) process has the PSD

$$\begin{eqnarray*}&&P(f)={\sigma }^{2}{\left|\sum _{j=0}^{q}{\beta }_{j}{\left(i2\pi f\right)}^{j}\right|}^{2}/{\left|\sum _{k=0}^{p}{\alpha }_{k}{\left(i2\pi f\right)}^{k}\right|}^{2}.\end{eqnarray*}$$

To illustrate the importance of fitting models with time-varying parameters, Kelly et al. (2014) simulated a light curve that switches from one CARMA process to another. More precisely, they constructed a nonstationary light curve by generating two CARMA processes of the same order (p = 5, q = 3), but with different parameters (see Kelly et al. 2014, Section 4.3):

$$\begin{eqnarray*}{\boldsymbol{\theta }}(t)=\left\{\begin{array}{ll}{{\boldsymbol{\theta }}}_{1} & {t}_{1}\leqslant t\lt {t}_{0}\\ {{\boldsymbol{\theta }}}_{2} & {t}_{0}\leqslant t\leqslant {t}_{N},\end{array}\right.\end{eqnarray*}$$

where ${\boldsymbol{\theta }}(t)={[{\alpha }_{1}(t),\,\ldots ,\,{\alpha }_{p}(t),{\beta }_{1}(t),\,\ldots ,\,{\beta }_{q}(t),{\sigma }^{2}(t)]}^{\top }$. The vector ${\boldsymbol{\theta }}\left(t\right)$ is a step-wise function that is constant before and after t₀. The approach based on piecewise-constant parameters is receiving growing interest in various areas of astrophysics. Wong et al. (2015) adopt a Poisson model for the photon counts. They define λ(t_j , w_i ) as the expected count per unit time and per unit wavelength averaged over the bin centered at (t_j , w_i ), and detect change points π such that {λ(t_j , w_i )∣t_j ≤ π} ≠ {λ(t_j , w_i )∣t_j > π}. Wong et al. (2015) estimate the number of change points and their values. Xu et al. (2021) develop a method for modeling a time series of images, and assume that the arrival times of the photons follow a Poisson process. They assume that all image stacks between any two adjacent change points (in the time domain) share the same unknown piecewise-constant function. Xu et al. (2021) estimate the number and the locations of all of the change points (in the time domain), as well as all of the unknown piecewise-constant functions between any pairs of the change points.

Instead of considering parameters that are piecewise-constant functions of time, in this paper we allow the parameters to be smooth functions of rescaled time, permitting the process to be locally stationary. The framework of local stationarity introduced by Dahlhaus (1997), where the parameter curves are defined in rescaled time u = t/T, provides a meaningful asymptotic theory. Locally stationary Y(t) means that if the functions m(u) and g(u) in Equation (1) are "smooth" and T is large, $m\left(\tfrac{t}{T}\right)\approx m\left(\tfrac{r}{T}\right)$ and $g\left(\tfrac{t}{T}\right)\approx g\left(\tfrac{r}{T}\right)$ for values of r close to t, that is, locally around t. More precisely, we assume that the functions m(x) and g(x) are Lipschitz continuous, that is, there exist constants C_m and C_g such that

$$\begin{eqnarray}&&\begin{array}{l}| m(z)-m(u)| \leqslant {C}_{m}| z-u| \qquad \text{and}\\ \ \ \times \,| g(z)-g(u)| \leqslant {C}_{g}| z-u| ,\quad {\rm{for}}\,{\rm{all}}\,u,z\in [0,1].\end{array}\end{eqnarray} \tag{ 2 }$$

To define m(u) and g(u) as functions of rescaled time let us consider, for each fixed u ∈ [0, 1] and increasing T, the sequence t = t_T = ⌊u T⌋, where ⌊x⌋ is the largest integer not exceeding x. Then we obtain the following uniform bound: ${\sup }_{t}| \tfrac{t}{T}-u| \lt \tfrac{1}{T}$. A time series is stationary if the moments of the underlying stochastic process, such as expectation and variance, are time-invariant. The idea behind local stationarity is to allow for time-varying parameters, in a way that locally the process behaves as stationary. Lipschitz continuity is a smoothness assumption that implies uniform continuity. The model in Equations (1)–(2) is a locally stationary process written in rescaled time in a way such that, as T grows, we observe more and more "observations" of the same type around u. That is, if m( · ) and g( · ) are smooth we have

$$\begin{eqnarray*}&&\begin{array}{l}| m\left(\displaystyle \frac{t}{T}\right)-m(u)| \leqslant {C}_{m}| \displaystyle \frac{t}{T}-u| \lt \displaystyle \frac{{C}_{m}}{T}\to 0\quad \text{and}\\ \ \ \times \,| g\left(\displaystyle \frac{t}{T}\right)-g(u)| \leqslant {C}_{g}| \displaystyle \frac{t}{T}-u| \lt \displaystyle \frac{{C}_{g}}{T}\to 0\quad {\rm{as}}\ T\to \infty ,\end{array}\end{eqnarray*}$$

for all t: = ⌊u T⌋. The locally stationary framework is important for handling, in a meaningful way, the asymptotic theory arising in statistics for processes with time-varying parameters. Suppose that we observe X_t = μ(t) + z_t , with ${z}_{t}\sim \mathrm{WN}(0,{\sigma }_{z}^{2})$ for t = 1, ..., T. Inference in this case means studying the properties of an estimator for the unknown function $\mu \left(t\right)$ on the grid {1, ..., T}. Given that μ changes over time, it is obvious that an asymptotic approach where T → ∞ is not suitable for describing a statistical method, since future "observations" {μ(t), t > T} do not necessarily contain any information on $\mu \left(t\right)$ on {1, ..., T}. To overcome these problems, Dahlhaus (1996) suggested to consider a triangular array of data. In analogy with nonparametric regression, it seems natural to set down the asymptotic theory in a way that we "observe" $\mu \left(t\right)$ on a finer grid (but on the same interval), i.e., that we observe the process ${Y}_{T}(t)=\mu \left(\tfrac{t}{T}\right)+{z}_{t}$, where Y_T is now a triangular array and μ is now rescaled to the interval [0, 1]. Working in rescaled time is often adopted also within the estimation framework of regression splines (see Zhou et al. 1998, among others).

Time series analysis of nonstationary sequences can be deterministic or stochastic. A popular example of stochastic nonstationarity is the well-known class of integrated processes, where the observed times series can be made stationary by differencing. Feigelson et al. (2018) apply autoregressive integrated moving average (ARIMA) models to light curves of several variable stars, discussing their effectiveness for different temporal characteristics. The process {X_t } is an ARIMA(p, d, q) process if Y_t = (1 − B)^d X_t , obtained by applying the operator 1 − B repeatedly d times, is a stationary ARMA(p, q) process. The most popular example of an ARIMA(p, d, q) process is the "random walk" X_t = X_t−1 + Z_t , which is an ARIMA with p = q = 0 and d = 1.

In the next section, we review the models proposed by Benkő et al. (2011) and Benkő (2018) for Blazhko light curves. Interestingly, our model in Equation (1) generalizes the models by Benkő et al. (2011) and Benkő (2018) in the sense that the modulating functions g( · ) are not confined to the class of parametric (sinusoidal or nonsinusoidal) functions.

The Blazhko effect is a periodic amplitude and phase variation in the light curves of RR Lyrae variable stars. In astronomy, the Blazhko effect is usually interpreted as a modulation phenomenon. Modulation is the process of transmitting a low-frequency signal into a high-frequency wave, called the carrier wave, by changing its amplitude, frequency, or phase angle through the modulating signal. In Appendix A we review the main mathematical definitions underlying the modulation phenomenon in astrophysics.

In this section, we review the models proposed by Benkő et al. (2011) and Benkő (2018), respectively, and we compare them with our novel model in Equation (1). To describe Blazhko light curves, Benkő et al. (2011) proposed to fit the following model:

$$\begin{eqnarray}\begin{array}{rcl}{\mu }^{* }(t) & = & {a}_{0}^{A}{a}_{0}+{a}_{0}{g}^{A}(t)\\ & & +\,\sum _{k=1}^{K}\left[{a}_{0}^{A}{a}_{k}+{a}_{k}{g}^{A}(t)\right]\sin [2\pi {{kf}}_{0}t+{\varphi }_{k}+{{kg}}^{F}(t)],\end{array}\end{eqnarray} \tag{ 3 }$$

where a_k and f₀ denote amplitude and frequency, respectively, and

$$\begin{eqnarray}\begin{array}{rcl}{g}^{M}(t) & = & \sum _{j=1}^{{{\ell }}^{M}}{a}_{j}^{M}\sin (2\pi {{jf}}_{m}t+{\varphi }_{j}^{M}),\\ M & = & A\ \mathrm{or}\ F.\end{array}\end{eqnarray} \tag{ 4 }$$

More recently, Benkő (2018) introduced a similar model:

$$\begin{eqnarray}\begin{array}{rcl}{\mu }^{* }(t) & = & {m}_{0}+\sum _{r=1}^{{\ell }}{b}_{r}\sin (2\pi {{rf}}_{m}t+{\varphi }_{r}^{b})\\ & & +\,\sum _{k=1}^{K}\left[{a}_{k}+{g}_{k}^{A}\left(t\right)\right]\sin [2\pi {{kf}}_{0}t+{\varphi }_{k}+{g}_{k}^{F}\left(t\right)],\end{array}\end{eqnarray} \tag{ 5 }$$

where b_r and ${\phi }_{r}^{b}$ denote amplitude and frequency of the modulating signal, respectively, and

$$\begin{eqnarray}&&{g}_{k}^{M}\left(t\right)=\sum _{j=1}^{{{\ell }}_{k}^{M}}{a}_{{kj}}^{M}\sin (2\pi {{jf}}_{m}t+{\varphi }_{{kj}}^{M}),\quad M=A\ \mathrm{or}\ F.\end{eqnarray} \tag{ 6 }$$

The functions ${g}^{M}\left(t\right)$ and ${g}_{k}^{M}\left(t\right)$ in Equation (4) and (6) are the modulating functions with subscripts M = A and M = F denoting amplitude and frequency modulation, respectively. The main pulsation frequency is denoted by f₀, whereas f_m is the modulating frequency. In this paper we improve the models in Equations (3)–(4) and (5)–(6) from two different viewpoints. From the modeling viewpoint, we relax the assumption of parametric amplitude and frequency modulations. Assuming parametric amplitude and frequency modulations results in a Fourier sum with time-invariant amplitudes and time-invariant frequencies, whereas our time-varying amplitudes and frequencies do not obey any particular form. From the estimation viewpoint, we do not rely on the nonlinear least-squares algorithms, such as the Levenberg−Marquardt algorithm, that are typically used to fit parametric nonlinear models. These methods require initial values close to the solution, which in some applications are difficult to find.

Both models proposed by Benkő et al. (2011) and Benkő (2018) and given by Equations (3) and (5), respectively, are a special case of our model defined by Equation (1). To see this, let us define

$$\begin{eqnarray}\begin{array}{rcl}v(t) & = & {a}_{0}^{A}{a}_{0}+{a}_{0}{g}^{A}(t),\\ {w}_{1,k}\left(t\right) & = & \left[{a}_{0}^{A}{a}_{k}+{a}_{k}{g}^{A}(t)\right]\sin [{\varphi }_{k}+{{kg}}^{F}(t)],\\ k & = & 1,\,\ldots ,\,K,\\ {w}_{2,k}\left(t\right) & = & \left[{a}_{0}^{A}{a}_{k}+{a}_{k}{g}^{A}(t)\right]\cos [{\varphi }_{k}+{{kg}}^{F}(t)],\\ k & = & 1,\,\ldots ,\,K,\end{array}\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}u(t) & = & {m}_{0}+\sum _{r=1}^{{\ell }}{b}_{r}\sin (2\pi {{rf}}_{m}t+{\varphi }_{r}^{b}),\\ {h}_{1,k}\left(t\right) & = & \left[{a}_{k}+{g}_{k}^{A}\left(t\right)\right]\sin [{\varphi }_{k}+{g}_{k}^{F}\left(t\right)],\\ k & = & 1,\,\ldots ,\,K,\\ {h}_{2,k}\left(t\right) & = & \left[{a}_{k}+{g}_{k}^{A}\left(t\right)\right]\cos [{\varphi }_{k}+{g}_{k}^{F}\left(t\right)],\\ k & = & 1,\,\ldots ,\,K.\end{array}\end{eqnarray} \tag{ 8 }$$

We now show how Equations (7) and (8) allow to us compare our model in Equation (1) with the models proposed by Benkő et al. (2011) and Benkő (2018), respectively. Comparing the models in Equations (1) and (3), the time-varying trend and amplitudes of the model in Equation (1) are expressed as

$$\begin{eqnarray}\begin{array}{rcl}m(t) & = & v(t),\\ {g}_{{\ell },k}\left(t\right) & = & {w}_{{\ell },k}\left(t\right),\quad {\ell }=1,2,\quad k=1,\,\ldots ,\,K.\end{array}\end{eqnarray} \tag{ 9 }$$

At the same time, comparing the model in Equation (1) with the model in Equation (5), the time-varying trend and amplitudes of the model in Equation (1) are

$$\begin{eqnarray}\begin{array}{rcl}m(t) & = & u(t),\\ {g}_{{\ell },k}\left(t\right) & = & {h}_{{\ell },k}\left(t\right),\quad {\ell }=1,2,\quad k=1,\,\ldots ,\,K,\end{array}\end{eqnarray} \tag{ 10 }$$

the ordinary frequency being f_k = kf₀.

As we can see in Equations (9) and (10), the functions m(t) and ${g}_{{\ell },k}\left(t\right)$ incorporate the amplitude and frequency modulation functions ${g}^{M}\left(t\right)$ and ${g}_{k}^{M}\left(t\right)$ in Equations (4) and (6). In this sense, the limitation of our approach is that it does not aim at identifying the amplitude and frequency modulating functions ${g}^{M}\left(t\right)$ and ${g}_{k}^{M}\left(t\right)$ in Equations (4) and (6). That said, the benefit of our approach from the estimation viewpoint is twofold. An important advantage of our model in Equation (1) over the models in Equations (3) and (5) is that the modulating frequency f_m does not need to be estimated. In other words, in order to describe statistically a Blazhko light curve using our model in Equation (1), we only need to estimate f₀. If the observed time series is indeed a Blazhko light curve, the modulating frequency f_m is included in the nonparametric trend m(t) and amplitude g_ℓ,k of our model in Equation (1). Moreover, assuming that the frequencies are known, for our model in Equation (1) we only need to estimate the functions m( · ) and g_ℓ,k ( · ), whereas for the model in Equations (3) and (5) the estimated parameters are the amplitudes ${a}_{0}^{A}$, a₀, m₀, a_k 's, b_r 's, and ${a}_{{}_{{kj}}}^{M}$'s, and the phases ${\varphi }_{r}^{b}$'s, φ_k 's, ${\varphi }_{j}^{M}$'s, and ${\varphi }_{{kj}}^{M}$'s.

In Section 4.1 we define estimators of the unknown trend m( · ), amplitudes {g_ℓ,k ( · ), ℓ = 1, 2, k = 1, ..., K}, and variance ${\sigma }_{z}^{2}$ of the model in Equation (1), and in Section 4.2 we explain how to select the tuning parameters associated with the B-splines and the penalization used in the estimation method. We denote by N the sample size, T = t_N − t₁ the time span, J the number of B-splines basis, d the degree of the B-splines, K the number of harmonics components, r the order of the penalty, and M the number of replications in Monte Carlo simulations.

We performed our calculations using the R Language for Statistical Computing (R Core Team 2021). Our codes combine existing functions (available as part of R packages) with our own development. The computations implemented in this paper are available as a GitHub public code repository. ⁶

4.1. Penalized Least Squares

As mentioned in Section 2, we use B-splines to estimate the trend and amplitudes of the model given by Equation (1). The smooth trend function m(t_i ) is modeled as a linear combination of B-splines basis

$$\begin{eqnarray*}&&m({t}_{i})=\sum _{j=1}^{J}{\alpha }_{j}{B}_{j}({t}_{i}),\quad i=1,\,\ldots ,\,N,\end{eqnarray*}$$

which can be written in matrix notation as

$$\begin{eqnarray*}&&{\boldsymbol{m}}={\bf{B}}{\boldsymbol{\alpha }},\end{eqnarray*}$$

where ${\boldsymbol{m}}={[m({t}_{1}),\,\ldots ,\,m({t}_{N})]}^{\top }$, B = [B_ij ] = [B_j (t_i )] is the N × J basis matrix (i = 1, ..., N, j = 1, ..., J) and ${\boldsymbol{\alpha }}={\left({\alpha }_{1},\,\ldots ,\,{\alpha }_{J}\right)}^{\top }$. The exact definition of B-splines is given in Appendix B.

The smooth amplitude functions, g_ℓ,k (t_i ), ℓ = 1, 2, are modeled in the same way:

$$\begin{eqnarray*}\begin{array}{rcl}{g}_{1,k}({t}_{i}) & = & \sum _{j=1}^{J}{\beta }_{k,j}{B}_{j}({t}_{i}),\\ {g}_{2,k}({t}_{i}) & = & \sum _{j=1}^{J}{\gamma }_{k,j}{B}_{j}({t}_{i}),\quad k=1,\,\ldots ,\,K.\end{array}\end{eqnarray*}$$

In matrix notation

$$\begin{eqnarray*}&&{{\boldsymbol{g}}}_{1,k}={\bf{B}}{{\boldsymbol{\beta }}}_{k}\,{\rm{a}}{\rm{n}}{\rm{d}}\,{{\boldsymbol{g}}}_{2,k}={\bf{B}}{{\boldsymbol{\gamma }}}_{k},\,k=1,\,\ldots ,\,K,\end{eqnarray*}$$

where ${{\boldsymbol{\beta }}}_{k}={\left({\beta }_{k,1},\,\ldots ,\,{\beta }_{k,J}\right)}^{\top }$, ${{\boldsymbol{\gamma }}}_{k}={\left({\gamma }_{k,1},\,\ldots ,\,{\gamma }_{k,J}\right)}^{\top }$, and ${{\boldsymbol{g}}}_{{\ell },k}\,={[{g}_{{\ell },k}({t}_{1}),\,\ldots ,\,{g}_{{\ell },k}({t}_{N})]}^{\top }$, ℓ = 1, 2, k = 1, ..., K. Thus, α , β _k , and γ _k , k = 1, ..., K, are vectors associated with the trend and amplitudes, respectively. We define the N × N matrices C_k and S_k as

$$\begin{eqnarray*}\begin{array}{r}\begin{array}{rcl}{{\bf{C}}}_{k} & = & {\rm{d}}{\rm{i}}{\rm{a}}{\rm{g}}\{\cos ({w}_{k}{t}_{1}),\,\ldots ,\,\cos ({w}_{k}{t}_{N})\}\,{\rm{a}}{\rm{n}}{\rm{d}}\,\\ {{\bf{S}}}_{k} & = & {\rm{d}}{\rm{i}}{\rm{a}}{\rm{g}}\{\sin ({w}_{k}{t}_{1}),\,\ldots ,\,\sin ({w}_{k}{t}_{N})\},\\ k & = & 1,\,\ldots ,\,K.\end{array}\end{array}\end{eqnarray*}$$

Thus, the model for the expected value of Y , in matrix notation, can be expressed as

$$\begin{eqnarray*}&&{\rm{E}}\left[{\boldsymbol{Y}}\right]={\boldsymbol{\mu }}={\boldsymbol{ \mathcal B }}{\boldsymbol{\theta }},\end{eqnarray*}$$

where ${\boldsymbol{ \mathcal B }}$ is the N × c design matrix given by

$$\begin{eqnarray*}&&{\boldsymbol{ \mathcal B }}=[{\bf{B}}| {{\bf{C}}}_{1}{\bf{B}}| ...| {{\bf{C}}}_{K}{\bf{B}}| {{\bf{S}}}_{1}{\bf{B}}| ...| {{\bf{S}}}_{K}{\bf{B}}],\end{eqnarray*}$$

with c = J(2K + 1), and

$$\begin{eqnarray*}&&{\boldsymbol{\theta }}={\left({{\boldsymbol{\alpha }}}^{\top },{{\boldsymbol{\beta }}}_{1}^{\top },\,\ldots ,\,{{\boldsymbol{\beta }}}_{K}^{\top },{{\boldsymbol{\gamma }}}_{1}^{\top },\,\ldots ,\,{{\boldsymbol{\gamma }}}_{K}^{\top }\right)}^{\top }\end{eqnarray*}$$

is the vector of regression coefficients of length c.

The ordinary least squares (OLS) estimator of θ is the vector ${\widehat{{\boldsymbol{\theta }}}}_{\mathrm{OLS}}$, which minimizes the sum of squares:

$$\begin{eqnarray*}&&{M}_{{\boldsymbol{\theta }}}=| | {\boldsymbol{Y}}-{\boldsymbol{ \mathcal B }}{\boldsymbol{\theta }}| {| }^{2}.\end{eqnarray*}$$

Equating to zero the partial derivatives with respect to each component of θ and assuming (as we shall) that ${{\boldsymbol{ \mathcal B }}}^{\top }{\boldsymbol{ \mathcal B }}$ is nonsingular, the estimator of θ is

$$\begin{eqnarray*}&&{\widehat{{\boldsymbol{\theta }}}}_{\mathrm{OLS}}={\left({{\boldsymbol{ \mathcal B }}}^{\top }{\boldsymbol{ \mathcal B }}\right)}^{-1}{{\boldsymbol{ \mathcal B }}}^{\top }{\boldsymbol{Y}}.\end{eqnarray*}$$

The OLS estimate also maximizes the likelihood of the observations when the errors z₁, ..., z_N are independent and identically distributed and Gaussian.

The size of the basis determines the amount of smoothing of the fitted curves. The larger the value of J, the bumpier the fitting will be. To avoid overfitting, Eilers & Marx (1996) proposed a penalty on the (high-order) finite differences of the coefficients:

$$\begin{eqnarray*}\begin{array}{cll}{M}_{{\boldsymbol{\theta }}}^{* } & = & | | {\boldsymbol{Y}}-{\boldsymbol{ \mathcal B }}{\boldsymbol{\theta }}| {| }^{2}+{\tau }_{1}| | {{\bf{D}}}_{r}{\boldsymbol{\alpha }}| {| }^{2}\\ & & +\,\displaystyle \sum _{k=1}^{K}\left\{{\tau }_{2k}| | {{\bf{D}}}_{r}{{\boldsymbol{\beta }}}_{k}| {| }^{2}+{\tau }_{2k+1}| | {{\bf{D}}}_{r}{{\boldsymbol{\gamma }}}_{k}| {| }^{2}\right\},\end{array}\end{eqnarray*}$$

where {τ_k , k = 1, ..., 2K + 1} are positive regularization parameters that control the smoothness of the curve, penalizing the coefficients that are far apart from one another. If τ_k = 0, k = 1, ..., 2K + 1, we have the standard normal equations of linear regression with a B-splines basis. The larger the value of τ_k , the closer the coefficient θ is to zero. When τ_k → ∞ we obtain a polynomial fit. The matrix D_r constructs rth-order differences of a vector η as

$$\begin{eqnarray*}&&{{\bf{D}}}_{r}{\boldsymbol{\eta }}={{\rm{\Delta }}}^{r}{\boldsymbol{\eta }}.\end{eqnarray*}$$

The first difference of η , Δ¹ η , is the vector with elements η _l+1 − η _l . Repeated differencing applied to Δ η results in higher differences, such as Δ² η and Δ³ η .

The penalties can be represented as θ ^⊤ P θ with the block-diagonal matrix ${\bf{P}}={\bf{T}}\otimes {{\bf{D}}}_{r}^{\top }{{\bf{D}}}_{r}$ and T = diag{τ₁, τ₂, τ₃, ..., τ_2K+1}. Then, minimizing

$$\begin{eqnarray*}&&{M}_{{\boldsymbol{\theta }}}^{* }=| | {\boldsymbol{Y}}-{\boldsymbol{ \mathcal B }}{\boldsymbol{\theta }}| {| }^{2}+{{\boldsymbol{\theta }}}^{\top }{\bf{P}}{\boldsymbol{\theta }}\end{eqnarray*}$$

with respect to θ , the penalized ordinary least squares estimator (POLS) of θ is

$$\begin{eqnarray}&&{\hat{{\boldsymbol{\theta }}}}_{\mathrm{POLS}}={\left({{\boldsymbol{ \mathcal B }}}^{\top }{\boldsymbol{ \mathcal B }}+{\bf{P}}\right)}^{-1}{{\boldsymbol{ \mathcal B }}}^{\top }{\boldsymbol{Y}}.\end{eqnarray} \tag{ 11 }$$

The prediction of Y at time t_i is given by

$$\begin{eqnarray}&&{\widehat{Y}}_{i}=\widehat{\mu }({t}_{i})={\boldsymbol{ \mathcal B }}{\left({t}_{i}\right)}^{\top }{\widehat{{\boldsymbol{\theta }}}}_{\mathrm{POLS}},\end{eqnarray} \tag{ 12 }$$

where ${\boldsymbol{ \mathcal B }}({t}_{i})$ is the ith row of ${\boldsymbol{ \mathcal B }}$, the residuals are $\widehat{{\boldsymbol{z}}}={\boldsymbol{Y}}-\widehat{{\boldsymbol{Y}}}$, with $\widehat{{\boldsymbol{Y}}}={\left({\widehat{Y}}_{1},\,\ldots ,\,{\widehat{Y}}_{N}\right)}^{\top }$, and the mean square error (MSE) is $\mathrm{MSE}={N}^{-1}{\sum }_{i=1}^{N}{\left({Y}_{i}-{\widehat{Y}}_{i}\right)}^{2}$.

The estimators of the trend m and amplitudes { g _ℓ,k , ℓ = 1, 2, k = 1, ..., K}, are

$$\begin{eqnarray}&&\hat{{\boldsymbol{m}}}={\bf{B}}\hat{{\boldsymbol{\alpha }}},\quad {\hat{{\boldsymbol{g}}}}_{1,k}={\bf{B}}{\hat{{\boldsymbol{\beta }}}}_{k},\quad {\hat{{\boldsymbol{g}}}}_{2,k}={\bf{B}}{\hat{{\boldsymbol{\gamma }}}}_{k}.\end{eqnarray} \tag{ 13 }$$

Another parameter of interest is the variance of the errors, ${\sigma }_{z}^{2}$, which can be estimated by

$$\begin{eqnarray*}&&{\hat{\sigma }}_{z}^{2}={\left[N-\mathrm{tr}\left(\hat{{\bf{S}}}\right)\right]}^{-1}\displaystyle \sum _{i=1}^{N}{\left\{{Y}_{i}-{\boldsymbol{ \mathcal B }}{\left({t}_{i}\right)}^{\top }{\hat{{\boldsymbol{\theta }}}}_{\mathrm{POLS}}\right\}}^{2}\end{eqnarray*}$$

where $\hat{{\bf{S}}}={\boldsymbol{ \mathcal B }}{\left({{\boldsymbol{ \mathcal B }}}^{\top }{\boldsymbol{ \mathcal B }}+{\bf{P}}\right)}^{-1}{{\boldsymbol{ \mathcal B }}}^{\top }.$

In addition to the point estimate, interval estimation for ${\widehat{Y}}_{i}$ is often of interest and is easy to construct. In Appendix C we derive parametric and nonparametric confidence intervals for ${\widehat{Y}}_{i}$.

4.2. Automatic Selection of the Tunable Parameters

Before calculating the estimator in Equation (11), it is necessary to select the tuning parameters ${\boldsymbol{\tau }}\,={\left({\tau }_{1},{\tau }_{2},\,\ldots ,\,{\tau }_{2K+1}\right)}^{\top }$. To choose the tuning parameters, we propose to use the Akaike information criterion (AIC).

The AIC penalizes the log-likelihood of a fitted model by considering the effective number of parameters. The definition of AIC given by Hastie et al. (2004) is

$$\begin{eqnarray*}&&\mathrm{AIC}({\boldsymbol{\tau }})=\overline{\mathrm{err}}({\boldsymbol{\tau }})+2\displaystyle \frac{\mathrm{df}}{N}{\widehat{\sigma }}_{0}^{2},\end{eqnarray*}$$

where $\overline{\mathrm{err}}({\boldsymbol{\tau }})$ corresponds to the mean square error in the case of Gaussian errors, df is the effective number of parameters, N is the number of observations used to fit the model, and ${\widehat{\sigma }}_{0}^{2}$ is given by the variance of the residuals from the ${\widehat{Y}}_{i}$ that are computed when τ = 0_2K+1.

The value for τ is chosen by minimizing the AIC, which is computed as

$$\begin{eqnarray}&&\mathrm{AIC}\left({\boldsymbol{\tau }}\right)=\frac{1}{N}\displaystyle \sum _{i=1}^{N}{\left\{{Y}_{i}-{\boldsymbol{ \mathcal B }}{\left({t}_{i}\right)}^{\top }{\hat{{\boldsymbol{\theta }}}}_{\mathrm{POLS}}\right\}}^{2}+2\frac{\mathrm{tr}\left(\hat{{\bf{S}}}\right)}{N}{\hat{\sigma }}_{0}^{2},\end{eqnarray} \tag{ 14 }$$

The AIC given by Equation (14) can also be used to select the number of B-splines J, the degree d of the B-spline, the order of penalty r, and the number of harmonic components K.

In Figure 1, we have generated N = 500 observations from the model described in Equation (1), with the Gaussian errors {z_i , 1 ≤ i ≤ 500} being simulated using the R function rnorm. We consider the following artificial signal:

$$\begin{eqnarray*}&&\begin{array}{l}\mu ({t}_{i})=-0.05{t}_{i}-(-0.0002{t}_{i}+0.0003{t}_{i}^{2})\\ \ \ \times \,\cos (0.2\pi {t}_{i})+(1-0.0005{t}_{i})\sin (0.2\pi {t}_{i}),\end{array}\end{eqnarray*}$$

with the errors following a Gaussian distribution with zero mean and variance ${\sigma }_{z}^{2}=1$. Time t is unequally spaced and was obtained from a uniform distribution U(θ₁, θ₂) with θ₁ = 0 and θ₂ = 55 using the R function runif. In the first plot of Figure 1 the observations Y are represented by the gray points, and the mean μ by the black curve. The orange, blue, and green curves illustrate three possible estimates for Y obtained using the method described in the Section 4.1 with increasing smoothing parameters. The orange line is the fit obtained with τ_j = 0, j = 1, 2, 3: the corresponding $\widehat{{\boldsymbol{Y}}}$ matches the data well, but fits the true μ poorly because it is wiggly. The blue curve is obtained using the smoothing parameters τ_j = 30, j = 1, 2, 3, and the green curve is obtained using τ_j = 200, j = 1, 2, 3. In the second plot of Figure 1, we observe that the optimal tuning parameters are τ_j = 30, j = 1, 2, 3, and as the values of τ increase the obtained curve fits the observed data less closely.

Zoom In Zoom Out Reset image size

A statistical model is an approximation to the true process that generates the observed data. After fitting the model given by Equation (1), it is necessary to check whether the residuals obtained from the fit behave like a white noise process. A significant departure from this assumption suggests the inadequacy of the assumed form of the model. Thus, it is important to assess whether the residuals follow a white noise process.

Detecting serial correlation becomes more challenging when the available observations are unequally spaced in time. If the observations are unequally spaced, so are the errors. In order to study the spectral density of the residuals obtained from the fitted model, in this section we derive the mathematical definition of the spectrum for unequally spaced time series.

Before presenting our approach, we briefly review the results given by Deeming (1975) about the relationships between the periodogram, the spectral density (PSD), and the autocorrelation function for continuous time series. Then, we extend the results given by Deeming (1975) to the case of discrete time series.

Let {ε_i } be a continuous, zero-mean stationary times series with spectral density

$$\begin{eqnarray*}&&{P}_{\varepsilon }(\lambda )={\int }_{-\infty }^{\infty }{r}_{\varepsilon }(h)\exp (i\lambda h){dh},\quad -\infty \lt \lambda \lt \infty ,\end{eqnarray*}$$

and autocovariance function

$$\begin{eqnarray*}&&{r}_{\varepsilon }(h)={\int }_{-\infty }^{\infty }{P}_{\varepsilon }(\lambda )\exp (-i\lambda h)d\lambda ,\quad h\in {\mathbb{R}}.\end{eqnarray*}$$

Consider a time series ε₁, ..., ε_N with spectrum P_ε ( · ) and autocovariance function r_ε ( · ), and assume that the observations ε₁, ..., ε_N are obtained at unequally spaced times t₁, ..., t_N , respectively. The periodogram of ${\boldsymbol{\varepsilon }}={\left({\varepsilon }_{1},\,\ldots ,\,{\varepsilon }_{N}\right)}^{\top }$ at frequency λ is defined as

$$\begin{eqnarray}&&{I}_{\,\varepsilon }(\lambda )=\sum _{k=1}^{N}\sum _{j=1}^{N}{\varepsilon }_{k}{\varepsilon }_{j}\exp (i\lambda [{t}_{k}-{t}_{j}]),\quad \lambda =2\pi f.\end{eqnarray} \tag{ 15 }$$

Deeming (1975) proved that the expectation of the periodogram of {ε_i } in Equation (15) is

$$\begin{eqnarray}&&{\rm{E}}\left[{I}_{\,\varepsilon }(\lambda )\right]={P}_{\varepsilon }(\lambda )\star {W}_{\varepsilon }(\lambda ),\end{eqnarray} \tag{ 16 }$$

where W_ε (λ) is the power spectral window given by

$$\begin{eqnarray*}&&{W}_{\varepsilon }(\lambda )=\sum _{j=1}^{N}\sum _{k=1}^{N}\exp (i\lambda [{t}_{k}-{t}_{j}]),\end{eqnarray*}$$

and P_ε (λ)⋆W_ε (λ) is the continuous convolution of P_ε (λ) with W_ε (λ) defined as

$$\begin{eqnarray*}&&{P}_{\varepsilon }(\lambda )\star {W}_{\varepsilon }(\lambda )={\int }_{-\infty }^{\infty }{P}_{\varepsilon }(\omega ){W}_{\varepsilon }(\lambda -\omega )d\omega .\end{eqnarray*}$$

The following lemma states that it is possible to extend the result in Equation (16) to the case of a discrete zero-mean stationary times series that is generated according to equally spaced times but observed at unequally spaced times. The lemma applies to unequally spaced time points t_i with index i belonging to a subset ${ \mathcal I }$ of the set ${\mathbb{N}}=\{1,\,2,\,...\}$ of positive integers.

Lemma 1. Let $\{{\varepsilon }_{i},\ {t}_{i}={t}_{0}+i{\rm{\Delta }},\ {\rm{\Delta }}\gt 0,\ i\in { \mathcal I }\subseteq {\mathbb{N}}\}$ be a zero-mean, stationary, discrete time series with spectral density

$$\begin{eqnarray}&&{P}_{\varepsilon }(\lambda )=\displaystyle \frac{1}{2\pi }\sum _{h=-\infty }^{\infty }\exp (i\lambda h{\rm{\Delta }}){r}_{\varepsilon }(h),\quad -\infty \leqslant \lambda \leqslant \infty ,\end{eqnarray} \tag{ 17 }$$

with autocovariance function defined as ${r}_{\varepsilon }(h)={\rm{E}}\left[{\varepsilon }_{k}{\varepsilon }_{j}\right]$, with $k=j+| h|$, $h\in {\mathbb{Z}}$, that can be expressed in term of the spectral density in Equation (17) as

$$\begin{eqnarray}&&{r}_{\varepsilon }(h)=\displaystyle \frac{2\pi }{{N}_{\,{ \mathcal I }}}\sum _{j=1}^{{N}_{\,{ \mathcal I }}}\exp (-i{\lambda }_{j}h{\rm{\Delta }}){P}_{\varepsilon }({\lambda }_{j}),\quad h\in {\mathbb{Z}},\end{eqnarray} \tag{ 18 }$$

where ${\lambda }_{j}=2\pi {f}_{j}$, with ${f}_{j}=j/({N}_{{ \mathcal I }}{\rm{\Delta }})$ and ${N}_{{ \mathcal I }}=\max {\{I\}}$. Then, the expectation of the periodogram in Equation (15) obtained from $\{{\varepsilon }_{i},i\in { \mathcal I }\}$ is

$$\begin{eqnarray}&&{\rm{E}}\left[{I}_{\,\varepsilon }(\lambda )\right]=\displaystyle \frac{2\pi }{{N}_{\,{ \mathcal I }}}\,{P}_{\varepsilon }(\lambda )* {W}_{\varepsilon }(\lambda ),\end{eqnarray} \tag{ 19 }$$

with power spectral window given by

$$\begin{eqnarray}&&{W}_{\varepsilon }(\lambda )=\sum _{k\in { \mathcal I }}\sum _{j\in { \mathcal I }}\exp (i\lambda [{t}_{k}-{t}_{j}])\end{eqnarray} \tag{ 20 }$$

and ${P}_{\varepsilon }(\lambda )* {W}_{\varepsilon }(\lambda )$ is the discrete convolution of ${P}_{\varepsilon }(\lambda )$ with ${W}_{\varepsilon }(\lambda )$ defined as

$$\begin{eqnarray*}\begin{array}{rcl}{P}_{\varepsilon }(\lambda )* {W}_{\varepsilon }(\lambda ) & = & \sum _{j=1}^{{N}_{\,{ \mathcal I }}}{P}_{\varepsilon }({\omega }_{j})* {W}_{\varepsilon }(\lambda -{\omega }_{j}),\\ {\omega }_{j} & = & 2\pi {f}_{j},\quad {f}_{j}=\displaystyle \frac{j}{{N}_{\,{ \mathcal I }}{\rm{\Delta }}}.\end{array}\end{eqnarray*}$$

Our result in Equation (19) differs from the result by Deeming (1975) in Equation (16). Deeming (1975) proved that the expectation of both discrete and continuous Fourier transforms of a continuous stochastic process f(t); (in the sense of Equations (31) and (32) in Deeming 1975) is equal to the continuous convolution of the spectral density of f(t) with a spectral window (see Equations (33) and (36) in Deeming 1975). In Lemma 1, instead, we prove that the expectation of the discrete Fourier transform of the discrete stochastic process ε_t is equal to the discrete convolution of the spectral density of ε_t with a spectral window (up to the constant $2\pi /{N}_{{ \mathcal I }}$).

When the time series is generated according to an equally spaced stochastic process and the observations are equally spaced, the periodogram is an (asymptotically) unbiased estimator of the spectral density (see Priestley 1981, page 418). However, when the observations are unequally spaced it does not make sense to estimate the spectral density in the same way. This is due to the power spectral window W_ε (λ) in Equations (19)–(20). Nevertheless, as we show in the following proposition, it is possible to disentangle the spectral density P_ε (λ) from the spectral window W_ε (λ).

Proposition 1. Let ${ \mathcal F }\{{g}_{j}\}[k]$ denote the Discrete Fourier Transform of the sequence of m numbers ${g}_{1},\,\ldots ,\,{g}_{m}$ into another sequence ${h}_{1},\,\ldots ,\,{h}_{m}$, that is,

$$\begin{eqnarray}\begin{array}{ccc}{h}_{k} & = & { \mathcal F }\{{g}_{j}\}[k]=\displaystyle \sum _{j=1}^{m}{g}_{j}\exp (-{ikj}2\pi \left./m\right),\\ k & = & 1,\,\ldots ,\,m.\end{array}\end{eqnarray} \tag{ 21 }$$

Accordingly, define ${{ \mathcal F }}^{-1}\{{h}_{k}\}[j]$ as the Inverse Discrete Fourier Transform of the sequence ${h}_{1},\,\ldots ,\,{h}_{m}$ into another sequence ${g}_{1},\,\ldots ,\,{g}_{m}$, that is,

$$\begin{eqnarray}\begin{array}{ccc}{g}_{j} & = & {{ \mathcal F }}^{-1}\{{h}_{k}\}[j]=\frac{1}{m}\displaystyle \sum _{k=1}^{m}{h}_{k}\exp ({ikj}2\pi /m),\\ j & = & 1,\,\ldots ,\,m.\end{array}\end{eqnarray} \tag{ 22 }$$

Assume that $\{{\varepsilon }_{i},i\in { \mathcal I }\}$ satisfy the same conditions as in Lemma 1. Then we can write the spectral density ${P}_{\varepsilon }({\lambda }_{j})$ in Equation (17) at frequency ${\lambda }_{j}=2\pi {f}_{j}$, with ${f}_{j}=j/({N}_{\,{ \mathcal I }}{\rm{\Delta }})$, as

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\varepsilon }({\lambda }_{j}) & = & \displaystyle \frac{{N}_{\,{ \mathcal I }}}{2\pi }{{ \mathcal F }}^{-1}\left\{\displaystyle \frac{{ \mathcal F }\{{\rm{E}}\left[{I}_{\varepsilon }({\lambda }_{j})\right]\}[k]}{{ \mathcal F }\{{W}_{\varepsilon }({\lambda }_{j})\}[k]}\right\}[j],\\ j & = & 1,\,\ldots ,\,{N}_{\,{ \mathcal I }}.\end{array}\end{eqnarray} \tag{ 23 }$$

The proofs of Lemma 1 and Proposition 1 are given in Appendix D. Equation (23) suggests that in order to estimate P_ε (λ), we need the value of Δ and W_ε (λ), as well as the estimate of ${\rm{E}}\left[{I}_{\,\varepsilon }(\lambda )\right]$. Notice that, in general, the time series denoted in this section as ${\boldsymbol{\varepsilon }}={\left({\varepsilon }_{1},\,\ldots ,\,{\varepsilon }_{N}\right)}^{\top }$ is possibly autocorrelated, whereas the errors ${\boldsymbol{z}}={\left({z}_{1},\,\ldots ,\,{z}_{N}\right)}^{\top }$ of our model in Equation (1) are assumed to be serially uncorrelated. The following algorithm is proposed to establish whether the unequally spaced residuals obtained when fitting our model in Equation (1) are uncorrelated.

• 1.  
  Obtain the residuals $\widehat{{\boldsymbol{z}}}={\boldsymbol{Y}}-\widehat{{\boldsymbol{Y}}}$.
• 2.  
  For each $j=1,\,\ldots ,\,{N}_{\,{ \mathcal I }}$, define ${\hat{I}}_{\,\widehat{z}}({\lambda }_{j})$ as the periodogram in Equation (15) computed upon the residuals $\widehat{{\boldsymbol{z}}}$ obtained in step (1), with ${\lambda }_{j}=\tfrac{2\pi j}{{N}_{\,{ \mathcal I }}{\rm{\Delta }}}$.
• 3.  
  For each $j=1,\,\ldots ,\,{N}_{\,{ \mathcal I }}$, define ${W}_{\widehat{z}}({\lambda }_{j})$ as the power spectral window in Equation (20) computed upon the residuals $\widehat{{\boldsymbol{z}}}$ obtained in step (1), with ${\lambda }_{j}=\tfrac{2\pi j}{{N}_{\,{ \mathcal I }}{\rm{\Delta }}}$.
• 4.  
  Smooth the periodogram obtained in step (2) over frequencies, and denote the smoothed periodogram by ${\tilde{I}}_{\,\widehat{z}}({\lambda }_{j})$.
• 5.  
  Calculate the Discrete Fourier Transform in Equation (21) of the power spectral window and the periodogram obtained in steps (3) and (4), respectively.
• 6.  
  For each frequency ${\lambda }_{j}=\tfrac{2\pi j}{{N}_{\,{ \mathcal I }}{\rm{\Delta }}}$, define the estimated spectral density of the errors z as

  $$\begin{eqnarray}\begin{array}{rcl}{\widehat{P}}_{z}({\lambda }_{j}) & = & \displaystyle \frac{{N}_{\,{ \mathcal I }}}{2\pi }{{ \mathcal F }}^{-1}\left\{\displaystyle \frac{{ \mathcal F }\{{\tilde{I}}_{\,\widehat{z}}({\lambda }_{j})\}[k]}{{ \mathcal F }\{{W}_{\widehat{z}}({\lambda }_{j})\}[k]}\right\}[j],\\ j & = & 1,\,\ldots ,\,{N}_{\,{ \mathcal I }},\end{array}\end{eqnarray} \tag{ 24 }$$

  where the inverse Fourier transform ${{ \mathcal F }}^{-1}$ is given by Equation (22).
• 7.  
  If the estimated spectral density obtained in step (6) does not vary significantly over frequencies, conclude that the errors are uncorrelated over time.

In this section we provide Monte Carlo simulations to illustrate the performance of the estimators $\widehat{\mu }\left(t\right)$, $\widehat{m}\left(t\right)$, ${\widehat{g}}_{{\ell },k}\left(t\right)$, ℓ = 1, 2, k = 1, ..., K, defined by Equations (12) and (13), and the estimator of the spectral density in Equation (23).

In Section 6.1 we simulate unequally spaced observations from the model in Equation (1), under two scenarios. In the first scenario both the trend and amplitudes are sinusoidal, whereas in the second scenario the trend and amplitudes are polynomial. In Section 6.2 we simulate a Blazhko light curve and fit the model in Equation (1). Finally in Section 6.3 we evaluate the performance of the estimator of the spectral density defined in Equation (23) of a discrete unequally spaced time series.

6.1. Simulating Our Novel Time-varying Model

In this section we generate the data according to the model described by Equation (1) with K = 2, N = 500, and time t is unequally spaced, obtained from a uniform distribution U(θ₁, θ₂) with θ₁ = 0 and θ₂ = 1. In order to illustrate the flexibility of our novel method, we consider two different scenarios for the trend and amplitudes. In the first scenario, we simulate a sinusoidal trend and amplitudes as $m(t)=\sin (2\pi t)$, ${g}_{\mathrm{1,1}}\left(t\right)=\cos (9\pi t)$, ${g}_{\mathrm{2,1}}(t)=\sin (6\pi t)$, ${g}_{\mathrm{1,2}}(t)=\cos (4\pi t)$, ${g}_{\mathrm{2,2}}(t)=\sin (7\pi t)$, with frequencies w₁ = 40π and w₂ = 100π. In the second scenario, we simulate a (global) polynomial trend and amplitudes as m(t) = 0.2t − 5t² + 5.5t³, ${g}_{\mathrm{1,1}}\left(t\right)=4{t}^{3}-5{t}^{2}$, g_2,1(t) = −0.5 − 0.5t + 2.5t² − 0.5t³, g_1,2(t) = −t + t² + 1.3t³, g_2,2(t) = 0.5 + 2t² − 3t³, with frequencies w₁ = 30π and w₂ = 40π. In both scenarios, we assume that the error terms {z_i , i = 1, ..., N} follow a Gaussian distribution with zero mean and variance ${\sigma }_{z}^{2}=2$.

In both scenarios, we simulate M = 200 realizations of the model in Equation (1). For each j = 1, ..., M, we compute the estimate ${\widehat{{\boldsymbol{\theta }}}}^{\left(j\right)}$ defined by Equation (11). In the first scenario, we select the smoothing parameter ${\boldsymbol{\tau }}={\left(\mathrm{50,}\,1,\,\mathrm{2,}\,\mathrm{10,}\,1\right)}^{\top }$, a total number of B-splines J = 33 of order d = 3, and an order penalty r = 2. In the second scenario, we choose the smoothing parameter ${\boldsymbol{\tau }}={\left(\mathrm{3,}\,3,\,\mathrm{3,}\,\mathrm{3,}\,3\right)}^{\top }$, a total number of B-splines J = 6 of order d = 3, and an order penalty r = 4. Figure 2 shows our estimates of μ , m , { g _ℓ,k , ℓ = 1, 2, k = 1, 2}, and their 95% confidence intervals. Figure 2 shows that our model in Equation (1) fits well the simulated data in both the sinusoidal and polynomial scenarios. That is, the trend and amplitudes are well fitted in both scenarios. The 95% confidence intervals are constructed in a nonparametric fashion using quantiles; see Appendix C.1 for more details.

Zoom In Zoom Out Reset image size

6.2. Simulating a Blazhko RR Lyrae Light Curve Characterized by Amplitude Modulation

We simulate a Blazhko RR Lyrae light curve with amplitude modulation according to Benkő et al. (2011) as

$$\begin{eqnarray}\begin{array}{rcl}{Y}_{i} & = & \mu ({t}_{i})+{z}_{i},\quad i=1,\,\ldots ,\,N=1000,\\ \mu ({t}_{i}) & = & \left[1+\displaystyle \frac{{U}_{m}({t}_{i})}{{U}_{c}}\right]\,c(t),\\ c(t) & = & {a}_{0}+\sum _{k=1}^{4}{a}_{k}\sin (2\pi {{kf}}_{0}{t}_{i}+{\varphi }_{k}),\\ {U}_{m}(t) & = & {a}_{m}\sin (2\pi {f}_{m}t+{\varphi }_{m}),\end{array}\end{eqnarray} \tag{ 25 }$$

where c(t) is the carrier wave with four harmonic components, ${U}_{m}\left(t\right)$ is the modulating signal, U_c = a_m /h is the amplitude of the unmodulated light curve, and {z_i , i = 1, ..., N} are the error terms. The values of the parameters used in Equation (25) and the time design are obtained from Benkő et al. (2011). In particular, a_m = 0.1 mag, h = 1.2, a₀ = 0.01 mag, f_m = 0.05 days⁻¹, φ_m = 270°, and the values {a_k , φ_k , 1 ≤ k ≤ 4} are presented in Table 1. We convert the Blazhko phase φ_m and the main phases {φ_k , 1 ≤ k ≤ 4} in Equation (25) from degrees to radians using the R function NISTdegTOradian (available in the R package NISTunits2016 by Gama 2016). The original time design {t_j , j = 1, ..., 28,799} in Benkő et al. (2011) is equally spaced. However, variable stars are often observed at irregular intervals. For this reason, in our simulation exercise we sample a subset of the original time points and use this subset to evaluate the performance of our method. We obtain the time design {t_i , i = 1, ..., 1000} in Equation(25) by sampling the original, equally spaced time design {t_j , j = 1, ..., 28,799}. We end up with N = 1000 unequally spaced observations ranging from t = 0.03819 days to t = 69.37847 days. The error terms are generated independently from a Gaussian distribution with zero mean and variance ${\sigma }_{z}^{2}=0.005$.

Table 1. Parameters of Simulated RR Lyrae Star

k	kf0	ak	φk
	(days−1)	(mag)	(°)
1	2	0.401	5.490
2	4	0.171	144.040
3	6	0.133	285.250
4	8	0.097	81.290

Note. Parameters (frequencies, amplitudes, and phases) obtained from Benkő et al. (2011), as explained in Section 6.2.

Download table as:  ASCIITypeset image

If we consider of our novel model in Equation (1), with time-varying trend and amplitudes specified as

$$\begin{eqnarray}\begin{array}{rcl}m({t}_{i}) & = & {a}_{0}\left[1+{U}_{m}({t}_{i})/{U}_{c}\right],\\ {g}_{1,k}({t}_{i}) & = & {a}_{k}\cos ({\varphi }_{k})[1+{U}_{m}({t}_{i})/{U}_{c}],\\ k & = & 1,\,\ldots ,\,4,\\ {g}_{2,k}({t}_{i}) & = & {a}_{k}\sin ({\varphi }_{k})[1+{U}_{m}({t}_{i})/{U}_{c}],\\ k & = & 1,\,\ldots ,\,4,\end{array}\end{eqnarray} \tag{ 26 }$$

we can rewrite the model in Equation (25) as a special case of our model given by Equation (1). The main advantage of fitting the model in Equation (1) instead of the model in Equation (25) is that one does not need to estimate the parameters a_m , h, a₀, φ_m , a_k , φ_k , k = 1, ..., 4, f_m . Moreover, we do not need to adopt any specific functional form for m( · ) and g( · ), such as those given by Equation (26), because they are well approximated by B-splines.

We fit the model in Equation (1) with K = 4 to the data generated according to the model in Equation (25). We assume that the frequencies f_k , k = 1, ..., 4, of each harmonic component are known; see Table 1. Also, we use a total of J = 18 B-splines of degree d = 3, an order penalty r = 1, and the smoothing parameters ${\boldsymbol{\tau }}={\left(\mathrm{5,}\,1,\,0.1,\,0.1,\,0.1,\,0.1\,,1,\,0.1,\,4\right)}^{\top }$.

We fit the model in Equation (1) to the simulated data obtained from the model in Equation (25), and present the results in Figure 3. The first row shows the simulation of the amplitude modulated RR Lyrae light curve given by Equation (25); (gray points), together with the true and fitted curve (solid red and solid black lines, respectively). The second row shows the residuals, and the third and fourth rows show the true trend and amplitudes (red lines) given by Equation (26) and their fits (black lines). We observe from Figure 3 that the model in Equation (1) fits well the simulated data. That is, the trend and amplitudes are well fitted, and the residuals satisfy the assumption of zero mean and constant variance. The 95% confidence intervals are constructed in a parametric fashion; see Appendix C.2 for more details.

Zoom In Zoom Out Reset image size

6.3. Estimating the Spectral Density of Unequally Spaced Time Series

In this section we estimate the spectral density of unequally spaced time series by means of our novel estimator in Equation (23). To this end, we simulate unequally spaced observations generated from the following AR(2) process:

$$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{i} & = & {\phi }_{1}{\varepsilon }_{i-1}+{\phi }_{2}{\varepsilon }_{i-2}+{z}_{i}\\ {z}_{i} & \sim & { \mathcal N }(0,{\sigma }_{z}^{2})\\ {t}_{i} & = & {t}_{0}+i{\rm{\Delta }},\end{array}\end{eqnarray} \tag{ 27 }$$

where i = 1, ..., N, with N = 500 equally spaced observations, starting time t₀ = 0.67, and Δ = 0.33. In order to simulate a realistic AR(2) process, we use the coefficients of the sunspot numbers in Example 3.2.9 of Brockwell & Davis (2016), where ϕ₁ = 1.318, ϕ₂ = −0.634, and ${\sigma }_{z}^{2}=289.2$. These ϕ coefficients ensure the existence of a causal solution

$$\begin{eqnarray}&&{\varepsilon }_{i}=\sum _{j=0}^{\infty }{\psi }_{j}{z}_{i-j}\end{eqnarray} \tag{ 28 }$$

of Equation (27). The time series in Equation (28) is causal in the sense that ε depends upon current and past (rather than future) values of the error term z. We simulate M = 500 times the AR(2) model given by Equation (27) obtaining the observations ${\varepsilon }_{1}^{\left(m\right)},\,\ldots ,\,{\varepsilon }_{N}^{\left(m\right)}$, m = 1, ..., M. Then, in order to obtain unequally spaced observations we use the following three steps.

• 1.  
  We divide time into 50 blocks, where each block has 10 observations, in a way to preserve the original time series structure.
• 2.  
  In order to preserve the autocorrelation between the observations, we select randomly 30 blocks and collect the time points corresponding to these blocks, obtaining a new set of time points $\{{t}_{i}^{* },i=1,\,\ldots ,\,n\}$, with n = 300 observations. In contrast to the simulation schemes of Sections 6.1 and 6.2 where time was sampled randomly, here data sets with uniformly sampled subsets are produced. While the former sampling is close to the data distribution of large ground-based surveys, the latter is the typical sampling of photometric space telescopes that are dedicated to high-cadence time series observations, such as Kepler (Koch et al. 2010).
• 3.  
  Finally, we collect the observations ${\varepsilon }_{i}^{\left(m\right)}$ corresponding to the new set of time points $\{{t}_{i}^{* },i=1,\,\ldots ,\,n\}$ and rename them as ${e}_{i}^{\left(m\right)}={\varepsilon }_{i}^{\left(m\right)}$, with ${e}_{i}^{\left(m\right)}$ being observed at time ${t}_{i}^{* }$, i = 1, ..., n.

Thus, we obtain the unequally spaced observations ${e}_{1}^{\left(m\right)},\,\ldots ,\,{e}_{n}^{\left(m\right)}$, which represent a subset of the equally spaced time series ${\varepsilon }_{1}^{\left(m\right)},\,\ldots ,\,{\varepsilon }_{N}^{\left(m\right)}$. For each m = 1, ..., M, and each fixed frequency λ_j = 2π f_j , f = j/(NΔ), j = 1, ..., N, we compute the periodogram of ${e}_{1}^{\left(m\right)},\,\ldots ,\,{e}_{n}^{\left(m\right)}$ as

$$\begin{eqnarray*}&&{I}_{\,e}^{\left(m\right)}({\lambda }_{j})=\sum _{k=1}^{n}\sum _{d=1}^{n}{e}_{k}^{\left(m\right)}{e}_{d}^{\left(m\right)}\exp (i{\lambda }_{j}[{t}_{k}^{* }-{t}_{d}^{* }]),\end{eqnarray*}$$

the average of the periodograms ${I}_{\,e}^{\left(m\right)}({\lambda }_{j})$ as

$$\begin{eqnarray*}&&{\overline{I}}_{\,e}({\lambda }_{j})=\displaystyle \frac{1}{M}\sum _{m=1}^{M}{I}_{\,e}^{\left(m\right)}({\lambda }_{j}),\end{eqnarray*}$$

and the power spectral window of ${e}_{1}^{\left(m\right)},\,\ldots ,\,{e}_{n}^{\left(m\right)}$ as

$$\begin{eqnarray*}&&{W}_{e}({\lambda }_{j})=\sum _{d=1}^{n}\sum _{k=1}^{n}\exp (i{\lambda }_{j}[{t}_{k}^{* }-{t}_{d}^{* }]).\end{eqnarray*}$$

For each frequency λ_j , j = 1, ..., N, replacing ${\rm{E}}\left[{I}_{\,e}({\lambda }_{j})\right]$ with ${\overline{I}}_{e}({\lambda }_{j})$ and substituting W_e (λ_j ) in Equation (23), the estimated spectral density of the unequally spaced time series ${e}_{1}^{\left(m\right)},\,\ldots ,\,{e}_{n}^{\left(m\right)}$ is given by

$$\begin{eqnarray}&&{\widehat{P}}_{e}({\lambda }_{j})=\displaystyle \frac{N}{2\pi }{{ \mathcal F }}^{-1}\left\{\displaystyle \frac{{ \mathcal F }\{{\overline{I}}_{e}({\lambda }_{j})\}[k]}{{ \mathcal F }\{{W}_{e}({\lambda }_{j})\}[k]}\right\}[j].\end{eqnarray} \tag{ 29 }$$

A smooth version of the estimated spectral density in Equation (29) is

$$\begin{eqnarray}&&{\widetilde{P}}_{e}({\lambda }_{j})=\tfrac{1}{{N}_{\lambda }}\sum _{i=1}^{{N}_{\lambda }}{K}_{h}({\lambda }_{j}-{\lambda }_{i}){\widehat{P}}_{e}({\lambda }_{j}).\end{eqnarray} \tag{ 30 }$$

The rescaled kernel function is defined as ${K}_{h}(x)=\tfrac{1}{h}K(x/h)$, where K is a second-order kernel and h is the bandwidth. For this application, we used the Gaussian kernel $K(y)=\tfrac{1}{\sqrt{2\pi }}\exp (-{y}^{2}/2)$ and a bandwidth h = 0.3.

Figure 4 compares the underlying spectral density P_ε (λ_j ) of the equally spaced time series {ε_i }, with the estimated spectral densities ${\widehat{P}}_{e}({\lambda }_{j})$ and ${\widetilde{P}}_{e}({\lambda }_{j})$ of the unequally spaced time series {e_i }. The underlying spectral density of the equally spaced time series {ε_i }, generated by the AR(2) process in Equation (27), is given by

$$\begin{eqnarray}&&\begin{array}{l}{P}_{\varepsilon }({\lambda }_{j})=\displaystyle \frac{{\sigma }_{z}^{2}}{2\pi }\left[1+{\phi }_{1}^{2}+{\phi }_{2}^{2}+2{\phi }_{2}\right.\\ \ +\,{\left.2({\phi }_{1}{\phi }_{2}-{\phi }_{1})\cos ({\lambda }_{j}{\rm{\Delta }})-4{\phi }_{2}{\cos }^{2}({\lambda }_{j}{\rm{\Delta }})\right]}^{-1}.\end{array}\end{eqnarray} \tag{ 31 }$$

The estimated spectral density ${\widehat{P}}_{e}({\lambda }_{j})$ of the unequally time series {e_i } is given in Equation (29), and its smooth version ${\tilde{P}}_{e}({\lambda }_{j})$ in Equation (30). Figure 4 shows that the estimated spectral density of the unequally time series, ${\widetilde{P}}_{e}({\lambda }_{j})$, fits very well the true spectral density P_ε (λ_j ).

Zoom In Zoom Out Reset image size

In this section, we fit our model in Equation (1) and the model proposed by Benkő (2018) in Equation (5) to the same light curve: the V783 Cyg, KIC 5559631. This time series has 61,351 unequally spaced observations and is available online from the Konkoly Observatory of the Hungarian Academy of Sciences webpage. ⁷ We choose this particular light curve for two reasons. First, the Blazhko effect of the V783 Cyg time series is known to be characterized by a sinusoidal amplitude and frequency modulation (Benkő et al. 2014). The light curve of V783 Cyg can be described by K = 15 significant harmonics with a sinusoidal amplitude and frequency modulations (Benkő et al. 2014), which makes V783 Cyg an ideal target for comparing the fits obtained with the models in Equations (1) and (5). Second, these two modulations are well captured and fitted by our novel model in Equation (1).

In order to reduce the computational time and to satisfy the condition t_i = t₀ + iΔ with Δ > 0 and $i\in { \mathcal I }\subseteq {\mathbb{N}}$ (which is required by Proposition 1), we analyze a ≈78 days segment of this light curve from t = 827.44 days to t = 904.9 days. For this segment, the time origin and the time spacing take the values t₀ = 827.42 days and Δ = 0.0204345 days, respectively, with a total of N = 2101 unequally spaced observations.

When fitting the models in Equations (1) and (5), the main pulsation and modulation frequencies are not estimated: they take the values f₀ = 1.611084 days⁻¹ and f_m = 0.036058 days⁻¹ (see Benkő et al. 2014), respectively. In addition to the K = 15 significant harmonics fitted by Benkő et al. (2014), we found, after pre-whitening and fitting our model in Equation (1), four significant frequencies taking the values ${f}_{11}^{{\prime} }=18.3254$ days⁻¹, ${f}_{12}^{{\prime} }=19.9365$ days⁻¹, ${f}_{13}^{{\prime} }=21.5476$ days⁻¹, and ${f}_{14}^{{\prime} }\,=23.1587$ days⁻¹. The values we obtain for $\{{f}_{j}^{{\prime} },11\leqslant j\leqslant 14\}$ demonstrate that these frequencies are not harmonics of the form k f₀, which might suggest that these four are independent frequencies. Interestingly, however, we find that the latter belong to a set of fourteen "reflection frequencies" of the form $\{{f}_{j}^{{\prime} }=2{f}_{N}-(30-j){f}_{0},1\leqslant j\leqslant 14\}$, where f_N = 24.46 days⁻¹ is the Nyquist frequency. Among these fourteen frequencies, only the last six $\{{f}_{j}^{{\prime} },9\leqslant j\leqslant 14\}$ exhibit significant peaks in the Lomb−Scargle periodogram (computed according to Lomb 1976). However, to avoid overfitting, we only consider the four frequencies $\{{f}_{j}^{{\prime} },11\leqslant j\leqslant 14\}$ corresponding to last four peaks of the estimated power spectrum (see the last row of Figure 5, bottom right panel). In summary, the only truly independent frequencies are f₀, f_m , and f_N , the other frequencies {f_k = kf₀, k = 1, ..., 15} and $\{{f}_{j}^{{\prime} }=2{f}_{N}-(30-j){f}_{0},1\leqslant j\leqslant 14\}$ being linear combinations (or harmonics) of those.

Zoom In Zoom Out Reset image size

The frequencies ${f}_{j}^{{\prime} }$ do not depend on the Blazkho frequency f_m , as the information regarding the Blazhko effect is captured by the time-varying trend m( · ) and amplitudes {g_ℓ,k ( · ), ℓ = 1, 2, k = 1, ..., K}; see Equations (9)−(10) and Figure 6. We use a different notation ($f^{\prime}$ rather than f) to avoid confusion, since in this case f_N < 16f₀ = 25.77 days⁻¹, and the four frequencies we are considering take values < 24 days⁻¹.

Zoom In Zoom Out Reset image size

After fitting the models in Equations (1) and (5), we compute their residuals and estimate their spectral densities using Equation (24). To estimate the spectral densities according to the procedure in Section 5, we adopt the Gaussian kernel $K(y)=\tfrac{1}{\sqrt{2\pi }}\exp (-{y}^{2}/2)$ with a bandwidth h = 7.2. We fitted both models with a PC having a 2.7 GHz 12 core Intel Xeon E5 processor and 64 GB of 1866 MHz DDR3 memory. Fitting our novel model in Equation (1) required 17 minutes and 13 s, whereas fitting the model by Benkő (2018) in Equation (5) required 12 minutes and 28 s.

The description provided so far applies to both fits of models in Equations (1) and (5). We now provide, separately, computational details about the estimation of these two models. Then in Sections 7.1 and 7.2 we compare and interpret the fits.

To fit our novel model in Equation (1), we apply the methodology described in Section 4. When fitting our model in Equation (1) we consider two sets of harmonic components. The first set is given by the harmonic components with frequencies {f_k = kf₀, k = 1, ..., 15} provided by Benkő et al. (2014), weighted by our amplitudes {g_ℓ,k (t_i ), ℓ = 1, 2, k = 1, ..., 15}. For the second set, the harmonic components are characterized by the four amplitudes $\{{g}_{{\ell },j}^{{\prime} }({t}_{i}),{\ell }=1,2,j=11,\,\ldots ,\,14\}$ weighting the corresponding four frequencies $\{{f}_{j}^{{\prime} },j=11,\,\ldots ,\,14\}$. That is, we fit the following extended version:

$$\begin{eqnarray*}\begin{array}{rcl}\mu ({t}_{i}) & = & m({t}_{i})+\sum _{k=1}^{15}\{{g}_{1,k}({t}_{i})\cos ({w}_{k}{t}_{i})\\ & & +\,{g}_{2,k}({t}_{i})\sin ({w}_{k}{t}_{i})\}\\ & & +\,\sum _{j=11}^{14}\{{g}_{1,j}^{{\prime} }({t}_{i})\cos ({w}_{j}^{{\prime} }{t}_{i})+{g}_{2,j}^{{\prime} }({t}_{i})\sin ({w}_{j}^{{\prime} }{t}_{i})\}\end{array}\end{eqnarray*}$$

of the model in Equation (1), with ω_k = 2π f_k and ${\omega }_{j}^{{\prime} }=2\pi {f}_{j}^{{\prime} }$. The resulting fitted model involves a total of 39 J parameters. Before fitting our model, we select the smoothing parameters τ and the number of B-splines J. These parameters are selected by the AIC criterion described in Section 4.2. To simplify the selection of the smoothing parameters τ , we consider the case τ₂ = ... = τ₁₁, τ₁₂ = ... = τ₂₁, τ₂₂ = ... = τ₃₁, and τ₃₂ = ... = τ₃₉. We pick the smoothing parameter τ₁ over the grid {0, 0.1, 10}, the parameters {τ_k , k = 2, ..., 39} over the grid {0, 0.1, 10, 100}, and the total number J of B-splines (of degree d = 3) over the grid {8, 13, 23, 33}. We apply the AIC formula in Equation (14). The lowest AIC value occurs for J = 33 B-splines, τ₁ = 0, τ_k = 0.1, k = 2, ..., 21, τ_k = 0, k = 22, ..., 31, and τ_k = 10, k = 32, ..., 39.

To fit the model in Equation (5), we implement the Levenberg−Marquardt algorithm using the R function nls.lm (available in the R package minpack.lm by Elzhov et al. 2016), with K = 15 and ${\ell }={{\ell }}_{k}^{A}={{\ell }}_{k}^{F}=1$, k = 1, ..., K, for a total of 93 parameters.

7.1. Comparing the Accuracy of the Fits

7.2. Comparing the Estimated Time-varying Parameters

In this article, we introduced a model for time series observations of variable stars that are modulated by smoothly time-varying mean magnitudes, amplitudes, and phases. Previous approaches assume that the underlying parameters are either time-invariant or piecewise-constant functions. From the modeling viewpoint, our approach is more flexible because it avoids assumptions about the functional form of the aforementioned time-dependent quantities. From the computational viewpoint, estimating our time-varying curves translates into the estimation of time-invariant parameters that can be performed by ordinary least squares.

An important challenge when dealing with astronomical time series is that observations are unequally spaced in time. In some cases, observations are unevenly spaced due to missing values. Missing values are sometimes handled via imputation, that is, the gap generated by the missing value is "filled in" by an estimated value. Our novel approach, which involves the classical periodogram, has the advantage of not relying on any imputation method.

We study the performance of our approach under several simulation scenarios. Finally, we apply our method to V783 Cyg (KIC 5559631), a well-known RR Lyrae star presenting the Blazhko effect. In this case, the effect is characterized by a sinusoidal amplitude and frequency modulation. When comparing the time-varying fit obtained with our novel model with the time-invariant fit obtained with the model proposed by Benkő (2018), we found that both amplitude and frequency modulations are well captured and fitted by our novel model, and also that our time-varying method outperforms the time-invariant fit. Indeed the estimation error obtained with our fit is significantly smaller than the error obtained with the time-invariant fit. In addition, the residuals obtained with our novel method are compatible with the assumption of stationary and uncorrelated errors, whereas the residuals obtained with the time-invariant model by Benkő (2018) exhibit a time-dependent trend and some significant spectral peaks.

In the future, we plan to extend our methodology in four important directions. First, we plan to apply our novel method to the study of a larger sample of Blazhko RR Lyrae stars. Second, our approach can be extended to the analysis of other classes of variable stars presenting long-term changes in their light-curve shapes. Third, our fitting method does not require the period(s), amplitude(s), and phase(s) of the Blazhko effect to be determined, as we obtain instead the empirical functions m( · ) and g_i,k ( · ). We are currently investigating what kind of (or how much more) information can be obtained from these empirical functions, as compared to conventional approaches. Finally, we aim to study Blazhko light curves characterized by more than one Blazhko frequency—V783 Cyg, which was addressed in some detail in this paper, is a special case, because this star does not show any additional Blazhko frequencies (Benkő et al. 2014).

The authors would like to thank two anonymous reviewers for helpful comments which led to significant improvements of the paper. D.S. was funded by the National Agency for Research and Development (ANID), Doctorado Nacional grant No. 2017-21171100. Support for G.M. and M.C. has been provided by ANID's Millennium Science Initiative through grant No. ICN12 120009, awarded to the Millennium Institute of Astrophysics (MAS). M.C. acknowledges additional support by Proyecto Basal AFB-170002 and FONDECYT grant #1171273. We thank all of the participants of the Astronomical Data Science Workshop organized by Texas A&M University on 2020 February 17–18, as well as the participants of the IISA 2021 conference organized by the University of Illinois Chicago on 2021 May 20–23. We thank the Astrostatistics group at the Center for Astrophysics of Harvard University for constructive criticism of this manuscript. Special thanks go to József Benkő for sharing the parameters we used to simulate the Blazhko star of Section 6.2, and to Gergely Hajdu for useful discussions.

Software: The Language R for Statistical Computing (R Core Team 2021), the R package NISTunits (Gama 2016), and the R package minpack.lm (Elzhov et al. 2016).

The Blazhko effect is a periodic amplitude and phase variation in the light curves of RR Lyrae variable stars. In astronomy, the Blazhko effect is usually interpreted as a modulation phenomenon. Modulation is the process of transmitting a low-frequency signal into a high-frequency wave, called the carrier wave, by changing its amplitude, frequency, and/or phase angle through the modulating signal. The function of the carrier wave is to carry the message or modulating signal from the transmitter to the receiver. The superposition of the signal and the carrier wave results in the so-called modulated signal.

In this Appendix we review two types of modulation, as given in Benkő et al. (2011): amplitude modulation and frequency modulation. This will be helpful for a comparison between our model (Equation (1)) and the models proposed by Benkő et al. (2011) and Benkő (2018), in the case of RR Lyrae stars presenting the Blazhko effect.

A.1. Amplitude Modulation

Amplitude modulation (AM) changes the amplitude of the carrier signal. Let the carrier wave c(t) be a sinusoidal signal of the form

$$\begin{eqnarray*}&&c(t)={U}_{c}\sin (2\pi {f}_{c}t+{\phi }_{c}),\end{eqnarray*}$$

where the constant parameters U_c , f_c , and ϕ_c are the amplitude, frequency, and phase of the carrier wave, respectively.

Let ${U}_{m}\left(t\right)$ represent a waveform that is the message to be transmitted, or the modulating signal. The transmitter uses the information signal ${U}_{m}\left(t\right)$ to vary the amplitude of the carrier U_c to produce the amplitude modulated signal U_AM:

$$\begin{eqnarray}\begin{array}{rcl}{U}_{\mathrm{AM}}(t) & = & [{U}_{c}+{U}_{m}(t)]\sin (2\pi {f}_{c}t+{\phi }_{c})\\ & = & [{U}_{c}+{U}_{m}(t)]\tfrac{c(t)}{{U}_{c}}=\left[1+\tfrac{{U}_{m}}{{U}_{c}}\right]c(t).\end{array}\end{eqnarray} \tag{ A1 }$$

In the simplest case, when the modulating signal is sinusoidal, that is,

$$\begin{eqnarray}&&{U}_{m}(t)={U}_{m}^{A}\sin (2\pi {f}_{m}t+{\phi }_{m}^{A}),\end{eqnarray} \tag{ A2 }$$

the amplitude modulated signal in Equation (A1) is

$$\begin{eqnarray}\begin{array}{rcl}{U}_{\mathrm{AM}}(t) & = & \left[{U}_{c}+{U}_{m}^{A}\sin (2\pi {f}_{m}+{\phi }_{m}^{A})\right]\\ & & \times \,\sin (2\pi {f}_{c}t+{\phi }_{c}).\end{array}\end{eqnarray} \tag{ A3 }$$

Clearly, a more complex example of amplitude modulation arises when K ≥ 1, where K denotes the number of harmonic components. Suppose the carrier wave c(t) is a linear combination of sine harmonics:

$$\begin{eqnarray*}&&c(t)={a}_{0}+\sum _{k=1}^{K}{a}_{k}\sin (2\pi {{kf}}_{0}t+{\phi }_{k}),\end{eqnarray*}$$

and the modulating signal is sinusoidal and given again by Equation (A2). Following the same idea as in Equation (A1), the amplitude modulated signal in Equation (A3) is

$$\begin{eqnarray}\begin{array}{rcl}{U}_{\mathrm{AM}}(t) & = & \left[1+\displaystyle \frac{{U}_{m}(t)}{{U}_{c}}\right]\,c(t)\\ & = & \left[1+\displaystyle \frac{{U}_{m}^{A}\sin (2\pi {f}_{m}t+{\phi }_{m}^{A})}{{U}_{c}}\right]\\ & & \times \,\left[{a}_{0}+\sum _{k=1}^{K}{a}_{k}\sin (2\pi {{kf}}_{0}t+{\phi }_{k})\right].\end{array}\end{eqnarray} \tag{ A4 }$$

If we call $h={U}_{m}^{A}/{U}_{c}$ and use the basic trigonometrical identities $\sin (a)\sin (b)=\tfrac{1}{2}[\cos (a-b)-\cos (a+b)]$ and $\sin (a)=\cos \left(a-\tfrac{\pi }{2}\right)$, Equation (A4) can be written as

$$\begin{eqnarray}\begin{array}{rcl}{U}_{\mathrm{AM}}(t) & = & {a}_{0}+\sum _{k=1}^{K}{a}_{k}\sin (2\pi {{kf}}_{0}t+{\phi }_{k})\\ & & +\,{a}_{0}h\sin (2\pi {f}_{m}t+{\phi }_{m}^{A})\\ & & +\,\sum _{k=1}^{K}\displaystyle \frac{{a}_{k}h}{2}\sin (2\pi ({{kf}}_{0}-{f}_{m})t+({\phi }_{k}-{\phi }_{m})+\pi /2)\\ & & -\,\sum _{k=1}^{K}\displaystyle \frac{{a}_{k}h}{2}\sin (2\pi ({{kf}}_{0}+{f}_{m})t+({\phi }_{k}+{\phi }_{m})+\pi /2).\end{array}\end{eqnarray} \tag{ A5 }$$

This example shows that when the time-varying amplitude ${U}_{m}\left(t\right)$ in Equation (A2) takes a sinusoidal form, the amplitude modulated model with time-varying amplitude in Equation (A4) can be written as a model with time-invariant parameters as in Equation (A5). This implies that, when frequencies and phases are known, the parameters {a_k , 0 ≤ k ≤ K} in Equation (A5) can be estimated by ordinary least squares.

A.2. Amplitude and Frequency Modulation

Frequency modulation (FM) changes the frequency of the carrier signal. We assume the sinusoidal carrier wave to be

$$\begin{eqnarray*}&&c(t)={U}_{c}\sin ({\rm{\Theta }}(t)),\end{eqnarray*}$$

where Θ(t) = 2π f_c t + ϕ_c is the angular part of the function. Suppose that the modulating signal is ${U}_{m}\left(t\right)$. Then the modulated angular part is given by

$$\begin{eqnarray*}&&{\rm{\Theta }}(t)=2\pi {f}_{c}t+2\pi {k}_{\mathrm{FM}}{\int }_{0}^{t}{U}_{m}^{F}(\tau ){\rm{d}}\tau +{\phi }_{c},\end{eqnarray*}$$

where k_FM is the frequency deviation, and the frequency modulated signal is expressed as

$$\begin{eqnarray}&&{U}_{\mathrm{FM}}(t)={U}_{c}\sin \left(2\pi {f}_{c}t+2\pi {k}_{\mathrm{FM}}{\int }_{0}^{t}{U}_{m}^{F}(\tau ){\rm{d}}\tau +{\phi }_{c}\right).\end{eqnarray} \tag{ A6 }$$

In the simplest case, when the modulating signal is represented by a sinusoidal wave with amplitude ${U}_{m}^{F}$ and frequency f_m , the integral of such a signal is

$$\begin{eqnarray*}&&{\int }_{0}^{t}{U}_{m}(\tau ){\rm{d}}\tau =\displaystyle \frac{{U}_{m}^{F}}{2\pi {f}_{m}}\sin (2\pi {f}_{m}t+{\phi }_{m}),\end{eqnarray*}$$

and the frequency modulated signal in Equation (A6) is

$$\begin{eqnarray}&&{U}_{\mathrm{FM}}(t)={U}_{c}\sin \left(2\pi {f}_{c}t+\displaystyle \frac{{k}_{\mathrm{FM}}{U}_{m}^{F}}{{f}_{m}}\sin (2\pi {f}_{m}t+{\phi }_{m})+{\phi }_{c}\right).\end{eqnarray} \tag{ A7 }$$

In practice, modulated signals can be a mixture of amplitude and frequency modulations, which can be used to described Blazhko RR Lyrae stars (Benkő et al. 2011). We review the simplest case when both AM and FM are sinusoidal. Combining the amplitude modulated signal in Equation (A3) and the frequency modulated signal in Equation (A7), the amplitude and frequency modulated signal is thus

$$\begin{eqnarray*}\begin{array}{rcl}{U}_{\mathrm{Comb}}(t) & = & [{U}_{c}+{U}_{m}^{A}\sin (2\pi {f}_{m}t+{\phi }_{m})]\\ & & \times \,\sin \left(2\pi {f}_{c}t+\displaystyle \frac{{k}_{\mathrm{FM}}{U}_{m}^{F}}{{f}_{m}}\sin (2\pi {f}_{m}t+{\phi }_{m})+{\phi }_{c}\right).\end{array}\end{eqnarray*}$$

A.3. Blazhko Modulation

Amplitude and frequency modulations have been observed in Blazhko RR Lyrae stars (e.g., Benkő et al. 2010; Chadid et al. 2010; Poretti et al. 2010; Sódor et al. 2012). Assuming that the observed data sets are precise and long enough, Benkő et al. (2011) proposed an amplitude and frequency modulation model for Blazhko RR Lyrae stars given by

$$\begin{eqnarray}&&{m}_{\mathrm{Comb}}^{* }(t)=\displaystyle \frac{{m}_{\mathrm{AM}}^{* }(t)}{{c}^{* }(t)}{m}_{\mathrm{FM}}^{* }(t),\end{eqnarray} \tag{ A8 }$$

where ${c}^{* }\left(t\right)$ is the carrier wave, and the functions ${m}_{\mathrm{AM}}^{* }\left(t\right)$ and ${m}_{\mathrm{FM}}^{* }\left(t\right)$ are the nonsinusoidal amplitude and frequency modulations, given respectively by

$$\begin{eqnarray}&&{m}_{\mathrm{AM}}^{* }(t)=\left[{a}_{0}^{A}+\sum _{p=1}^{q}{a}_{p}^{A}\sin (2\pi {{pf}}_{m}t+{\varphi }_{p}^{A})\right]{c}^{* }(t),\end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray}\begin{array}{rcl}{m}_{\mathrm{FM}}^{* }(t) & = & {a}_{0}+\sum _{k=1}^{K}{a}_{k}\sin \left[2\pi {{kf}}_{0}t+{{ka}}_{0}^{F}\right.\\ & & \left.+\,k\sum _{p=1}^{q}{a}_{p}^{F}\sin (2\pi {{pf}}_{m}t+{\varphi }_{p}^{F})+{\varphi }_{k}\right].\end{array}\end{eqnarray} \tag{ A10 }$$

Here, the modulating signal used in the amplitude and frequency modulation is an arbitrary periodic signal represented by a Fourier sum with a constant frequency f_m . Superscripts A and F denote the amplitude modulation and frequency modulation parameters, respectively, and f₀ and f_m are the main pulsation and modulation frequencies, respectively.

Substituting Equations (A9) and (A10) in Equation (A8), we finally have

$$\begin{eqnarray*}\begin{array}{rcl}{m}_{\mathrm{Comb}}^{* }(t) & = & \left[{a}_{0}^{A}+\sum _{p=1}^{q}{a}_{p}^{A}\sin (2\pi {{pf}}_{m}t+{\varphi }_{p}^{A})\right]\\ & & \times \,\left\{{a}_{0}+\sum _{k=1}^{K}{a}_{k}\sin \left[2\pi {{kf}}_{0}t+{{ka}}_{0}^{F}+k\sum _{p=1}^{q}\right.\right.\\ & & \times \,\left.\left.{a}_{p}^{F}\sin (2\pi {{pf}}_{m}t+{\varphi }_{p}^{F})+{\varphi }_{k}\right]\right\}.\end{array}\end{eqnarray*}$$

In this Appendix we define B-splines and give some details about the estimation method that we used in this manuscript. For more details we refer the reader to the book by de Boor (1978).

A B-spline curve f(t) of degree d is defined as

$$\begin{eqnarray}&&f(t)=\sum _{j=1}^{J}{P}_{j}{B}_{j,d}(t),\end{eqnarray} \tag{ B1 }$$

where P_j are the control points and ${B}_{j,d}\left(t\right)$ are the B-spline basis functions. Let ${t}_{\min }$ and ${t}_{\max }$ be, respectively, the lower and upper bounds of the domain of interest. In order to build the B-spline basis of degree d, we first divide the domain into n intervals, with n being a positive integer, obtaining the n + 1 knots ξ_d , ξ_d+1, ..., ξ_d+n . Each knot satisfies ξ_j < ξ_j+1, for all j. Second, we define 2d additional knots ξ ₀, ξ ₁, ..., ξ_d−1, ξ_n+d+1, ..., ξ_n+2d−1, ξ_n+2d . Then, the jth B-spline basis, ${B}_{j,d}\left(t\right)$, can be defined recursively as

$$\begin{eqnarray}\begin{array}{rcl}{B}_{j,d}(t) & = & \displaystyle \frac{t-{\xi }_{\,j-1}}{{\xi }_{\,j+d-1}-{\xi }_{\,j-1}}{B}_{j,d-1}(t)\\ & & +\,\displaystyle \frac{{\xi }_{\,j+d}-t}{{\xi }_{\,j+d}-{\xi }_{\,j}}{B}_{j+1,d-1}(t),\quad j=1,\,\ldots ,\,J,\end{array}\end{eqnarray} \tag{ B2 }$$

with

$$\begin{eqnarray}{B}_{j,0}(t)=\left\{\begin{array}{ll}1 & t\in [{\xi }_{\,j-1},{\xi }_{\,j}),\\ 0 & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ B3 }$$

being used to initialize the recursion. Thus, to build the B-spline curve given by Equation (B1), we need n + 2d + 1 knots, and the total number of B-splines basis functions is J = n + d.

To illustrate how to construct a B-spline basis, consider the case of degree d = 2 and assume that the domain $[{t}_{\min },{t}_{\max }]$ has been divided into n = 3 intervals, obtaining the knots ξ ₂, ..., ξ ₅. In this instance, the 2d = 4 additional knots are defined as ξ ₀, ξ ₁, ξ ₆, and ξ ₇. Using Equation (B2), we obtain

$$\begin{eqnarray*}\begin{array}{rcl}{B}_{\mathrm{1,2}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,0}}{{\xi }_{\,2}-{\xi }_{\,0}}{B}_{\mathrm{1,1}}\left(t\right)+\displaystyle \frac{{\xi }_{\,3}-t}{{\xi }_{\,3}-{\xi }_{\,1}}{B}_{\mathrm{2,1}}(t),\\ {B}_{\mathrm{2,2}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,1}}{{\xi }_{\,3}-{\xi }_{\,1}}{B}_{\mathrm{2,1}}\left(t\right)+\displaystyle \frac{{\xi }_{\,4}-t}{{\xi }_{\,4}-{\xi }_{\,2}}{B}_{\mathrm{3,1}}(t),\\ {B}_{\mathrm{3,2}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,2}}{{\xi }_{\,4}-{\xi }_{\,2}}{B}_{\mathrm{3,1}}\left(t\right)+\displaystyle \frac{{\xi }_{\,5}-t}{{\xi }_{\,5}-{\xi }_{\,3}}{B}_{\mathrm{4,1}}(t),\\ {B}_{\mathrm{4,2}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,3}}{{\xi }_{\,5}-{\xi }_{\,3}}{B}_{\mathrm{4,1}}\left(t\right)+\displaystyle \frac{{\xi }_{\,6}-t}{{\xi }_{\,6}-{\xi }_{\,4}}{B}_{\mathrm{5,1}}(t),\\ {B}_{\mathrm{5,2}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,4}}{{\xi }_{\,6}-{\xi }_{\,4}}{B}_{\mathrm{5,1}}\left(t\right)+\displaystyle \frac{{\xi }_{\,7}-t}{{\xi }_{\,7}-{\xi }_{\,5}}{B}_{\mathrm{6,1}}(t),\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{B}_{\mathrm{1,1}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,0}}{{\xi }_{\,1}-{\xi }_{\,0}}{B}_{\mathrm{1,0}}\left(t\right)+\displaystyle \frac{{\xi }_{\,2}-t}{{\xi }_{\,2}-{\xi }_{\,1}}{B}_{\mathrm{2,0}}(t),\\ {B}_{\mathrm{2,1}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,1}}{{\xi }_{\,2}-{\xi }_{\,1}}{B}_{\mathrm{2,0}}\left(t\right)+\displaystyle \frac{{\xi }_{\,3}-t}{{\xi }_{\,3}-{\xi }_{\,2}}{B}_{\mathrm{3,0}}(t),\\ {B}_{\mathrm{3,1}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,2}}{{\xi }_{\,3}-{\xi }_{\,2}}{B}_{\mathrm{3,0}}\left(t\right)+\displaystyle \frac{{\xi }_{\,4}-t}{{\xi }_{\,4}-{\xi }_{\,3}}{B}_{\mathrm{4,0}}(t),\\ {B}_{\mathrm{4,1}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,3}}{{\xi }_{\,4}-{\xi }_{\,3}}{B}_{\mathrm{4,0}}\left(t\right)+\displaystyle \frac{{\xi }_{\,5}-t}{{\xi }_{\,5}-{\xi }_{\,4}}{B}_{\mathrm{5,0}}(t),\\ {B}_{\mathrm{5,1}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,4}}{{\xi }_{\,5}-{\xi }_{\,4}}{B}_{\mathrm{5,0}}\left(t\right)+\displaystyle \frac{{\xi }_{\,6}-t}{{\xi }_{\,6}-{\xi }_{\,5}}{B}_{\mathrm{6,0}}(t),\\ {B}_{\mathrm{6,1}}\left(t\right) & = & \displaystyle \frac{t-{\xi }_{\,5}}{{\xi }_{\,6}-{\xi }_{\,5}}{B}_{\mathrm{6,0}}\left(t\right)+\displaystyle \frac{{\xi }_{\,7}-t}{{\xi }_{\,7}-{\xi }_{\,6}}{B}_{\mathrm{7,0}}(t),\end{array}\end{eqnarray*}$$

and the coefficients $\{{B}_{j,0}\left(t\right)$, j = 1, ..., 7} are defined in Equation (B3).

Suppose we have N observations {t₁, ..., t_N }, that might be either equally or unequally spaced, with ${t}_{i}\in [{t}_{\min },{t}_{\max }]$ for all i = 1, ..., N. The B-splines basis matrix evaluated at time {t₁, ..., t_N }, denoted by B, is the N × J matrix with entries {B_j,d (t_i ), i = 1, ..., N, j = 1, ..., J}, in a way that each row contains a B-spline basis. The jth B-spline basis function satisfies the following properties:

$$\begin{eqnarray*}\left\{\begin{array}{ll}{B}_{j,d}(t)\gt 0 & {\xi }_{\,j-1}\lt t\lt {\xi }_{\,j+d},\\ {B}_{j,d}(t)=0 & {\xi }_{\,0}\leqslant t\leqslant {\xi }_{\,j-1}\,{\rm{and}}\,{\xi }_{\,j+d}\leqslant t\leqslant {\xi }_{\,n+2d},\\ \sum _{j=1}^{J}{B}_{j,d}(t)=1 & {\xi }_{\,d}\lt t\lt {\xi }_{\,n+d},\\ \tfrac{{\partial }^{k}{B}_{j,d}(t)}{\partial {t}^{k}}{| }_{t={\xi }_{{\ell }}} & 0\leqslant k\leqslant d-1\,{\rm{are}}\,{\rm{continuous}}.\end{array}\right.\end{eqnarray*}$$

For ease of notation we use, throughout our manuscript, ${B}_{j}\left(t\right)$ instead of ${B}_{j,d}\left(t\right)$. Let us now consider the example of estimating the mean function $\mu \left(t\right)$ of model Y_i = μ(t_i ) + ε_i using B-splines. Let ${\boldsymbol{Y}}=({Y}_{1},\,\ldots ,\,{Y}_{6})^{\prime}$ be the available N = 6 responses observed, respectively, at time {t₁, ..., t₆}, with ${t}_{\min }={t}_{1}$ and ${t}_{\max }={t}_{6}$. Then assume that $\mu (t)={\sum }_{j=1}^{J}{P}_{j}{B}_{j}\left(t\right)$, for all t ∈ [t₁, t₆]. We use here B-splines basis functions of degree d = 2; in order to construct them, we divide the domain [t₁, t₆] into n = 3 intervals. Hence, the total number of knots ξ ₀, ..., ξ ₇ is n + 2d + 1 = 8, and the total number of B-splines basis functions is J = n + d = 5. The 6 × 5 design matrix B has entries B_ij = B_j (t_i ),with i = 1, ..., 6 and j = 1, ..., 5, which permits estimating the coefficients {P_j , j = 1, ..., 5} by ordinary least squares. Indeed, if ${\boldsymbol{Y}}={\left({Y}_{1},\,\ldots ,\,{Y}_{6}\right)}^{\top }$ denotes the response vector and ${\boldsymbol{\theta }}={\left({P}_{1},\,\ldots ,\,{P}_{5}\right)}^{\top }$ the parameter vector, we can rewrite the model as Y = B θ + z , where ${\boldsymbol{z}}={\left({z}_{1},\,\ldots ,\,{z}_{{}_{6}}\right)}^{\top }$ is the error vector. The estimated parameters are defined as $\hat{{\boldsymbol{\theta }}}\,={\left({\hat{P}}_{1},\,\ldots ,\,{\hat{P}}_{{}_{5}}\right)}^{\top }={({{\bf{B}}}^{\top }{\bf{B}})}^{-1}{{\bf{B}}}^{\top }{\bf{Y}}$, and the estimated mean as $\widehat{\mu }(t)={\sum }_{j=1}^{5}{\widehat{P}}_{j}{B}_{j}\left(t\right)$, for all t ∈ [t₁, t₆].

C.1. Nonparametric Quantiles

C.2. Parametric Quantiles

We use the parametric quantiles to construct the confidence intervals for our simulation in Section 6.2 and our application in Section 7. Assuming that the error terms {z_i , i = 1, ..., N} follow a Gaussian distribution with zero mean and variance ${\sigma }_{z}^{2}$, the (1 − α) × 100% prediction interval for ${\rm{E}}\left[{\widehat{Y}}_{i}\right]={{\boldsymbol{ \mathcal B }}}^{\top }({t}_{i}){\rm{E}}\left[\widehat{{\boldsymbol{\theta }}}\right]$, with i = 1, ..., N, is

$$\begin{eqnarray*}&&{{\boldsymbol{ \mathcal B }}}^{\top }({t}_{i})\widehat{{\boldsymbol{\theta }}}\pm z(1-\alpha /2)\sqrt{{{\boldsymbol{ \mathcal B }}}^{\top }({t}_{i})\mathrm{Var}\left[\widehat{{\boldsymbol{\theta }}}\right]{\boldsymbol{ \mathcal B }}({t}_{i})},\end{eqnarray*}$$

where z(1 − α/2) denotes the (1 − α/2) quantile of the standard Gaussian distribution, and

$$\begin{eqnarray*}&&\mathrm{Var}\left[\hat{{\boldsymbol{\theta }}}\right]={\sigma }_{z}^{2}{\left({{\boldsymbol{ \mathcal B }}}^{\top }{\boldsymbol{ \mathcal B }}+{\bf{P}}\right)}^{-1}{{\boldsymbol{ \mathcal B }}}^{\top }{\boldsymbol{ \mathcal B }}{\left({{\boldsymbol{ \mathcal B }}}^{\top }{\boldsymbol{ \mathcal B }}+{\bf{P}}\right)}^{-1}.\end{eqnarray*}$$

The (1 − α) × 100% confidence interval for the trend ${\rm{E}}\left[\hat{m}({t}_{i})\right]={\mathscr{B}}{\left({t}_{i}\right)}^{{\rm{\top }}}{{\boldsymbol{Q}}}_{m}{\rm{E}}\left[\hat{{\boldsymbol{\theta }}}\right]$ is

$$\begin{eqnarray*}&&\begin{array}{c}\begin{array}{c}{\mathscr{B}}{\left({t}_{i}\right)}^{{\rm{\top }}}{{\boldsymbol{Q}}}_{m}\hat{{\boldsymbol{\theta }}}\pm z(1-\alpha /2)\\ \,\,\,\times \,\sqrt{{\mathscr{B}}{\left({t}_{i}\right)}^{{\rm{\top }}}{{\boldsymbol{Q}}}_{m}{\rm{V}}{\rm{a}}{\rm{r}}\left[\hat{{\boldsymbol{\theta }}}\right]{{\boldsymbol{Q}}}_{m}^{{\rm{\top }}}{\mathscr{B}}({t}_{i})},\end{array}\end{array}\end{eqnarray*}$$

and the (1 − α) × 100% confidence intervals for the amplitudes ${\rm{E}}\left[{\hat{g}}_{\ell ,k}({t}_{i})\right]={\mathscr{B}}{\left({t}_{i}\right)}^{{\rm{\top }}}{{\boldsymbol{Q}}}_{g}(\ell ,k){\rm{E}}\left[\hat{{\boldsymbol{\theta }}}\right]$, ℓ = 1, 2, k = 1, ..., K, are

$$\begin{eqnarray*}&&\begin{array}{c}\begin{array}{l}{\mathscr{B}}{\left({t}_{i}\right)}^{{\rm{\top }}}{{\boldsymbol{Q}}}_{g}(\ell ,k)\hat{{\boldsymbol{\theta }}}\pm z(1-\alpha /2)\\ \,\,\,\times \,\sqrt{{\mathscr{B}}{\left({t}_{i}\right)}^{{\rm{\top }}}{{\boldsymbol{Q}}}_{g}(\ell ,k){\rm{V}}{\rm{a}}{\rm{r}}\left[\hat{{\boldsymbol{\theta }}}\right]{{\boldsymbol{Q}}}_{g}{\left(\ell ,k\right)}^{{\rm{\top }}}{\mathscr{B}}({t}_{i})},\end{array}\end{array}\end{eqnarray*}$$

where ${\mathscr{B}}({t}_{i})$ is the ith row of the matrix ${\mathscr{B}}$, and ${\mathscr{B}}=[{\bf{B}}| {\bf{B}}| ...| {\bf{B}}]$ is a matrix of dimension N × c. The c × c matrices Q_m , {Q_g (ℓ, k), ℓ = 1, 2, k = 1, ..., K} satisfy ${{\bf{Q}}}_{m}{\boldsymbol{\theta }}={\left({{\boldsymbol{\alpha }}}^{\top },{{\bf{0}}}_{J}^{\top },\,\ldots ,\,{{\bf{0}}}_{J}^{\top }\right)}^{\top }$ and

$$\begin{eqnarray*}\begin{array}{cc}{{\bf{Q}}}_{g}(1,1){\boldsymbol{\theta }}={\left({{\bf{0}}}_{J}^{\top },{{\boldsymbol{\beta }}}_{1}^{\top },{{\bf{0}}}_{J}^{\top },\,\ldots ,\,{{\bf{0}}}_{J}^{\top }\right)}^{\top }, & {{\bf{Q}}}_{g}(2,1){\boldsymbol{\theta }}={\left({{\bf{0}}}_{J}^{\top },\,\ldots ,\,{{\bf{0}}}_{J}^{\top },{{\boldsymbol{\gamma }}}_{1}^{\top },{{\bf{0}}}_{J}^{\top },\,\ldots ,\,{{\bf{0}}}_{J}^{\top }\right)}^{\top },\\ {{\bf{Q}}}_{g}(1,2){\boldsymbol{\theta }}={\left({{\bf{0}}}_{J}^{\top },{{\bf{0}}}_{J}^{\top },{{\boldsymbol{\beta }}}_{2}^{\top },{{\bf{0}}}_{J}^{\top },\,\ldots ,\,{{\bf{0}}}_{J}^{\top }\right)}^{\top }, & {{\bf{Q}}}_{g}(2,2){\boldsymbol{\theta }}={\left({{\bf{0}}}_{J}^{\top },\,\ldots ,\,{{\bf{0}}}_{J}^{\top },{{\bf{0}}}_{J}^{\top },{{\boldsymbol{\gamma }}}_{2}^{\top },{{\bf{0}}}_{J}^{\top },\,\ldots ,\,{{\bf{0}}}_{J}^{\top }\right)}^{\top },\\ & \vdots \vdots \\ {{\bf{Q}}}_{g}(1,K){\boldsymbol{\theta }}={\left({{\bf{0}}}_{J}^{\top },\,\ldots ,\,{{\bf{0}}}_{J}^{\top },{{\boldsymbol{\beta }}}_{K}^{\top },{{\bf{0}}}_{J}^{\top },\,\ldots ,\,{{\bf{0}}}_{J}^{\top }\right)}^{\top }, & {{\bf{Q}}}_{g}(2,K){\boldsymbol{\theta }}={\left({{\bf{0}}}_{J}^{\top },\,\ldots ,\,{{\bf{0}}}_{J}^{\top },{{\boldsymbol{\gamma }}}_{K}^{\top }\right)}^{\top }.\end{array}\end{eqnarray*}$$

In this Appendix we prove the results in Lemma 1 and Proposition 1 (see Section 5).

D.1. Proof of Lemma 1

Notice that the expectation of the periodogram in Equation (15) of the observations $\{{\varepsilon }_{i},\ i\in { \mathcal I }\}$ is

$$\begin{eqnarray*}&&{\rm{E}}\left[{I}_{\,\varepsilon }(\lambda )\right]=\sum _{k\in { \mathcal I }}\sum _{j\in { \mathcal I }}{\rm{E}}\left[{\varepsilon }_{k}{\varepsilon }_{j}\right]\exp (i\lambda [{t}_{k}-{t}_{j}]).\end{eqnarray*}$$

If we replace the expectation ${\rm{E}}\left[{\varepsilon }_{k}{\varepsilon }_{j}\right]$ with the right-hand side of Equation (18), we obtain

$$\begin{eqnarray}&&\begin{array}{l}{\rm{E}}\left[{I}_{\,\varepsilon }(\lambda )\right]=\displaystyle \frac{2\pi }{{N}_{\,{ \mathcal I }}}\sum _{j=1}^{{N}_{\,{ \mathcal I }}}{P}_{\varepsilon }({\omega }_{j})\sum _{k\in { \mathcal I }}\sum _{j\in { \mathcal I }}\exp (i[\lambda -{\omega }_{j}]\\ \ \ \times \,[{t}_{k}-{t}_{j}])=\displaystyle \frac{2\pi }{{N}_{\,{ \mathcal I }}}\sum _{j=1}^{{N}_{\,{ \mathcal I }}}{P}_{\varepsilon }({\omega }_{j}){W}_{\varepsilon }(\lambda -{\omega }_{j})\\ \ \ \ =\,\displaystyle \frac{2\pi }{{N}_{\,{ \mathcal I }}}{P}_{\varepsilon }(\lambda )\,* \,{W}_{\varepsilon }(\lambda ).\end{array}\end{eqnarray} \tag{ D1 }$$

D.2. Proof of Proposition 1

Let ${ \mathcal F }\{{g}_{j}\}[k]$ denote the discrete Fourier transform (DFT) of the sequence of m numbers g₁, ..., g_m into another sequence h₁, ..., h_m , that is,

$$\begin{eqnarray*}&&\begin{array}{l}{h}_{k}={ \mathcal F }\{{g}_{j}\}[k]=\sum _{j\,=\,1}^{m}{g}_{j}\exp (-{ikj}2\pi /m),\\ \ \ \ k=1,\,\ldots ,\,m,\end{array}\end{eqnarray*}$$

and ${{ \mathcal F }}^{-1}\{{h}_{k}\}[j]$ denote the inverse DFT of the sequence h₁, ..., h_m into another sequence g₁, ..., g_m , that is,

$$\begin{eqnarray*}\begin{array}{ccc}{g}_{j} & = & {{ \mathcal F }}^{-1}\{{h}_{k}\}[j]=\frac{1}{m}\displaystyle \sum _{k=1}^{m}{h}_{k}\exp ({ikj}2\pi /m),\\ j & = & 1,\,\ldots ,\,m.\end{array}\end{eqnarray*}$$

Let ${ \mathcal F }\{{g}_{j}\}[k]$ and ${ \mathcal F }\{{{\ell }}_{j}\}[k]$ be, respectively, the DFTs of the sequences {g_j } and {ℓ_j } into the sequences {h_k } and {m_k }. Then, the convolution theorem states that

$$\begin{eqnarray}&&{ \mathcal F }\{{g}_{j}\,* \,{{\ell }}_{j}\}[k]={ \mathcal F }\{{g}_{j}\}[k]{ \mathcal F }\{{{\ell }}_{j}\}[k].\end{eqnarray} \tag{ D2 }$$

Applying the convolution theorem in Equation (D2) to Equation (D1), we obtain

$$\begin{eqnarray*}&&{ \mathcal F }\{{P}_{\varepsilon }({\lambda }_{j})\}[k]=\displaystyle \frac{{N}_{\,{ \mathcal I }}}{2\pi }\displaystyle \frac{{ \mathcal F }\{{\rm{E}}\left[{I}_{\,\varepsilon }({\lambda }_{j})\right]\}[k]}{{ \mathcal F }\{{W}_{\varepsilon }({\lambda }_{j})\}[k]}.\end{eqnarray*}$$

The inverse DFT of the last equation gives the result in Equation (23).