Visualization and Assessment of Spatio-temporal Covariance Properties

TL;DR

This paper introduces a functional data analysis framework for visualizing and testing the properties of spatio-temporal covariance functions, specifically separability and symmetry, using real and simulated data.

Contribution

It proposes a novel approach combining functional boxplots and nonparametric tests to assess covariance properties from spatio-temporal data.

Findings

Effective visualization of covariance properties using functional boxplots.

Nonparametric tests successfully detect non-separability and asymmetry.

Method performs well on simulated and real-world datasets.

Abstract

Spatio-temporal covariances are important for describing the spatio-temporal variability of underlying random processes in geostatistical data. For second-order stationary processes, there exist subclasses of covariance functions that assume a simpler spatio-temporal dependence structure with separability and full symmetry. However, it is challenging to visualize and assess separability and full symmetry from spatio-temporal observations. In this work, we propose a functional data analysis approach that constructs test functions using the cross-covariances from time series observed at each pair of spatial locations. These test functions of temporal lags summarize the properties of separability or symmetry for the given spatial pairs. We use functional boxplots to visualize the functional median and the variability of the test functions, where the extent of departure from zero at all temporal lags indicates the degree of non-separability or asymmetry. We also develop a rank-based nonparametric testing procedure for assessing the significance of the non-separability or asymmetry. The performances of the proposed methods are examined by simulations with various commonly used spatio-temporal covariance models. To illustrate our methods in practical applications, we apply it to real datasets, including weather station data and climate model outputs.

Tables

Table 1. Table 1: Notations of the chosen covariance models.

Model	Notation	Formula
Gneiting	$C^G(𝐡,u;\beta )$	$\frac{1}{0.5\left|u\right|+1}\text{exp}\left\{-\frac{\Vert 𝐡\Vert }{{\left(0.5\left|u\right|+1\right)}^{\beta /2}}\right\}$
Cressie and Huang separable	$C_{sep}^{CH}(𝐡,u)$	$\frac{1}{{\left(0.5\left|u\right|+1\right)}^2{\left({\Vert 𝐡\Vert }^2+1\right)}^{3/2}}$
Cressie and Huang non-separable	$C_{nsep}^{CH}(𝐡,u)$	$\frac{\left(0.5\left|u\right|+1\right)}{{\left\{{\left(0.5\left|u\right|+1\right)}^2+{\Vert 𝐡\Vert }^2\right\}}^{3/2}}$
Cesare	$C^C(𝐡,u;k)$	$\frac{1}{3}\left\{\exp \left(-𝐡\right)\exp \left(-u\right)+\exp \left(-u\right)+k\exp \left(-𝐡\right)\right\}$

Table 2. Table 2: Type I error and power for simulated data from the Gneiting model C G ​ ( 𝐡 , u ; β ) superscript 𝐶 𝐺 𝐡 𝑢 𝛽 C^{G}(\mathbf{h},u;\beta) with different values of β 𝛽 \beta in 100 100 100 experiments.

	Nominal Level	Type I error	Power
	$\alpha$	$\beta =0$	$\beta =0.25$	$\beta =0.5$	$\beta =0.75$	$\beta =1$
Proportion of	$5\%$	0.04	0.26	0.63	0.96	1.00
rejections	$1\%$	0.02	0.10	0.45	0.85	0.98

Table 3. Table 3: Proportion of rejections for simulated data from the Cressie-Huang separable and non-separable models in 100 100 100 experiments.

	Nominal Level	Model
	$\alpha$	$C_{sep}^{CH}$	$C_{nsep}^{CH}$
Proportion of	$5\%$	0.04	0.97
rejections	$1\%$	0.00	0.87

Table 4. Table 4: Type I error and power for simulated data from the Gneiting model C G ′ ​ ( 𝐡 , u ; λ ) superscript 𝐶 superscript 𝐺 ′ 𝐡 𝑢 𝜆 C^{G^{\prime}}(\mathbf{h},u;\lambda) with different values of λ 𝜆 \lambda in 100 100 100 experiments.

	Nominal Level	Type I error	Power
	$\alpha$	$\lambda =0$	$\lambda =0.025$	$\lambda =0.05$	$\lambda =0.075$	$\lambda =0.1$
