1 Introduction

Multivariate spatial field data are increasingly common and whose modeling typically relies on building cross-covariance functions to describe cross-process relationships. An alternative viewpoint is to model the matrix of spectral measures. We develop the notions of coherence, phase and gain for multidimensional stationary processes. Coherence, as a function of frequency, can be seen to be a measure of linear relationship between two spatial processes at that frequency band. We use the coherence function to illustrate fundamental limitations on a number of previously proposed constructions for multivariate processes, suggesting these options are not viable for real data. We also give natural interpretations to cross-covariance parameters of the Matérn class, where the smoothness indexes dependence at low frequencies while the range parameter can imply dependence at low or high frequencies. Estimation follows from smoothed multivariate periodogram matrices. We illustrate the estimation and interpretation of these functions on two datasets, forecast and reanalysis sea level pressure and geopotential heights over the equatorial region. Examining these functions lends insight that would otherwise be difficult to detect and model using standard cross-covariance formulations.

Keywords: coherency; gain; multivariate random field; periodogram; phase; reanalysis; squared coherence; spectral density

Coherence for Random Fields

William Kleiber111Department of Applied Mathematics, University of Colorado, Boulder, CO. Author e-mail: william.kleiber@colorado.edu

July 5, 2019

1 Introduction

The theory of univariate continuous stochastic processes has become well developed over nearly a century of research. The past quarter century or so has seen an increasing interest and development of models for multivariate spatial processes. The recent review by Genton and Kleiber (2015) gives a relatively comprehensive treatment of the basic approaches that have been explored to build stochastic spatial models. In the discussion, Bevilacqua et al. (2015) pose the question, given the recent deluge of multivariate constructions, “which parametric model is more flexible?” Indeed, the relative strengths and weaknesses of multivariate models have been explored only empirically, that is, by testing a battery of different models on particular datasets, and comparing performance either by likelihood values or by predictive cross-validation (in the multivariate context, spatial prediction is known as co-kriging). Thus, a fundamental open question is: (1) to what extent can the flexibilities of model constructions be compared theoretically? Additionally, most models are motivated in the covariance domain, and the natural follow-up question is: (2) are there other approaches than covariance to measure and quantify spatial dependence?

We introduce the notion of spectral coherence, phase and gain for multidimensional and multivariate spatial random fields. We propose that these functions allow for natural partial answers to the critical open questions (1) and (2). We show that a number of previously proposed models lack sufficient practical flexibility in terms of prediction, and we suggest insights into well established models such as the multivariate Matérn, where parameters such as the cross-covariance smoothness and range have had elusive direct interpretations that relate to process behavior.

Let us illustrate the ideas developed in this manuscript by considering the time series case first. Suppose is a bivariate complex-valued weakly stationary time series on with covariance functions and cross-covariance functions for . The corresponding spectral densities are for . Spectral modeling in time series is well developed; for example, traditional autoregressive moving average models imply processes with rational spectral densities. The function

is known as the squared coherence function, and can be interpreted as a quantification of the linear relationship between and at frequency (Brockwell and Davis, 2009).

We entertain two example datasets from the atmospheric sciences. Both are reforecast and reanalysis data products over the equatorial region based on a well established numerical weather prediction (NWP) model. Reanalysis forecasts are from a fixed version of a NWP model that are run retrospectively to generate a large database of model forecasts and analyses (in this context, an analysis can be considered a best estimate of the current state of the atmosphere). First we look at forecasted surfaces of sea level pressure at daily forecast horizons between 24 and 192 hours. We show that coherence can be used as a diagnostic to assess forecast quality, and additionally illustrate frequency bands at which forecasts improve over time. The second dataset involves geopotential heights at differing pressure levels. We show that coherence and phase extract and highlight qualities of the spatial relationship between different pressure levels that are difficult to model using extant multivariate covariance constructions, and indeed illustrate some fundamental limitations of existing popular constructions.

2 Spectra for Multivariate Random Fields

Suppose is a -variate weakly stationary random field on admitting a matrix-valued covariance function where . For simplicity of exposition we suppose is a mean zero process. For complex-valued stationary processes, , so that . The main obstacle to multivariate process modeling is developing flexible classes of matrix-valued covariance functions that are nonnegative definite. We say is nonnegative definite if, for any choices of and locations for and we have

Note this reduces to the usual definition of nonnegative definiteness for a univariate covariance, .

For univariate processes, Bochner’s Theorem states that is a valid (i.e., nonnegative definite) function if and only if it can be written

where is a positive finite measure (Stein, 1999). If admits a density with respect to the Lebesgue measure on , we call it the spectral density for . The multivariate extension of Bochner’s fundamental result is given by Cramér (1940), and is contained in the following theorem specialized to covariances admitting spectral densities.

Theorem 1 (Cramér 1940).

A matrix-valued function is nonnegative definite if and only if

for such that the matrix is nonnegative definite for all .

The functions are the spectral and cross-spectral densities for the marginal and cross-covariance functions . Note that . When the spectral density exists, it can be solved for as the Fourier transform of the covariance function,

