{'fraction_non_alphanumeric': 0.04249691674467674, 'fraction_numerical': 0.017087921068663723, 'mean_word_length': 4.355295615500915, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Levy-Dynamic Spatio-Temporal Process: Towards Big Data Analysis

In this era of big data, all scientific disciplines are evolving fast to cope up with the enormity of the available information. So is statistics, the queen of science. Big data are particularly relevant to spatio-temporal statistics, thanks to much-improved technology in satellite based remote sensing and Geographical Information Systems. However, none of the existing approaches seem to meet the simultaneous demand of reality emulation and cheap computation. In this article, with the Levy random fields as the starting point, e construct a new Bayesian nonparametric, nonstationary and nonseparable dynamic spatio- temporal model with the additional realistic property that the lagged spatio-temporal correlations converge to zero as the lag tends to infinity. Although our Bayesian model seems to be intricately structured and is variable-dimensional with respect to each time index, we are able to devise a fast and efficient parallel Markov Chain Monte Carlo (MCMC) algorithm for Bayesian inference. Our simulation experiment brings out quite encouraging performance from our Bayesian Levy-dynamic approach. We finally apply our Bayesian Levy-dynamic model and methods to a sea surface temperature dataset consisting of 139,300 data points in space and time. Although not big data in the true sense, this is a large and highly structured data by any standard. Even for this large and complex data, our parallel MCMC algorithm, implemented on 80 processors, generated 110,000 MCMC realizations from the Levy-dynamic posterior within a single day, and the resultant Bayesian posterior predictive analysis turned out to be encouraging. Thus, it is not unreasonable to expect that with significantly more computing resources, it is feasible to analyse terabytes of spatio-temporal data with our new model and methods.


Introduction
This is the era of "big data" and the scientific community, including the statistical community, is mesmerized by the sheer charm of the phrase! So much so that editors and reviewers of so-called reputed journals keep rejecting papers submitted to those journals on account of their perception of inapplicability of the papers' contribution to big data (in our experiences)! Such papers are often related to genetics and spatio-temporal statistics where complex dependence structures play the key roles. However, big data refers to at least one terabyte of data, and from that perspective, no sensible modeling approach to account for complex dependence in the underlying data-generating process can be feasible without supercomputing resources. Even with much less amount of such structured data, very powerful and well-maintained computing facilities are necessary, which are usually unavailable to individual researchers. The response of the statistical community to the big data challenge (and indeed the advices of the editors and reviewers!) is to simplify the model by ignoring most of the dependence structures, assuming linearity, and so on. One reviewer also very kindly advised not to develop sophisticated theories and methods, since in his opinion, the linearity assumption is almost always sufficient! Breaking up the data into as many parts as possible despite its highly dependent structure, implementing the model on broken-up sub-datasets in the embarrassingly parallel way treating the sub-datasets as independent, and then combining the results in some manner, is another response being adopted by the statistical community. Unfortunately, despite the existence of such magic tricks of the trade, much to the likings of the editors and reviewers, there does not seem to exist any model-based statistical work that analyzes any structured big data of the order of terabytes.
In this article, we confine ourselves to the area spatio-temporal statistics, and propose and develop a new, highly structured Bayesian nonparametric and dynamic spatio-temporal model based on Lévy random fields for analyzing reasonably large spatio-temporal data. We refer to the new underlying process as Lévy-dynamic spatio-temporal process. With the very limited computing facilities of Indian Statistical Institute, which provides us access to only 80 VMWare cores, analyzing terabytes of spatio-temporal data is still infeasible. But we are able to analyze a highly structured spatio-temporal sea surface temperature data consisting of 139, 300 observations (we chose 300 spatial locations, each with 398 time points as training points and set aside another 50 spatial locations, each with 398 time points for prediction) in less than 24 hours. It is important to remark in this context that although our Markov Chain Monte Carlo (MCMC) algorithm is highly intricate, it consists of parallelizable structures exploiting which we have parallelized the algorithm over the available cores with the C language and the Message Passing Interface (MPI) protocol, leading to significant computational savings. We point out that our our approach is flexible enough to model space-time data where the dynamic structure is less pronounced, or even purely spatial data.
In general, and definitely in the big data context, researchers do not concern themselves with the theoretical properties of their spatio-temporal models. For instance, most real datasets are expected to arise from nonstationary, non-Gaussian stochastic processes, but it is common practice for the sake of convenience to assume stationary Gaussian processes, usually with isotropic covariance structures. Assumption of separability of the covariance structure with respect to spatial and temporal structures is also very much common in the literature. The drawbacks of such simplistic approaches did motivate researchers to develop nonstationary, nonparametric and nonseparable approaches to modeling space-time data. However, a common limitation in all such approaches is the failure to account for the realistic property that the lagged spatio-temporal correlations converge to zero as the lags tend to infinity, despite nonstationarity, non-Gaussianity and nonseparability. Das and Bhattacharya (2020) provide an example of such lagged correlation property in the case of a PM10 pollution dataset, which they also established to be nonstationary and non-Gaussian. In this article, we demonstrate the same properties in the case of the sea surface temperature data. Thus, all the published approaches seem to be inadequate for modeling realistic spatial/spatiotemporal data. A comprehensive account of the strengths and limitations of the existing spatio-temporal approaches is detailed in Das and Bhattacharya (2020). In this endeavor, we show that our Lévy-dynamic space-time process possesses all the aforementioned realistic properties. Thus, our new model harnesses powerful parallel computing ability with desirable realistic properties for analyzing large datasets.

Other approaches with the desirable spatio-temporal properties
It is important to point out that our Lévy-dynamic approach is not the first one to consist of the desirable spatio-temporal properties. Indeed, the spatio-temporal process of Das and Bhattacharya (2020) is a very flexible process in this regard. The process, which results from an appropriate kernel convolution of order-based dependent Dirichlet processes (Griffin and Steel (2006)), is nonstationary, nonparametric, nonseparable, and possesses the property that the lagged spatio-temporal correlations converge to zero as the lags tend to infinity. The continuity and smoothness properties are also accounted for. However, the temporal part of the process does not have the dynamic structure, and considers time as an argument of the functional form of the stochastic process. Note that such a strategy is very appropriate for various datasets where the numbers of time points vary significantly with the spatial locations, with many locations having only a few time points. Pollution datasets on PM10 and PM2.5, for instance, are of this nature, and have been analyzed by Das and Bhattacharya (2020). In such cases, temporal dynamics are inappropriate. However, for other cases, incorporation of temporal dynamics is important. Computationally, the method is not too demanding, but analysis of very large data in reasonable time still seems to be infeasible. Importantly, several aspects of the MCMC algorithm can be parallelized, which might make analyses of many large spatio-temporal datasets feasible. Guha and Bhattacharya (2017) proposed a nonstationary, nonparametric, nonseparable dynamic state space spatio-temporal model, based on compositions of Gaussian processes in both the observational and the latent evolutionary levels. Under suitable conditions, the lagged spatio-temporal correlations also converge to zero. Continuity and smoothness properties of the process are investigated as well. But as it stands, the computational aspects seem to be too demanding to allow analysis of very large space-time datasets within reasonable time. However, we do have ideas to significantly improve the computational method, along with suitable parallelization.
To our knowledge, other than our Lévy-dynamic process, the approaches of Das and Bhattacharya (2020) and Guha and Bhattacharya (2017) are the only available ones that realistically account for nonstationarity, non-Gaussianity, nonseparability and convergence of the lagged correlations to zero. That all these properties are to be expected of real data, has been aptly demonstrated in Das and Bhattacharya (2020) with the pollution data, as already referred to. In this article, we shall demonstrate all these properties in detail, with respect to the sea surface temperature data that we analyze.

Existing approaches for large spatio-temporal data analyses
Although our intention is to provide an overview of the existing models and methods for large space-time data, most of the relevant existing literature seems to be exclusively concentrated on spatial data. Hence, we shall include mostly spatial methods in this brief review.
Banerjee (2017) reviews methods based on Gaussian processes for large spatio-temporal data, with focus on low-rank models and methods based on sparse covariance matrices associated with Gaussian processes. The essence of low-rank models is to represent the underlying (Gaussian) process in terms of realizations of some latent process with a relatively small number of points, so that dimension is effectively reduced. There are various approaches in this regard based on kernel convolutions and posterior expectations of the original process given the process values at a small set of points. Sparsity in covariance matrices is induced by specifying that the spatio-temporal distance between two points in space and time is zero beyond some specified threshold. Various issues related to the basic methods are discussed, with references to computational gains in large datasets of size of the order 10 5 . Guhaniyogi et al. (2011) analyse a forest biomass spatial dataset consisting of about 6000 observations, using Gaussian predictive process model, based on 25, 000 MCMC iterations. In another work, Guhaniyogi and Banerjee (2018) divide up the available spatial data into several sub datasets, fit Bayesian Gaussian process model to each sub dataset in parallel, and combine the results using geometric mean of the posteriors given the sub datasets. The procedure allowed them to analyse a spatial sea surface temperature dataset consisting of 120, 000 spatial observations (they used 117, 600 observations as training data points and set aside the rest for prediction). The approach does not have a temporal component and that is the reason that they were forced to consider only spatial analysis of the data, confining attention the same month (October) across the domain. It would have been useful if the MCMC details and computing time for this dataset were also reported in the paper. Heaton et al. (2018) presents a new flavour by not only reviewing Gaussian process based methods for large spatial data, but also reporting the details of a competition among various research groups on the basis of their preferred methods for analysing given simulated and real datasets. Both the simulated and real datasets consisted of 105, 569 spatial training observations, while the test data sets consisted of 44, 431 and 42, 740 observations, respectively.
The major concern in Gaussian process models is the large matrix-based computations which are necessarily inefficient, and much effort of the existing works has been directed towards simplification of such matrix computations, by various means. Clearly, far efficient computational algorithms can be achieved for approaches that are matrix-free. The Whittle likelihood approach (Whittle (1954), Guyon (1995)) associated with the spectral domain is also matrix-free and hence amenable to fast computation, but in reality the approach has limited application (see Banerjee (2017)). In the realm of classical spatio-temporal linear dynamic Gaussian state-space models where the spatial points are on a lattice grid, efficient computational strategies, that are essentially matrix-free, can be designed; see Dutta and Mondal (2015), for instance. However, for other setups available in the literature, and particularly in the Bayesian paradigm, appropriate matrix-free methods are difficult to devise. Furthermore, issues such as nonstationarity, non-Gaussianity, nonseparability and properties of lagged correlation structures do not seem to find importance in the existing works related to large data. The relatively small sizes of the datasets and much smaller MCMC sample sizes employed in the Bayesian spatial/spatio-temporal literature also leaves much to be desired.
The above issues provide the motivation for introduction of our Lévy-dynamic process. Indeed, our Lévy-dynamic approach is completely matrix-free, and hence, needless to mention, is a right candidate for analyzing large datasets. In addition, our model encapsulates all the realistic properties of spatial/spatio-temporal processes that are overlooked by the existing methods.
The rest of our article is structured as follows. In Section 2 we provide a brief overview of Lévy random fields, and in Section 3 introduce our Lévy-dynamic spatio-temporal process. The properties of its covariance structure, as well as continuity and smoothness properties, are investigated in Section 4. Specifications of relevant stochastic processes driving the Lévy-dynamic spatio-temporal process, are provided in Section 5. In Section 6, we introduce spatio-temporal random effects in our Bayesian model to account for finer details of the underlying real phenomenon, and in Section 7, we provide the hierarchical form of our complete Bayesian Lévy-dynamic model, along with the prior specifications. An overview of our parallel MCMC algorithm for implementing the Bayesian model is provided in Section 8. In Section 9, we provide details of our simulation experiment for assessing the performance of our proposed model and methodologies. Details of our analysis of the large sea surface temperature dataset are provided in Section 10. Finally, we summarize our contributions and make concluding remarks in Section 11.
Proofs of our results, the forms of the joint posterior and the full conditionals, and the complete parallel MCMC algorithm for Bayesian Lévy-dynamic inference are provided in the supplement, whose sections, equations and algorithms have the prefix "S-" when referred to in this paper.
Then, L(·) is a random signed measure such that for disjoint Borel sets A i , L(A i ) are independent, infinitely-divisible random variables. This random measure is referred to as the Lévy random measure, which is endowed with the following form of characteristic function (see, for example, Wolpert et al. (2011)) where ν(dµ, dβ) = λπ(dµ, dβ), is referred to as the Lévy measure. The Lévy measure is not required to be finite, provided that (2) is well-defined for all ζ ∈ R. For details regarding integrability in the case of infinite Lévy measure, see Wolpert et al. (2011) and Applebaum (2004). For our purpose, we shall consider only finite Lévy measure. When (2) is well-defined, it is possible to construct (1) using integrals with respect to Poisson random measures. That is, let N (dµ, dβ) ∼ P oisson(ν(dµ, dβ)) be the Poisson random measure, so that for disjoint Borel sets C i ⊆ R p+1 , N (C i ) ∼ P oisson (ν(C i )) independently. Then for any Borel set A with compact closure, given J = N (R p+1 ), where given J, (µ j , β j ) : j = 1, . . . , J is the random set of support points of the Poisson random measure. For examples of the Lévy measure and the corresponding Lévy random measures, see Wolpert et al. (2011) andApplebaum (2004).

Lévy random field
Consider any real-valued measurable function g on R p . Then, again using integration with respect to Poisson random measure as in (3), consider the following representation (see Wolpert et al. (2011)): The representation (4) constitutes the Lévy random field. This is well-defined for bounded measurable functions g when the Lévy measure is finite. Now, extending g(µ) to g(x, µ), where x ∈ R p , we obtain from (4): The extension (5), the similar form of which is provided in (Wolpert et al. (2011)) will play an important role in our spatio-temporal modeling strategy. Next, we introduce our proposed idea of nonstationary, nonseparable, dynamic spatiotemporal process that uses aspects of Lévy random fields of the form (5) as building blocks.
Since π, the stochastic process for (µ jt , β jt ) : t = 1, 2, . . . , is stationary, the marginal distributions of (µ jt , β jt ) are the same for all t and j. It follows that when M(s) = s, the marginal distribution f (s, t) reduces to the same form as that of Wolpert et al. (2011), when Σ j = Σ and τ j = τ in their case. Now, there might arise the question regarding independence of Σ, τ and ξ of j and t. Note that, since J t depends upon t, it is not possible to assume that Σ = Σ j , τ = τ j and ξ = ξ j , as this would imply that for different t, the dimensions of {(Σ j , τ j , ξ j ) : j = 1, . . . , J t } are different, which would not make sense. The assumptions Σ = Σ jt , τ = τ jt and ξ = ξ jt are sensible in this regard, but since µ jt and β jt are already time-dependent, time-dependence of Σ, τ and ξ might lead to temporal bias emerging from too many temporal dependence structures. Moreover, τ and ξ are anyway associated with the temporal part of the kernel, which is a function of the time index t.
There might also arise the question on not allowing J t to depend upon s. Again, by the same argument as above, this does not make sense unless µ jt and β jt are also made sdependent, which might lead to spatial bias in this case since the kernel is already spatially dependent. Moreover, making µ jt and β jt spatially dependent is expected to bring in much computational burden, while no inferential gain is expected.
Also note that since unlike space, time is dynamic in nature, postulation of dynamic structures for µ jt and β jt is indispensable for imparting a dynamic structure to our model, but for spatial dependence no more structure is necessary given spatial dependence of the kernel.
There is also an important issue regarding the temporal dynamics. In many spatiotemporal datasets, the numbers of temporal data points for the spatial locations are not only very different, but a large number of locations usually contain only a few temporal data points, even as small as just 2 or even 1. These are common issues in pollutant datasets such as PM10 or PM2.5 (see, for example, Das and Bhattacharya (2020)). To such data, application of our Lévy-dynamic process with its temporal dynamics, would not be sensible. However, in these situations we can set J t = J, µ jt = µ j and β jt = β j . We may then also set Σ = Σ j , τ = τ j and ξ = ξ j . That is, in the aforementioned situation, we propose modification of (7) to The temporal component of the spatio-temporal process will be taken care of in the t − τ j part of the kernel in (8). This modified model is of course very much applicable to the above kinds of spatio-temporal data. With the exponential kernel form K(s, t|Σ, ξ) = exp − 1 2 s T Σs − ξ|t| , the covariance structure will be separable in this case, but if desired, nonseparability can be easily enforced by slightly modifying the kernel toK(s, t|Σ) = exp − 1 2s T Σs , wheres = (s T , t) T , and where Σ is a positive definite matrix with nonzero off-diagonal elements.
Note that our ideas are applicable to purely spatial context as well, by simply replacing 4 Theoretical properties of the Lévy-dynamic process As iterated several times, in keeping with reality, it is important to ensure nonstationarity of the proposed spatio-temporal process. It is also desirable that the spatio-temporal correlations converge to zero as the lag tends to infinity, which is expected of the underlying real phenomenon. For the Lévy-dynamic process, we establish both the properties in Section 4.1. Nonseparability of the covariance structure of our process with respect to space and time is evident from the covariance form that we provide with respect to these results, and hence we do not provide any separate result on nonseparability. Furthermore, we investigate continuity and smoothness properties of our process in Sections 4.2 and 4.3, respectively.

Covariance properties
The purpose of this subsection is two-fold. First, we show that the covariance structure of the Lévy-dynamic process is nonstationary in the sense that it does not depend upon the locations and times only through their differences (Theorem 1). Thus, the process is not even weakly stationary, implying strong nonstationarity. Then we show that the lagged covariance structure, in spite of nonstationarity, converges to zero as the spatio-temporal lag tends to infinity (Theorem 2). Theorems 1 and 2 thus establish the realistically desirable properties of the Lévy-dynamic process.
Theorem 1. The covariance structure with respect to the function f (·, ·) given by (7) is nonstationary, satisfying the following properties: (i) Given any t, Cov(f (s 1 , t), f (s 2 , t)) does not depend upon s 1 and s 2 only through s 1 −s 2 .
(ii) Given any s ∈ R p , Cov(f (s, t 1 ), f (s, t 2 )) does not depend upon t 1 and t 2 only through t 1 − t 2 .

Continuity properties
If the process is almost surely continuous for every x 0 ∈ R p then the process is said to have continuous realizations.
First, it is clear that f (s, t) = 0≤j<Jt K(M(s) − µ jt , t − τ |Σ, ξ)β jt can not be almost surely continuous in (s, t), since f (s, t) is a jump process with respect to t (and J t ). However, for fixed t, f (s, t) can be almost surely continuous with respect to s, as the following result shows.
Theorem 5. Assume that M(s) is almost surely continuous in s and that K(x−µ, t−τ |Σ, ξ) is continuous in x. Then f (s, t) is almost surely continuous in s.
Theorem 6. Assume that M(s) is almost surely continuous in s and that K(x−µ, t−τ |Σ, ξ) is continuous in x. Also assume that K(·, ·|·, ·) is uniformly bounded. Then f (s, t) is L 1continuous in s.
The following two results show that f (s, t) is not L 2 -continuous even with respect to s.
Theorem 7. Assume that M(s) is almost surely continuous in s and that K(x, t) is continuous in (x, t). Also assume that K(·, ·) is uniformly bounded. Even then f (s, t) is not L 2 -continuous with respect to s for any fixed t.
The next theorem shows that if only convergence in expectation in considered, then f (s, t) converges to f (s 0 , t 0 ) in expectation as (s, t) → (s 0 , t 0 ).
Theorem 8. Assume that M(s) is almost surely continuous in s and that K(x, t) is continuous in (x, t). Also assume that K(·, ·) is uniformly bounded. Then, as (s, t) → (s 0 , t 0 ),

Smoothness properties
Now we examine differentiability of our spatio-temporal process.
If the process is almost surely differentiable at all x 0 ∈ R p , then it is said to be differentiable almost surely.
Since f (s, t) is not even continuous in (s, t), it is certainly not differentiable. However, for fixed t, differentiability of our process is given by the following result.
Theorem 10. Assume that for almost all paths, all the partial derivatives of the elements of M(s) with respect to the elements the s exist and are continuous. Also assume that all the partial derivatives of K(x − µ, t − τ |Σ, ξ) with respect to the elements of x exist and are continuous. Then f (s, t) is almost surely differentiable with respect to s.
Then f (s, t) is L r -differentiable with respect to s. 5 Choice of increasing stochastic processes for M and stationary processes for µ t and β t

Smooth increasing stochastic processes for the components of M
As valid increasing stochastic processes, the subordinators, which are almost surely increasing Lévy processes, merit serious consideration. Examples of such increasing processes are Poisson processes, α-stable subordinators, the Lévy subordinator, inverse Gaussian subordinators, Gamma subordinators, etc. See Applebaum (2004) for details on subordinators.
Since these are Lévy processes, they have stationary and independent increments. However, increasing stochastic processes with independent increments must be jump processes; see Ferguson and Klass (1972). Hence, although subordinators qualify as models for the components of M, they fail to satisfy smoothness, or even continuity of f , as required by Theorems 5, 6, 8, 10, 12. This requires us to create new increasing processes that are also smooth. Details follow. For each = 1, . . . , p, let us consider the following stochastic process for M : for s 1 , s 2 ∈ R such that s 1 > s 2 , where C is some positive constant and X is a positive random variable independent of s 1 and s 2 . We set r ≥ 1 in (11). Almost sure continuity and differentiability of M is achieved even with r = 1. Note that under (11), although M have stationary increments, the increments are not independent, due to the presence of X . For = 1, . . . , p let s We shall also set X = |Z |, where Z ∼ N (ν , ω 2 ), where ν , ω 2 will be treated as unknown.
The positive constant C will also be treated as unknown in our setup. We shall set M s r , whereC > 0 will be treated as unknown. We shall set r = 2 in our applications.

Stationary stochastic process models for µ t and β t
In Section 3 we have assumed that µ t and β t are stationary stochastic processes. This assumption was important in proving our theoretical results results. In practice, particularly, for large spatio-temporal data analysis, it is important to keep the forms of the stationary stochastic processes as simple as possible. Thus, in practice, it is useful to model the components of µ t and β t as independent stationary AR(1) processes.
In spite of such simplicity, the actual time series f (·, t) is nonstationary and has rich enough temporal covariance structure, borne out by the kernel, which is rendered timedependent through µ jt as well as t, and even through β jt . The relevant results on temporal nonstationarity are provided by (ii) and (iii) of Theorem 1, and Theorem 2 shows that the covariance structure has the desirable asymptotic property.
Of course, if necessary, we can consider any desired stationary stochastic process models for µ t and β t , with dependence among the components of µ t and β t . Theorems 1 and 2 would continue to hold in all such situations, and for any almost surely monotonically increasing stochastic processes M ; = 1, . . . , p, bringing out the flexibility and generality of our strategies.
In particular, if t k − t k−1 = 1 for all k, then the stationary, first order autoregressive model may be the default choice. For irregularly spaced time series, we recommend the irregular autoregressive (IAR) model introduced by Eyheramendy et al. (2018), which we briefly review below.

Irregular autoregressive model
For an increasing sequence of observation times {t k : k ≥ 1}, an IAR process {ζ t k : k ≥ 1} is defined by Eyheramendy et al. (2018) as the following: where t k are iid random variables with zero mean and unit variance. It is assumed that ζ t 1 is a zero-mean random variable with variance σ 2 ζ . It can be seen that E(ζ t k ) = 0 and V ar(ζ t k ) = σ 2 ζ , for k ≥ 1. For k 1 ≥ k 2 , the covariance between ζ t k 1 and ζ t k 2 is given by implying that for any t > s, the autocovariance function can be defined as γ(t − s) = E(ζ t ζ s ) = σ 2 ζ ρ t−s ζ , signifying a second-order weakly stationary process. However, under some mild conditions, strict stationarity and ergodicity can be ensured, as shown in the following result of Eyheramendy et al. (2018): Theorem 13 (Eyheramendy et al. (2018)). Consider the IAR process defined by (13), and let 0 < ρ ζ < 1. Assume that t k − t k−n ≥ C log n, as n → ∞, where C is a positive constant satisfying C log ρ 2 ζ < −1. Then, there exists a solution to the IAR process and the sequence {ζ t k : k ≥ 1} is stationary and ergodic.
It is noted in Eyheramendy et al. (2018) that the case t k − t k−1 = 1 for all k ≥ 1 corresponds to regular AR(1), which satisfies the conditions of Theorem 13, since t k − t k−n = n > log n and ρ 2 ζ < 1 is a part of the assumptions regarding stationary AR(1) processes. Observe that regular AR(1) allows −1 < ρ ζ < 1 for stationarity, while IAR requires 0 < ρ ζ < 1. Indeed, from (14) it is clear that non-integer positive real values t k 1 − t k 2 , ρ σ 2 ζ will be denoted by ρ and σ 2 , respectively. In the case of β t , we shall denote ρ ζ and σ 2 ζ by ρ β and σ 2 β , respectively. Although derivation of the important statistical properties of the IAR (or AR(1)) process does not require the normality assumption of the errors, for our applications, we shall assume normality.
6 Incorporation of random effects in the Lévy-dynamic spatio-temporal model In practice, the functional form f (s i , t k ) driven by specific choices of the kernel K need not be always sufficient to explain the underlying spatio-temporal structure in precise details. Hence, we shall attempt to further enhance inference by considering spatio-temporal random effects in our model. In other words, we shall consider the following model for data analysis: where α is the overall effect and φ(s i , t k ) are the spatio-temporal random effects. We assume independently for i = 1, . . . , n and k = 1, . . . , m. In the above, In the case of prediction of y(s,t) at locations and time pointt, where at least one ofs ort is not in the training dataset, then, assuming that for N ≥ 1, {(i * r , k * r ) : r = 1, . . . , N } = arg min{ s − s i 2 + (t − t k ) 2 : i = 1, . . . , n; j = 1, . . . , m}, we define φ 0 (s,t) = N r=1 y(s i * r , t k * r )/N . Thus, although the spatio-temporal dependence structure is encapsulated in the dependence among f (s i , t k ), the random effects φ(s i , t k ), along with the overall effect α, are introduced to capture the finer details of the specific spatial location and time point associated with the data and enhance inference. In particular, precisions of the predictions at given locations and time points where data are not observed, are likely to be sharper with these random effects. However, given the basic spatio-temporal structure offered by f (s i , t k ) and the way φ 0 (s i , t k ) are constructed, α and φ(s i , t k ) are not expected to have significant variabilities. As such, we shall consider the priors for σ 2 α and σ 2 φ to reflect the opinion that with relatively high certainty they are not much different from zero. Note that α may be viewed as the average of all the y(s i , t k ) and can be set to zero when y(s i , t k ) are standardized to have mean zero and variance one.
Observe that although there are nm random effects in our model, these can be integrated out from (15), so that under the marginalized model, y(s i , t k ) admits the representation where˜ ik ∼ N (0, σ 2 + σ 2 φ ), independently. As argued above, there are reasons to consider a prior for σ 2 φ that concentrates around zero. Thus, it would make sense to deterministically set σ 2 φ ≈ 0, which would also make σ 2 identifiable in the variance σ 2 + σ 2 φ of˜ ik . This entire exercise certainly leads to huge computational savings compared to the original, nonmarginalized version. Importantly, with σ 2 φ = 0, we shall demonstrate with our simulation experiment that although the non-marginalized version has better MCMC mixing properties, the final Bayesian prediction results are remarkably similar for the two versions. Thus, we shall also consider the marginalized version with σ 2 φ = 0 in the real data scenario. It is worth mentioning that ours is not the first spatio-temporal work to incorporate random effects. Random effects in spatial and spatio-temporal setups have also been considered in Kang and Cressie (2011) and Wu et al. (2016); see also Section 4.4.1 of Wikle et al. (2019).

Hierarchical form of the Lévy-dynamic model with prior details
For the sake of generality, we assume the time points to be of the form {t k : k ≥ 1}. Then our Bayesian Lévy-dynamic spatio-temporal model admits the following hierarchical form: for i = 1, . . . , n and k = 1, . . . , m, X ∼ N ν , ω 2 ; = 1, . . . , p; We shall further set Σ to be a diagonal matrix with unknown positive diagonal elements σ 2 1 , . . . ,σ 2 p . The specific prior forms for (29) and (30) would be the following in our applica-15 tions:σ 2 ∼ IG(aσ2 , bσ2 ); = 1, . . . , p; In the above, G(a, b) denotes the Gamma distribution with positive parameters a and b having density g 1 (x) ∝ x a−1 exp(−bx), for x > 0. IG(a, b) denotes the inverse Gamma distribution with positive parameters a and b with density g 2 (x) ∝ x −a−1 exp(−b/x), for x > 0. We recommend a = 2.01 and b = 1.01 for the inverse gamma priors except those for σ 2 α , σ 2 φ and σ 2 . Thus the means and variances are b/(a−1) = 1 and b 2 /((a−1) 2 (a−2)) = 100 in these cases. For σ 2 α and σ 2 φ we recommend a = 10 4 and b = 1, so that the means and variances are close to zero, to reflect the opinion that the spatial and temporal random effects do not have significant variabilities in the presence of f (s i , t k ). Moreover, in our applications, we fit our Bayesian model after standardizing the space-time datasets, and hence set α = 0 for model implementation. We finally convert our Bayesian predictions to the original locations and scales, for the reporting purpose.
For σ 2 , we again recommend a = 10 4 and b = 1. Again, this encapsulates our opinion that σ 2 is not much different from zero. The reason for this opinion about σ 2 is that when n and m are even reasonably large, with very high probability, the overall variability of the spatio-temporal data is expected to be very large, which would drastically increase the posterior mean and variance of σ 2 , unless its prior means and variances are fixed to be very small. Note that large posterior mean and variance of σ 2 would render predictions at desired spatial locations and time points highly unreliable. Consequently, we set such small values of the mean and variance to obtain reasonable predictions. For λ, we set b λ = 0.001 and a λ = 10b λ in our applications, so that the prior mean and variance of λ are 10 and 10 4 , respectively.
As regards the zero-mean normal priors for ν , log ρ β 1−ρ β and log ρ 1−ρ , we set the variances to be 100 in our applications.
Recall that when t k − t k−1 = 1, the IAR model boils down to the regular AR(1) model, so that in such a situation we set |ρ | < 1 for = 1, . . . , p and |ρ β | < 1 for stationarity. For the priors on ρ and ρ β , we then set ρ = −1 + 2 exp(ρ ) 1+exp(ρ ) and To ensure stationarity of the regular AR(1) setup we also set β jt 1 ∼ N 0, Now note that the exponential kernel K(M(s) − µ, t − τ |σ 2 1 , . . . ,σ 2 p , ξ) that we use for our purpose will have negligible values for large values of τ , ξ,σ 2 1 , . . . ,σ 2 p . We reparameterize these non-negative parameters generically by exp(ϕ), where −∞ < ϕ < ∞, and obtain the prior distribution for ϕ corresponding to the priors for the original non-negative parameters. We then truncate ϕ on the interval [−20, 5]. We adopt the same strategy for the nonnegative parameters σ 2 , C andC as well. Thus, although we allow small positive values of the original non-negative parameters, large positive values are ruled out. We also truncate X and µ ( ) jt k on [−10, 10]. Further, the mixing behaviour of our MCMC is improved by truncating the normal distributions associated with ρ and ρ β on [−10, 10] and adopting the same aforementioned reparameterization and truncation strategy for σ 2 β .

An overview of our parallel MCMC algorithm
The form of the joint posterior distribution is provided in Section S-20, using which the forms of the full conditional distributions of the parameters are detailed in Section S-21. In our MCMC method, we use Gibbs sampling steps to simulate from most of the standard full conditional distributions. We refer to the relevant set of parameters updated using Gibbs steps by ζ. Although the full conditionals of σ 2 1 , . . . , σ 2 p and σ 2 β are also available in closed forms, instead of using Gibbs steps for these, we update these parameters simultaneously in a single block consisting of θ = (X 1 , . . . , X p ,C 1 , . . . ,C p , C 1 , . . . , C p ,σ 2 1 , . . . ,σ 2 p , τ, ξ, ρ 1 , . . . , ρ p , σ 2 1 , . . . , σ 2 p , ρ β , σ 2 β ) using Transformation based Markov Chain Monte Carlo (TMCMC) introduced by Dutta and Bhattacharya (2014). The essence of TMCMC is to simultaneously update many parameters in a single block using simple deterministic transformations of some low-dimensional (usually, one-dimensional) random variable. As can be anticipated, this drastic dimension reduction leads to great improvement in acceptance rates and faster convergence compared to traditional MCMC methods. For details on TMCMC, see Dutta and Bhattacharya (2014), Dey and Bhattacharya (2016), Dey and Bhattacharya (2017), Dey and Bhattacharya (2019). In our case, we consider a mixture of additive and multiplicative TMCMC; such a mixture outperforms both additive and multiplicative TMCMC by combining the localised moves of additive TMCMC and the non-local moves of multiplicative TMCMC (see Dey and Bhattacharya (2016) for details). The multiplicative transformation has been referred to as "random dive" by Dutta (2012).
After implementing the mixture TMCMC, we supplement this with another deterministic move type consisting of additive and multiplicative transformations to further enhance mixing properties of our methodology. Now note that due to the Markov property, (U k , β k , J t k ), for all the odd values of k ∈ {1, . . . , m} can be updated simultaneously in parallel processors. Once (U k , β k , J t k ) are updated for odd values of k, those for the even values of k can then be updated simultaneously in parallel processors.
Since J t k is a random variable, this makes the dimensions of U k and β k random, rendering the updating problem of (U k , β k , J t k ) a variable-dimensional problem, for every k = 1, . . . , m.
Since reversible jump MCMC introduced in Green (1995) is well-known to be a very inefficient method for handling variable-dimensional problems, Das and Bhattacharya (2019) came up with a novel and efficient alternative to solving variable-dimensional cases, using appropriate deterministic transformations of fixed and low-dimensional random variables. The method, referred to as Transdimensional Transformation based Markov Chain Monte Carlo (TTMCMC), is an extension of TMCMC for fixed-dimensional setups to general variabledimensional problems.
In our case, for each value of k, we update (U k , β k , J t k ) using TTMCMC in separate parallel processors. Also, given all other unknowns, we update φ(s i , t k ) simultaneously in separate parallel processors by sampling from their full conditional distributions. Note that in order to sample from the full conditional distributions of λ, σ 2 φ , σ 2 and α, computations of the sums ) 2 are required. We compute these by splitting the sums into the available parallel processors, each processor computing only a small part of each sum. The final sum is aggregated into a single processor, which then updates the relevant parameters. The TMCMC update required in step (S-78) requires computing n i=1 m k=1 log y(s i , t k )|X 1 , . . . , X p , C 1 , . . . , C p ,C 1 , . . . ,C p , U k , β k , J t k ,σ 2 1 , . . . ,σ 2 p , τ, ξ, α, φ(s i , t k ), σ 2 , which we again compute by splitting the sum into the available parallel processors, finally aggregating the result into a single processor where TMCMC is applied. The complete algorithm is provided as Algorithm S-1 in the supplement. Note that in the marginalized model where the random effects are integrated out, the algorithm is simplified.