Proportion of	$5\%$	0.04	0.07	0.20	0.47	0.70
rejections	$1\%$	0.00	0.01	0.07	0.25	0.39

Equations

$$C(\mathbf{h},u)=\frac{\sigma^{2}}{(a|u|^{2\alpha}+1)^{\beta}}\hbox{exp}\Big{\{}-\frac{c\|\mathbf{h}\|^{2\gamma}}{(a|u|^{2\alpha}+1)^{\beta\gamma}}\Big{\}},$$

$$C(\mathbf{h},u)=\frac{{\sigma }^2(a|u|+1)}{\{(a|u|+1{)}^2+b^2\Vert \mathbf{h}{\Vert }^2{\}}^{3/2}},$$

$$C(\mathbf{h},u)=\frac{{\sigma }^2}{(a|u|+1{)}^2(b^2\Vert \mathbf{h}{\Vert }^2+1{)}^{3/2}}.$$

$$C(\mathbf{h},u)=k_1{C}_s(\mathbf{h})C_t(u)+k_2{C}_s(\mathbf{h})+k_3{C}_t(u),$$

$$f_{\mathbf{h}}(u)=C(\mathbf{h},u)/C(\mathbf{h},0)-C(0,u)/C(0,0).$$

$${\hat{f}}_{\mathbf{h}}(u)=\hat{C}(\mathbf{h},u)/\hat{C}(\mathbf{h},0)-\hat{C}(0,u)/\hat{C}(0,0),u=0,\dots ,p-1,$$

$$\widehat{C}(\mathbf{h},u)=\dfrac{\displaystyle\sum_{i=1}^{p-u}\big{\{}Z(\mathbf{s}_{1},i)-\frac{1}{p-u}\sum_{j=1}^{p-u}Z(\mathbf{s}_{1},j)\big{\}}\big{\{}Z(\mathbf{s}_{1}+\mathbf{h},i+u)-\frac{1}{p-u}\sum_{j=1}^{p-u}Z(\mathbf{s}_{1}+\mathbf{h},j+u)\big{\}}}{p-u},$$

$$C(\mathbf{h};u)=\frac{(1-\lambda )}{a|u{|}^{2\alpha }+1}exp\left\{-\frac{c\Vert \mathbf{h}\Vert }{(a|u{|}^{2\alpha }+1{)}^{\beta /2}}\right\}+\lambda {\left(1-\frac{1}{2v}|h_1-vu|\right)}_{+},$$

$$C^{G^{\prime }}(\mathbf{h},u;\lambda )=\frac{1-\lambda }{0.2|u|+1}exp(-\Vert \mathbf{h}\Vert )+\lambda {\left(1-\frac{1}{2}|h_1-0.2u|\right)}_{+}.$$

$$g_{\mathbf{h}}(u)=C(\mathbf{h},u)-C(\mathbf{h},-u).$$

Full text

Visualization and Assessment of Spatio-temporal Covariance Properties

Huang Huang111 CEMSE Division, King Abdullah University of Science and Technology, Thuwal 23955-6900 Saudi Arabia. E-mail: huang.huang@kaust.edu.sa; ying.sun@kaust.edu.sa and Ying Sun1

Abstract

Spatio-temporal covariances are important for describing the spatio-temporal variability of underlying random processes in geostatistical data. For second-order stationary processes, there exist subclasses of covariance functions that assume a simpler spatio-temporal dependence structure with separability and full symmetry. However, it is challenging to visualize and assess separability and full symmetry from spatio-temporal observations. In this work, we propose a functional data analysis approach that constructs test functions using the cross-covariances from time series observed at each pair of spatial locations. These test functions of temporal lags summarize the properties of separability or symmetry for the given spatial pairs. We use functional boxplots to visualize the functional median and the variability of the test functions, where the extent of departure from zero at all temporal lags indicates the degree of non-separability or asymmetry. We also develop a rank-based nonparametric testing procedure for assessing the significance of the non-separability or asymmetry. The performances of the proposed methods are examined by simulations with various commonly used spatio-temporal covariance models. To illustrate our methods in practical applications, we apply it to real datasets, including weather station data and climate model outputs.

Some key words: full symmetry, functional boxplot, functional data ranking, rank-based test, separability, spatio-temporal covariance

1 Introduction