Theorem 1 has primarily been used in practice to build multivariate covariance models, by specifying matrices of spectral densities that are nonnegative definite for all frequencies.

2.1 Coherence

In time series, the notion of frequency coherence is well developed, and can be used, for instance, to assess whether one time series is related to another by a time invariant linear filter. These notions carry over to the spatial case, and form the point of entry for our analyses.

If form a matrix-valued covariance function with associated spectral densities , then define the coherence function (or coherency function)

We might assume for all for , but can define if . The coherence function can be complex-valued, so in practice we examine the absolute coherence function, . The real-valued function is the squared coherence function, and by Theorem 1, for all . Values of near unity indicate a linear relationship between and at particular frequency bands.

The following theorem relates the coherence to optimal prediction of a random process based on another process. The predictive estimator is based on a kernel smoothed process which is a natural predictor given the interpretation of the univariate kriging weights as a kernel function (Kleiber and Nychka, 2015).

Theorem 2.

Suppose is a complex-valued mean zero weakly stationary bivariate field with matrix-valued covariance admitting a spectral density matrix that is everywhere nonzero. Then the continuous square integrable function that minimizes is


Additionally, the spectral density of the predictor is


for all .

In particular, the relationship (1) implies that the optimal weighting function is modulated by the coherence between the two processes, and indeed has greater spectral weight on frequencies with high coherence. An immediate corollary to Theorem 2 is

Thus, the coherence has an attractive interpretation as the amount of variability that can be attributed to a linear relationship between two processes at a particular frequency. In the following development, we use the coherence function to illuminate fundamental limitations on some popular multivariate covariance constructions.

The coherence function can be used as a tool to compare proposed multivariate models, as an indicator of the amount of flexibility of bivariate relationships at differing frequencies. For example, a rather classic approach to specifying covariances is separability, setting where is a univariate covariance function and is a positive definite matrix (Mardia and Goodall, 1993; Helterbrand and Cressie, 1994; Bhat et al., 2010). This approach has been empirically shown to be insufficiently flexible, and the following lemma contributes to the empirical results.

Lemma 3.

If where is a covariance function and is a positive definite matrix with th entry , then the squared coherence between the th and th process is constant, in particular .

A more sophisticated method of generating multivariate covariance structures is to convolve univariate square integrable functions (Gaspari and Cohn, 1999; Oliver, 2003; Gaspari et al., 2006; Majumdar and Gelfand, 2007). In particular, if are square integrable functions for then is a valid matrix-covariance function where denotes the convolution operator. This is sometimes known as covariance convolution (especially when are positive definite functions to begin with). The following proposition suggests this approach to model building is overly-restrictive, and indeed implies that the resulting coherence is necessarily constant over all frequencies.

Proposition 4.

If and are square integrable functions on and a matrix-valued covariance is defined via for where denotes convolution, then for all such that the Fourier transforms of and are nonzero.

Multivariate processes can sometimes be modeled as being related by local averaging. For example, the relationship between column integrated ozone observations and local ozone might be plausibly modeled as observations being locally averaged over the true underlying field (Cressie and Johannesson, 2008). Wind observations are often time averaged over moving windows to produce smoother and more stable observation series (Hering et al., 2015). The following proposition characterizes the coherence in such situations, and serves to illustrate the intimate link between process relationship and coherence.

Proposition 5.

If is a weakly stationary stochastic process and for some continuous square integrable kernel function that is symmetric, then for all such that the Fourier transform of is nonzero.

According to Proposition 5, estimated coherences near unity over all frequency bands may be indicative of a linear or local averaged relationship between processes, and this result may serve as the theoretical basis for testing such a hypothesis. Fuentes (2006) uses a similar notion to develop a test for separability of space-time processes.

The kernel convolution method, introduced by Ver Hoef and Barry (1998) and Ver Hoef et al. (2004), originally involved representing a process as a moving average against a white noise process. In simple cases this yields the covariance convolution model. This can be generalized to


where is a mean zero stationary process with covariance , and are square integrable symmetric kernel functions for . As with each previous construction, this approach also yields constant coherence.

Proposition 6.

If and are constructed as in (3), then is constant for all .

For all of these models, separable, covariance convolution and kernel convolution, the resulting multivariate structure is restricted to constant coherence. In light of Theorem 2, this suggests that none of these models can attain optimal prediction for any multivariate processes exhibiting nontrivial coherences. Indeed the examples below in Section 4 exhibit nonconstant coherence, and call for more flexible modeling frameworks.

2.2 Phase and Gain

Similar notions to frequency coherence can be motivated by examining spectral density matrices. If is a stationary random vector with spectral density matrix then define . Note that is possibly complex-valued. We define the gain function which is sometimes referred to as the gain of on in time series (Brockwell and Davis, 2009). Additionally define the phase function at frequency as , the complex argument of . Note that the phase function satisfies and .

The interpretations of gain and phase are most clear when considering processes built by the relationship for some and , i.e., is a shifted and rescaled version of . Then it is straightforward to show that the phase function is

