Engineering
Computer Science-2021.gz:143459
Computer Science
8e80d6af979234c4235ea3578f02616155b119f9
Bayesian autoregressive spectral estimation

Autoregressive (AR) time series models are widely used in parametric spectral estimation (SE), where the power spectral density (PSD) of the time series is approximated by that of the \emph{best-fit} AR model, which is available in closed form. Since AR parameters are usually found via maximum-likelihood, least squares or the method of moments, AR-based SE fails to account for the uncertainty of the approximate PSD, and thus only yields point estimates. We propose to handle the uncertainty related to the AR approximation by finding the full posterior distribution of the AR parameters to then propagate this uncertainty to the PSD approximation by \emph{integrating out the AR parameters}; we implement this concept by assuming two different priors over the model noise. Through practical experiments, we show that the proposed Bayesian autoregressive spectral estimation (BASE) provides point estimates that follow closely those of standard autoregressive spectral estimation (ASE), while also providing error bars. BASE is validated against ASE and the Periodogram on both synthetic and real-world signals.


I. INTRODUCTION
In Signal Processing, spectral estimation (SE) refers to the determination of the energy that is contributed by each component in a time series. This work focuses, in particular, on the spectral representation given by the Fourier power spectral density (PSD). In practice, the challenge is to estimate the PSD from a set of noisy observations of the time series; this is seen in applications such as biomedical engineering [1] and telecommunications [2]. Classical methods for computing the PSD can be roughly divided in two categories: parametric and nonparametric. Parametric methods impose a generative model on the data such as autoregressive moving average (ARMA) models [3], sums of sinusoids (Lomb-Scargle [4]), Bernstein polynomials [5] or Gaussian mixtures [6]. On the other hand, nonparametric methods for spectral estimation do not assume any structure on the data and are thus more flexible, yet they do not discriminate the signal from the noise, which critically affects the spectral estimates. Typical examples of the nonparametric SE methods are the classic Periodogram [7] which can be considered as a histogram of frequencies as well as recent Bayesian nonparametric approaches [8], [9], [10].
Even though the randomness conveyed by time series observations is evident, dealing with the uncertainty related to the data-driven computation of the PSD is not a standard practice within SE. In parametric approaches, the model fit to the data through an optimisation procedure, e.g. using least squares or maximum likelihood. Therefore, the PSD is approximated deterministically, by a point estimate, and all the information about the mismatch between the chosen model and the data, which could lead to probabilistic PSD estimates, is unfortunately usually neglected.
To address this void in the SE literature, which is particularly detrimental in noisy data and model mismatch, we aim to cater for uncertainty in PSD estimates finding a distribution over PSDs, rather than point estimates. This can be achieved by equipping the existing framework for parametric spectral estimation with a Bayesian treatment. Instead of committing to a point-estimate of parameters of the time-series model, we can find their posterior distribution, thus constructing a posterior distribution over models for time series. Then, we can transport this distribution via the Fourier transform to find the posterior distribution over PSDs, where posterior in this case is in the sense of conditional to the time series observations. This rationale becomes of particular importance for model misspecification, since, strictly speaking, a parametric model may always be considered to be misspecified for real data and therefore the distribution of models allows us to account for the limited capacity of the parametric representation.
In this work, we apply this generic idea to the particular class of autoregressive time series models, thus our proposed method is referred to as Bayesian autoregressive spectral estimation (BASE). In particular, we propose two models, termed Model I and Model II which difference stems from the noise prior considered: Model I uses a half-Normal prior while Model II uses an Inverse-Gamma (and thus conjugate) prior. Both models assume a normal prior over the AR weights and Model I is implemented using Markov chain Monte Carlo (MCMC) [11] due to the nonconjugate prior. Via numerical simulations, we validate the proposed BASE method and the two models proposed on two synthetic signals and a real world example, these experiments show the advantage of the proposed BASE in different settings against its deterministic counterparts.
where p ∈ N, [a 1 , . . . , a p ] ∈ R p are the (autoregressive) parameters and { t } t∈N is a noise process. We consider the noise to be Gaussian independent and identically distributed random variables with zero mean and variance σ 2 .
A. Power spectral density of AR(p) processes AR(p) processes are ubiquitous in parametric spectral estimation as they have closed-form Fourier PSDs [12] given explicitly by the AR parameters. Specifically, using the Ztransform [3] the PSD of AR(p) can be expressed by The modes (or peaks) of the PSD in eq. (2), i.e., the frequencies where the AR process convey more power, are given by the roots of the denominator. These roots are the poles of the dynamical system that generates the AR(p) process {x t } t∈N . Since we can regard the model order p as the degrees of freedom of the parametric PSD for the AR model, the larger the p the more flexible the PSD. This is in line with the interpretation in the temporal domain, where a model with more lags can cater for more complex time series and thus represents a more general model.
The fact that the PSD above is directly parametrised by {a 1 , . . . , a p , σ n } is perhaps the main motivation for using the AR model for spectral estimation. However, even though there have been Bayesian approaches to fit the AR parameters they have not been applied to the SE problem. In practice, the literature shows that the use of AR parameters in SE has been mainly deterministic, e.g., by computing the PSD in eq. (2) with a plug-in approximation of the AR parameters such as those found by least squares, the method of moments, or maximum likelihood. Such an approach, does not allow for modelling the uncertainty related to the peaks of the resulting PSD, due to noisy data, or to account for model misspecification (i.e., incorrect model order) which can lead to either peak splitting (over-estimation) or onto fewer peaks (under-estimation).