Spatio-temporal covariance functions play an important role in parametrically modeling geostatistical data, particularly for Gaussian random fields. Spatio-temporal models often have complex structures and rely on various assumptions to simplify the model and reduce the computational burden (Gneiting et al., 2006), such as full symmetry and separability (Cressie and Huang, 1999; Kyriakidis and Journel, 1999; Gneiting, 2002; Stein, 2005). Specifically, let $Z(\mathbf{s},t)$ be a second-order stationary random process at spatial location $\mathbf{s}\in\mathbb{R}^{d}$, and temporal point $t\in\mathbb{R}.$ The covariance function $C(\mathbf{h},u)=\hbox{Cov}\{Z(\mathbf{s},t),Z(\mathbf{s}+\mathbf{h},t+u)\}$ is positive definite and only depends on the spatial lag $\mathbf{h}$ and temporal lag $u$. Then, the covariance function is called fully symmetric if $C(\mathbf{h},u)=C(-\mathbf{h},u)=C(\mathbf{h},-u)$, and separable if $C(\mathbf{h},u)=C(\mathbf{h},0)C(\mathbf{0},u)/C(\mathbf{0},0)$, for any spatial lag $\mathbf{h}$ and temporal lag $u$. One can show that separability implies full symmetry. The separable covariance model is a product of purely spatial and temporal covariances, ignoring the interaction between space and time; the symmetric covariance model acknowledges the interaction but assumes that the spatio-temporal cross-covariances depend on neither the direction of temporal lags nor the direction of spatial lags. Although these models with simplified structures are easier to fit, they have limitations and may not be realistic in practical applications. In separable models, even small spatial changes can lead to large changes in the correlations, causing a lack of smoothness (Stein, 2005), and full symmetry is often violated when data are influenced by dynamic processes with a prevailing flow direction (Gneiting, 2002).

Therefore, it has been an active research area to develop more flexible spatio-temporal covariance models motivated by real applications. Jones and Zhang (1997) developed families of spectral densities that lead to non-separable covariance models without a closed form. Cressie and Huang (1999) proposed families of non-separable covariance functions based on the Fourier transform of non-negative, finite measures with explicit expressions, but these are limited to classes with known analytical solutions of Fourier integrals. Later, Gneiting (2002) extended this approach to more general classes with Fourier-free implementations. All these models are non-separable but fully symmetric. Stein (2005) generated covariance functions that are isotropic but not fully symmetric by taking derivatives of fully symmetric models. Gneiting et al. (2006) proposed anisotropic, asymmetric models by adding a compactly supported Lagrangian correlation function to a fully symmetric covariance model, where the coefficient of the Lagrangian component controls the extent of asymmetry.

In real data analysis, before choosing a covariance function from the existing models, it is necessary to assess the properties of separability and symmetry either by exploratory data analysis or through formal hypothesis tests. However, it is challenging to visualize and assess such properties from spatio-temporal observations. Brown et al. (2001) determined the separability from the closeness to zero of the maximum likelihood estimates of blurring parameters associated with separability, without considering the uncertainty. Shitan and Brockwell (1995) developed an asymptotic test for separability, but it is restricted to spatial autoregressive processes. Fuentes (2006) presented a test of separability in the spectral domain using a simple two-factor analysis of variance. Regarding likelihood-based tests, Mitchell et al. (2006) proposed a likelihood ratio test of separability for replicated multivariate data that examines whether the covariance matrix is a Kronecker product of two matrices of smaller dimensions. Mitchell et al. (2005) extended this approach to test the separability of spatio-temporal covariances by partitioning the observations into approximate replicates. After Scaccia and Martin (2005) presented tests of symmetry and separability for lattice processes using the periodogram, Li et al. (2007) developed a unified testing framework for both separability and symmetry by constructing contrasts of covariances at selected spatio-temporal lags.

Existing methods mostly build separability and symmetry tests based on given spatial and temporal lag sets, and work on visualization is sparse. In this paper, we propose formulating separability and symmetry test functions as functional data, and provide a visualization tool to illustrate these properties. These test functions are functions of temporal lags and constructed from the cross-covariances of each pair of the time series observed at spatial locations. With the obtained functional data of test functions, we develop rank-based testing procedures that are model-free with data depth-based functional data ranking techniques. Because the test statistics are rank-based, the results are more robust to outliers, which may come from erroneous measurements of the variables under study or poor estimates of the covariances due to limited sample size.