Li and Zhang (2011) develop an approach to modeling this type of asymmetric cross-covariance behavior. This result shows that their construction will have phase shift function that depends on the angle . This result may be used as an exploratory data approach or as the basis for a statistical test to suggest whether a pair of spatial processes exhibit an asymmetric relationship; Li and Zhang (2011) use the empirical cross-correlation function to visually assess such asymmetric behavior. The gain function in this case is simply ; all frequency components of are exaggerated by an amount for .

Below we consider some multivariate constructions that are particular to real-valued processes having real-valued spectral matrices. Any model with real-valued cross-spectral density has , but a possibly non-trivial gain function. Thus, testing for can be viewed as a test for a real-valued cross-spectral density, which seems very relevant given most multivariate models are developed under this assumption.

2.3 Revisiting the Multivariate Matérn

The multivariate Matérn is a model for matrix-valued covariance functions such that each process is marginally described by a Matérn covariance function, and all cross-covariance functions also fall in the Matérn class (Gneiting et al., 2010; Apanasovich et al., 2012). Specifically, the multivariate Matérn imposes for and for . Here, where is a modified Bessel function of the second kind of order . Gneiting et al. (2010) and Apanasovich et al. (2012) discuss restrictions on the parameters and that result in a valid model.

The Matérn class is popular due to the smoothness parameters which continuously index smoothnesses of the sample paths of the process. In particular, sample paths are times differentiable if and only if , and there is an additional relationship between and the fractal dimension in that sample paths have dimension (Goff and Jordan, 1988; Handcock and Stein, 1993). These interpretations and implications also hold in the multivariate case, where indexes the smoothness of the th component process . The parameters act as range parameters, and control the rate of decay of spatial correlation away from the origin.

A standing issue with the multivariate Matérn is that the cross-covariance parameters, and for , do not have straightforward interpretations that are analogous to the marginal smoothness and range interpretations, and indeed nowhere in the literature have these parameters been linked directly to process behavior. We find that these parameters have direct interpretations when considering the coherence function between two processes.

The coherence function for a bivariate process with multivariate Matérn correlation is

We explore in detail two simplified versions of this coherency function. First, consider the case where , all covariance and cross-covariance functions share a common range. Then

Figure 1: Coherence functions for various bivariate Matérns with , varying and .

The multivariate Matérn is a valid model only if , and thus by inspecting (4) we note some modeling implications resulting from these restrictions. First, if , the coherency is constant across all frequencies, implying a constant linear relationship between two processes at all frequency bands. Second, if , we have greater coherency at low frequencies, with as . This analysis seems to suggest a natural interpretation of the cross-covariance smoothness in that it controls the amount of cross-process dependence at various frequencies, but can only imply greater or equal coherence at low frequencies versus high frequencies. Figure 1 serves to illustrate these points, showing various coherence functions for common length scale parameters , varying and .

Perhaps surprisingly, a similar analysis suggests the cross-covariance range parameter induces potentially greater flexibility in the coherence function than the smoothness parameter . In particular, if , the coherence function is


Examining (5), we see that, depending on whether or the coherence will have distinct behavior. For simplicity, set . If , the coherence will be greater for small frequencies than high frequencies, similar to the behavior implied by (4) when . However, here as , implying non-negligible coherence between processes at high frequencies, unlike that in (4). If , we have the complementary result that for , that is, two processes share behavior at high frequencies rather than low. Figure 2 illustrates these scenarios, and seems to suggest that, at least as far as coherence is concerned, the cross-covariance range parameter yields potentially greater flexibility than the cross-smoothness.

Figure 2: Coherence functions for various bivariate Matérns with , varying and .

The so-called parsimonious Matérn model is defined by imposing common range parameters as well as (Gneiting et al., 2010). This model has been empirically shown to produce inferior model fits to datasets as compared to more general versions of the multivariate Matérn as well as other multivariate classes (Gneiting et al., 2010; Apanasovich et al., 2012). The coherence function for a bivariate parsimonious Matérn model is constant, , which suggests an inflexible model for the spectral behavior of spatial processes.

We close this section with an empirical illustration of the implications of cross-covariance parameter choice on random field realizations and the associated low and high frequency behavior. We simulate two bivariate Matérn models on an equally-spaced grid of in . In both cases we low-pass and high-pass filter the resulting bivariate field. The low-pass filter is a matrix of zeros except for a grid of , while the high-pass filter is similar with a grid of on the edge and in the center. The effect of the filters is to remove high frequency behavior (low-pass filtering) or low frequency behavior (high-pass filtering).

Figure 3 shows low-pass filtered realizations of two bivariate Gaussian processes with bivariate Matérn covariances. The top row is the case with equal ranges but with a greater cross-smoothness . According to Figure 1, we should expect the realizations to show similar low frequency behavior, but dissimilar high frequency behavior. Indeed, panels (a) and (b) show similar low frequency behavior, while the pairwise scatterplot of high-passed values, panel (c), suggests little correlation at high frequencies. Complementary, the second row is the case with equal smoothnesses, but a greater cross-range parameter, . Figure 2 suggest we should expect less coherence at low frequencies while having greater correlation at high frequencies. Again, these theoretical results are reinforced: panels (d) and (e) are not suggestive of strong low frequency coherence, while panel (f) exhibits positively correlated high frequency characteristics (and indeed has an empirical correlation coefficient of approximately ).

