{'fraction_non_alphanumeric': 0.038661878417289464, 'fraction_numerical': 0.016952991038208937, 'mean_word_length': 4.405736271423574, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonparametric Spectral Methods for Multivariate Spatial and Spatial-Temporal Data

We propose computationally efficient methods for estimating stationary multivariate spatial and spatial-temporal spectra from incomplete gridded data. The methods are iterative and rely on successive imputation of data and updating of model estimates. Imputations are done according to a periodic model on an expanded domain. The periodicity of the imputations is a key feature that reduces edge effects in the periodogram and is facilitated by efficient circulant embedding techniques. In addition, we describe efficient methods for decomposing the estimated cross spectral density function into a linear model of coregionalization plus a residual process. The methods are applied to two storm datasets, one of which is from Hurricane Florence, which struck the souteastern United States in September 2018. The application demonstrates how fitted models from different datasets can be compared, and how the methods are computationally feasible on datasets with more than 200,000 total observations.


Introduction
There are several exciting scientific campaigns to produce and observe multivariate data that vary over space and time. For example, large climate centers such as the National Center for Atmospheric Research (NCAR) produce high resolution simulations of the Earth system. These models include dozens of variables evolving in concert over space and time. The National Aeronautics and Space Administration (NASA) and the National Oceanic and Atmospheric Administration (NOAA) have deployed numerous satellites that collect observations of the Earth surface and atmosphere. NASA and NOAA recently launched a pair of satellites in its ongoing geostationary operational environmental satellite (GOES) program, the GOES-16 and GOES-17 spacecraft, that sit in geostationary orbit, continually monitoring light reflected by the Earth surface and atmosphere in 16 separate wavelength bands. In turn, the raw data are processed to produce dozens of physically relevant variables, such as atmospheric water vapor content and land surface temperature.
A statistical framework for analyzing these data should, at a minimum, include (1) sufficiently flexible statistical models capable of capturing complex multivariate dependencies, and (2) computationally efficient tools for estimating the models from data. Genton and Kleiber (2015) provide a thorough review of existing modeling and estimation frameworks. Kleiber (2017) conducted a theoretical analysis of a number of multivariate spatial models from the literature and concluded that many of them impose oversimplistic restrictions on the coherence between pairs of variables. For example, separable and kernel convolution multivariate spatial models (Majumdar and Gelfand, 2007) have constant coherence. In the linear model of coregionalization (LMC) (Banerjee et al., 2014), when component processes have spectral densities that decay with different rates, the coherence always converges to a non-zero constant as the frequency increases.
It is desirable to have models that do not impose such restrictions, for example, to have a model class that allows the coherence to decay to zero as frequency increases. Kleiber (2017) showed that the multivariate Matérn models (Gneiting et al., 2010;Apanasovich et al., 2012) do possess such flexiblity. However, likelihood-based estimation of parameters in multivariate Matérn models is not computationally feasible for the massive datasets mentioned in our opening paragraph. Computational issues come from two sources. The first is that the data size prohibits formation and factoring of the covariance matrix, necessary operations for evaluating the likelihood function. There exist promising approximations that could in principle apply to the multivariate case, such as Vecchia's likelihood approximation (Vecchia, 1988;Datta et al., 2016;Guinness, 2018;Katzfuss and Guinness, 2017) and hierarchical matrix approximations (Saibaba et al., 2012;Ambikasaran et al., 2016;Litvinenko et al., 2017), but these approximations are untested for multivariate models, and it is not yet clear how to optimally implement them, and how their performance will compare to univariate cases. Such work would be a welcome development to the literature. However, a second and perhaps more serious computational problem, even for fast approximate methods, is that multivariate spatial models contain many more parameters, usually on the order of the square of the number of components. The large number of parameters poses a serious issue for optimization of the (approximate) likelihood function, which typically requires many iterations until convergence over a large parameter space. We demonstrate this problem with a simple example in Section 4 Nonparametric spectral methods offer the potential of addressing both the modeling and computational demands for multivariate spatial and spatial-temporal data. In this framework, a discrete Fourier transform (DFT) is applied to each multivariate component, and then the periodogram vectors are smoothed over frequency to generate a nonparametric estimate of the cross spectral density matrices. This approach is reasonably flexible from a modeling standpoint, since the only assumption is that the cross spectral density matrices vary smoothly with frequency, and it is computationally inexpensive, since we can evaluate the DFTs quickly with FFT algorithms. However, nonparametric spectral methods have limited applicability since they typically apply only to complete, gridded data on rectangular domains. Moreover, even if we have such data, edge effects can introduce severe biases in the spatial and especially the spatial-temporal cases (Guyon, 1982). Tapering (Dahlhaus and Künsch, 1987) and differencing (Lim and Stein, 2008) have the potential to reduce edge effects, but differencing is not possible when there are missing values, and Guinness (2017) showed that tapering can be ineffective under certain missingness scenarios and suboptimal compared to imputation-based methods.
This paper provides methods that addresses edge effects and extend the applicability of nonparametric spectral methods to incomplete gridded multivariate spatial-temporal data. This is achieved by leveraging the simple, yet powerful technique of periodic imputation introduced for the univariate case in Guinness (2017). We also introduce a method for decomposing the estimated cross spectral density function into an LMC cross spectral density function plus a residual cross spectral density function, which can be used as tool for exploring the spatial-temporal variation in the data. We apply the methods to compare the multivariate spatial-temporal covariances of data from an ordinary storm and from Hurricane Florence, demonstrating that these methods can computationally feasibly estimate multivariate spatial-temporal models with more than 200,000 obervations, and how the estimates can be interpreted to distinguish among models.

Model and Notation
Let x ∈ Z d be a location on the d-dimensional integer lattice, and let Y (x) = (Y 1 (x), . . . , Y p (x) T ) ∈ R p be a p-variate mean-zero random process on Z d . The process Y is stationary if  (20,20). Dark gray represents locations for observations in U , light gray and white represent locations for observations in V . Note that locations of missing data on X a need not be the same in each component.
as the covariance between the process at pairs of locations separated by lag h = (h 1 , . . . , h d ). Cramer's Theorem (Wackernagel, 1998, cf.) states that K has the representation where i = √ −1, and ω · h = ω 1 h 1 + · · · + ω d h d . The matrix-valued function f (ω) ∈ C p×p has (j, k) entry f jk (ω), and must be a Hermitian positive definite matrix in order that K(h) is a real and positive definite function. We call f (ω) the cross-spectral density (CSD) function.
Let j 1 , . . . , j n ∈ {1, . . . , p} and x 1 , . . . , x n be a sequence of locations and U = (Y j 1 (x 1 ), . . . , Y jn (x n )) be the vector of observations. In this notation, n is the total number of individual observations, and we may observe between 0 and p of the components at any location. Because the process is stationary, we can assume without loss of generality that each location x k falls "northeast" of the location (1, . . . , 1) and "southwest" of some location a = (a 1 , . . . , a d ). Formally we say that x k falls in the observation lattice X a , where X a = {(r 1 , . . . , r d ) : r k ∈ {1, . . . , a k } for every k}. ( The methods in this paper involve imputing the multivariate process onto a domain larger than the observation lattice X a . To this end, define b = (b 1 , . . . , b d ) with b j > a j for each j, and the corresponding embedding lattice X b . Let (U, V ) denote a complete vector of observations on X b , that is (U, V ) contains each of the p elements of Y (x) at each location x ∈ X b , ordered with the observations U first, then the missing values V second. Figure 1 provides a visual description of these definitions.
In the next subsection, we describe methods for iteratively imputing the missing values V on X b . However, we first describe how we obtain a CSD estimate given a complete set of observations (U, V ) on X b . We propose the following routine for obtaining a CSD estimate f [(U, V )]: The symbol F b refers to the set of Fourier frequencies on a grid of size b, and * is complex conjugate. Each computation written in terms of ω is computed for every ω ∈ F b . If the equation is written in terms of j, it is computed for each j = 1, . . . , p, and likewise, equations written in terms of j and k are computed also for each k = 1, . . . , p. In (3), we compute the DFT of the complete dataset on X b . In (4), we maximize Whittle's loglikelihood approximation for some parametric spectral density f θ . In Sections 4 and 5 we use the quasi Matérn spectral density proposed in Guinness and Fuentes (2017). In (5), we smooth the sample covariances of the normalized DFT entries with kernel α, and multiply by the square roots of the estimated parametric spectral densities. Guinness (2017) suggested the parametric normalization in (5) as a method for reducing smoothing bias, which allows for the use of wider smoother kernels, resulting in estimates with smaller variance. If instead f θ is assumed to be constant, the routine for obtaining a CSD estimate reduces to a standard routine described, for example, in Brockwell and Davis (2006) for multivariate time series. Nonparametric or semiparametric CSD estimation methods are desirable because the estimated CSD function f (ω) is automatically positive definite at every frequency if the support of α includes at least p Fourier frequencies in F b . However, when the dimension d of the field is larger than 1, substantial bias is imparted by edge effects (Guyon, 1982;Lim and Stein, 2008). In the next subsection, we describe iterative imputation methods for reducing edge effect bias.

Iterative Periodic Imputation Methods
Lee and Zhu (2009) proposed a method for estimating a univariate time series spectrum from incomplete data by iteratively imputing the missing data on the observation domain. Guinness (2017) showed that when d > 1, univariate spectral density estimates can be improved dramatically when data are iteratively imputed onto embedding lattice X b with a conditional simulation from a periodic model on X b . Here, we propose an extension of the periodic imputation method to the estimation of multivariate cross spectral densities. At iteration , the missing values V are imputed according to a Gaussian process model with a covariance function that is periodic on X b , where f ( ) is the CSD function at the th iteration, and m = b 1 · · · b d is the total number of locations in X b . With the complete dataset (U, V ), we update the spectrum estimate with and iterate. After a burn-in period of B iterations, we update with a weighted average of the previous and current spectrum, and we monitor changes in f ( ) for convergence. The full estimation algorithm is as follows: Estimation Algorithm: 1. Simulate V given U from Gaussian process with covariance R ( ) .

Given
For > B 1. Simulate V given U from Gaussian process with covariance R ( ) .

Given
Drawing Gaussian V given U is the main computational challenge for the iterative methods. All of the other operations are either pointwise multiplications, divisions, DFTs, or convolutions, which can be computed with a DFT. Drawing V given U involves an unconditional simulation from R ( ) on X b , which can be done with DFTs and pointwise multiplications. It also involves solving a linear system with the covariance matrix for U under R ( ) . This is the most demanding step. We use preconditioned conjugate gradient (PCG), which is an iterative method for solving positive definite linear systems (Saad, 2003). The forward multiplications in PCG can be done in O(pm log m + p 2 m) with circulant embedding techniques (Wood and Chan, 1994). We use a preconditioner based on Vecchia's approximation (Vecchia, 1988) that can be computed in O(pm) time.

Factor Decompositions
The CSD f (ω) has an interpretation as the covariance matrix for Y(ω), the variation in Y at frequency ω. This is admittedly difficult to communicate to a diverse audience. One solution is to simply report the estimated cross covariance function, the inverse DFT of f . To further facilitate interpretability of a p-dimensional multivariate model, it can be useful to consider lower-dimensional representations for the CSD functions. To this end, we propose a factor model consisting of the linear model of coregionalization (LMC) plus a residual process Z(x), where A j ∈ R p , W 1 , . . . , W J are univariate process components of the LMC, independent of each other and of Z, a p-variate spatial process with correlated components. The CSD function for this process is where g j (ω) > 0 is the univariate spectral density for W j , and h is the CSD function for Z.
There is a clear lack of identifiability in assigning power from f (ω) to either A j A T j g j (ω) or h(ω); for example, any model written as (8) can be rewritten as For this reason, we argue that it is not meaningful to ask for an "optimal" or "true" value of J. We take the position here that since the purpose of the factor model is to explain variation in the response with a lower-dimensional model, we should assign as much power as possible to the factor term over a range of values of J. Then the amount of variation explained by the factor term can be used as one way to distinguish among models.

Profiled Optimization
After the CSD function f has been estimated using the proposed periodic imputation methods, we aim to select A 1 , . . . , A J and g 1 , . . . , g J in order to minimize the amount of power assigned to the residual process Z. To quantify the power, we consider the sum of the variances of the residual processes, which can be written in terms of the sum of traces of the residual CSD function, subject to g j (ω) > 0 and f (ω) − A j A T j g j (ω) nonnegative definite for each ω. Even for J = 1, minimizing the loss function over all g j (ω) is demanding because there are potentially thousands or millions individual frequencies ω ∈ F b . This intractability has motivated a theoretical study to find closed-form minimizers over g j (ω) for fixed A = (A 1 , . . . , A J ). This allows us to solve the global optimization problem by profiling out g j (ω) and numerically minimizing the sum-of-variance criterion over A. Further, we can assume without loss of generality that each A j lies on the unit sphere because the radius of A j can be absorbed into g j (ω), so the numerical optimization is over a J(p − 1) space.
The sum-of-variance criterion can be rewritten as The first term on the right does not depend on g j or A j , and so minimizing the criterion corresponds to maximizing the second term subject to g j (ω) > 0 and h(ω) nonnegative definite for each ω.
For fixed A with k A 2 jk = A j 2 = 1, this corresponds to maximizing ω j g j (ω). Since the nonnegative definiteness criterion applies separately to each ω, we can maximize j g j (ω) separately for each ω.
The theorem gives a result for finding the minimizers over g 1 (ω), . . . , g J (ω) for fixed A. Following the theorem are two corollaries that give specific results for the cases J = 1 and J = 2. To state the theorem, parameterize nonnegative g j (ω) as exp(c j ), and define the J × J matrices Further, define the Lagrangian function L(c 1 , . . . , c J , λ) = e c 1 + · · · e c J + λ det(C −1 − B).
Theorem 1. For fixed A = (A 1 , . . . , A J ), the minimizer of with respect to c 1 , . . . , c J , subject to h(ω) nonnegative definite, is either a solution to ∇L = 0 or for some k ∈ {1, . . . , J}, the minimizer of The theorem states that the minimizers either solve the gradient of the Lagrangian or are the solution to a subproblem. This arises because we must check the boundaries of the parameter space when using the method of Lagrange multipliers. We must reparameterize g j (ω) = exp(c j ) because it is possible that the maximizer of j g j (ω) subject to h(ω) nonnegative definite includes one or more negative g j (ω). The following two corollaries give concrete solutions for the cases J = 1 and J = 2.
According to Theorem 1, the solution is either the above, or the solution to one of the one-factor problems, establishing the corollary.

Conditional Expectation of Component Processes
The conditional expectations of the jth factor component given the data is The conditional expectation of Y (ω) given U is where W j is the DFT of W j . Therefore, in order to compute the conditional expectation of W j , we simply compute E(V |U ), then take the DFT of (U, E(V |U )) to obtain E(Y(ω)|U ), and then take the inverse DFT of We use the expected factors in Section 5 to explore the decomposition of the estimated spectrum from a thunderstorm dataset.

Simulations
The purpose of the simulation study is to demonstrate some of the computational issues with maximum likelihood estimation of the multivariate Matérn parameters, and show how our proposed methods are capable of producing fast and accurate estimates of the CSD function, even when we do not assume knowledge of the parametric form of the model.
M(r ; ν) is the Matérn function, written in terms of K ν , a modified Bessel function of the second kind, with parameter ν. The parameters σ 2 jj , α jj , and ν jj represent the marginal variance, the inverse range, and the smoothness of Y j . These marginal parameters must all be positive. To ensure positive definiteness of the multivariate Matérn, Gneiting et al. (2010) give conditions on the parameters in the bivariate Matérn model, and Apanasovich et al. (2012) give conditions on the parameters in a more general p-variate setting.
In the simulation study, we consider p = 2, 3, and 4. The inverse range and smoothness parameters are α jk = 0.25 and ν jk = 0.5 + 0.5(j + k − 2)/(2p − 2). To define the variance parameters, let β be a p × p matrix with β jk = 0.8 |j−k| . The variance parameters are σ jk = jk Γ(ν jk ) Γ(ν jk + 1) Γ(ν jj + 1)Γ(ν kk + 1) This particular model parameterization is from the parsimonious multivariate Matérn family (Gneiting et al., 2010). Table 1 presents the specific parameter values used when p = 3. Due to the symmetry constraints σ jk = σ kj , α jk = α kj , and ν jk = ν kj , the full multivariate Matérn has 3p(p + 1)/2 unique parameters. For p = 2, this is 6 parameters; for p = 3, 18 parameters; for p = 4, 30 parameters. A mentioned in Section 1, the large number of free parameters is one aspect of the computational difficulty of working with the multivariate Matérn model. We specify grid size (16, 16) and simulate 70 datasets for each value of p ∈ {2, 3, 4} with no missing values on the grid. To estimate multivariate Matérn parameters, we use the exact Gaussian likelihood function, maximized using the optim function in R with the default Nelder-Mead algorithm. The estimation procedure makes no assumptions that the true model is a member of the parsimonious multivariate Matérn family. We run the Nelder-Mead algorithm for 6000 iterations, stopping at each 1000 iterations to report estimation progress and elapsed time.
For our proposed periodic imputation methods, we take three choices of expansion parameter τ ∈ {1.00, 1.25, 1.50}, and we use the parametric variant of the iterative algorithm, run for 50 burn-in iterations, and an averaging tolerance of ε = 0.01. We use a Gaussian smoothing kernel, The parameter σ 2 can be profiled out, so at each iteration, Whittle's likelihood is optimized over the two-dimensional parameter space (α, ν) for each of the p components.
To evaluate the estimators, we consider a spectral norm criterion where λ max (A) is the largest absolute eigenvalue of matrix A, f is the true CSD function, and f is the estimate.
The results of the simulation study are given in Table 2. When no periodic embedding is used (τ = 1.0), the estimates are computed extremely fast but are poor relative to the other estimators. The poor performance when τ = 1.0 highlights the fact that spectral estimators based on the periodogram of the observed data can suffer from severe edge effects. The periodically embedded estimators are nearly as accurate or more accurate than maximum likelihood, even after 6000 Nelder-Mead iterations. When p = 2 the periodic imputation estimators are slightly worse in terms of spectral norm, whereas the periodic imputation estimators are better when p = 3 and p = 4, even though no knowledge of the true CSD functions is assumed. Further, the periodic imputation methods are much faster. They converge on average in less than 25 seconds in every case, whereas the maximum likelihood estimates converge in several thousand iterations and take on the order of one to three hours.
The simulations were conducted on a single node of the Cheyenne supercomputer managed by the National Center for Atmospheric Research. Each node has 36 2.3-GHz Intel Xeon E5-2697V4 processors, and each processor has 2 threads. One processor was used to manage the parallel tasks, and each of the remaining 35 processors analyzed two datasets in parallel on the two threads. The computing times in Table 2 should be understood in the context that only a single thread was used to analyze each dataset. We can expect up to, for example, a four-fold speedup for each task when running the code on a two core machine with four total threads.

Multivariate Spatial-Temporal Analysis of Storms
Launched in November 2016, the GOES-16 satellite now sits in geostationary orbit, allowing its advanced baseline imager (ABI) to capture images of the Earth's atmosphere in 16 wavelength bands at up to 500m resolution in space and up to 30 seconds in time. GOES 16 produces terabytes of data per day and provides a wealth of information about atmospheric processes. In this section, we perform a multivariate spatial-temporal analysis of images from two separate storms, one of which is an ordinary convective storm formed over Florida from July 22, 2018, while the other is from Hurricane Florence on September 14, 2018, the day it made landfall with the southeastern U.S. coastline. The analysis is demonstrative of the sort of comparisons that are possible with our proposed methods.
For each storm, we analyze 60 images from 4 wavelength bands separated by 1 minute. Band 1 has a native resolution of 1km, while the resolution of the other bands is 2km. We locally average the Band 1 data to coincide with the resolution of the other bands. Figure 3 contains images from Band 1 at the beginning, middle, and end of the hour. A stationary model is unrealistic for the entire scene, and so we subset the data to the black polygons in Figure 3, where the stationary assumption is more tenable. Figure 4 contains the subsetted images from the middle time point at the four wavelength bands. Table 3 contains information about the bands selected for the analyses. We choose one band in the visible spectrum since the visible spectrum bands are all highly correlated due to the grayscale color of clouds; Band 6 contains information about particle sizes; Band 7 can be used for ice detection; and Band 9 contains information about water vapor content in the mid-troposphere (NOAA, 2018) The ordinary storm has observation lattice size a = (55, 57, 60) in longitude, latitude, time, with 253,504 total observations; Florence has a = (47, 31, 60), with 167,804 total observations. For both storms, we specify expansion factor τ = 1.25, extending the lattice by 25% in each of the three dimensions, and we again use the quasi Matérn parametric spectral density. Before estimation, we subtract off the sample mean µ k from each component. We use B = 20 burn-in iterations, and we Wavelength Nickname Band 1 0.47 microns Blue Band Band 6 2.2 microns Cloud Particle Size Bands Band 7 3.9 microns Shortwave Window Band Band 9 6.9 microns Mid-Level Tropospheric Water Vapor Band Table 3: Information about bands used in multivariate spatial-temporal analysis. specify convergence tolerance parameter ε = 0.005. The ordinary storm CSD estimates converge after 6 averaging iterations, which took 3.40 hours total for burn-in and convergence. Florence converged after 15 iterations and took 1.43 hours. All computations are in the R programming language and run on an Intel Core i5-7200 CPU (2 cores, 4 threads at 2.50GHz) with 8GB memory. Vecchia's preconditioner is implemented in C++ with the Rcpp package (Eddelbuettel and François, 2011).
To visualize the fits, we include plots of the estimated real part of the coherences for the two storms in Figure 5. We see that for the ordinary storm, the coherence between Bands 1,6, and 7 persists throughout the frequency domain, while Band 9 has a slight negative coherence with the other bands throughout the domain. The situation is quite different for the Hurricane Florence data; Band 9 has a small amount of positive coherence with Bands 1 and 7 and is negatively coherent with Band 6. The coherences weaken at the highest frequencies. Bands 1, 6, and 7 are weakly coherent.
The total variation differs widely between the four bands; for example, in the ordinary storm, the estimated standard deviations are 41.15, 0.70, 0.01 and 0.20. These differences in variation are not necessarily a reflection of the relative importance of the bands, so instead of decomposing the estimated CSD functions f directly-which would undoubtedly return vectors A j attempting to Hurricane Florence (bottom row) at various frequencies ω = (ω 1 , ω 2 , ω 3 ). "Low" refers to ω j = 0, "mid" to ω j = 1/4, and "high" to ω j = 1/2. In each 3×3 group of coherence matrices, rows are corresond to latitudinal frequencies, columns to longitudinal frequencies.
Ordinary Storm  Table 4: Summary of factor decompositions for the two storm datasets.
explain variation in Band 1-we decompose the normalized CSD function where C jj (0) is estimated covariance function at lag zero, that is, the estimated variance. Table  1 provides a summary of the decompositions. For the ordinary storm, the J = 1 decomposition puts most of its weight on Bands 1, 6, and 7. When J = 2, A 1 is nearly unchanged, and A 2 explains variation in Band 9. For the Hurricane Florence data, when J = 1, A 1 points mostly in the direction of Bands 1, 6, and 9. When J = 2, A 1 changes little, and A 2 explains common variation in Bands 1 and 6. Lastly, in Figure 6, we plot the original ordinary storm data from t = 30, and the expected J = 1 and J = 2 factor representations of each band k, We can see that the one-factor representation ignores variation in Band 9 but captures much of the variation in Bands 1, 6, and 7. This is expected since these three bands are strongly coherent across frequencies. The two-factor captures some of the variation in Band 9 but still oversmooths.

Discussion
We have introduced simple, flexible, and computationally efficient methods for estimating stationary multivariate spatial-temporal models from incomplete gridded data. The methods rely on successive imputation of data onto an expanded lattice under a model that is periodic on the expanded lattice.
The simulation studies demonstrate that the periodic domain expansion is crucial for addressing edge effects; when no expansion is performed, the estimates of the spectrum are poor. The new estimates are competitive with maximum likelihood-though much faster-when there are p = 2 multivariate components, and the estimates are faster and more accurate than maximum likelihood with p > 2 components. We have argued that this arises both from the computational demand of evaluating the likelihood function, and of maximizing over the large number of parameters in the multivariate Matérn model. The paper describes a method for decomposing the estimated spectrum into a linear model of coregionalization plus a residual multivariate process. The decomposition is meant as an exploratory tool for understanding the variation in the data. Finding an optimal such decomposition proved to be an interesting topic, and we have provided some theoretical results that make the numerical search over possible decompositions feasible. The factor decomposition was applied to two storm datasets, where we found that two-component decompositions explained roughly half of the variation in the data.
Only one-and two-factor decompositions were pursued here. Higher-order decompositions are in principle computationally feasible using Theorem 1 and some tedious algebra, but we leave this problem for future work. We have also not explored the predictive capabilities of the fitted models. Prediction is an interesting and important topic in multivariate spatial-temporal models, see for example Zhang and Cai (2015), but we have decided to focus instead on exploring the model fits and how they differ for two types of storms.   Table 5 contains simulation results for the periodic imputation methods for all four bandwidths. Figure 7 contains estimated imaginary parts of the coherences for both storm datasets. Figure 8 shows Hurricane Florence data and estimated one-and two-factor representations of the data. Florence (bottom row) at various frequencies ω = (ω 1 , ω 2 , ω 3 ). "Low" refers to ω j = 0, "mid" to ω j = 1/4, and "high" to ω j = 1/2. In each 3 × 3 group of coherence matrices, rows are corresond to latitudinal frequencies, columns to longitudinal frequencies.

A Supplementary Tables and Figures
Figure 8: Hurricane Florence data from time t = 30 (top row), one factor representation (middle row), and two factor representation (bottom row). Factor representations use the expected values of W j given the data.