Abstract
Multivariate spatial field data are increasingly common and whose modeling typically relies on building crosscovariance functions to describe crossprocess 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 crosscovariance 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 crosscovariance formulations.
Keywords: coherency; gain; multivariate random field; periodogram; phase; reanalysis; squared coherence; spectral density
Coherence for Random Fields
William Kleiber^{1}^{1}1Department of Applied Mathematics, University of Colorado, Boulder, CO. Author email: 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 crossvalidation (in the multivariate context, spatial prediction is known as cokriging). 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 followup 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 crosscovariance 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 complexvalued weakly stationary time series on with covariance functions and crosscovariance 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 matrixvalued covariance function where . For simplicity of exposition we suppose is a mean zero process. For complexvalued stationary processes, , so that . The main obstacle to multivariate process modeling is developing flexible classes of matrixvalued 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 matrixvalued function is nonnegative definite if and only if
for such that the matrix is nonnegative definite for all .
The functions are the spectral and crossspectral densities for the marginal and crosscovariance 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 matrixvalued 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 complexvalued, so in practice we examine the absolute coherence function, . The realvalued 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 complexvalued mean zero weakly stationary bivariate field with matrixvalued covariance admitting a spectral density matrix that is everywhere nonzero. Then the continuous square integrable function that minimizes is
(1) 
Additionally, the spectral density of the predictor is
(2) 
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 matrixcovariance 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 overlyrestrictive, and indeed implies that the resulting coherence is necessarily constant over all frequencies.
Proposition 4.
If and are square integrable functions on and a matrixvalued 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 spacetime 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
(3) 
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 complexvalued. 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 crosscovariance 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 crosscorrelation 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 realvalued processes having realvalued spectral matrices. Any model with realvalued crossspectral density has , but a possibly nontrivial gain function. Thus, testing for can be viewed as a test for a realvalued crossspectral 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 matrixvalued covariance functions such that each process is marginally described by a Matérn covariance function, and all crosscovariance 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 crosscovariance 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 crosscovariance functions share a common range. Then
(4) 
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 crosscovariance smoothness in that it controls the amount of crossprocess 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 crosscovariance range parameter induces potentially greater flexibility in the coherence function than the smoothness parameter . In particular, if , the coherence function is
(5) 
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 nonnegligible 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 crosscovariance range parameter yields potentially greater flexibility than the crosssmoothness.
The socalled 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 crosscovariance parameter choice on random field realizations and the associated low and high frequency behavior. We simulate two bivariate Matérn models on an equallyspaced grid of in . In both cases we lowpass and highpass filter the resulting bivariate field. The lowpass filter is a matrix of zeros except for a grid of , while the highpass 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 (lowpass filtering) or low frequency behavior (highpass filtering).
Figure 3 shows lowpass 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 crosssmoothness . 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 highpassed values, panel (c), suggests little correlation at high frequencies. Complementary, the second row is the case with equal smoothnesses, but a greater crossrange 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 ).
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,
(6) 
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 realvalued 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
(7) 
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 everincreasing domain in all axial directions, and typically asymptotic results here echo those in time series. The complementary version is infill asymptotics (sometimes called fixeddomain 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 matrixvalued periodogram. We use the same assumptions to Fuentes (2002), generalized to the multivariate setting.

The true spectral densities decay as as , .

The marginal and crosscovariances satisfy .

, and for all such that .
Theorem 8.
Under the assumptions A1A3, we have

,

and

for .
According to Theorem 8, the matrixvalued 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 matrixvalued periodogram is
(8) 
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 filtervariational 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 shortterm forecasts to share higher coherence with the reanalyses than the longterm 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 lowpass 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 crossperiodograms for all available days allowing for forecast validation (e.g., the , 24h horizon, has 89 available days, ). The crossperiodograms are smoothed using the same lowpass filter as the marginals. If denotes the smoothed crossperiodograms, we estimate the squared coherence function as
for , that is, the average over all available daily smoothed crossperiodograms.
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 northsouth direction for sea level pressure in this region. As the forecast horizon decreases the coherence begins building between lowtomid 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 QQ 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 . Crosscovariance 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 contains the crosscovariance 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 crosscovariance smoothness will decay to the marginal value after approximately 16 days, whereas the crosscovariance scale meets the marginal value between 34 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 timevarying NadarayaWatson 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 lowpass filter, and calculate smoothed empirical crossperiodograms yielding . Then the squared coherence is estimated as the arithmetic average of each day’s empirical squared coherence estimate,
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 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 realvalued crossspectral 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 crosscovariance function. Phase and gain also yield straightforward interpretations as a physical spaceshift 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 crossperiodograms, 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 crosscovariance 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 spacetime 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).
Appendix
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 crosscovariance 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 vectorvalued process with matrixvalued 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 complexvalued 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 . ∎
Acknowledgements
The author thanks Michael Scheuerer for many helpful discussions during the development of this research. This research was supported by National Science Foundation grants DMS1417724 and DMS1406536.
References
 Apanasovich et al. (2012) Apanasovich, T. V., Genton, M. G., and Sun, Y. (2012), “A valid Matérn class of crosscovariance 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 spacetime 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, SpringerVerlag, 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 spatialtemporal 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), “Crosscovariance 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, SpringerVerlag, Berlin.
 Gneiting et al. (2010) Gneiting, T., Kleiber, W., and Schlather, M. (2010), “Matérn crosscovariance 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 secondorder 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 crossvariogram 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 ensemblebased 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 secondgeneration global mediumrange 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 cokriging under intrinsic coregionalization,” Mathematical Geology, 26, 205–226.
 Hering et al. (2015) Hering, A. S., Kazor, K., and Kleiber, W. (2015), “A Markovswitching 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 southeastern 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 crossperiodograms using fixeddomain 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), “Spatialtemporal 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 crosscovariance,” 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), “Fixeddomain 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: SpringerVerlag.
 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: SpringerVerlag, third edn.
 Wang and Xie (1996) Wang, B. and Xie, X. (1996), “Lowfrequency 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), “Lowfrequency 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.