Frequency-Domain Stochastic Modeling of Stationary Bivariate or Complex-Valued Signals

There are three equivalent ways of representing two jointly observed real-valued signals: as a bivariate vector signal, as a single complex-valued signal, or as two analytic signals known as the rotary components. Each representation has unique advantages depending on the system of interest and the application goals. In this paper we provide a joint framework for all three representations in the context of frequency-domain stochastic modeling. This framework allows us to extend many established statistical procedures for bivariate vector time series to complex-valued and rotary representations. These include procedures for parametrically modeling signal coherence, estimating model parameters using the Whittle likelihood, performing semi-parametric modeling, and choosing between classes of nested models using model choice. We also provide a new method of testing for impropriety in complex-valued signals, which tests for noncircular or anisotropic second-order statistical structure when the signal is represented in the complex plane. Finally, we demonstrate the usefulness of our methodology in capturing the anisotropic structure of signals observed from fluid dynamic simulations of turbulence.



There are no comments yet.


page 2

page 21


Learning Representations Using Complex-Valued Nets

Complex-valued neural networks (CVNNs) are an emerging field of research...

Stratified Graph Spectra

In classic graph signal processing, given a real-valued graph signal, it...

Classification of simulated radio signals using Wide Residual Networks for use in the search for extra-terrestrial intelligence

We describe a new approach and algorithm for the detection of artificial...

The Widely Linear Complex Ornstein-Uhlenbeck Process with Application to Polar Motion

Complex-valued and widely linear modelling of time series signals are wi...

Associative Memories Using Complex-Valued Hopfield Networks Based on Spin-Torque Oscillator Arrays

Simulations of complex-valued Hopfield networks based on spin-torque osc...

Tikhonov Regularization of Circle-Valued Signals

It is common to have to process signals or images whose values are cycli...

A Probabilistic Spectral Analysis of Multivariate Real-Valued Nonstationary Signals

A class of multivariate spectral representations for real-valued nonstat...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

This paper develops new modelling and estimation procedures for complex-valued time series. We derive the interpretation of such models, and efficient methods for their estimation based on sampled time series. Complex-valued series are a useful representation of bivariate series, and arise in a significant number of important application areas, such as functional Magnetic Resonance Imaging (Rowe, 2005), Doppler ultrasound signals of bloodflow (Brands et al., 1997), optical tomography (Erkmen and Shapiro, 2006), seismic readings (Samson, 1977)

, neural networks 

(Mandic and Goh, 2009), meteorology and oceanography (Gonella, 1972) and blind source separation (Eriksson and Koivunen, 2006). We refer to Section 1.9 in Schreier and Scharf (2010) for a more extensive list.

There have also been a number of recent advances in methodological statistics for the analysis of complex-valued vectors, as well as for time-series samples (Walden and Rubin-Delanchy, 2009; Schreier and Scharf, 2010; Walden, 2013). A common need in applications is the ability to specify a simple parametric model for the signal structure, and then to estimate these parameters from a finite sample. However, approaches to inference for complex-valued series have mainly been non-parametric, and little focus has been placed on methods for efficient parametric inference. For reasons that will be come clear, our approach to parametric modelling will be in the frequency-domain. To this end, a main contribution of this paper is extending the Whittle likelihood (Whittle, 1953b) to complex-valued time series, addressing new ways to deal with finite sample effects in the process.

Figure 1: (Top left) trajectories from the Global Drifter Program. (Top right) trajectory of a North Atlantic drifter ID #44000. (Bottom left) velocities of this drifter over time in each Cartesian direction. (Bottom right) the estimated power spectrum of the complex-valued velocities for the first 40 days only.