Simulation study
For our simulation experiment, we generate data from the so-called general quadratic nonlinear (GQN) model (Wikle and Hooten (2010), Cressie and Wikle (2011)), to which we fit our Lévy-dynamic spatio-temporal model with random effects and make predictions at various locations and time points. Specifically, our data-generating GQN model is of the following form: for i = 1, . . . , n and k = 1, . . . , m, We assume that independently, φ 1t k (·), φ 2t k (·), η t k (·), t k (·), β 0 (·) ∼ GP (0, c(·, ·)), a zero-mean Gaussian process with covariance function c(s 1 , s 2 ) = exp (− s 1 − s 2 ), for any s 1 , s 2 ∈ R 2 , where · denotes the Euclidean norm. We further assume that independently, for i = 1, . . . , n, j = 1, . . . , n, l = 1, . . . , n, a ij ∼ N (0, 0.001 2 ) and b ijl ∼ N (0, 0.001 2 ). We set g(β t k −1 (s l )) = β 2 t k −1 (s l ) and for i = 1, . . . , n, independently simulate s i ∼ U (0, 1) × U (0, 1). We generate y(s i , t k ) for n = 120 spatial locations and m = 50 time points from the above model defined by (46) and (47) and the associated distributions. For the purpose of prediction, we set aside 20 spatial locations and the associated 50 times points for each of the locations. Our main goal in this simulation experiment is to make reliable predictions at the set-aside locations and time points. As mentioned earlier, we first standardize the dataset, to which we apply our model and methods, and make Bayesian predictions. Finally, we transform the predictions to the original location and scale for reporting. Thus, for our purpose, we set α = 0 and consequently, updating α and σ 2 α are not required.