III. BAYESIAN SPECTRAL ESTIMATION OF AUTOREGRESSIVE SIGNALS
Our aim is to compute the posterior distribution of the PSD for the AR case, thus providing a natural account for model uncertainty, which arises from i) the consideration of a finite dataset, ii) model mismatch, iii) noisy observations. From now on, we will refer to the standard (deterministic) autoregressive spectral estimation approach as ASE and to the proposed Bayesian counterpart as BASE.

A. Generic Bayesian parametric spectral estimation
The posterior distribution of S x , the PSD of a stationary stochastic process {x} t∈N , conditional to a set of observations x = [x 1 , x 2 , . . . , x n ], can be computed by integrating out the model of the time series. Denoting the space of all possible time series models by M, the posterior PSD is given by where the last expression only takes the (reasonable) assumption that the PSD, S x , and the observations, x, are conditionally independent given the model M . This follows from the fact that the data cannot say more about the PSD than the model itself.
Notice that the above integral is defined over the entire, possibly infinite-dimensional, model space M. To calculate this integral we need to assume a structure for the time-series models to be considered; this can be achieved by choosing a finite-dimensional prior over the models, namely p(M ). A straightforward way to choose this prior is to first assume a parametrisation over models, say M = M θ , θ ∈ R d , and then choose a prior over the parameters p(θ). This construction is know as the push-forward measure [13] that transports a distribution over the finite-parameter, θ ∈ R d , towards the space of models through the mapping θ → M θ .
With this parametrisation, the posterior in eq. (3) can be expressed with respect to the model parameters, θ, via the change of variable theorem as where the integration is now performed over the (finite) parameter space, that is, Under the assumption that the process {x t } t∈N is stationary and given by a model M θ , its PSD is uniquely defined by the model's parameters θ and thus the distribution p(S x |θ) is a Dirac measure supported on a single points on the space of PSDs. Therefore, sampling from the posterior over PSDs, p(S x |x), is straightforward: simply sample a parameter from the posterior distribution over parameters θ * ∼ p(θ|x), and then map this parameter sample to PSDs according to where PSD(M ) denotes the PSD corresponding to the model M .
In order to perform this procedure, we need to choose: i) a model space M θ , ii) a prior over the parameters p(θ), and iii) a sampling procedure. We refer to these in the remaining part of this section.