The rest of the paper is organized as follows. In §2, we introduce separability and symmetry test functions, which are used to visualize the covariance property of separability and symmetry, and provide the detailed rank-based testing procedure. In §3, we illustrate the visualization tools for simulated data and demonstrate the rank-based test for various spatio-temporal covariance models. In §4, we apply our methods to two real datasets, wind speed from monitoring stations, and surface temperatures and wind speed from numerical model outputs, to study their covariance structures.

2 Methodology

2.1 Non-separable covariance functions

Gneiting (2002) introduced the following rich family of non-separable spatio-temporal covariance functions, which includes separable ones as a special case,

[TABLE]

where $\sigma^{2}$ is the variance, $a,c\geq 0$ determine the temporal and spatial range, $\alpha,\gamma\in[0,1]$ control the temporal and spatial smoothness, and $\beta\in[0,1]$ is the spatio-temporal interaction. The covariance function is separable when $\beta=0$, and a larger $\beta$ is associated with a more non-separable model.

Cressie and Huang (1999) also proposed a non-separable model,

[TABLE]

and a possible corresponding separable model was discussed in Mitchell et al. (2005) as,

[TABLE]

Later, Rodrigues and Diggle (2010) further defined the properties of positive and negative non-separability and showed that the two aforementioned models are both positively non-separable. To investigate models of negative non-separability, we consider the one from Cesare et al. (2001),

[TABLE]

where $k_{1}>0$, $k_{2},k_{3}\geq 0$, and $C_{s}(\cdot)$ and $C_{t}(\cdot)$ are valid spatial and temporal covariance functions, respectively.

Table 1 summarizes the specific models used to illustrate our methods, after we plugged in the appropriate parameter values, chose suitable building covariance functions, and applied some trivial transformations, as necessary.

2.2 Visualization for separability

Noting that separable covariance functions require $C(\mathbf{h},u)=C(\mathbf{h},0)C(\mathbf{0},u)/C(\mathbf{0},0)$, for any spatial lag $\mathbf{h}$ and temporal lag $u$, or equivalently, $C(\mathbf{h},u)/C(\mathbf{h},0)$ remains a constant $C(\mathbf{0},u)/C(\mathbf{0},0)$ for any given $\mathbf{h}$. Therefore, we introduce the separability test function defined below.

Definition 1**.**

Given a valid spatio-temporal covariance function $C(\mathbf{h},u)$, the separability test function $f_{\mathbf{h}}(u)$ is a function of temporal lag $u$ for any spatial lag $\mathbf{h}$, defined as

[TABLE]

By Definition 1, $f_{\mathbf{h}}(u)$ is [math] for any $u$ and $\mathbf{h}$ if $C(\mathbf{h},u)$ is separable, and moves away from zero for non-separable models. Suppose $f_{\mathbf{h}_{1}}(u),\ldots,f_{\mathbf{h}_{n}}(u)$ are a set of test functions for pairs of locations with spatial lags $\mathbf{h}_{1},\ldots,\mathbf{h}_{n}$. We use functional boxplots (Sun and Genton, 2011) to visualize these separability test functions. The motivation is that functional boxplots can characterize the most representative functional realization, or the functional median, and remove the outliers, which are often involved in the datasets, especially when we are going to use sample covariances to estimate the test functions. Functional boxplots order the functional observations by their functional data depths, identify the observation with the largest depth value as the functional median and construct the $50\%$ central region that covers half of the data with largest depth values. Although depth values in the functional boxplot are calculated by the modified band depth proposed by (López-Pintado and Romo, 2009) by default, it is also possible to choose other depth notions developed for functional data, such as integrated data depth (Fraiman and Muniz, 2001), half-region depth (López-Pintado and Romo, 2011), and extremal depth (Narisetty and Nair, 2016). Outliers are detected by the $1.5$ times the $50\%$ central region empirical rule. The factor $1.5$ can be adjusted for outlier detection proposes in an adjusted functional boxplot, but the adjustment is not necessary for visualizing functional data (Sun and Genton (2012a)). Therefore, it is worth pointing out that we use the functional boxplot for visualizing covariance properties in this paper rather than detecting outliers. Alternatively, contour plots are the traditional way to visualize the pattern of covariances with respect to different spatial and temporal lags. We show next that the functional boxplots provide a much better visualization and interpretation of the separability than the contour plots.