Implementation details
We implement MCMC Algorithm S-1 (wth p = 2) for the non-marginalized Bayesian Lévydynamic model, written in C in conjunction with the message passing interface (MPI) protocol, on 25 parallel processors in our 2 TB memory VMWare (each core has about 2.8 GHz CPU speed). Generation of 11 × 10 4 MCMC realizations took 1 hour 27 minutes in our VMWare. However, on 80 parallel processors (the maximum number of cores on our VMWare), the time taken is 1 hour 47 minutes. This is a consequence of slower communications among much larger number of processors with a relatively small amount of data. That is, the computational overhead with 80 processors is less than the communication overhead among the 80 processors. Also note that the number of time points is only 50 and only 25 processors can be used at a time to update the associated parameters corresponding to even or odd indices, making the other processors redundant for these parameters. More experimentations led to the conclusion that 25 processors provide the most efficiency for this problem.
To prevent extreme propositions in the multiplicative moves of our algorithm brought about by dividing the current realization by too close to zero, we set U (−1, 1) subject to | | > 0.01. This is of course a theoretically valid step and irreducibility of our algorithm is preserved by the additive transformation. For details, see Dey and Bhattacharya (2016). We discarded the first 10 4 realizations as burn-in and stored every 10-th realization in the next 10 5 iterations to obtain 10 4 MCMC realizations for Bayesian inference. Figure 1 exhibit the trace plots of J k for different k (for which we display and provide the analyses of the corresponding spatial predictions) and Figure 2 provides the trace plots of some other parameters of our Lévy-dynamic Bayesian model. All the trace plots vindicate excellent convergence. Panel (i) of Figure 2 shows that even though the prior for σ 2 has mean and variance close to zero, the posterior distribution still takes on high values occasionally.