B. Model I: Half-normal noise
We shall choose the model space as that of AR models owing to their appealing properties for spectral estimation outlined in Section II-A. Furthermore, we will assume that all parameters in the AR model are independent, that a 1:p , p ∈ N, are Normally distributed and the noise variance σ 2 follows a half-Normal prior. Thus, the priors are where the half-Normal distribution corresponds to the multiplication between a Normal distribution and an indicator function 1 {σ 2 ≥0} , and then properly renormalised (i.e., multiplied by 2). The choice of normality for the autoregressive parameters follows from the fact that the Gaussian prior is conjugate to the (conditional) likelihood as verified in the next section, and thus the posterior is also normally distributed.
Regarding the definition of the likelihood, given a set of observations {x 1:T }, T ∈ N consider the conditional likelihood of the AR(p), that is, the expression p(x p+1:T |x 1:p , θ), where θ = [a 1 , . . . , a p , σ 2 ] denotes all model parameters and we consider {x 1:p } to be fixed and not part of the generative model for simplicity of presentation.
Following eq. (1), the conditional likelihood of the AR(p) model is given by where we have adopted the notation a = [a 1 , . . . , a p ] and The conditional likelihood is thus Gaussian on the AR parameter vector a, meaning that the Gaussian prior on a established on the previous section is conjugate to the conditional likelihood and results on the following Gaussian posterior Remark. Since the complete (rather than conditional) likelihood is given by it requires calculating the distribution p(x 1:p |θ); this is challenging since it involves the infinite-order moving-average representation of the AR(p) process. Therefore, under the assumption that the first p observations are fixed they become independent from θ, and the density p(x 1:p |θ) becomes a constant for theta in eq. (10). A similar reasoning applies for the marginal distribution p(x p+1:T , x 1:p ). Notice that the discrepancy between the conditional and complete likelihoods becomes negligible for increasing amounts of data (i.e., T 1). For this reason, we will consider only the conditional likelihood.
With this choice, generating a sample, θ * , from the posterior requires the evaluation of p(x p+1:T |x 1:p , θ)p(a)p(σ 2 ). Then, this sample can be used to compute the PSD explicitly, according to eq. (2). This procedure avoids explicit computation of the model M θ , since the PSD is computed directly from the parameters.
Notice that, by construction, all the posterior PSD samples will have the form as in eq. (2) which means that no probability will be assigned to expressions that do not follow this form. This advantage of the parametric approach to SE results in that, using a set of samples of PSDs constructed as explained above, we can numerically compute the mean PSD and error bars conditional to a set of signal values x.

C. Model II: Inverse-Gamma noise
In order to derive a closed-form expression for the posterior over the AR parameters, we can assume the noise variance σ 2 follows an Inverse-Gamma prior distribution, rather than the Half-Normal one in the previous section. This posterior is known as the Normal-Inverse-Gamma distribution.
Specifically, let us consider the following priors where I p denotes the p−dimensional identity matrix. Then the posterior distribution for θ = (a 1:p , σ 2 ) is proportional to a Normal-Inverse-Gamma: where, using the notation X p+1: From this Norma-Inverse-Gamma posterior, we can sample the parameters (a 1:p , σ 2 ) simply by specifying the hyperparameters (µ 0 , λ, α, β) and using standard sampling techniques. Our conjecture is that avoiding the MCMC stage required above can result in faster posterior computations.
The MAP estimators for the parameters (a 1:p , σ 2 ) are given by where recall our notation X := X p+1:T and x := x p+1:T and we have assumed that we have m observations of the signal (though in practice usually only m = 1 observation is available). In this model, the hyperparameters can be chosen using a Grid search and cross-validation with a 5-fold split: for each set of hyperparameters, compute the MAP estimate on the train split, to then score it on the test split using the likelihood given that MAP for a 1:p and σ 2 as the score function. The final score assigned to the hyperparameters is the average of the scores obtained on each fold. Lastly, the set of hyperparameters with the best score will be used in the posterior.