Firstly, we use functional boxplots of the separability test functions and contour plots of the covariances to visualize the Gneiting model $C^{G}(\mathbf{h},u;\beta)$ for different values of $\beta$. The spatial locations are on a $4\times 4$ regular grid within the unit square $[0,1]\times[0,1]$. The top panels in Figure 1 show the functional boxplots applied to the separability test functions for the Gneiting model $C^{G}(\mathbf{h},u;\beta)$ with different values of $\beta$. Each curve is associated with a specific pair of spatial locations. The central region (magenta area) and the median (black curve) move away from zero as $\beta$ increases. The bottom panels in Figure 1 show the contours of the covariances for different spatial and temporal lags. In the separable case, the covariances at non-zero temporal lags should decay more rapidly as the spatial lag increases. We do see this trend in the contour plots; however, the separability is much more obvious in the functional boxplot.

Similarly, we draw the same plots to compare the Cressie and Huang separable model $C^{CH}_{sep}(\mathbf{h},u)$ and the Cressie and Huang non-separable model $C^{CH}_{nsep}(\mathbf{h},u)$, and the Cesare model $C^{C}(\mathbf{h},u;k)$ for different values of $k$. The results are shown in Figures 2 and 3. The functional boxplots clearly show the separability or non-separability of each test function, and, if non-separable, the positive or negative non-separability. While the contour plots show different patterns and the separable covariance decays much faster than non-separable covariances at non-zero temporal lags along an increasing spatial lag, the indication for separability is still unclear.

2.3 Visualization and rank-based testing procedure for real data

In §2.1–2.2, we discussed the visualization of true spatio-temporal covariance functions, and we saw that the functional boxplots produce intuitive indications of separability in the test functions. However, the true covariance model is unknown in practice, and we must use sample covariance estimators to estimate the test functions without assuming any parametric model to fit. For each pair of spatial locations, we use the two time series to estimate all the components in the separability test function $f_{\mathbf{h}}(u)$ by sample covariances. In other words, if we have observations $Z(\mathbf{s},t)$ at $p$ time points $1,\ldots,p$, at two locations $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ with a spatial lag $\mathbf{h}=\mathbf{s}_{2}-\mathbf{s}_{1}$, then the test function is estimated by

[TABLE]

where

[TABLE]

and other quantities are calculated similarly.

With these estimated separability test functions, we also develop a rank-based test for the separability hypothesis. López-Pintado and Romo (2009) proposed a test for functional data to decide whether two sets of samples are from the same distribution. It can be viewed as a functional data generalization of the univariate Wilcoxon rank test and has many applications. For example, Sun and Genton (2012b) proposed a robust functional analysis of variance (ANOVA) model for testing whether the differences in functional data observed from multiple groups are significant. The main idea is to see whether the positions of the samples from the two sets are similar, with respect to a reference distribution. We follow a similar procedure to test the null hypothesis, separability, by comparing the given observations to a simulated dataset under $H_{0}$: separability. In this rank-based test, two additional datasets are required. One follows the distribution under the null hypothesis $H_{0}$, and the other is a reference dataset following the distribution under the null hypothesis $H_{0}$ or the same distribution as the observations. We propose to generate these two datasets under $H_{0}$ by simulating data from a constructed separable covariance function $\widehat{C}^{H_{0}}$, for which we estimate the covariance at observed spatial and temporal lags by sample covariances, respectively, and let $\widehat{C}^{H_{0}}(\mathbf{h},u)=\widehat{C}(\mathbf{h},0)\widehat{C}(0,u)$. After we get all estimated separability test functions, the hypothesis test can be performed to compare the two populations. More details can be found in the paper by López-Pintado and Romo (2009), and we briefly explain the procedure as follows:

1. Step 1

Estimate the test functions for the observations, denoted by $\widehat{f}_{\mathbf{h}_{i}}(u)$, $i=1,\ldots,n$. 2. Step 2