Figure 3: Low-pass filtered bivariate Matérn with and for the (a) first process, (b) second process. Panel (c) is the pairwise scatterplot of high-pass filtered values of the two processes with contour levels for comparison. Panels (d), (e) and (f) are analogous plots for a bivariate Matérn simulation with and .

2.4 The Linear Model of Coregionalization

The linear model of coregionalization (LMC) is a competing framework for multivariate modeling, and is built by decomposing a multivariate process as linear combinations of uncorrelated, univariate processes (Goulard and Voltz, 1992; Royle and Berliner, 1999; Wackernagel, 2003; Schmidt and Gelfand, 2003). In particular, we entertain the following version,


The matrix is known as the coregionalization matrix, and controls the strength of dependencies on the latent uncorrelated processes . For the following, suppose and are uncorrelated processes with spectral densities and , respectively.

Given the number of parameters in the LMC, it is often useful to impose restrictions on the coregionalization matrix , such as setting (Berrocal et al., 2010). The following Lemma is the unsurprising result that the LMC yields multivariate processes that are exactly coherent when the coregionalization matrix has zero determinant.

Lemma 7.

In the linear model of coregionalization (6), if then the coherence function between and is unity if and only if .

Note that this result simply states that, under the LMC, two processes are exactly coherent when they differ only by a scalar multiplier.

Under the same working assumptions, , we have the gain function of on is

If, as is common in using the LMC, we set , we have the gain function is simply ; that is, there is constant gain at all frequencies by the amount of coregionalization, . The complementary case where yields the gain

that is, the relative contribution of component to the combined spectrum of . As mentioned previously, if the latent processes have real-valued spectral densities, the phase function is exactly zero at all frequencies.

3 Estimation of Spectra

Suppose is a -variate process that has been observed at a regular grid of points , of marginal dimensions where . If grid spacing in the th dimension is , define . Then the spatial periodogram matrix with th entry is defined as where


is available at Fourier frequencies where for . Note that .

Whereas in the time series literature it is natural to consider asymptotics as time , resulting in effectively uncorrelated blocks of a process, in the spatial realm there are two competing asymptotic frameworks. Increasing domain asymptotics is similar to the time series case where samples are taken on an ever-increasing domain in all axial directions, and typically asymptotic results here echo those in time series. The complementary version is infill asymptotics (sometimes called fixed-domain asymptotics) in which the domain boundary is fixed and points are sampled at an ever finer resolution within the domain (Zhang and Zimmerman, 2005). Depending on the asymptotic framework under consideration, the large sample properties of the periodogram (7) change.

Using infill asymptotics, Lim and Stein (2008) show that the raw multivariate periodogram can exhibit bias at low frequencies, and suggest prewhitening the process to overcome this inadequacy (in the univariate case Stein (1995) gives a simulated example where the bias is quite substantial). However, under a mixture of infill and increasing domain asymptotics, Fuentes (2002) showed (for univariate processes) the analogous result to the time series case that the periodogram is asymptotically unbiased and is uncorrelated at differing Fourier frequencies. Additionally, in this latter case it is not a consistent estimator, but must be smoothed to gain consistency.

Analogous to the time series and univariate spatial field case, under certain assumptions the nonparametric periodogram (7) is asymptotically unbiased, and generates asymptotically uncorrelated random variables between distinct Fourier frequencies. The following theorem illustrates this feature of the matrix-valued periodogram. We use the same assumptions to Fuentes (2002), generalized to the multivariate setting.

  • The true spectral densities decay as as , .

  • The marginal and cross-covariances satisfy .

  • , and for all such that .

Theorem 8.

Under the assumptions A1-A3, we have

  • ,

  • and

  • for .

The proof for Theorem 8 follows directly from Fuentes (2002) and is not included here.

According to Theorem 8, the matrix-valued periodogram is not an asymptotically consistent estimator. To produce a consistent estimator of the spectral density at a particular frequency , in practice we locally smooth adjacent periodogram values and appeal to property (iii) of Theorem 8. In particular, the smoothed matrix-valued periodogram is


where is the empirical cumulative distribution function of Fourier frequencies . Here, is some kernel function with bandwidth , where, as we have it written, the same kernel is applied to each process. Naturally, different kernels may be used for different processes if the scientific context calls for such an approach.

Note that we can’t directly use the nonparametric periodogram fraction to estimate the coherence as at all Fourier frequencies. Thus, we estimate the coherence functions by using the smoothed periodograms,

for .

4 Illustrations

We examine two datasets from the atmospheric sciences, gridded reforecasts and reanalyses of sea level pressure and geopotential heights over the equatorial region. Reforecast data are produced retrospectively from a fixed version of a numerical weather prediction model, in this case the 2nd generation National Oceanic and Atmospheric Administration’s (NOAA) Global Ensemble Forecast System Reforecast (Hamill et al., 2013). Forecasts are generated at 3 hour increments from to hours, with the 0h forecast being a reanalysis, that is, an estimate of the current state of the atmosphere. For the data below, the control initial conditions were produced using a hybrid ensemble Kalman filter-variational analysis system (Hamill et al., 2011).