To show the relevance of our work, we plot measurements from a large array of freely-drifting oceanographic instruments, created and distributed by the Global Drifter Program (GDP, in the first panel of Figure 1. The GDP is a database of global circulation measurements important for our understanding of decadal to century-scale climate change. Each series in this data set can be represented as a complex-valued series. For example, we show a 150-day trajectory from Global Drifter Program ID #44000 in the top right of Figure 1. The velocities corresponding to the eastward and northward displacements, together with the periodogram (the magnitude square of the Fourier transform) of the complex-valued velocity for the first 40 days of these two velocity signals, are displayed on the bottom row of Figure 1. We have plotted the negative and positive frequency axes overlaid, and it is clear that certain features are only present on one side of the spectrum, indicating oscillatory behaviour with a preferred direction of spin. Such structure is significantly easier to model in the frequency domain than in the time domain. The first key contribution of this paper is the creation of appropriate frequency-domain models for complex-valued time series like the surface drifter velocity records.

The second key contribution is the development of an improved inference method. Specifically, we practically and efficiently adapt Whittle likelihood for real and complex-valued time series. For asymptotic understanding, well-established results have been the consequence of the initially derived theory, e.g. Whittle (1953b); Dzhaparidze and Yaglom (1983). For finite however, more subtle effects are known to be in evidence (Dahlhaus, 1988). Most observed time series are nonstationary (as in Figure 1), such that stationarity can only be reasonably assumed over small windows, for example by assuming local stationarity or segmenting the series. This motivates the development of parsimonious models, and indicates the need for bias-corrected estimation procedures for finite-length samples.

To achieve accurate inference with small sample sizes, we here propose the use of the blurred Whittle likelihood. This novel approach accounts for sampling effects such as leakage and aliasing, while maintaining the computational efficiency of the standard Whittle likelihood. We prove consistency and discuss properties of our likelihood, and analyse performance through simulations against state-of-the-art time domain approximations (Anitescu et al., 2012). While our focus is on complex-valued series, the advantages of the blurred Whittle likelihood are just as relevant for applications to real-valued series.

In addition, we extend a number of related statistical procedures to complex-valued time series. We propose semi-parametric procedures for models that are misspecified at certain frequencies, relating this to Robinson (1995). We construct a parametric test for impropriety, which tests for anisotropy if the time series track spatial positions. Finally, we provide model choice criteria for selecting between nested models.

At first glance, it may seem that formulating theory and methods for complex-valued time series is just a special case of the study of bivariate real-valued series. In this paper we will show that this is not the case. We introduce flexible models for complex-valued signal structure that incorporate the possibility of anisotropy or impropriety. A focus here is on the Matérn process, a natural choice for modelling multiscale variability such as that seen in the surface drifter trajectories. We shall illustrate such concepts with practical examples.

The paper is organised as follows. In Section 2 we provide some necessary preliminaries for real- and complex-valued time series. In Section 3 we introduce our models for Cartesian and rotary structure in complex-valued time series. Section 4 introduces the Whittle likelihood for complex-valued time series and also the blurred likelihood which correctly accounts for finite-sample bias in both real- and complex-valued time series. We then prove consistency of the blurred Whittle likelihood in Section 5. In Section 6 we discuss semi-parametric inference, likelihood ratio tests and model choice. Finally, in Section 7 we demonstrate our modelling and inference procedures, including a detailed analysis of oceanographic flow and seismic data. Concluding remarks can be found in Section 8. Readers only interested in developments for real-valued processes should only refer to Sections 2.1, 4.1, 4.2, 5 and 7.1.

2 Preliminaries

2.1 Real-valued Time Series

The sequence is defined from a corresponding zero-mean continuous-time Gaussian process , by , where is the sampling period and . Assuming that the process is second-order stationary, one may define the autocovariance sequence for , where denotes expectation. The power spectral density of forms a Fourier pair with the autocovariance sequence:


The angular frequency is given in radians. For a sample , a simple, but statistically inconsistent estimate of is the periodogram—denoted —which is the squared absolute value of the Discrete Fourier Transform (DFT) defined as


Using the spectral representation theorem we can write in terms of orthogonal increments


where . Note that because is a discrete sequence, has already been aliased. Thus there may be departures between and the continuous-time process’ spectral density, which we denote as .

A direct consequence of equations (2.2) and (2.3) is that by calculating the DFT from a finite sample size, we have automatically “blurred” the spectral information, due to leakage and aliasing effects. Leakage is associated with the truncation of the sequence over a finite time interval, whereas aliasing results from the discrete sampling of the continuous-time process to generate an infinite sequence; both effects are a consequence of analysing a sample of a continuous-time process. As it does not account for these effects, the periodogram is considered a naïve estimator of the continuous time process’ spectral density. The expected value of the periodogram is equal to the convolution of the true spectrum with the Fejér kernel (the magnitude-squared of the Dirichlet kernel), which can badly bias spectral estimates (e.g. (Percival and Walden, 1993)). Leakage effects are often tackled by using tapering (e.g. Velasco and Robinson (2000)

) which reduces bias and variance at each frequency, at the expense of increased correlation between neighbouring Fourier frequencies; we will discuss this trade-off in more detail in this paper. The issue of aliasing results in energy outside of the Nyquist frequency (the highest observable frequency,

) being folded back into the observable spectrum, which is often ignored when estimating spectra of continuous sampled processes (Percival and Walden, 1993).

The classical approach to frequency-domain inference, the Whittle likelihood, makes use of the periodogram and therefore relies on asymptotic behaviours for very large sample sizes to minimise bias effects. In Section 4, we propose a new modification of the Whittle likelihood that directly incorporates finite sampling effects from both leakage and aliasing correctly, which can be used with either the periodogram or tapered spectral estimates.

2.2 Complex-valued Time Series

As is already done in my applications areas, it is often convenient to express bivariate time series as complex-valued. The real and imaginary parts often represent displacements or velocities in orthogonal spatial directions, and will be henceforth referred to as the Cartesian components. A more useful decomposition of the variability, however, is to separate the process into rotary components, that is, into contributions from positively-rotating and negatively-rotating complex exponentials. Statistical understanding of complex-valued time series, in particular in the frequency domain, has been developed over the past 20 years (see e.g. Walker (1993) and Walden (2013)).

To form a complex-valued series, we start from a pair of real-valued sequences, and . These processes can be combined to form the complex-valued sequence via . We shall assume that and are jointly Gaussian and zero mean, such that any collection of the processes at a number of sampled points are Gaussian. Therefore to fully characterise their properties we need to specify the joint covariance as well as and . If the process is stationary then we have and , as well as for the cross-covariance sequence; the latter then forms a Fourier pair with , as in equation (2.1). Alternatively, in terms of the stationary complex-valued Gaussian process , the process may be fully characterised by specifying the covariance and relation sequence


with denoting the complex conjugate of . Clearly specifying and for the complex-valued process is equivalent to specifying , , and for the bivariate real-valued process. It is directly apparent from equation (2.4) that (Hermitian symmetry) and (symmetry). The spectrum and complementary spectrum are defined from and via a Fourier transform

where and . The spectrum and complementary spectrum together describe the second order properties of the Fourier transform of the complex-valued process . We note that is not the relation sequence of the Fourier-domain orthogonal increment process , but rather is the Fourier transform of the relation sequence of . If the relation sequence, or the complementary spectrum, is zero everywhere then the process is said to be proper, and otherwise it is improper (Schreier and Scharf, 2010).

A key motivation behind this paper is not to interpret complex-valued time series as bivariate real-valued time series, but instead as processes that have two complex-valued counter rotating components. This decomposes a process according to the direction of spin. Specifically, we write


where and are both analytic, that is, supported only on positive frequencies. Under the assumption of stationarity, the spectra of the two rotary processes are given by and , while the cross-spectrum is ; all of these are defined only for , and vanish otherwise.

The relationships between three different representations of a complex-valued process – Cartesian (), complex () and rotary () – are presented in Table 1. The derivations can be found in Section A of the appendix. Note that relationships involving do not hold at frequency , as the energy at this frequency is divided between the two rotary components – these specific relationships can also be found in the appendix. The rotary components are defined only for positive frequencies as they are analytic and hence equal to zero for negative frequencies. The relationships between the Cartesian and complex/rotary representations are nontrivial, and highlight some of the differences between working with bivariate and complex-valued data. We will refer to Table 1 in future sections to provide more intuition and simplify calculations. Similar expressions relating the covariance (and relation) sequences of each representation can be found from Table 1 using inverse Fourier transforms.

Cartesian () Complex () Rotary () Cartesian 3-5 3-5 Complex 3-5 Rotary 3-5 3-5

Table 1: Relationships between spectra in three different representations of a complex-valued process: Cartesian (left column), complex-valued (middle column), and rotary (right column). All three elements in each row are equal, as shown in Section A of the appendix.

Sometimes the decomposition is made slightly differently, such that is analytic and the second component is defined to be anti-analytic, i.e. have a spectrum only defined on negative frequencies. We prefer our representation, and consider it more natural as we then only need to specify a covariance structure between and , as the relation sequence is zero everywhere, in contrast to Walden and Rubin-Delanchy (2009) for example.

At first the division in equation (2.5) may seem artificial and opaque. One might argue that representing a bivariate time series as complex-valued is unnecessarily complicated, making the inclusion of two extra complex-valued processes seem yet more redundant. We shall see that splitting up a process into two components in this way is natural, as it divides the process into forward rotating and reverse rotating phasors, as is clear from equation (2.5). The rotary decomposition has been used as a modelling tool for many years and across many applications, as detailed in Section 1.9 in Schreier and Scharf (2010).

3 Modelling techniques for complex-valued time series

In this section we discuss key aspects of the specification of stochastic models for complex-valued time series, including the possibility of coherence between time series components. Specifying rotary coherence is different from specifying Cartesian coherence. As a practical example, we focus on models based on the Matérn process.

3.1 Cartesian coherency

As real-world time series may be characterised by significant coherence between components, we must develop stochastic models that are capable of capturing such behaviour. A common way of defining covariance structure in bivariate series is in the time domain, but in this paper we do this in the frequency domain by modelling the frequency-domain coherency. We represent the cross-spectrum between and as


where is the coherency of and ; is the coherence, and is the group delay, quantifying whether or are leading or lagging in time at frequency-cycle . We refer to as the Cartesian coherency. The group delay and coherence must satisfy


on account of the Hermitian symmetry of .

From a modelling standpoint, we may also regard equation (3.1) as specifying the Cartesian cross-spectrum for a given choice of the individual spectra and , together with a suitable choice of the coherency . It can be shown that choosing with and does indeed lead to a valid spectral matrix that is always non-negative definite. Realistic functions specifying would include the logistic function


as we note directly from equation (3.2) that must be an even function. Furthermore it is reasonable to have decaying polynomials such that , which will be the case in equation (3.3) provided the highest-order polynomial coefficient is positive, i.e. . Similarly, as

is odd, we may propose the form


The simplest choice for the orders and is , which imposes a simple structure of fixed correlation across all frequencies. Note that while a constant phase is not valid, we may have , where is the signum function; or we may have phase that is locally constant within particular frequency bands (see Hamon and Hannan (1974, p. 137)).

Next we ask what constraint is imposed on the phase function if we require the time series to be proper. It follows from Table 1 that we require and in order for the time series to be proper. For the second condition we require that at each positive frequency , either or and , such that and are strictly “out of phase.” The latter possibility follows from the fact that a time series is zero-mean Gaussian and proper if and only if it is circular (i.e. a series whose properties are invariant under rotation), see Schreier and Scharf (2010). An example of such would be a deterministic time series which follows an exact circle, which corresponds to and .

For spatially generated processes such as the oceanographic time series we shall study subsequently, propriety is found to be related to the condition of spatial isotropy, with impropriety in turn implying anisotropy. Thus we can test for anisotropy in the spatial process that generates the sampled time series, and so if the complex-valued time series is generated from a spatial process, we to refer to propriety and isotropy interchangeably.

3.2 Rotary coherency

Alternatively one may employ a model which specifies the coherence between rotary components and . As we have chosen the two to be supported only on positive frequencies however, there can be no covariance between and . Therefore in addition to specifying the variance of and , corresponding to specifying the rotary spectra, we only need to model their covariance at the same frequency. The rotary coherency (Gonella, 1972) is defined as


where again . Taking equation (3.5) as a model for the rotary cross spectrum, we find he flexibility for to vary across frequencies particularly useful. In oceanographic data, for example, it is unlikely that clockwise and anti-clockwise oscillations will be correlated identically at all frequencies—but at different frequencies we would expect a different mode of correlation, see e.g. Davis et al. (1981). Another desirable aspect of the rotary model specification is that we can account for the possibility of different slopes governing the spectral decay across positive and negative frequencies. As we shall see later such would be extremely challenging to incorporate into the model for and , but considerably easier to incorporate for the complex-valued rotary decomposition.

As in the Cartesian case, one may construct models for the coherence and group delay comparable to equations (3.3) and (3.4), with the difference that these functions now only need to be supported over positive frequencies. If we desire to model proper complex-valued series, we see from Table 1 that we require that for . Thus, a proper complex-valued time series has no rotary coherency at non-zero frequencies, which is an important interpretation of propriety. As an aside, we point out that if one desires to model a non-zero-mean proper time series, we must additionally ensure that vanishes; this is because only for . Modelling the rotary coherency can be used to formulate a test for impropriety, as is detailed later in Section 6.2.

3.3 The complex-valued Matérn process

To illustrate the specification of a model with frequency-domain coherence, we will employ a particular type of random process that is in increasingly widespread use. The Matérn process (Gneiting et al., 2010) is a continuous stationary Gaussian process that has received much attention in the spatial and spatiotemporal modelling literature. After carrying some theoretical extensions in the remainder of this section, we will also employ the Matérn process in the applications of Section 7.

The univariate Matérn process is a three-parameter, stationary Gaussian process with covariance given by


where and is the modified Bessel function of the second kind. The smoothness parameter defines the Hausdorff dimension (equal to max) and the degree of differentiability of the process, the range parameter is a timescale parameter, where can be referred to as the correlation timescale, and defines the magnitude of the variability of the process. The Matérn is a useful process in time series analysis, as the limiting behaviour over short timescales (or as ) is like nonstationary fractional Brownian motion, whereas the limiting behaviour over long timescales (or as

) is white noise – the Matérn hence provides a continuum between these two regimes over different timescales. We can therefore view the process as truly multiscale. The process has the potential to capture complicated structure while simultaneously being parsimonious, yielding useful modelling potential for structured processes with which we wish to infer physically meaningful parameters from relatively few observations.

3.3.1 Multivariate Matérn process

There has been much interest in constructing valid multivariate Matérn processes in the spatial statistics literature. In particular, Gneiting et al. (2010) and Apanasovich et al. (2012) have constructed multivariate Matérn processes such that all marginal covariances and cross-covariances are themselves also of Matérn type, i.e. for dimensions and (where ) of the -variate process

where is the correlation coefficient between and (which is trivially set to 1 if ), and is the univariate Matern covariance function defined in equation (3.6). Gneiting et al. (2010) specifies necessary and sufficient conditions for the bivariate Matérn process () to have cross-covariances that are of Matérn type and for the covariance matrix to be nonnegative definite. In this paper, however, we choose to specify the cross-covariance structure in the frequency domain. Our rationale for this choice is twofold. First, for the oceanography application we study, it is easier to specify the covariance of the flow process in terms of specially chosen complex-valued processes from physical considerations. Second, it is very easy to check whether a specified choice corresponds to a valid covariance (we simply require in equation (3.5)), unlike the torturous calculations sometimes necessary in the time domain. The spectrum of the continuous-time Matérn process has a succinct and interpretable representation:


where indicates the spectral density of the continuous-time process. This relates to the aliased spectral density of a sequence —as defined in equation (2.1)—via the relationship for .

We can specify that and are Matérn processes with spectra that follow the form given in equation (3.7), for . We denote these parameters as and , which are permitted to differ between the processes. If we specify a model for the cross-spectra using equation (3.5), then the overall process is valid as long as . By “valid” we mean it possesses an invertible covariance/spectral matrix, but not necessarily a cross-covariance/cross-spectrum that is also Matérn, unlike the definition of Gneiting et al. (2010). In this sense our specification is more general, particularly as we can propose frequency-dependent coherence. If we do wish to formulate cross-spectra that are of Matérn form, we can set and as real-valued and constant across frequency, but still allow and . It is easy to show using equation (3.5), that this specifies a valid Matérn as per the Gneiting et al. (2010) definition, however for comparison, we verify that this parameterisation also meets the conditions required by Gneiting et al. (2010), in Section B of the appendix.

3.3.2 Parametric rotary models

Finally, it is interesting to discuss the Matérn family in the context of a transformation framework, which we will make use of later in the applications. We define the deformation operator , where is an amplitude scaling transformation and is a temporal dilation or compression. We take , , assuming that and are independent realisations of the same stochastic process. This deforms the original stochastic processes and also introduces correlations. It follows that , and we see that the processes remain jointly stationary, but stretched or compressed in time by the factor . If we look at the cross-spectrum then it transpires that it takes the form

We find that if takes a Matérn form, then so does , but with altered parameters. Both and change from such deformations but not . This indicates how a jointly Matérn process can be generated by a simple transformation of an original Matérn.

Alternatively we may combine a Matérn process in each of the two Cartesian or rotary components, with the parametric specification of coherency discussed in Sections 3.1 and 3.2. Choosing the rotary representation, and setting the damping and slope parameters of the positive and negative rotary components to be equal, the form of the spectral matrix becomes


where the coherence function remains to be specified. This type of bivariate model would regulate the same scaling behaviour between and , but the energy of each process being individually regulated by and . We use this model to capture both isotropic and anisotropic structure in oceanographic flow processes in Section 7.2.

The usefulness of the Matérn family is that multiscale behaviour is achievable from one single family. A subtle issue, however, is the relationship between the rotary and Cartesian decompositions. Specifying two Matérns with unequal slope and/or range parameters in one representation will, in general, force the other to posses a mixture of Matérns; this can be inferred from the identities in Table 1. In general, A Matérn specified in and only specifies a Matérn in and if and , and if and are constant across frequencies333However, we note that one could artificially construct an unrealistic tailored form for such that the cross-spectra are coincidentally of Matérn form even when or .; the reverse it also true. In such cases the cross-spectra and are in general complex-valued and out-of-phase. This means that a simple model in one decomposition does not lead to a simple model in the other, and the practitioner must therefore understand the data and the processes driving the data to find the simplest decomposition in which to model.

4 The Whittle Likelihood

In this section we turn to the problem of estimating the parameters of a chosen model from observations. Maximum likelihood is the preferred approach but has two drawbacks. First, its large computational cost renders it unwieldy to apply in practice; this is due to the computational cost associated with inverting the full time-domain covariance matrix. Second, we wish to parameterise the structure in terms of properties of and , which is not conveniently implemented in the time-domain. The method we shall choose, however, will closely approximate the time-domain log-likelihood, preserving its desirable qualities.

Henceforth we let and indicate length  samples of the corresponding random processes arranged as column vectors. In the univariate case of a sample from a zero-mean real-valued Gaussian process, after choosing a particular parametric model characterised by a vector of parameters , the time domain log-likelihood can be written as


where we note that the subscript “” is a label standing for “time”, the superscript “” denotes the transpose, and the superscript “” is the matrix inverse. Here is an theoretical covariance matrix that would take, under the assumption that it is a sample drawn from the proposed model. denotes the determinant of . To find the log-likelihood appropriate for the bivariate process consisting of and we concatenate the two time domain samples into a single vector: . The theoretical covariance matrix for this vector under the assumed model is denoted by . The log-likelihood of the vector can then be written as


in direct parallel with equation (4.1). The best choice of parameters for our chosen model to characterise the observed data is found by maximising the likelihood function

Because the time-domain maximum likelihood is known to have optimal properties, any other estimator will be compared with the properties of this quantity.

A standard technique to approximating equation (4.1) is in the frequency domain, from the seminal work of Whittle (1953b), which results in vastly increased computational speed. This approach approximates using a Fourier representation, and and utilises the special properties of Toeplitz matrices. For a single real-valued time series the Whittle likelihood is given by


where the subscript stands for “Whittle” and the superscript “” denotes the Hermitian transpose. is the parametric form of the theoretical spectrum of the continuous-time process as a function of frequency, and is the set of discrete Fourier frequencies: . The simplification of equation (4.4) follows as and are scalar-valued at each frequency; however, we keep the form of equation (4.3) for comparison with later models in which the spectrum will be matrix-valued at each frequency. Note that may be chosen to be smaller that the full range of Fourier frequencies, in which case the Whittle likelihood can be considered as a semi-parametric fitting procedure (see Section 6.1 for more discussion).

The Whittle likelihood approximates the time-domain likelihood when , i.e. , and this statement can be made precise (Dzhaparidze and Yaglom, 1983). However its computational cost is versus for the time-domain likelihood, a vast increase in speed that renders it more suitable for applications. There have been recent advances however in constructing time-domain approximations to maximum likelihood (see Anitescu et al. (2012)) and we compare performance in a simulation study in Section 7.2.

4.1 The Blurred Whittle Likelihood

The Whittle likelihood is based on the periodogram spectral estimate , which is known to be a problematic quantity. The periodogram is in general an inconsistent and biased measure of the continuous time process’ spectral density, due to the blurring caused by leakage and aliasing effects. It is important to note that the desirable properties of the Whittle likelihood rely on the asymptotic behaviour of the periodogram for very large sample sizes. In order to prevent problems with the periodogram for finite sample sizes from adversely affecting our inference process, we define an alternative likelihood function to and based on , written as


where the subscript “” stands for “frequency.” Here is a version of the theoretical modelled spectrum that has been blurred by a convolution with the Fejér kernel . It can readily be shown that if is drawn from a process with a given theoretical spectrum, the expected value of its periodogram will blur the theoretical spectrum in precisely the same way. Therefore the bias in the estimation process caused by leakage and aliasing arising from finite sample effects has now been correctly accounted for. We call this the blurred Whittle likelihood and prove consistency of estimators in Section 5, as well as document its performance in a Monte Carlo study in Section 7.1.

The usual form of the Whittle likelihood can be computed efficiently by using the Fast Fourier Transform (FFT) to calculate the DFT, which is an operation. The blurred Whittle likelihood however, is slightly more complicated to compute as must now be calculated. The best way to do this is by working with the Fourier pair of the blurred theoretical spectrum , which is the expected biased sampled autocovariance sequence of length , denoted for ,


Observe that the truncation to finite length has the effect of multiplying the true autocovariance sequence with the triangle kernel. In fact, both aliasing and leakage effects are automatically accounted for in equation (4.6); the effect of aliasing is accounting for by sampling the theoretical autocovariance function at discrete times, while the effect of leakage due to the truncation of the sample to finite length is accounted for by the triangle kernel. By forming the blurred and aliased spectrum in this way from a time domain multiplication, we avoid an expensive summation over the aliased frequencies, and the algorithm may still be performed in operations.

Note that if the closed form of the autocovariance sequence is unknown, then we can inverse Fourier transform the modelled spectrum of the process to approximate the autocovariance sequence numerically (see equation (2.1)). Subsequently we would take the Fourier transform of this sequence, multiplied with the triangle kernel (as in equation (4.6)), to recover an estimate of , still in operations. As an aside, we note that computing the Whittle likelihood with the aliased, but not blurred, spectrum —defined in equation (2.1)—can be complicated, as this seldom has an analytic form and must be instead approximated by either arithmetically wrapping in contributions from frequencies higher than the Nyquist, or via an approximation to the Fourier transform in equation (2.1). This is in contrast to the blurred Whittle likelihood, where the blurred spectrum can be computed exactly in one operation, as in equation (4.6). It is therefore easier to incorporate both leakage and aliasing effects together, rather than separately.

4.2 Parametric Tapered Whittle Inference

So far we have encountered the time-domain likelihood, the standard Whittle likelihood, and the blurred Whittle likelihood. There is a fourth possible likelihood that could be defined by extending the Whittle likelihood to tapered discrete Fourier transforms (Velasco and Robinson, 2000). This corresponds to replacing the direct spectral estimator formed from in equation (2.2) with one using the taper

where is real-valued. Setting recovers the periodogram estimate of equation (4.3). The tapered Fourier transform is used to reduce bias in non-parametric estimation from leakage (Percival and Walden, 1993). In the parametric setting we take


Velasco and Robinson (2000) demonstrated that for certain discrete processes it is beneficial to use this estimator, rather than the standard Whittle likelihood, for parameter estimation—particular when the spectrum exhibits a large dynamic range. Nevertheless, tapering will not remove all leakage effects, and the issue of aliasing for continuous sampled processes remains. We therefore define a fifth possible type of likelihood, which accounts for finite sampling effects from tapering


We call this the tapered blurred Whittle likelihood, and it can be computed using a similar calculation to equation (4.6) to find

Note that the time-domain kernel must be pre-computed (unless it has an analytic form) which requires operations, but this can be stored for each taper and data-length , and does not have to be recomputed each time when performing numerical optimisation.

With the modifications to the Whittle likelihood we have proposed, the practitioner is free to select between a tapered or periodogram spectral estimate in the same way as before, but can now account for finite sample bias effects by using our methods. Both the blurred tapered and blurred periodogram likelihoods have advantages and disadvantages. Tapering reduces broadband leakage, but increases variance (by effectively throwing away statistical information) and also increases correlations between local frequencies. The periodogram may exhibit greater global bias and correlation with other frequencies, on account of the implicit convolution with the Fejér kernel; the magnitude of this effect depends on the degree of dynamic range of the observed process in frequency. We explore the comparison between these various methods, and demonstrate the importance of blurring the likelihood, through Monte Carlo simulations in Section 7.1.

4.3 The Whittle Likelihood for Complex-Valued Time Series

In the previous section, modified versions of the Whittle likelihood appropriate for a single real-valued time series were developed. Next we construct the Whittle likelihood for complex-valued time series. We start from the bivariate Whittle likelihood (Whittle, 1953a) and then reformulate this for rotary () and complex () components using the relationships given in Section 2.2. While the Whittle Likelihood is normally derived with continuous-time spectra (see equation (4.3)), for notational simplicity in this subsection we formulate in terms of aliased spectra for sampled sequences, denoted , as defined in equation (2.1); the representations however do not change for continuous-time spectra. We first define the Cartesian spectral matrix and the Cartesian Fourier transform vector as

Using these definitions, the standard bivariate likelihood is (Whittle, 1953a)


where again is the set of discrete Fourier frequencies: . We note that this likelihood is the bivariate form of equation (4.3), and it approximates the bivariate likelihood of given in equation (4.2). The bivariate likelihood (4.9) can be equivalently expressed in term of the rotary components, as will now be shown. We start by defining the rotary spectral matrix and rotary Fourier vector:

The complex-valued length- sequences and are defined implicitly by equation (2.5) for . There is a simple linear relationship between and


as can be seen directly from together with equation (2.5). By inverting equation (4.10) we have that


We then substitute equation (4.11) into equation (4.9)


By expanding and rearranging, and using at , we see that

which from Table 1 is simply equal to , meaning that equation (4.12) simplifies to