Generate two sets of spatio-temporal datasets from the constructed separable covariance function $C^{H_{0}}(\mathbf{h},u)$ and estimate the test functions, denoted by $\widehat{f}^{0}_{\mathbf{h}_{j}}(u)$, $j=1,\ldots,m$ and $\widehat{f}^{R}_{\mathbf{h}_{l}}(u)$, $l=1,\ldots,r$. 3. Step 3

For each $i$, combine $\widehat{f}_{\mathbf{h}_{i}}(u)$ and all $\widehat{f}^{R}_{\mathbf{h}_{l}}(u)$ for $l=1,\ldots,r$, and compute their functional band depths and modified band depths (López-Pintado and Romo, 2009). Order these functions by their band depth values, and use modified band depths to separate ties. Let $r_{i}$ be the proportion of functions $\widehat{f}^{R}_{\mathbf{h}_{l}}(u)$ with smaller depth values than $\widehat{f}_{\mathbf{h}_{i}}(u)$. 4. Step 4

Repeat Step 3 for each $\widehat{f}^{0}_{\mathbf{h}_{j}}(u)$, and let $r^{\prime}_{j}$ be the proportion of functions $\widehat{f}^{R}_{\mathbf{h}_{l}}(u)$ with smaller depth values than $\widehat{f}_{\mathbf{h}_{i}}(u)$. 5. Step 5

Order $r_{1},\ldots,r_{n},r^{\prime}_{1},\ldots,r^{\prime}_{m}$ from smallest to largest, and assume the ranks of $r_{1},\ldots,r_{n}$ are $q_{1},\ldots,q_{n}$. Then the test statistics is $W=\sum_{i=1}^{n}q_{i}$.

Critical values can be determined from the asymptotic distribution of the test statistics, which is the sum of $n$ numbers that are randomly chosen from $1,\ldots,n+m$. In practice, however, we found that the distribution with a limited sample size is not close enough to the asymptotic one. Hence, we use the bootstrap to calculate the critical values. Under $H_{0}$, we generate $b$ sets of spatio-temporal datasets from $\widehat{C}^{H_{0}}(\mathbf{h},u)$. We repeat Steps 1–5 and obtain the test statistics for each dataset, which provides an approximate distribution under the null hypothesis. We examine the performance of this technique through simulations in §3, where we use $b=100$.

2.4 Symmetric covariance functions and visualization

Gneiting et al. (2006) proposed the following family of asymmetric covariance models that include symmetric ones as a special case,

[TABLE]

where $(\cdot)_{+}=\max(0,\cdot)$, $h_{1}$ is the first component of the spatial lag $\mathbf{h}$, $\lambda$ determines the symmetry, $v$ is the velocity along the direction of $h_{1}$, $a$ and $c$ correspond to the temporal and spatial ranges, $\alpha$ corresponds to smoothness, and $\beta$ corresponds to the separability. After plugging in appropriate parameter values, we get the following covariance model to study symmetry, referred as the Gneiting model $C^{G^{\prime}}(\mathbf{h},u;\lambda)$.

[TABLE]

We know that a fully symmetric covariance function satisfies $C(\mathbf{h},u)=C(\mathbf{h},-u)$. Analogously to Definition 1, we introduce the symmetry test functions with definition below.

Definition 2**.**

Given a valid spatio-temporal covariance model $C(\mathbf{h},u)$, the symmetry test function $g_{\mathbf{h}}(u)$ is a function of temporal lag $u$ for any spatial lag $\mathbf{h}$, defined as

[TABLE]

It is easy to see that $g_{\mathbf{h}}(u)$ remains zero when the underlying covariance model is fully symmetric. An illustration is given in Figure 4 for the Gneiting model $C^{G^{\prime}}$ with different values of $\lambda$, where we see the symmetry test function $g_{\mathbf{h}}(u)$ moves away from zero as $\lambda$ increases. Contour plots of the covariances are also shown in the bottom panels of Figure 4. Since the covariance relies on the temporal lag $u$ and both the first and second components $h_{1},h_{2}$ of the spatial lag $\mathbf{h}$, we fix $h_{2}=0$, and draw the contours for different values of $u$ and $h_{1}$. It is easy to observe that the covariances are clearly asymmetric when $\lambda>0$. For real data, the estimated symmetry test functions are used in the functional boxplots, and a similar rank-based testing procedure described in §2.3 is developed to test the significance of the asymmetry.

3 Simulation

3.1 Simulation design