4.1 Sea Level Pressure

The first dataset we consider is a set of reforecast sea level pressures (SLP) over the equatorial region. Sea level pressures in this region are approximately stationary, and we compare forecast horizons in 24 hour increments from 0h to 192h (8 days out). The data consist of gridded reforecasts from the first 90 days of 2014 at increments over 360 longitude and 47 latitude bands between to defining the equatorial region.

One approach to examining the quality of forecasts is the coherence between the forecast with the corresponding reanalysis. For example, we might compare the 24h forecast of SLP generated on January 1, 2014 to the 0h reanalysis generated on January 2, 2014. It is well known that forecast skill decays with horizon, and we expect the short-term forecasts to share higher coherence with the reanalyses than the long-term forecasts.

We begin by standardizing each analysis and forecast horizon grid cell by subtracting the temporal average and diving by empirical standard standard deviation to produce forecast anomalies. Denote these anomalies by for forecast horizons corresponding to forecast horizons hours, spatial locations in the equatorial region on days , i.e., the first 90 days of 2014.

Each day’s marginal process empirical periodogram (7) is calculated for all forecast horizons , yielding . The smoothed periodogram is a convolution with a simple low-pass filter, a matrix of zeros with a constant block of . Interest focuses on comparing various forecast horizons with the reanalysis at , so we calculate empirical cross-periodograms for all available days allowing for forecast validation (e.g., the , 24h horizon, has 89 available days, ). The cross-periodograms are smoothed using the same low-pass filter as the marginals. If denotes the smoothed cross-periodograms, we estimate the squared coherence function as

for , that is, the average over all available daily smoothed cross-periodograms.

Figure 4: Estimated absolute coherence functions for the GEFS sea level pressure reforecast data, comparing the 168, 96 and 24 hour forecast horizons with the zero hour analysis. The vertical axis spans .

Figure 4 shows estimated absolute coherence functions for horizons 168, 92 and 24h with the 0h analysis. Even at long lead lead times there is substantial coherence, which increases by a substantial margin at very low longitudinal frequencies. For any given longitude frequency band, the coherence appears to be relatively constant across latitudes, which is sensible given that there appears to be greater variability in the equatorial direction than in the north-south direction for sea level pressure in this region. As the forecast horizon decreases the coherence begins building between low-to-mid frequency bands in the latitudinal direction, suggesting that the statistical characteristics of short term forecasts are more similar to observed sea level pressure than the longer term forecasts. However, note that at the highest frequencies there is not a substantial improvement in forecast skill, bordering on no improvement, which suggests that small scale events are difficult to forecast even at one day out.

To quantify the differences in forecast horizon skill, we estimate a set of bivariate Matérn models. As each forecast and analysis arises from the same physical model, we assume that the marginal spectra follow the same statistical behavior, that is, we suppose the marginal spectral densities at equal at all horizons, . Suggested by Figure 4 and exploratory analysis for the marginal spectra, we additionally suppose the spectrum is constant across latitude frequencies. A Gaussianity assumption does not appear to be justifiable, based on empirical Q-Q plots. Thus, we follow the suggestion of Fuentes (2002) and estimate marginal Matérn parameters by minimizing the squared difference between theoretical log spectral density and log average periodogram, having averaged over all days, forecast horizons and latitude bands. The resulting least squares estimates are and . Cross-covariance parameters are estimated by minimizing least squares distance to the average empirical coherence functions, estimated by averaging over days and latitude bands.

Forecast horizon 24 48 72 96 120 144 168 192
0.96 0.94 0.91 0.88 0.85 0.81 0.76 0.73
0.079 0.076 0.075 0.073 0.071 0.070 0.068 0.067
1.06 1.06 1.05 1.04 1.03 1.02 1.01 1.01
Table 1: Cross-covariance Matérn family parameter estimates for GEFS sea level pressure reforecast data.

Table 1 contains the cross-covariance parameter estimates corresponding to all forecast horizons. As forecast horizon increases, all parameters decay; in fact, the decay is almost exactly linear for each variable beyond the 24h horizon. Fitting a linear model to the parameters as a function of forecast horizon suggests the cross-covariance smoothness will decay to the marginal value after approximately 16 days, whereas the cross-covariance scale meets the marginal value between 3-4 days. Although these ideas can be used to generate scientific hypotheses, there is substantial extrapolation for that strongly depends on a linear assumption.

A word of caution is in order; the estimates in Table 1 do not always imply a valid multivariate covariance structure. However, on any given grid the estimated parameters may yield a valid model, it just cannot be guaranteed for all grids. One possibility is that the bivariate Matérn is not sufficiently flexible to describe the stochastic structure of these fields, which is a call for further research in this area.

4.2 Geopotential Height