MCMC convergence
That the acceptance rates for the TTMCMC and TMCMC rates are adequate, are evident from the trace plots. Specifically, the average TTMCMC acceptance rates of the birth, death and no-change moves over 50 time points are approximately 0.09, 0.726 and 0.619, respectively. The average overall TTMCMC acceptance rate is 0.415. The fixeddimensional parameters that are updated using TMCMC have acceptance rate 0.921, and the corresponding acceptance rate associated with the mixing-enhancement step is 0.656. All these acceptance rates are calculated with respect to the entire set of 11 × 10 4 MCMC realizations of our algorithm. Figure 3 depicts the densities of the temporal predictions at various spatial locations as color plots. In the color plots progressively intense colors correspond to higher densities associated with 16 quantiles dividing the density support (that is, the minimum and the maximum quantiles undertaken for the color plots are 1/16 = 0.0625 and 15/16 = 0.9375). The true time series at the spatial locations are denoted by the thick, black line. Notice that in almost all the cases, the entire true time series is included in the high density regions of our Bayesian-predicted time series using our Lévy-dynamic spatio-temporal process. However, in a few cases, where the actual temporal data points take on very highly positive or negative values, our predictions have not been adequate (not shown). Figure 4 shows the spatial surface plots at various time points. The middle surface is the actual spatial surface (supplemented by spline-based interpolations), while the lower and upper surfaces (again, supplemented by spline-based interpolations) are the lower and upper bounds, respectively, of the 0.875 credible region. In other words, the lower and upper surfaces correspond to 1/16-th and 15/16-th quantiles of the respective posterior predictive distributions. The detailed spatial posterior predictive densities are shown in Figure 5, as color plots akin to the temporal color plots of Figure 3, with spatial indices replacing the time indices. Unlike the time points, no ordering is intended with respect to the spatial indices. That is, we simply refer to i = 1, . . . , 20, as the spatial indices associated with s i . Thus, it is evident from the figures that almost all the actual spatial data fall within the high-density regions of the corresponding posterior predictive distributions. However, again for some time points, there are a few spatial data points that are highly positive or negative, and our predictions failed in such cases (not shown).     Figure 4 are shown as colour plots with progressively higher densities depicted by progressively intense colours.