Spatio-temporal data collected from monitoring sites are often sparse in space but dense in time. In this simulation, we consider data from $4\times 4=16$ regularly spaced locations in the unit square $[0,1]\times[0,1]$, and $p=2,000$ (a comparable number to many real applications) equally spaced time points.

Noting that the correlation between observations with temporal lags greater than $100$ is negligible in each model, we generate our spatio-temporal data sequentially in time. Let $\mathbf{Z}_{i}$ be a matrix of size $16\times 100$, indicating the spatio-temporal observations from time point $(i-1)\times 100+1$ to $i\times 100$ at the $16$ locations, for $i=1,\ldots,20$. We first generate $\mathbf{Z}_{1}$, and then generate $\mathbf{Z}_{i}$, $i=2,\ldots,20$, only conditioning on $\mathbf{Z}_{i-1}$. Thus the entire generated spatio-temporal data is $\mathbf{Z}=(\mathbf{Z}_{1},\mathbf{Z}_{2},\ldots,\mathbf{Z}_{20})$.

3.2 Visualization and assessment for separability

The Gneiting model $C^{G}(\mathbf{h},u;\beta)$ with five different values of $\beta=0,0.25$, $0.5$, $0.75$, and $1$ is used to generate observations $\mathbf{Z}$. Three examples of the estimated separability test functions are shown in Figure 5. We see that, in the functional boxplot, the distance between the median (black line) and zero (green dashes) indicates the degree of separability. As $\beta$ increases, the majority of the separability test functions move away from zero. This gives us a clear visualization of the separability.

To formally test the separability and examine the performance of the proposed rank-based testing procedure in §2.3, we show the Type I error and power for simulated data over $100$ simulations. The results are shown in Table 2. We see the Type I error is close to the nominal level of the hypothesis test, and the power for non-separable models increases with an increasing $\beta$. The power is quite close to one when $\beta$ approaches $1$.

In addition, we also take a look at the Cressie-Huang models $C^{CH}_{nsep}$ and $C^{CH}_{sep}$. Estimated separability test functions from one simulated dataset are shown in Figure 6. We see the majority of the separability test functions from the non-separable model move away from zero. Similarly, a rank-based test is used to examine the performance; the results are shown in Table 3.

3.3 Visualization and assessment for symmetry

We choose five different values of $\lambda=0$, $0.025$, $0.05$, $0.075$, and $0.1$ in the Gneiting model $C^{G^{\prime}}(\mathbf{h},u;\lambda)$ to generate observations $\mathbf{Z}$. The reason for choosing relatively small values for $\lambda$ is that the power of the rank-based test is already strong enough for small values of $\lambda$. Three examples of the estimated symmetry test functions are shown in Figure 7. Since the chosen values of $\lambda$ are very small, the deviation of the symmetry test functions from zero in our most extreme case when $\lambda=0.1$ is still subtle; see Figure 4. Although the functional medians of the estimated test functions do not show obvious changes for the different values of $\lambda$, we can still see that the majority of the symmetry test functions move away from zero as $\lambda$ increases. Furthermore, the results of the rank-based test in Table 4 show that the power is very strong, even when $\lambda$ is as small as $0.01$, and the Type I errors are all close to the nominal levels.

4 Applications

4.1 Wind speed

Wind speed is an important atmospheric variable in weather forecasting and many other environmental applications. One source of wind speed data is from monitoring stations. In this section, we use our proposed method to visualize and assess properties of covariances in the wind speed from ten monitoring stations in the northwestern U.S.A. The locations are shown in Figure 8, where we divide them into two regions. Coastal sites are in green, and inland sites are in orange. We analyze the hourly data observed at each station from 2012 to 2013 and visualize the spatio-temporal covariance structure from two seasons, summer (June to August) and winter (December to February).

The plots of the estimated separability test functions are shown in Figure 9, and the p-value results of the rank-based test are given in the titles. We see the separability test functions are away from zero in all cases, and all the p-values for testing separability are much smaller than $0.05$, indicating rejection of separability at $\alpha=0.05$. When computing the p-value, we choose $10,000$ replicates from the bootstrap to approximate the distribution of the test statistics under the null hypothesis. Since the covariance functions for all the four cases are significantly non-separable, we next investigate their properties of symmetry. The symmetry test functions as well as the p-values of the test for different cases are shown in Figure 10. The functional boxplot for the coastal region in winter suggests the strongest symmetry, and its corresponding p-value is indeed the largest, while p-values for the other cases are much smaller than $0.05$. One possible reason might be that there is no prevailing wind direction during the winter in the coastal region.