Our second example is on the same spatial domain, but whose values are geopotential heights. Geopotential height is the height (in meters) above sea level at which the atmospheric pressure is a certain level. In the atmospheric sciences, it is common to examine geopotential heights as indicators of climatic regimes, for instance Knapp and Yin (1996) examines the relationship between heights and temperature anomalies over a portion of the United States.

Three of the most common geopotential height maps are the 850hPa, 500hPa and 300hPa maps. The first, 850hPa, approximately defines the planetary boundary layer, which is the lowest level of the atmosphere that interacts with the surface of the Earth (note 1000hPa is approximately sea level). The 300hPa level is at the core of the jet stream, while the 500hPa approximately divides the atmosphere in half, and whose anomalies are used in part to assess climatological temperature variations. The vertical structure of geopotential heights is a focus of some interest within atmospheric sciences (Blackmon et al., 1979).

We examine geopotential height reanalysis anomalies for representing the 850hPa, 500hPa and 300hPa pressure levels on days , the first 6 months of 2014. The anomalies are differences between the reanalysis height and a time-varying Nadaraya-Watson kernel smoothed estimate of the mean with a bandwidth of 5 days. Experiments suggest the results below are qualitatively robust against choices of the bandwidth and smoothing kernel.

Similar to the previous section, we smooth the marginal process periodogram (7) using a low-pass filter, and calculate smoothed empirical cross-periodograms yielding . Then the squared coherence is estimated as the arithmetic average of each day’s empirical squared coherence estimate,

Figure 5: Estimated absolute coherence functions for the GEFS geopotential height data between (a) 850hPa and 500hPa, (b) 850hPa and 300hPa and (c) 500hPa and 300hPa.

Figure 5 contains the three estimated pairwise absolute coherence functions. There is high coherence between the lower pressure levels at low frequencies, and some evidence of moderate coherence between all levels at low frequencies. We also note some strikingly different behavior than for the sea level pressure example. First, there is an apparent ridge in coherence at low frequencies (approximately ) which may be indicative of equatorial planetary waves (Wang and Xie, 1996; Xie and Wang, 1996; Kiladis et al., 2009). Planetary waves can play crucial roles in the formation of tropical cyclones (Molinari et al., 2007). Additionally, note that there is high coherence at the highest Fourier frequencies for all coherence functions (capping out at approximately and ). This is evidence of a nonseparable relationship in the frequency domain, and we are unaware of any current multivariate models that can adequately capture such behavior. One possible explanation for this high coherence at high frequencies is artifacts in the data assimilation scheme, in particular aberrant observational data leading to unusually large anomalies in geopotential height. Indeed, variables such as sea level pressure are well constrained by a wealth of observational data, while geopotential heights are less constrained, usually being observed by sparsely released weather balloons.

Figure 6: Estimated phase functions for the GEFS geopotential height data between (a) 850hPa and 500hPa, (b) 850hPa and 300hPa and (c) 500hPa and 300hPa.

Figure 6 shows pairwise plots for each pair of geopotential height anomalies. In particular there is strong evidence of a phase shift at a low frequency band between the 850hPa and both lower pressure heights. These frequencies indicate wavelengths of approximately km, which is a typical wavelength for planetary equatorial waves, or Rossby waves (Wang and Xie, 1996; Xie and Wang, 1996; Kiladis et al., 2009). On the other hand, there is not substantial evidence of phase shift between the pairs of heights at other frequencies. Most extant multivariate models utilize real-valued cross-spectral densities, and thus are insufficiently flexible to capture this type of phase shifted behavior at specific spectra.

5 Discussion

The notion of coherence, phase and gain are common in the time series literature, but have been yet unexplored for multivariate spatial processes. We casted these functions for multidimensional processes. The coherence between two variables can be interpreted as a measure of linear relationship at particular frequency bands, resulting in a complementary framework for comparing processes than the usual cross-covariance function. Phase and gain also yield straightforward interpretations as a physical space-shift and relative amplitude of frequency dependence when comparing two processes. We developed these ideas for stationary processes, and future research may be directed toward the analogous cases for nonstationary processes, perhaps extending the work of Fuentes (2002).

Coherence, phase and gain can be estimated using smoothed cross-periodograms, and in our examples we showed that, as exploratory tools, these functions can be very useful in detecting structure that may not be readily captured using extant multivariate models. We additionally illustrated that the coherence function gives a natural interpretation to the multivariate Matérn cross-covariance parameters that have otherwise been uninterpretable, lending insight into an outstanding problem.

A number of future research directions may be considered, including the adaptation of coherence to multivariate space-time processes. This work can also be seen as a call to develop more flexible multivariate models, perhaps working directly in the spectral domain rather than the covariance domain, where most previous work has fallen, echoing the call of Simpson et al. (2015).


The following Lemma is useful in proving some results of the manuscript.

Lemma 9.

Suppose is a stationary processes on with covariance function having spectral density and where is continuous, symmetric and square integrable with Fourier transform . Then has covariance function with associated spectral density . Additionally, the cross-covariance function between and is with spectral density .

The proof of this Lemma involves straightforward calculations involving convolutions and is not included here.