Faster implementation after integrating out the random effects
Note that simulation of φ, although parallelised, requires many parallel processors for efficiency if either n or m is even moderately large. As already reported in Section 9.1, this constitutes significant communication overhead among the processors and slows down implementation. Redundancy of the cores is also a consequence in the case of moderately large datasets. On the other hand, for significantly large datasets with a large number of time points, a large number of parallel processors are necessary for efficiency. However, such large number of processors are usually not available.
To get around this problem, for i = 1, . . . , n and k = 1, . . . , m, we invoke the marginalized model (17) and the discussion thereafter and set φ(s i , t k ) = φ 0 (s i , t k ), provided the posterior uncertainties in φ(s i , t k ) are negligible. Since the prior distributions of φ(s i , t k ) are concentrated around φ 0 (s i , t k ) (recall that σ 2 φ is concentrated around zero a priori), our experiments reveal that this is indeed the case (not shown for brevity).
To justify our standpoint, we conduct a further experiment by setting φ(s i , t k ) = φ 0 (s i , t k ) in our model. Not only does the time taken reduce to just 24 minutes from 1 hour 27 minutes with 25 cores, but the prediction results with this setup (see Figures 6, 7 and 8 are remarkably similar to those in the previous non-marginalized implementation (Figures 3, 4 and 5).
In this case, the average TTMCMC acceptance rates of the birth, death and no-change moves over 50 time points are approximately 0.103, 0.673 and 0.617, respectively, and the average overall TTMCMC acceptance rate is 0.420. The TMCMC step has acceptance rate 0.017, and that of the mixing-enhancement step is 0.152. Thus, compared to the nonmarginalized model, here the acceptance rates of the TMCMC and the mixing-enhancement steps have significantly decreased, but this is not a strong enough reason for concern since both the models yield remarkably similar performances with respect to our main goal, Bayesian prediction.

Analysis of sea surface temperature data
We now consider analysis of a real, sea surface temperature dataset available at http: //iridl.ldeo.columbia.edu/SOURCES/.CAC/ in the netCDF file format. The data pertains to monthly sea surface temperatures during January 1970 -December 2003 at the tropical Pacific Ocean region covering 124 • E − 70 • W and 30 • S − 30 • N , gridded at a 2 • by 2 • resolution. Analysis of a somewhat similar (and much smaller) dataset has been reported in Cressie and Wikle (2011), on the basis of some simple linear and non-linear dynamic state-space models, but that data represented monthly temperature anomalies from the normal, rather than the actual temperatures. The models of Cressie and Wikle (2011) are not intended to cover enough grounds like ours, namely, weak and strong nonstationarity, non-separability, nonparametric non-Gaussianity and convergence of the lagged correlations to zero. Nevertheless, the simplicity of their models enabled them to perform simple Gibbs sampling based Bayesian analysis, for both of their linear and nonlinear dynamic models. However, the samplers are run for only 6000 MCMC iterations (the first 1000 discarded as burn-in). Wikle et al. (2019) also consider the anomalies dataset for some simplistic    spatio-temporal analyses. Our dataset consists of space-time data at n = 2520 spatial locations, for each of m = 398 time points. That is, the size of our dataset is nm = 10, 02, 960. Figure 9 displays the sea surface temperature plots during January 1989January , 1993January and 1998 to exhibit the effects of La Niña (colder than normal temperatures), normal temperatures and El Niño (warmer than normal temperatures), respectively. Wikle et al. (2019) provide similar plots in essence with their anomalies data using different colouring schemes in the R package, as opposed to ours in Python, in the context of the actual temperatures.

Spatio-temporal nonstationarity of the sea surface temperature data
In real life situations, the assumption of stationarity, or even covariance stationarity of the underlying spatio-temporal process, are usually too naive to be realistic. In our case, Figure  9 alone vindicates that for different time points, the spatial distributions of temperature are different, thus empirically ruling out stationarity. Moreover, for different spatial locations, the distributions of the time series are also different, as can be observed from the thick black lines in Figure 13. For more formal conclusion regarding nonstationarity of the sea surface temperature data, we adopt the recursive Bayesian theory and methods developed by Roy and Bhattacharya (2020) in this regard. In a nutshell, their key idea is to consider the Kolmogorov-Smirnov distance between distributions of data associated with local and global space-times. Associated with the j-th local space-time region is an unknown probability p j of the event that the underlying process is stationarity when the observed data corresponds to the j-th local region and the Kolmogorov-Smirnov distance falls below c j , where c j is any non-negative sequence tending to zero as j tends to infinity. With suitable priors for p j , Roy and Bhattacharya (2020) constructed recursive posterior distributions for p j and proved that the underlying process is stationary if and only if for sufficiently large number of observations in the j-th region, the posterior of p j converges to one as j → ∞. Nonstationarity is the case if and only if the posterior of p j converges to zero as j → ∞.
Covariance stationarity has been treated by Roy and Bhattacharya (2020) using similar principles, replacing the local and global distributions by local and global covariances of spatio-temporal lag h , where h is the difference between the spatial-temporal co-ordinates, and · is the Euclidean distance. The process is covariance stationarity if and only if for sufficiently large number of observations in the j-th region, the posterior of p j converges to one as j → ∞, for all h > 0. On the other hand, the process is covariance nonstationary if and only if there exists h > 0 such that the posterior of p j converges to zero as j → ∞.
In our implementation of the ideas of Roy and Bhattacharya (2020), we set the j-th local region to be the entire time series for the spatial location s j , for j = 1, . . . , 2520. Thus, the size of each local region is 398, which is sufficiently large for our purpose. The number of regions, 2520, is also large enough for the Bayesian recursive theories to be applicable.
To check stationarity, we choose c j to be of the same nonparametric, dynamic and adaptive form as detailed in Roy and Bhattacharya (2020). The dynamic form requires an initial value for the sequence. It is important to remark here that in practice, the choice of the  initial value has significant effect on the convergence of the posteriors of p j , and hence such a choice must be carefully made. However, in our case, for such a large dataset we expect initial values even close enough to zero to suffice for inferring stationarity if the underlying phenomenon is indeed stationary. As it turned out, for all initial values less than or equal to 0.26, the recursive Bayesian procedure led to the conclusion of nonstationarity of the underlying spatio-temporal process. We implemented the idea with our parallelised C code on 2 parallel processors of our ordinary laptop; the time taken is just 2 seconds. For the initial value 0.26, Figure 10 displays the means of the posteriors of p j ; j = 1, . . . , 2520, showing clear convergence to zero. The respective posterior variances are negligibly small and hence not shown. Thus, as already anticipated, the spatio-temporal process that generated the sea surface temperature data, can be safely regarded as nonstationary. Figure 11 shows the results of our investigation of covariance stationarity. Panels (b), (c), (e) and (f) show convergence of the posterior means of p j in this context to zero for different partitioned intervals of h associated with sufficient data such that the covariances are well-defined. The posterior variances are again negligibly small as before. Thus, covariance nonstationarity of the underlying spatio-temporal process is also clearly indicated. In these cases, we chose the initial values of c j to be 0.05. The time taken for our parallelised C code implementation on our dual-core laptop is only 5 seconds for each h . Hence, the sea surface temperature phenomenon is not even weakly stationary.

Convergence of lagged spatio-temporal correlations to zero
Recall that two major purposes of our Lévy-dynamic spatio-temporal model is to account for nonstationarity of most real-life data and to emulate the property of most real datasets that the lagged spatio-temporal correlations tend to zero as the spatio-temporal lag h tends to infinity, in spite of nonstationarity. For the sea surface temperature dataset, we have already confirmed strict nonstationarity as well as covariance nonstationarity. It now remains to address the issue of the lagged correlations.  For our purpose, we randomly select 25 spatial locations and consider the entire time series associated with each of them. We further augment the data with the first 50 time points of another random location, thus yielding a dataset of size 10, 000. We compute the lagged correlations on 80 parallel processors on our VMWare, each processor computing the correlation for a partitioned interval of lag h such that the interval is associated with sufficient data making the correlation well-defined. The time taken for this exercise is about 3 minutes. Figure 12 shows that with respect to this real dataset, our expectation that the lagged spatio-temporal correlations converge to zero in spite of nonstationarity, is not unfounded. Further experiments with larger data sizes, but with much longer implementation times, corroborated this result.

Non-Gaussianity of the sea surface temperature data
Simple quantile-quantile plots (not shown for brevity) revealed that the distributions of the time series data at the spatial locations, distributions of the spatial data at the time points, and the overall distribution of the entire dataset, are far from normal. Thus, traditional Gaussian process based models of the underlying spatio-temporal process are ruled out. Since the temporal distributions at the spatial locations and the spatial distributions at different time points are also much different, it does not appear feasible to consider parametric stochastic process models for the data. In this regard as well, relevance of our nonparametric Lévy-dynamic process is quite pronounced.
In other words, we have validated that the underlying spatio-temporal process that generated the sea surface temperature data is non-Gaussian, strictly and weakly nonstationary, and the lagged correlations converge to zero as the lags tend to infinity. Moreover, there is no reason to traditionally assume separability of the spatio-temporal covariance structure. Since our nonparametric Lévy-dynamic process is endowed with all the aforementioned characteristics, it seems to be a very appropriate candidate for analysing the data.

Bayesian Lévy-dynamic model implementation and results
The current computational resources at Indian Statistical Institute are certainly not adequate for analysing the entire sea surface temperature dataset within a reasonable time frame. Hence, we randomly chose 300 spatial locations and the entire time series associated with each of them. Thus, our selected subsample consists of 119, 400 spatio-temporal observations, which is not a small dataset with respect to our sophisticated Bayesian hierarchical modeling framework with complex dependence structures. Indeed, we are not aware of application of realistically sophisticated Bayesian hierarchical models to spatio-temporal datasets as large. We further randomly choose another set of 50 locations from the remaining set of locations and the corresponding time series data of size 398 for each location for evaluation of the predictive performance of our model. As before, we standardize the dataset and report the prediction results after transforming them back to the original location and scale.
Simplification of the structure induced by the spatio-temporal random effects by setting φ(s i , t k ) = φ 0 (s i , t k ) for i = 1, . . . , n = 300 and k = 1, . . . , m = 398 (for training) and φ(s i , t k ) = φ 0 (s i , t k ) for i = 1, . . . , n = 50 and k = 1, . . . , m = 398 (for prediction) in our marginalized Lévy-dyanmic model (17) and the following discussion brought down the implementation time (on 80 parallel processors) from more than 3 days (estimated) to less than a single day. Here the average TTMCMC acceptance rates of the birth, death and nochange moves over 398 time points are approximately 0.095, 0.709 and 0.661, respectively, and the average overall TTMCMC acceptance rate is 0.439. The acceptance rates for TMCMC and the mixing enhancement steps are 0.033 and 0.03, respectively. Note that these acceptance rates are broadly similar to those reported in the context of the simulation experiment with the marginalized random effects model (Section 9.4). Thus, the acceptance rates seem to exhibit a tendency of robustness with respect to different datasets of varying sizes. Figure 13 shows the time series predictions at a few of the 50 spatial locations set aside for prediction. Observe that for each spatial location, the entire true time series falls well within the associated 0.875 posterior predictive density region. Although in panels (c) and in particularly panel (f), the time series do not pass through the highest posterior predictive density regions depicted by the most intense colours, the trends in these cases seem to be still well-captured by our Bayesian model and methods. Figure 14 displays the spatial predictions at a few time points, in the forms of 95% lower and upper Bayesian spatial prediction surfaces, while the middle, true spatial surface corresponds to the 50 spatial locations meant for prediction. The reason for choosing 95% surfaces, rather than 0.875, as in the simulation experiments is that for many time points, the 0.875% surfaces failed to capture several true spatial data points. Even the 95% surfaces failed to satisfactorily contain either one or two spatial data points, for several time points (not shown for brevity). Setting σ 2 φ to some non-negligible positive value might have solved the issue, but would have increased the Bayesian prediction intervals for most of the other data points which are already well-captured. Figure 15 exhibits the density based colour plots of the spatial predictions corresponding to Figure 14. In this case we use 40 percentiles to make the figures correspond to 95% credible regions, for comparability with Figure 14.

Summary and conclusion
The Gaussian process is overused in the spatial/spatio-temporal literature, particularly in large data scenarios. The reality issues such as non-Gaussianity, nonstationarity, nonseparability and properties of the lagged correlations are often relegated to the background in favour of convenience, even in general data-fitting scenarios. These seem to induce some inflexibility in the current state-of-the-art spatial/spatio-temporal statistics, as even a plethora of existing Gaussian process based methods and a competition among them failed to yield analysis of any big data in the order of terabytes. The key impediment in the implementation of all such methods is matrix-based computation, which may be ameliorated, but can not be avoided.
In the zest for simplifying computations in Gaussian processes, the realistic issues are often forgotten, as mentioned above. Thus, there is a need to develop realistic spatial and spatio-temporal models and methods that satisfy the realistic properties, and are also amenable to fast and efficient matrix-free computation to meet the challenges of large data. In this regard, we introduce our Bayesian Lévy-dynamic spatio-temporal model based upon Lévy random fields, and show that it satisfies the desirable realistic properties. As we have shown, the approach is flexible enough for modeling space-time data with weak temporal dynamics or even purely spatial data.
For capturing micro-scale spatio-temporal variations, we introduce spatio-temporal random effects, which are amenable to marginalization that enormously simplify computations. The model is completely matrix-free, but is variable-dimensional with respect to each of the time indices. We handle the variable-dimensional parameters using TTMCMC and the fixeddimensional parameters using TMCMC, all embedded in a novel parallel MCMC algorithm, which we code in C in the MPI paradigm for parallelism. The model structure allows us to update the variable-dimensional parameters for all the even (odd) time indices in parallel, followed by updating those for all odd (even) time indices. Even for fixed-dimensional updates, we compute the acceptance ratios and several other quantities in parallel. Thus, in conjunction with integrating out the random effects, our parallel MCMC algorithm leads to huge computational savings. However, the mixing properties are enhanced when the random effects are not integrated out. Despite this, the Bayesian predictions are almost unaffected by the issue of marginalization of the random effects as borne out by our simulation experiment, providing the green signal to consider the marginalized model for analysis of large datasets.
We indeed analyse a relatively large sea surface temperature dataset consisting of 139, 300 space-time observations using our marginalized Bayesian Lévy-dynamic model; the time taken being less than 24 hours on our VMWare with 80 cores. The encouraging results suggest that with more powerful and well-maintained computing facilities, we may ambitiously begin analyzing "big data" in the order of terabytes, without any compromise whatsoever on the theoretical properties with respect to data realism.

Supplementary Material
S-12 Proof of Theorem 1 Proof. Let us first prove (i). Note that given M(s 1 ) and M(s 2 ), the covariance between f (s 1 , t) and f (s 2 , t) is given by

(S-50)
Let us now consider the term (S-49). Note that for any s ∈ R p , (S-51) Hence, (S-49) is given by for some function g(·). Now, taking expectation of both sides of (S-51) yields for any s ∈ R p . Hence, At least the second term of (S-55) is not a function of s 1 − s 2 . Hence, it follows from (S-53), (S-54) and (S-55) that Cov (f (s 1 , t), f (s 2 , t)) = E Cov f (s 1 , t), f (s 2 , t) M(s 1 ), M(s 2 ) does not depend upon s 1 and s 2 only through s 1 − s 2 . This proves (i).
To prove (ii), note that for any s ∈ R p , for any t. Hence, the covariance term (S-57) is In (S-56), the covariance is given by Hence, 1 T Jt 1 Cov(X, Y M(s), J t 1 , J t 2 )1 Jt 2 simplifies to Consequently, (S-56) is given by Because of (S-58), Cov f (s, t 1 ), f (s, t 2 )) M(s) is also the same as the above expression for (S-56). Since at least E [min{J t 1 , J t 2 }] does not depend upon t 1 and t 2 through t 1 − t 2 , it is clear that the unconditional covariance structure Cov (f (s, t 1 ), f (s, t 2 )) does not depend upon t 1 and t 2 only through t 1 − t 2 . The proof of (iii) is similar to that of (ii) with s corresponding to t 1 and t 2 replaced by s 1 and s 2 , respectively.
The treatment of the term (S-60) is as follows.
Proof. A sufficient condition for differentiability of (multivariate) functions with multiple