IV. SIMULATIONS
We validated the proposed Bayesian autoregressive spectral estimation method, BASE, over three case studies, two for synthetic data and one for real world data. The first one considered a synthetic AR signal, where the true and approximate model orders were different, this simulation is aimed to show that BASE is robust to model misspecification due to the chosen prior. The second experiment illuminates the unbiasedness and concentration properties of the approximate posterior over PSDs by applying BASE on a continuoustime signal generated by a Gaussian process (GP) [14] with known PSD. Finally, the third experiment validates the ability of the proposed BASE model to detect periodicities from a sub-sampled real world astronomical time series. We also compare, over all the case studies, the performance of the method in both versions: the Model I in Sec. III-B based on a half-Normal noise prior (using MCMC) and the Model II in Sec. III-B based on the conjugate Inverse-Gamma noise prior with closed-form posterior. In the experiments using MCMC we used a NUTS sampler [15] in the PyMC3 Python toolbox [16].

A. Misspecified model order: AR-generated time series
We implemented a stable AR(4) time-series to generate data that was later processed by BASE (Model I) with the assumption of a AR(10) model. Fig. 1 shows the true values of the regression coefficients and their approximate posterior using MCMC; notice that these posteriors contain the true AR(4) coefficients, and they also converge to narrow posteriors centred around zero for higher orders coefficients. This validates the self-regularisation property, or robustness to misspecified model order, within the proposed BASEframework which will be key when computing the sample PSDs. Next, using the approximate posterior over the AR parameters (coefficients and variance), we generated PSD samples according to eq. (2). Fig. 2 shows the density over PSDs using the proposed BASE (95% error bars), the standard autoregressive spectral estimation (ASE), the true PSD and the Periodogram. Notice that the proposed BASE succeeded in estimating the true PSD while remaining unbiased and concentrated around the true value.

B. Continuous-time signal: Gaussian process with Laplace covariance
Let us consider a synthetic signal generated from a Gaussian process (GP) [14] with Laplace covariance function given by where σ 2 denotes the marginal covariance and l is known as the process lenghtscale. Furthermore, we considered zeromean Gaussian noise added to the GP sample. Recall that a sample from a GP can be regarded as a signal generated by an AR(∞) model. The proposed BASE (Model I) was implemented on this synthetic signal with an AR(4) model, the order of which was chosen from the rate of decay of the Laplace kernel. Then, the PSD estimate of the BASE model was compared against those of the standard ASE of the same order and the Periodogram. Fig. 3 shows the PSD estimates (95% error bars for BASE) as well as the true PSD of the signal, given by the Fourier Transform of the Laplace kernel. Notice that the PSD posterior within BASE was concentrated around the true PSD and was able to capture the mean of its deterministic counterpart.