We recall the spectral representation for a stationary vector-valued process with matrix-valued covariance function having spectral measures defined on the Borel -algebra on . There is a set of complex random measures on such that if are disjoint, and for . Then has the spectral representation

see Gihman and Skorohod (1974) for details. If all admit associated spectral densities , then in shorthand we write .

Proof of Theorem 2.

The spectral representation implies

for complex-valued random measures , . Then if has Fourier transform ,

by a change of variables. Then, using that and

we have

The integrand is minimized for each if

That the density of is now follows by the convolution theorem for Fourier transforms. ∎

Proof of Proposition 4.

If is the Fourier transform of , the result immediately follows as the spectral density of is . ∎

Proof of Proposition 5.

This result follows directly from Lemma 9. ∎

Proof of Proposition 6.

This result follows from Lemma 9 and that the spectral density for white noise is constant over all frequencies. ∎


The author thanks Michael Scheuerer for many helpful discussions during the development of this research. This research was supported by National Science Foundation grants DMS-1417724 and DMS-1406536.


  • Apanasovich et al. (2012) Apanasovich, T. V., Genton, M. G., and Sun, Y. (2012), “A valid Matérn class of cross-covariance functions for multivariate random fields with any number of components,” Journal of the American Statistical Association, 107, 180–193.
  • Berrocal et al. (2010) Berrocal, V. J., Gelfand, A. E., and Holland, D. M. (2010), “A bivariate space-time downscaler under space and time misalignment,” Annals of Applied Statistics, 4, 1942–1975.
  • Bevilacqua et al. (2015) Bevilacqua, M., Hering, A. S., and Porcu, E. (2015), “On the flexibility of multivariate covariance models,” Statistical Science, in press.
  • Bhat et al. (2010) Bhat, K. S., Haran, M., and Goes, M. (2010), “Computer model calibration with multivariate spatial output: a case study in climate parameter learning,” in Frontiers of Statistical Decision Making and Bayesian Analysis, eds. M. H. Chen, P. Müller, D. Sun, K. Ye, and D. K. Dey, pp. 401–408, Springer-Verlag, New York.
  • Blackmon et al. (1979) Blackmon, M. L., Madden, R. A., Wallace, J. M., and Gutzler, D. S. (1979), “Geographical variations in the vertical structure of geopotential height fluctuations,” Journal of the Atmospheric Sciences, 36, 2450–2466.
  • Brockwell and Davis (2009) Brockwell, P. J. and Davis, R. A. (2009), Time series: theory and methods, Springer Science & Business Media.
  • Cramér (1940) Cramér, H. (1940), “On the theory of stationary random processes,” Annals of Mathematics, 41, 215–230.
  • Cressie and Johannesson (2008) Cressie, N. and Johannesson, G. (2008), “Fixed rank kriging for very large spatial data sets,” Journal of the Royal Statistical Society, Series B, 70, 209–226.
  • Fuentes (2002) Fuentes, M. (2002), “Spectral methods for nonstationary spatial processes,” Biometrika, 89, 197–210.
  • Fuentes (2006) Fuentes, M. (2006), “Testing for separability of spatial-temporal covariance functions,” Journal of Statistical Planning and Inference, 136, 447–466.
  • Gaspari and Cohn (1999) Gaspari, G. and Cohn, S. E. (1999), “Construction of correlation functions in two and three dimensions,” Quarterly Journal of the Royal Meteorological Society, 125, 723–757.
  • Gaspari et al. (2006) Gaspari, G., Cohn, S. E., Guo, J., and Pawson, S. (2006), “Construction and application of covariance functions with variable length fields,” Quarterly Journal of the Royal Meteorological Society, 132, 1815–1838.
  • Genton and Kleiber (2015) Genton, M. G. and Kleiber, W. (2015), “Cross-covariance functions for multivariate geostatistics,” Statistical Science, in press.
  • Gihman and Skorohod (1974) Gihman, I. I. and Skorohod, A. V. (1974), The Theory of Stochastic Processes, Vol. 1, Springer-Verlag, Berlin.
  • Gneiting et al. (2010) Gneiting, T., Kleiber, W., and Schlather, M. (2010), “Matérn cross-covariance functions for multivariate random fields,” Journal of the American Statistical Association, 105, 1167–1177.
  • Goff and Jordan (1988) Goff, J. A. and Jordan, T. H. (1988), “Stochastic modeling of seafloor morphology: Inversion of sea beam data for second-order statistics,” Journal of Geophysical Research, 93, 13589–13608.
  • Goulard and Voltz (1992) Goulard, M. and Voltz, M. (1992), “Linear coregionalization model: Tools for estimation and choice of cross-variogram matrix,” Mathematical Geology, 24, 269–282.
  • Hamill et al. (2011) Hamill, T. M., Whitaker, J. S., Kleist, D. T., Fiorino, M., and Benjamin, S. G. (2011), “Predictions of 2010’s tropical cyclones using the GFS and ensemble-based data assimilation methods,” Monthly Weather Review, 139, 3243–3247.
  • Hamill et al. (2013) Hamill, T. M., Bates, G. T., Whitaker, J. S., Murray, D. R., Fiorino, M., Galarneau Jr., T. J., Zhu, Y., and Lapenta, W. (2013), “NOAA’s second-generation global medium-range ensemble reforecast dataset,” Bulletin of the American Meteorological Society, 94, 1553–1565.
  • Handcock and Stein (1993) Handcock, M. S. and Stein, M. L. (1993), “A Bayesian analysis of kriging,” Technometrics, 35, 403–410.
  • Helterbrand and Cressie (1994) Helterbrand, J. D. and Cressie, N. (1994), “Universal co-kriging under intrinsic coregionalization,” Mathematical Geology, 26, 205–226.
  • Hering et al. (2015) Hering, A. S., Kazor, K., and Kleiber, W. (2015), “A Markov-switching vector autoregressive stochastic wind generator for multiple spatial and temporal scales,” Resources, 4, 70–92.
  • Kiladis et al. (2009) Kiladis, G. N., Wheeler, M. C., Haertel, P. T., Straub, K. H., and Roundy, P. E. (2009), “Convectively coupled equatorial waves,” Reviews of Geophysics, 47.
  • Kleiber and Nychka (2015) Kleiber, W. and Nychka, D. W. (2015), “Equivalent kriging,” Spatial Statistics, 12, 31–49.
  • Knapp and Yin (1996) Knapp, P. A. and Yin, Z. (1996), “Relationships between geopotential heights and temperature in the south-eastern US during wintertime warming and cooling periods,” International Journal of Climatology, 16, 195–211.
  • Li and Zhang (2011) Li, B. and Zhang, H. (2011), “An approach to modeling asymmetric multivariate spatial covariance structures,” Journal of Multivariate Analysis, 102, 1445–1453.
  • Lim and Stein (2008) Lim, C. Y. and Stein, M. (2008), “Properties of spatial cross-periodograms using fixed-domain asymptotics,” Journal of Multivariate Analysis, 99, 1962–1984.
  • Majumdar and Gelfand (2007) Majumdar, A. and Gelfand, A. E. (2007), “Multivariate spatial modeling for geostatistical data using convolved covariance functions,” Mathematical Geology, 39, 225–245.
  • Mardia and Goodall (1993) Mardia, K. and Goodall, C. (1993), “Spatial-temporal analysis of multivariate environmental monitoring data,” in Multivariate Environmental Statistics, eds. G. P. Patil and C. R. Rao, pp. 347–386, Amsterdam: North Holland.
  • Molinari et al. (2007) Molinari, J., Lombardo, K., and Vollaro, D. (2007), “Tropical cyclogenesis within an equatorial Rossby wave packet,” Journal of the Atmospheric Sciences, 64, 1301–1317.
  • Oliver (2003) Oliver, D. S. (2003), “Gaussian cosimulation: Modelling of the cross-covariance,” Mathematical Geology, 35, 681–698.
  • Royle and Berliner (1999) Royle, A. and Berliner, L. M. (1999), “A hierarchical approach to multivariate spatial modeling and prediction,” Journal of Agricultural, Biological and Environmental Statistics, 4, 1–28.
  • Schmidt and Gelfand (2003) Schmidt, A. M. and Gelfand, A. E. (2003), “A Bayesian coregionalization approach for multivariate pollutant data,” Journal of Geophysical Research – Atmospheres, 108.
  • Simpson et al. (2015) Simpson, D., Lindgren, F., and Rue, H. (2015), “Beyond the valley of the covariance function,” Statistical Science, in press.
  • Stein (1995) Stein, M. L. (1995), “Fixed-domain asymptotics for spatial periodograms,” Journal of the American Statistical Association, 90, 1277–1288.
  • Stein (1999) Stein, M. L. (1999), Interpolation of Spatial Data: Some Theory for Kriging, New York: Springer-Verlag.
  • Ver Hoef and Barry (1998) Ver Hoef, J. M. and Barry, R. P. (1998), “Constructing and fitting models for cokriging and multivariable spatial prediction,” Journal of Statistical Planning and Inference, 69, 275–294.
  • Ver Hoef et al. (2004) Ver Hoef, J. M., Cressie, N., and Barry, R. P. (2004), “Flexible spatial models for kriging and cokriging using moving averages and the Fast Fourier Transform (FFT),” Journal of Computational and Graphical Statistics, 13, 265–282.
  • Wackernagel (2003) Wackernagel, H. (2003), Multivariate Geostatistics, Berlin: Springer-Verlag, third edn.
  • Wang and Xie (1996) Wang, B. and Xie, X. (1996), “Low-frequency equatorial waves in vertically sheared zonal flow. Part I: stable waves,” Journal of the Atmospheric Sciences, 53, 449–467.
  • Xie and Wang (1996) Xie, X. and Wang, B. (1996), “Low-frequency equatorial waves in vertically sheared zonal flow. Part II: unstable waves,” Journal of the Atmospheric Sciences, 53, 3589–3605.
  • Zhang and Zimmerman (2005) Zhang, H. and Zimmerman, D. L. (2005), “Towards reconciling two asymptotic frameworks in spatial statistics,” Biometrika, 92, 921–936.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description