4.2 Climate model outputs

In climate studies, numerical models are used to simulate many important variables by solving a set of dynamic equations. Two types of models are available: General Circulation Model (GCM) is used to describe the dynamic system of the entire earth, with a global coverage but a comparatively low resolution; however, Regional Climate Model (RCM) with GCM outputs as their boundary conditions is used to simulate locally high-resolution data. Different combinations of GCMs and RCMs may lead to different results. In this section, we apply our visualization and assessment tools to data based on four combinations of two GCMs, GFDL and HADCM3, and two RCMs, ECP2 and HRM3, which are provided by the North American Regional Climate Change Assessment Program (NARCCAP) (Mearns et al., 2012). More details about these numerical models can be found at http://www.narccap.ucar.edu/. We focus on a small region which consists of $16$ spatial locations in the North Atlantic Ocean (as shown in Figure 11) and study the daily surface temperature and wind speed from 1976 to 1979.

We first look at the mean of the surface temperature across the four years at each location, which is shown in Figure 12. We see different GCMs as boundary conditions give different patterns of surface temperatures, while different RCMs show similar results. We then apply our visualization and assessment tools to the covariance structures of these four GCM and RCM combinations. To eliminate seasonality, we remove the monthly mean from the daily temperatures. The estimated separability test functions and the p-values of the rank-based test for separability are shown in Figure 13. With all the p-values greater than $0.05$, there is no significant evidence that the covariances of these four cases are non-separable at the $5\%$ significance level. However, the output from GFDL and HRM3 shows stronger non-separability and the separability is rejected at the $10\%$ significance level with a p-value $0.09$.

Next, we look at the daily wind speed in the same area. The mean wind speeds at these locations are all close to zero and do not show much difference. Again, we remove the monthly mean from the daily wind speed and then apply our visualization and assessment methods. The estimated separability test functions are shown in Figure 14, where the medians are away from zero. These test functions show large variability in the functional boxplot, and all the p-values from the rank-based tests are as small as zero, indicating strong separability. To further investigate the symmetry, we plot the estimated symmetry test functions; the illustration is given in Figure 15. All the symmetry test functions are also very far from zero, and we find that the p-values of the test for symmetry are zero. Therefore, we conclude that the covariance of wind speed is neither separable nor symmetric.

In this application, the covariance of the daily surface temperatures tends to be separable and thus also symmetric; for daily wind speed, the covariance is non-separable and asymmetric. Our analysis suggests that, for the region we have considered, the daily surface temperatures produced by the climate models do not show significant space-time interaction, whereas the daily wind speeds clearly show asymmetric space-time interaction.

5 Discussion

We presented a functional data analysis approach to visualizing and assessing spatio-temporal covariance properties. The proposed method is suitable for visualizing the separability and the symmetry of a stationary spatio-temporal process. The p-value calculated in the rank-based test can be used to measure the degree of separability or symmetry. We illustrated our methods using various classes of spatio-temporal covariance models, and our simulations demonstrated the good performance of our proposed methods. In the practical applications, we applied our method to temperature and wind speed data from either monitoring sites or climate model outputs, and illustrated how to interpret the separability and symmetry in these different scenarios.

However, we have only considered spatio-temporal datasets collected from a small number of spatial locations. Because we need to estimate the test functions for each pair of locations, computation will become an issue as the number of locations grows. One possible solution is to divide the region of interest into smaller subregions, and estimate all the test functions separately for each of these subregions. Then, the obtained test functions can be combined for the visualization and the rank-based test. This would significantly reduce the number of test functions, and each small subregion is more likely to be stationary. However, future research is required to assess the performance of this approach.

Acknowledgement

We wish to thank the North American Regional Climate Change Assessment Program (NARCCAP) for providing the data used in this paper. NARCCAP is funded by the National Science Foundation (NSF), the U.S. Department of Energy (DoE), the National Oceanic and Atmospheric Administration (NOAA), and the U.S. Environmental Protection Agency Office of Research and Development (EPA).