C. Finding the periodicity of an astronomical time series
For the third experiment, we considered a real-world time series which may not be an AR process and therefore its PSD will not be precisely given by the expression in eq. (2). This experiment was constructed to validate the Bayesian approach to handle model misspecification by averaging over the posterior models. The signal considered was the sunspots dataset [17], known for having a period of 11 years (frequency ≈ 0.0909[years −1 ]), we implemented BASE (Model I) to detect this periodicity only using 1/6 of the available data.
We first computed both the autocorrelation (ACF) and partial autocorrelation (PCF) functions of the sunspots series. Fig. 4 shows the ACF and PCF, and show that due to the slow decay of PCF, the order of the autoregressive component of the signal cannot be specified with any certainty; this supports the evidence that the sunspot series is not an AR process. We implemented BASE (Model I) with an AR(9) model to then compare it against the ASE and the Periodogram, where ASE and the Periodogram used all the available data. Fig. 5 shows the peak estimates of the PSDs using BASE, where the main peak (colour-coded in red) is 0.09181278[years −1 ] and thus in line with the expected frequency peak. Fig. 6 shows the PSD estimates for the three models, with the proposed BASE model exhibiting its peak near the correct known period of the sunspots signal, as it was expected from Fig. 5.
Notice that the errors bars of BASE are now larger than those of the previous experiments, this is because the signal considered in this experiment is not coming from an AR process and thus we are in the scenario of model mismatch. In this sense, due to the fact that BASE averages over possible AR models, its error bars are larger but they contain other peaks of the Periodogram. Finally, observe that despite the fact that an order p = 9 for the model was chosen, this did not result in multiple peaks which validates BASE's robustness to overfitting arising from its Bayesian nature.
D. Performance of the closed-form posterior BASE (Model II) was also implemented in the same experiments described above for comparison. Since this model has a closed form posterior, we refer to it as "CF" in the remaining part of this section. Following the experiment in Sec. IV-A, the proposed CF method was implemented on the AR(4) synthetic time series, where the real parameters are known. Fig. 7 shows the PSD estimate (95% error bars) along with the ASE estimation, the true PSD and the Periodogram. By comparing these results with those shown in Fig. 2, note that the estimation of the mean is only slightly better, likewise, the variance of the estimation is slightly reduced.
Next, we implemented CF on the same synthetic Gaussian process signal generated in IV-B. Fig. 8 shows the PSD estimates (95% error bars for CF) as well as the true PSD of the signal, given by the Fourier Transform of the Laplace kernel. The estimates provided by the CF method (BASE Model II) are very similar to those of BASE Model I, which uses MCMC, as they also concentrates around the true PSD and its mean captures the estimation of the ASE. Notice that the error bars are slightly better for CF than those given by BASE Model I. This suggest that the posterior distribution is more concentrated on the true model. Lastly, following Sec. IV-C we validated the CF model on real data by testing it on the sunspots data-set [17]. Fig. 9 shows the results with the CF approach (BASE Model II), which again exhibits slightly less variance than BASE Model I using MCMC. Although the peak frequency found by CF has shifted to the left, the posterior PSD is still precise and concentrated on the true peak (found by computing the Periodogram over the entire dataset). This suggest that the CF could also recover the fundamental frequency of the data from a limited number of observations.

V. DISCUSSION AND FUTURE WORK
We have proposed a novel framework termed BASE, a Bayesian approach to autoregressive spectral estimation. BASE exploits the closed-form properties of the AR model to compute power spectral densities (PSD), and also complements it with a probabilistic treatment, thus allowing us to quantify uncertainty in the estimates of the PSDs through the use of confidence intervals.
The description of the proposed method follows the idea that we can sample from the posterior distribution over PSDs (given a set of observations of the time series) by (i) sampling the posterior parameters of an AR model, to then (ii) compute the corresponding PSD of such sample parameter. In this setting, we proposed two models: both assume a multivariate normal prior on the AR parameters by they differ in that one assumes a half-Normal prior for the noise and the other one adopts an Inverse-Gamma; we dubbed them Model I and Model II respectively. BASE Model I uses MCMC to find the posterior samples while Model II features a conjugate prior and this results in a closed-form posterior.
Both BASE models have been implemented on three experiments using: i) data generated by an AR model, ii) data generated by a Gaussian process with Laplace kernel, and iii) a real-world astronomical time-series. Through these simulations, we have validated BASE in terms of robustness to model misspecification, unbiasedness, accuracy, and periodicity detection. Regarding the comparison between Model I and Model II, we can say that, Model II gives posteriors that are tighter than those of Model I, arguably due to its use of a closed from posterior. Furthermore, the main feature of Model II is the reduction on the computation time due to avoiding the use of MCMC. The main setback of Model II, however, is that it requires the use of a specific prior. Future work includes: i) heuristics to determine the hyperparameters µ 0 , λ, α, β, ii) a quantitative assessment of the reduction in computational complexity provided by Model II, iii) a development of a hierarchical prior to simultaneously identify the AR model order (e.g., a Dirichlet prior), and iv) extensions to ARMA models to cater for moving-average spectral components and spectra with zeros, iv) comparisons against nonparametric models such as [18].