%\documentclass{statsoc}

\documentclass[11pt,a4paper]{article}
%\usepackage{authordate1-4}
% Pacotes utilizados
%\usepackage{graphicx} % Incluir graficos
\usepackage{graphics}
\usepackage[dvipdfm]{graphicx} % Incluir graficos
%%\usepackage[brazil]{babel} % traduz para o portugues Capitulo, resumo, etc...
%%\usepackage[latin1]{inputenc} % `traduz' os acentos do portugues
\usepackage{color}
\usepackage{natbib}

\def\E{\mbox{E}}
\def\Ex{\mbox{Exp}}
\def\U{\mbox{U}}
\def\V{\mbox{V}}
\def\Vec{\mbox{Vec}}
\def\Cov{\mbox{Cov}}
\def\Var{\mbox{Var}}
\def\N{\mbox{N}}
\def\Ga{\mbox{Ga}}
\def\Di{\mbox{Di}}
\def\St{\mbox{St}}
\def\Be{\mbox{Be}}
\def\Bi{\mbox{Bi}}
\def\I{\mbox{I}}
\def\diag{\mbox{diag}}
\def\PT{\mbox{PT}}
\def\Wi{\mbox{Wi}}
\def\GP{\mbox{GP}}
\def\DP{\mbox{DP}}
\def\Pn{\mbox{Pn}}
\def\Pg{\mbox{Pg}}
\def\d{\mbox{d}}


\setlength{\oddsidemargin}{-0.1in}
\setlength{\evensidemargin}{-0.1in}
\setlength{\topmargin}{-0.5in}
\setlength{\parindent}{0cm}
\setlength{\parskip}{2ex}
\textheight=9.0in
\textwidth=6.5in
\renewcommand{\baselinestretch}{1.5}

%%%%%%%%%%%%%%%%
\usepackage[T1]{fontenc} % codifica�o da fonte em 8-bits
\hyphenation{bi-blio-gra-fia}
\sloppy % preferencia a underfull (em vez de overfull)
%%%%%%%%%%%%%%%%%
\newtheorem{teo}{Theorem}[section]
\newtheorem{defin}{Definition}[section]
\newtheorem{prop}{Proposition}[section]
\newtheorem{proof}{Proof}[section]
\newtheorem{cor}{Corollary}[section]
\newtheorem{ape}{Appendix}

\newcommand{\appsection}[1]{\let\oldthesection\thesection
  \renewcommand{\thesection}{Appendix \oldthesection}
  \section{#1}\let\thesection\oldthesection}


\title{Non-Gaussian Bayesian spatiotemporal modelling}

\author{Tha\'{\i}s C.O.~da Fonseca and Mark F.J.~Steel\\
Department of Statistics, University of Warwick, U.K.}


\begin{document}

\maketitle

\abstract{The aim of this work is to construct non-Gaussian and nonseparable covariance functions for processes that vary continuously in space and time. Stochastic modelling of phenomena over space and time is important in many areas of application. But choice of an appropriate model can be difficult as one must take care to use valid covariance structures. We start from a general and flexible way of constructing valid nonseparable covariance functions derived through mixing over separable Gaussian covariance functions. We then generalize the resulting models by allowing for individual outliers as well as regions with larger variances. We induce this through scale mixing with separate positive-valued processes. Smooth mixing processes are applied to the underlying correlated Gaussian processes in space and in time, thus leading to regions in space and time of increased spread. We also apply a separate uncorrelated mixing process to the nugget effect to generate individual outliers. We consider posterior and predictive Bayesian inference with these models. We implement this through a Markov chain Monte Carlo  sampler and apply our modelling approach to temperature data in the Basque country.}

{\it Key words: Bayesian Inference; Flexible tails; Mixtures; Outliers; Temperature data;  Spatiotemporal modelling.}

\section{Introduction}

The motivation of this work is to develop and study non-Gaussian models for processes that vary continuously in space and time. This is a problem of interest in many fields of science such as geology, hydrology and meteorology. Consider the problem of modelling a phenomenon of interest over space and time as a random process
\begin{equation}\label{ProcessZ}
\{Z(s,t);(s,t)\in D\times T\},
\end{equation}
where $(s,t)\in D\times T$, $D\subseteq\Re^d,\;T\subseteq\Re$ are space-time coordinates that vary continuously in $D\times T$. We usually observe a realization of this process at locations $s_i, i=1,\ldots ,I$ and time points $t_j, j=1,\ldots ,J$. The usual assumption for the finite dimensional distributions implied by this process is that, for these spatiotemporal coordinates, the random vector $Z=(Z(s_1,t_1),\ldots,Z(s_I,t_1),\ldots,Z(s_1,t_J),\ldots,Z(s_I,t_J))'$ has a multivariate normal distribution with covariance matrix $\Sigma$ with elements $\Sigma_{kk'}={\rm Cov}(Z_k,Z_{k'})$, $k,k'=1,2,\ldots,N=IJ$. This class is mathematically very convenient, %since the finite dimensional distributions are characterized by the mean vector and covariance matrix. And thus, second-order stationarity is equivalent to full stationarity. On the other hand,
but it is a very restrictive assumption %since the fit might not be satisfactory if
and the data may well present non-Gaussian characteristics. For instance, if there are aberrant observations in the data set it would be useful to consider heavy tailed distributions in order to accommodate these observations. In recent years some models describing departures from Gaussianity for spatial processes were presented in the literature. \citet{Cressie80} discuss robust estimation of the variogram when the distribution has heavier tails than the normal in spatial models. %The fourth root of $[Z(t+h)-Z(t)]^2$ is used in the estimation of the variogram and an almost completely unbiased estimator is proposed. After a simulation study using contaminated data sets the $M$-estimators were chosen because of their robustness in the estimation of the variogram for data sets with outliers.
\citet{DeOliveira97} used nonlinear transformations of random fields in order to accommodate moderate departures from Gaussianity. For instance, their proposal includes the Gaussian and the lognormal models as sampling distributions. %The latter can accommodate skewness and heavier tails than the normal distribution.
\citet{PSteel06} proposed a geostatistical model that accommodates non-Gaussian tail behavior. Their proposal allows the Gaussian case as a limiting case. The proposed class of processes is based on scale mixing a Gaussian process which allows for modelling regions with larger observational variance. Here we consider similar mixing in the nugget effect component allowing for individual outliers. In addition, the ideas in \citet{PSteel06} are here extended to processes in space and time, while avoiding the restrictive assumption of separability between space and time.

In the context of larger observational variance, \citet{Gut01} considered temporally independent samples from a spatiotemporal process $Z(s,t)=\sqrt{\nu(s)}Y_1(s)+Y_2(s,t)$, where $Y_1(s)$ is a spatial process and $Y_2(s,t)$ accounts for the nugget effect. Their model addresses the problem of anisotropy through deformation of the spatial coordinates and  uses Bayesian semi-parametric modelling of the deformation function. The general model potentially also accounts for different variances in space but they adopted the simplifying assumption of constant variances $\nu(s)=\nu,\forall s$. %They comment that in their application to the French precipitation data, the assumption of constant variances in space is not ideal and that they would address this issue in future work.
\citet{Gut03} consider the complete model that incorporates spatial heterogeneity by modelling $\nu(s)$ as latent variables with a log-Gaussian distribution. %The posterior distribution of these latent variables are obtained using a MCMC sampler.

To deal with heterogeneity in time, \citet{Stein08} proposed a model that can account for occasional bursts of increased variability in time. This was done by considering the transformed spatiotemporal process divided by a function of time which was estimated by computing the sample standard deviation at each time point and then smoothed by cubic splines. Notice that this approach does not allow for predictions in time. %The author also assumed that for this data the changes in the variance occur almost simultaneously at all sites. But if this was not the case the model would not be good even for interpolation in space. It is pointed out in the paper that the adjusted process still is not able to capture extreme deviations in the tails.

The model we propose here is able to capture heterogeneous variability both in time and space, as well as outliers in space through a mixed nugget effect. In addition, the covariance structure is nonseparable between space and time. We present an application to maximum temperature data in the Spanish Basque Country, in which the model massively outperforms Gaussian modelling, both in terms of within-sample data support and out-of-sample predictive fit. Moreover, the model easily allows for prediction in space and in time, since we can also predict the mixing processes. We use a Bayesian inferential framework with mildly informative priors. The flexibility of the proposed model does not substantially complicate posterior and predictive inference since conditional on the mixing processes the finite dimensional distributions are all Gaussian distributions.

%We also explore the performance of the proposed model under departures from Gaussianity using data sets contaminated by outliers in several ways.
%Finally, we assess the predictive ability of the proposed models by leaving some testing data out of the estimation procedure and comparing the predictive distributions at these points. This was done using log predictive scores and interval scores that can be applied direct to the $95\%$ posterior credible intervals for the predictions.

\section{Spatiotemporal modelling\label{SecST}}

Building adequate models for processes observed over space and time is not an easy task. Many features have to be considered like stationarity, separability, isotropy and Gaussianity, to mention a few. Adequate specification of the sampling distribution plays an important role in this context since misspecification can lead to poor forecasts/interpolations in space and time. In particular, models based on Gaussianity will not perform well if the data are contaminated by outliers or if there are regions in space or time with larger observational variance. For this reason, we propose a general model able to capture individual outliers as well as regions with larger variance. We use the idea of scale mixing in order to construct processes that imply finite dimensional distribution with heavier tail than the normal distribution.

We consider nonseparable models in space and time generated as proposed in \citet{Fons08}. This construction takes advantage of the models proposed for spatial and temporal processes separately and combines them by using a continuous mixture of separable covariance functions. Let $(U,V)$ be a bivariate nonnegative random vector with distribution $\mu(u,v)$ and independent of $\{Z_1(s);s \in D\}$ and $\{Z_2(t);t \in T\}$ which are purely spatial and temporal random processes, respectively, taken to be independent. Define the process
\begin{equation}\label{mixST}
Z(s,t)=Z_1(s;U)Z_2(t;V),
\end{equation}
where $Z_1(s;u)$ is a purely spatial random process for every $u\in \Re_+$ with covariance function $C_1(s;u)=\sigma_1 \exp\{-\gamma_1(s)u\}$ which is a stationary covariance for $s\in D$ and every $u\in \Re_+$ and a measurable function of $u \in \Re_+$ for every $s\in D$. $Z_2(t;v)$  is a purely temporal random process for every $v\in \Re_+$ with covariance function $C_2(t;v)=\sigma_2 \exp\{-\gamma_2(t)v\}$ which is a stationary covariance for $t\in T$ and every $v\in \Re_+$ and a measurable function of $v \in \Re_+$ for every $t\in T$. $\gamma_1(s)$ is a purely spatial variogram on $D$ and $\gamma_2(t)$ is a purely temporal variogram on $T$. Then the corresponding covariance function of $Z(s,t)$ is a convex combination of separable covariance functions. It is valid \citep[see][]{Ma02a,Ma03b} and generally nonseparable, and is given by
\begin{equation} \label{MixtureCov2}
C(s,t)=\int C_1(s;u)C_2(t;v)d\mu(u,v)
\end{equation}
In particular, if we define $U=X_0+X_1$ and $V=X_0+X_2$ where $X_0$, $X_1$ and $X_2$ are independent nonnegative random variables with finite moment generating functions $M_{0}$, $M_{1}$ and $M_{2}$, respectively, then the resulting covariance function is given by
\begin{equation}\label{CovNonSep}
C(s,t)=\sigma^2M_0(-\gamma_1(s)-\gamma_2(t))\;M_1(-\gamma_1(s))\;M_2(-\gamma_2(t)),
\end{equation}
where $\sigma^2=\sigma_1\sigma_2$. For some interesting classes generated by this approach see \citet{Fons08}. We now consider a more general process $\{\tilde Z(s,t);(s,t)\in D\times T\}$ defined by
\begin{equation}\label{TildeProcess}
\tilde Z(s,t)=\tilde Z_1(s;U)\tilde Z_2(t;V),
\end{equation}
where
\begin{equation}\label{GeneralModSpace}
\tilde Z_1(s;U) = \sqrt{1-\tau^2}\frac{Z_1(s;U)}{\sqrt{\lambda_1(s)}}+\tau \frac{\epsilon(s)}{\sqrt{h(s)}}
\end{equation}
where $\{\lambda_1(s);s\in D\}$ is a mixing process which is uncorrelated to $\epsilon(s)$ and $Z_1(s;u)$. $\{\epsilon(s);s\in D\}$ denotes an uncorrelated Gaussian process with zero mean and unitary variance which introduces a nugget effect parameterised by $\tau$. $\{h(s);s\in D\}$ is also an uncorrelated process in $\Re_+$ with distribution $P_{h}$. The mixing process $\lambda_1(s)$ is spatially correlated and responsible for allowing regions in space with larger variances while the process $h(s)$ can create traditional outliers, {\it i.e.}~observations with unusually large nugget effects. We also want to allow larger observational variances in time, thus we consider the following process in time
\begin{equation}\label{GeneralModTime}
\tilde Z_2(t;V) = \frac{Z_2(t;V)}{\sqrt{\lambda_2(t)}}
\end{equation}
where $\{\lambda_2(t);t\in T\}$ is a mixing process which is uncorrelated with $Z_2(t;v)$. The covariance function for the process $\{\tilde Z(s,t);(s,t)\in D\times T\}$ is given by
\begin{equation}\label{GeneralCov2}
\tilde C(s,t)={\rm Cov}(\tilde Z(s_0,t_0),\tilde Z(s_0+s,t_0+t))=\E_{U,V}\{\tilde C_1(s;U)\tilde C_2(t;V)\},
\end{equation}
where $s,s+s_0\in D$ and $t,t+t_0\in T$, $\tilde C_1(s;u)={\rm Cov}(\tilde Z_1(s_0;u),\tilde Z_1(s_0+s;u))$ and $\tilde C_2(t;v)={\rm Cov}(\tilde Z_2(t_0;v),\tilde Z_2(t_0+t;v))$.



\section{Scale mixing in space \label{ScaleSpace}}

In this section we consider scale mixing in the space dimension. This will account for individual outliers (through the process $h(s)$) and regions in space with larger observational variance (through the process $\lambda_1(s)$). The latter is quite common in meteorological applications where outliers are often associated with severe weather events such as tornados and hurricanes. \citet{Chen07} pointed out that these events do not usually happen in a single location but cover an extended region.

Initially we consider the case where $\lambda_2(t)=1,\;\forall t\in T$. The results in \citet{PSteel06} are directly applicable to the process $\tilde Z_1(s;u)$, that is a purely spatial process. Firstly, the mixing process $\lambda_1(s)$ needs to be correlated to induce mean square continuity in the process $\tilde Z_1(s;u)$ for $\tau=0$, that is, we need to correlate the mixing variables so that locations that are close together will have very similar values of $\lambda_1(s)$. The simplest way to do this is to consider a common mixing variable $\lambda_1(s)=\lambda\sim P_{\lambda}$. Then we have that
\begin{equation}
\tilde C_1(s;u)=(1-\tau^2)\E\left[\lambda^{-1}\right ]C_1(s;u)+\tau^2\E\left [h^{-1}\right ]I_{(s=0)}
\end{equation}
$\E[h^{-1}]=\E[h(s)^{-1}]$, $\forall s$, where $h\sim P_h$. Solving the integral in (\ref{GeneralCov2}) and assuming throughout that $\sigma^2=\sigma_2$ and $\sigma_1=1$ we obtain
\begin{equation}\label{covSpace}
\tilde C(s,t)=\sigma^2M_0(-\gamma_1(s)-\gamma_2(t))\;\tilde M_1(-\gamma_1(s))\;M_2(-\gamma_2(t))
\end{equation}
where
$\tilde M_1(-\gamma_1)=(1-\tau^2)\E[\lambda^{-1}]M_1(-\gamma_1(s))+\tau^2\E[h^{-1}]I_{(s=0)}$. Therefore the correlation structure for $s>0$ is given by
\begin{equation}
\tilde \rho(s,t)=\frac{\tilde C(s,t)}{\tilde C(0,0)}=\frac{\E[\lambda^{-1}]}{\E[\lambda^{-1}]+w^2\E[h(s)^{-1}]}M_0(-\gamma_1(s)-\gamma_2(t))\;M_1(-\gamma_1(s))\;M_2(-\gamma_2(t))
\end{equation}
where $w^2={\tau^2}/{(1-\tau^2)}$. When $\tau^2=0$ (no nugget effect) the mixing does not affect the correlation structure, that is, $\tilde\rho(s,t)=\rho(s,t)$, where $\rho(s,t)=C(s,t)/C(0,0)$ is the correlation function of $\{Z(s,t),(s,t)\in D\times T\}$. For instance, if we take $P_{\lambda}={\rm Ga}(\nu/2,\nu/2)$ then the unconditional distribution of $\tilde Z=({\tilde Z}(s_1,t_1),\ldots,{\tilde Z}(s_I,t_J))'$ is $IJ$-variate student-$t$ with $\nu$ degrees of freedom. \citet{Omre06} presents some characteristics of student-$t$ random fields.

But setting $\lambda_1(s)=\lambda\sim P_{\lambda}$ is an extreme situation and we would like to have individual mixing variables that account for spatial heterogeneity. Thus, we consider a process $\lambda_1^{-1/2}(s)$ that is mean square continuous which implies that $\tilde Z_1(s;u)$ is also mean square continuous. This means that we need to satisfy the condition $\E\left[\left(\lambda_1^{-1/2}(s_i)-\lambda_1^{-1/2}(s_{i'})\right)^2\right]\to 0$ as $s_i\to s_{i'}$. %This condition is equivalent to $\E[(\lambda_1^{-1/2}(s_i)\lambda_1^{-1/2}(s_{i'}))^2]\to \E[\lambda_1^{-1}(s_i)]$ as $s_i\to s_{i'}$.
The latter is satisfied for the log-Gaussian process proposed in \citet{PSteel06} where $\{\ln(\lambda_1(s));s\in D\}$ is a Gaussian process with mean $-\frac{\nu}{2}$ and covariance structure $\nu C^*_1(.)$. This implies a lognormal distribution for each $\lambda_1(s_i)$ with mean one and $\Var[\lambda_1(s_i)]=\exp(\nu)-1$, so that the marginal distribution becomes more spread out as $\nu$ increases. For large $\nu$ the distribution also becomes more right-skewed with the mode shifting towards zero, allowing for substantial variance inflation for some spatial regions. The spatial covariance function is then given by
\begin{equation}\label{C1loggaus}
\tilde C_1(s;u)=(1-\tau^2)\exp\left\{ \frac{\nu}{4}\left [C^*_1(s)-1\right ]+\nu \right\}C_1(s;u)+\tau^2\E\left [h^{-1}\right ]I_{(s=0)},
\end{equation}
and the resulting covariance function for the spatiotemporal process obtained by solving (\ref{GeneralCov2}) is given by (\ref{covSpace})
%\begin{equation} \label{covSpace}
%\tilde C(s,t)=\sigma^2 M_0(-\gamma_1(s)-\gamma_2(t))\;\tilde %M_1(-\gamma_1(s))\;M_2(-\gamma_2(t))
%\end{equation}
with %$\sigma^2=\sigma_2$, $\sigma_1=1$ and
$\tilde M_1(-\gamma_1)=(1-\tau^2)\exp\left\{ \frac{\nu}{4}\left [C^*_1(s)-1\right ]+\nu \right\}M_1(-\gamma_1)+\tau^2\E[h^{-1}]I_{(s=0)}$. Therefore the correlation structure for $s>0$ is given by
\begin{equation}
\tilde \rho(s,t)=\frac{\exp\left\{ \frac{\nu}{4}\left [C^*_1(s)-1\right ]+\nu \right\}}{\exp\{\nu\}+w^2\E[h^{-1}]}M_0(-\gamma_1(s)-\gamma_2(t))\;M_1(-\gamma_1(s))\;M_2(-\gamma_2(t)).
\end{equation}
Throughout, we will use $C^*_1(s)=M_1(-\gamma_1(s))$. If, in addition, we take $h(s)=1$ we have the same model presented in \citet{PSteel06} for the space dimension. As commented in \citet{PSteel06}, we could use a different covariance function $C^*_1(.)$ for the mixing process but then we would need to estimate the parameters in $C^*_1(.)$ and this might not be easy. This is still a subject for further research as well as the study of the effect of this assumption in the estimation of the separability parameter, for example.

The smoothness properties presented in \citet{PSteel06} and \citet{Fons08} extend to this framework. In the following, $f^{(q)}(x)$ will denote the $q^{th}$ derivative of a function $f(x)$ with respect to $x$.
\begin{prop}\label{Prop1}
In the case without nugget effect and with $C^*_1(s)=M_1(-\gamma_1(s))$, the purely spatial process $\{\tilde Z(s,t_0),s\in D\}$ at a fixed time point $t_0\in T$ is $m$ times mean square differentiable if and only if $M_1^{(2m)}(r)$ and $\gamma_1^{(2m)}(s)$ exist and are finite at $0$.
\end{prop}



\subsection{Process $h(s)$}

We define the process $\{h(s);s\in D\}$ as an uncorrelated mixing process that allows for larger nugget effects. This would mean a traditional outlier observation originating from e.g.~measurement error. The process is uncorrelated with $\epsilon(s)$, $Z_1(s;U)$ and $Z_2(t;V)$ and the finite dimensional distribution is $P_{h}$.

Aberrant observations are common in time series analysis and might also be encountered in processes observed in space. Therefore it is essential to consider this possibility when modelling phenomena over space-time. It is important to understand the effect of outliers in the estimation of the parameters in the correlation structure as this will affect directly the predictions. %Most of the proposed methods for detecting outliers in the time series literature assume ARIMA type models, for reference see \citet{Tsay00} where multivariate time series are considered. \citet{Pena01} considered multivariate outlier detection and robustness of the covariance matrix estimation. The latter paper proposes the analysis of projections of the data onto a set of directions obtained minimizing and maximizing the kurtosis coefficient of the projections. After identifying the outliers, the mean and covariance matrix are estimated using the observations that are not labeled as outliers. In an applied setting, \citet{Davis06} developed a detection procedure for multivariate outliers using a control chart based on the statistic $T^2=(Z-\mu)'\Sigma^{-1}(Z-\mu)$. First the data in space and time is modeled considering seasonality, heteroscedasticity and variogram fitting. After that, the residuals are considered for outlier detection. It is pointed out that under this procedure multiple out-of-control points close in time suggest the need for model adaptation or re-estimation of the parameters.

We consider the detection of outliers jointly in the estimation procedure. Therefore the variables $h_i=h(s_i)$, $i=1,\ldots, I$ are considered latent variables and the posterior distribution provides an indication of outlying observations. If the marginal posterior distribution of $h_i$ has a lot of mass close to 0, this indicates inflation of the scale $\tau^2$ and therefore an outlying observation. We consider the following distributions for $h_i$, $i=1,\ldots,I$ where $\nu_h>0$:
\begin{enumerate}
\item $\ln( h_i)\sim {\rm N}(-\nu_{h}/2,\nu_{h})$, that is, $\E[h_i]=1$ and $\Var(h_i)=\exp(\nu_h)-1$. If $\nu_h$ is close to 0 then the distribution of $h_i$ is very tight around 1.
\item $h_i\sim {\rm Ga}(1/\nu_h,1/\nu_h)$, where ${\rm Ga}(a,b)$ denotes the Gamma distribution with density function $f(x)=b^a/\Gamma(a)x^{a-1}\exp(-bx)$, so that $\E[h_i]=1$ and $\Var(h_i)=\nu_h$. If $\nu_h$ is close to 0 then again the distribution of $h_i$ is very tight around 1.
\end{enumerate}

In order to evaluate the tail behavior of the finite dimensional distribution of the proposed process we find the kurtosis coefficient that is given by ${\E[\tilde Z^{4}_{ij}]}/{\E^2[\tilde Z^{2}_{ij}]}$, where $\tilde Z_{ij}=\tilde Z(s_i,t_j)$, $i=1,\ldots,I$ and $j=1,\ldots,J$. The kurtosis coefficient of the marginal finite dimensional distributions implied by the process defined in (\ref{TildeProcess})-(\ref{GeneralModSpace}) in combination with a log-Gaussian $\lambda_1(s)$ is given by
\begin{equation}
{\rm kurt}[\tilde Z_{ij}] = \frac{3\{\exp(3\nu)+2w^2\exp(\nu)\E[h_i^{-1}]+w^4\E[h_i^{-2}]\}}{\exp(2\nu)+2w^2\exp(\nu)\E[h_i^{-1}]+w^4\E^2[h_i^{-1}]}.
\end{equation}
Notice that when $\tau^2=0$, that is, if there is no nugget effect then the kurtosis is given by $3\exp(\nu)$ as in \citet{PSteel06}. In the case of $\ln( h_i)\sim {\rm N}(-\nu_{h}/2,\nu_{h})$, $\E[h_i^{-1}]=\exp(\nu_h)$ and $\E[h_i^{-2}]=\exp(3\nu_h)$, for any $\nu_h>0$. If instead $h_i\sim {\rm Ga}(1/\nu_h,1/\nu_h)$ then $\E[h_i^{-1}]={1}/{1-\nu_h}$ and $\E[h_i^{-2}]={1}/{(1-\nu_h)(1-2\nu_h)}$. The latter case requires that $\nu_h<0.5$. Figure \ref{fg1} shows the implied kurtosis for several values of $w^2$ for both models when $\nu=0.5$. The Gamma distribution for $P_h$ looks less flexible since it gives kurtosis very close to $3\exp (\nu )$ for almost all values of $\nu_h$. We need to go very close to $\nu_h=0.5$ to get larger values for the kurtosis coefficient. For $\nu=0$ the kurtosis is an increasing function of $\nu_h$ whereas for $\nu>0$ it is not monotonous, as illustrated in Figure \ref{fg1}.

\begin{figure}[htb]
\center
\begin{tabular}{cc}
%\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/KurtosisLogNormal.pdf} &
%\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/KurtosisGamma.pdf} \\
% (a) $h_i$ Lognormal and $\nu=0$. & (b) $h_i$ Gamma and $\nu=0$. \\
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/KurtosisLogNormalc.pdf} &
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/KurtosisGammad.pdf} \\
 (a) $h_i$ Lognormal & (b) $h_i$ Gamma  \\
\end{tabular}
\caption{Kurtosis for different values of $w^2$ and $\nu=0.5$.}
\label{fg1}
\end{figure}

%\clearpage
%\newpage

\subsection{Parameterisation\label{secTheta}}

In what follows we consider some particular choices for the variogram functions $\gamma_1(s)$, $\gamma_2(t)$ and variables $X_0$, $X_1$ and $X_2$ introduced in Section \ref{SecST}. Consider $\gamma_1(s) = ||s/a||^{\alpha}$, $\gamma_2(t) = |t/b|^{\beta}$, $a,b>0$, $\alpha\in (0, 2]$ and $\beta\in (0, 2]$. Let $X_0\sim {\rm Ga}(\eta_0,1)$, $X_2\sim {\rm Ga}(\eta_2,1)$, which results in Cauchy covariance functions, that is, $M_0(x) = (1-x)^{-\eta_0}$ and  $M_2(x) = (1-x)^{-\eta_2}$, $\eta_0,\eta_2>0$. For $X_1$ we consider the following distributions
\begin{enumerate}
\item $X_1\sim {\rm Ga}(\eta_1,1)$ resulting in the Cauchy function $M_1(x) = (1-x)^{-\eta_1}$, $\eta_1>0$;
\item $X_1\sim {\rm InvGa}(\eta_1,1)$ resulting in the Mat\'ern function $M_1(x) = \frac{2(-x)^{\eta_1/2}}{2^{\eta_1-1}\Gamma(\eta_1)}K_{\eta_1}(2\sqrt{-x})$, $\eta_1>0$. See \citet{Stein99} for details of this class of covariance functions;
\item $X_1\sim {\rm GIG}(\eta_1,\delta,\delta)$ resulting in the Generalized Mat\'ern function $M_1(x) = (1-\frac{x}{\delta})^{-\eta_1/2}{K_{\eta_1}(2\delta\sqrt{1-\frac{x}{\delta}})}/{K_{\eta_1}(2\delta)}$, $\eta_1\in \Re$, $\delta>0$. See \citet{Shka68} for details of this class of covariance functions.
\end{enumerate}
Here ${\rm InvGa}(\eta_1,1)$ denotes the Inverse Gamma distribution with density function $f(x)=1/\Gamma(\eta_1)x^{-\eta_1-1}\exp(-1/x)$ and ${\rm GIG}(\eta_1,\delta,\delta)$ denotes the Generalized Inverse Gaussian distribution with density function $f(x)=1/(2K_{\eta_1}(2\delta))x^{\eta_1-1}\exp\{-(\delta x+\delta/x)\}$. We use the correlation between the variables $U$ and $V$ as an indication of interaction between space and time components. This correlation is given by
\begin{equation}\label{correl}
c=\frac{\eta_0}{\sqrt{(\eta_0+V_1)(\eta_0+\eta_2)}},
\end{equation}
where $V_1={\Var}(X_1)$. Thus, $0\leq c\leq 1$ could be used as a measure of space-time interaction, with $c=0$ indicating separability and $c=1$ meaning high dependence between space and time. Notice that in the case $X_1\sim {\rm InvGa}(\eta_1,1)$ the variance of $X_1$ does not exist (unless $\eta_1>2$ is imposed through the prior) and the dependence between space and time is then measured by
\begin{equation}
\tilde c=\frac{\eta_0}{\sqrt{(\eta_0+\tilde V_1(\eta_1))(\eta_0+\eta_2)}},
\end{equation}
where $\tilde V_1(\eta_1)=(Q(0.75;\eta_1)-Q(0.25;\eta_1))^2$ and $Q(x;\eta_1)$ is the quantile of $X_1$ corresponding to $100 x\%$.

\subsection{Inference\label{SecInf}}

Suppose we observe realizations $\tilde z_{ij}$ of $\tilde Z_{ij}\equiv\tilde Z(s_i,t_j)$ at locations $s_i$, $i=1,\ldots,I$ and time points $t_j$, $j=1,\ldots,J$. Defining $\lambda_1=(\lambda_1(s_1),\dots,\lambda_1(s_I)), h=(h(s_1),\dots,h(s_I))$ and $\mu=(\mu(s_1,t_1,\dots\mu(s_I,t_J))'$ a location function, the likelihood function is given by
\begin{equation}\label{likelihoodSpace}
L_1(\theta,\lambda_1,h;\tilde z)=f_N({\tilde z}\mid \mu,\Sigma),
\end{equation}
where $\theta=(\eta_0,\eta_1,\eta_2,\alpha,\beta,a,b,\sigma^2,\tau^2)$ or $\theta=(\delta,\eta_0,\eta_1,\eta_2,\alpha,\beta,a,b,\sigma^2,\tau^2)$ depending on which distribution we choose for $X_1$, $f_N(.|\mu,\Sigma)$ denotes the multivariate Gaussian density function with  mean $\mu$ and covariance matrix $\Sigma$ and
$\Sigma_{kk'}={\rm Cov}[(\tilde Z)_k,(\tilde Z)_{k'}]$ where
$${\rm Cov}(\tilde Z_{ij},\tilde Z_{i'j'})=\sigma^2M_0(-\gamma_1-\gamma_2)\left [(1-\tau^2)\frac{M_1(-\gamma_1)}{\sqrt{\lambda_{1i}\lambda_{1i'}}}+\tau^2\frac{I({s_i=s_{i'})}}{\sqrt{h_ih_{i'}}}\right ]M_2(-\gamma_2),$$
$i,i'=1,\ldots, I$, $j,j'=1,\ldots, J$, $\gamma_1=\gamma_1(s_i-s_{i'})$, $\gamma_2=\gamma_2(t_j-t_{j'})$, $\lambda_{1i}=\lambda_1(s_i)$ and $(\tilde Z)_k$ is the $k^{th}$ element of the $IJ$-dimensional vector $\tilde Z$.

%$$\Sigma_{kk'}=\sigma^2M_0(-\gamma-\gamma_2)\left [(1-\tau^2)\frac{M_1(-\gamma_1)}{\sqrt{\lambda_{1k}^*\lambda_{1k'}^*}}+\tau^2\frac{I({s_k^*=s_{k'}^*)}}{\sqrt{h_k^*h_{k'}^*}}\right ]M_2(-\gamma_2),k,k'=1,\ldots, IJ$$
%$\lambda_1^*=1_J\otimes \lambda_1$, $\lambda_1=(\lambda_1(s_1),\ldots ,\lambda_1(s_I))^t$ and $h^*=1_J\otimes h$, $h=(h(s_1),\ldots ,h(s_I))^t$. $\gamma_1=\gamma_1(s_k^*-s_{k'}^*)$, $s^*=1_J\otimes s$ and $\gamma_2=\gamma_2(t_k^*-t_{k'}^*)$, $t^*=t \otimes 1_I$.
Note that Gaussian behaviour is only assumed given the mixing variables $\lambda_1$ and $h$. Integrating out with respect to these mixing variables leads to non-Gaussian distributions. We augment with the latent variables $\lambda_1$ and $h$ in order to identify possible regions with larger observational variance and/or traditional outliers. The vector $\ln(\lambda_1)$ is multivariate normally distributed with mean $-\nu/2$ and covariance matrix with elements $\nu M_1(-\gamma_1(s_i-s_{i'}))$. And $\ln(h_i),\;i=1,\ldots, I$ is either normal with mean $-\nu_h/2$ and variance $\nu_h$ or $h_i\sim\Ga(1/\nu_h,1/\nu_h)$.


We use stochastic simulation via MCMC to obtain an approximation of the posterior distribution of $(\theta,\lambda_1,h)$. We obtain samples from the target distribution $p(\theta,\lambda_1,h|z)$ by successive generations from the full conditional distributions. More specifically, we use a hybrid Gibbs sampler scheme with Metropolis-Hastings steps. We use random walk proposals to generate values of $\lambda_1$ and $h$. We also consider groups in space in order to define blocks in the sampler. For a more elaborate algorithm see \citet{PSteel06}.

Model comparison is conducted on the basis of Bayes factors. These are computed from the MCMC output using methods to approximate the marginal predictive density of $z$. In previous simulated studies \citep{Fons08} we noticed that the estimator $p_4$ of \citet{NewR94} (with their $d$ as small as $0.01$), the  optimal bridge sampling approach of \citet{meng}, and the shifted Gamma estimator proposed by \citet{Raft07} (with values of their $\lambda_1$ close to one) give essentially the same results, especially the last two. %As the latter are relatively easily calculated. %Since we perform many simulations and the calculation of the bridge sampling estimator is more time demanding we use only the estimator $p_4$ and the shifted gamma estimator in this work.

\subsection{Prediction and interpolation}

{\it Gaussian model}

Suppose we are interested in predicting $Z_p$ at a location $s_p$ at time point $t_p$, where $s_p$ is not necessarily included within the sampling design. Under a Bayesian approach, the prediction of $Z_p$ is based on the posterior distribution $P(Z_p|Z_o)$ where $Z_o$ are the available observations of the process $Z(s,t)$. We have that
\begin{equation}
P(Z_p|Z_o)=\int P(Z_p|Z_o,\theta) P(\theta|Z_o) d\theta
\end{equation}
Under the Gaussian model, $(Z_o,Z_p|\theta)$ has a multivariate Gaussian distribution. Using the properties of conditional distributions in the Gaussian family we obtain the mean and the variance for $Z_p|Z_o,\theta$ given by
\begin{equation}
\bar \mu=\mu_p+\Sigma_{po}\Sigma_{oo}^{-1}(z_{o}-\mu_o)
\end{equation}
and
\begin{equation}
\bar \Sigma=\Sigma_{pp}-\Sigma_{po}\Sigma_{oo}^{-1}\Sigma_{op}
\end{equation}
Suppose we have $\theta^{(1)},\ldots,\theta^{(M)}$ generated from the posterior distribution $\theta|Z_o$ (by the MCMC sampler). Then we approximate the posterior distribution of $Z_{p}$ by generating $Z_p^{(k)}$ from $P(Z_p|Z_o,\theta^{(k)})$, for $k=1,\ldots,M$.

{\it Non-Gaussian model}

Suppose we have samples from the distribution of $\theta,\lambda_{1o}|Z_o$. We need to sample $\lambda_{1p}|\theta,\lambda_{1o},Z_o$ and then use the Gaussian sampler just described to sample $Z_p|Z_o,\theta,\lambda_{1p},\lambda_{1o}$. Notice that $p(\lambda_{1p}|\theta,\lambda_{1o},Z_o)\equiv p(\lambda_{1p}|\theta,\lambda_{1o})$, where
\begin{equation}
(\ln\lambda_{1o},\ln\lambda_{1p})'|\theta\sim N(-\nu/2,\nu M_1),
\end{equation}
with $M_{1ij}=M_1(-\gamma_1(s_i-s_j))$. Thus, given $\lambda_{1o}$ and $\theta$ we can easily generate values of $\lambda_{1p}$.

{\it Predictive model comparison}

In order to check the predictive accuracy of each model we use a predictive scoring rule. Scoring rules provide summaries for the evaluation of probabilistic forecasts by assigning a number associated with the prediction and the actual value observed for the process. For more details about scoring rules see \citet{GneitRaf07}. %Reversing the sign of the scoring rule,
In particular, we use the log predictive score (LPS) based on the predictive distribution $p$ (which can be multivariate) and on the observed value $x$,
\begin{equation}\label{lps}
LPS(p;x)=-\ln(p(x)).
\end{equation}
The smaller LPS  is, the better the model does in forecasting $Z_p$.

%In most spatiotemporal applications the pratitioner will be interested in reporting prediction intervals. If the forecaster quotes the $(1-\xi)100\%$ interval $(q_1,q_2)$ where $q_i=quantile(p,\xi_i)$ and $x$ materializes, then he or she will be rewarded by the score $S_p(q_1,q_2;x)$. Setting $\xi_1=\xi/2$ and $\xi_2=1-\xi/2$ and reversing the sign of the scoring rule, we get the negatively oriented interval score,
%\begin{equation} \label{eqIS}
%IS(q_1,q_2;x)=(q_2-q_1)+\frac{2}{\xi}(q_1-x)I(x<q_1)+\frac{2}{\xi}(x-q_2)I(x>q_2).
%\end{equation}
%This score rewards narrow prediction intervals and penalizes (depending on $\xi$) if the observation misses the interval. Usually, we use $\xi=0.05$ resulting in a $95\%$ credible interval.



\section{Scale mixing in time \label{ScaleTime}}

We also want to consider the case where $\lambda_2(t)\not =1$. Let $\{\lambda_2(t);t\in T\}$ be a mixing process in time. %If we have that $\lambda_1(s) =1,\;\forall s\in D$ then the results in section \ref{ScaleSpace} apply directly to this situation. Otherwise,
Notice that the process $\{Y(s,t);(s,t)\in D\times T\}$ where
\begin{equation}\label{ProcY}
Y(s,t)=\lambda_2(t)^{1/2}\tilde Z(s,t)
 \end{equation}
is exactly the process proposed in section \ref{ScaleSpace}. We derive properties and inference by using this fact. We will not consider the case with temporal nugget effect but the model presented here can be easily extended to allow for a nugget effect in time as done in the spatial dimension.

Defining  $\lambda_2=(\lambda_2(t_1),\dots,\lambda_2(t_J))$ and using (\ref{likelihoodSpace}) and (\ref{ProcY}) the resulting likelihood function is given by
\begin{equation}\label{likelihoodTime}
L(\theta,\lambda_1,h,\lambda_2 ;\tilde z)=L_1(\theta,\lambda_1,h ;\Lambda_2 \tilde z)\prod_{j=1}^{J} \lambda_{2j}^{I/2},
\end{equation}
where $\Lambda_2=\diag((\lambda_{21}^{1/2},\ldots ,\lambda_{2J}^{1/2})'\otimes 1_I)$, $\lambda_{2j}=\lambda_2(t_j)$  and $1_I$ is a vector of ones of size I.

We use in the time dimension the same kind of mixing process used in the space, where $\{\ln(\lambda_2(t)):t\in T\}$ is a Gaussian process with mean $-{\nu_2}/{2}$ and covariance structure $\nu_2 C^*_2(.)$. The temporal covariance function is given by
\begin{equation}\label{C2loggaus}
\tilde C_2(t;v)=\exp\left\{ \frac{\nu_2}{4}\left [C^*_2(t)-1\right ]+\nu_2 \right\}C_2(t;v)
\end{equation}
And the resulting covariance function for the spatiotemporal process obtained by solving (\ref{GeneralCov2}) is given by
\begin{equation}
C(s,t)=M_0(-\gamma_1(s)-\gamma_2(t))\;\tilde M_1(-\gamma_1(s))\;\tilde M_2(-\gamma_2(t))
\end{equation}
where $\tilde M_1(-\gamma_1(s))=(1-\tau^2)\exp\left\{ \frac{\nu}{4}\left [C^*_1(s)-1\right ]+\nu \right\}M_1(-\gamma_1(s))+\tau^2\E[h(s)^{-1}]I_{(s=0)}$ and $\tilde M_2(-\gamma_2(t))=\exp\left\{ \frac{\nu_2}{4}\left [C^*_2(t)-1\right ]+\nu_2 \right\}M_2(-\gamma_2(t))$. Therefore, the correlation structure for $s,t>0$ is given by
\begin{equation}
\tilde \rho(s,t)=\frac{\exp\left\{ \frac{\nu}{4}\left [C^*_1(s)-1\right ]+\nu+\frac{\nu_2}{4}\left [C^*_2(t)-1\right ]+\nu_2 \right\}}{[\exp\{\nu\}+w^2\E[h(s)^{-1}]]\exp\{\nu_2\}}M_0(-\gamma_1(s)-\gamma_2(t))\;M_1(-\gamma_1(s))\;M_2(-\gamma_2(t)),
\end{equation}
%where $w^2=\frac{\tau^2}{(1-\tau^2)}$.
where we set $C^*_1(s)=M_1(-\gamma_1(s))$ and $C^*_2(t)=M_2(-\gamma_2(t))$.

This scale mixing in time will capture periods in time with larger observational variance, which can be seen as a way to address the issue of volatility clustering, which is quite a common occurrence in {\it e.g.}~financial time series data. This aspect of our model is reminiscent of a stochastic volatility model used in this literature.

To summarize, the full model, with scale mixing through $\lambda_1, \lambda_2$ and $h$ will be able to accommodate smooth spatial heterogeneity in the variance (through $\lambda_1$), gradual temporal changes in the variance (through $\lambda_2$) and spatial outliers (through the fat-tailed distribution of the nugget effect, induced by mixing with $h$).

%\section{Interaction in the mixing variables?}

%\section{Other distributions for $\lambda$?}

%\section{Mixing using other ST models presented in the literature}


\section{Empirical Results \label{secSim}}

Throughout, the prior distribution for the hyperparameters $\nu, \nu_h$ and $\nu_2$ **is this true?** are ${\rm Ga}(1,5)$, assigning a large probability mass to values close to zero which indicate simpler cases (without scale mixing, {i.e.}~no spatial or temporal heterogeneity or fat tails). The prior distribution for the remaining parameters is mildly informative as discussed in \citet{Fons08}.


\subsection{Simulation results}

We have analysed a substantial number of generated datasets, with and without perturbations with respect to Gaussianity. This has illustrated that the priors are reasonable and the inference methods are reliable and  efficient.  In addition, we can successfully identify outliers and regions of increased variance in time and space. In fact, we are able to separate these effects even when they are all present in the data. The use of Bayes factors leads to sensible model choices, in line with the way the data were generated.

\subsection{Application to temperature data}

We now present an application to the maximum temperatures recorded daily in July of 2006 ($J=31$) in 70 locations within the Spanish Basque country. We consider $I=67$ of these locations for estimation of the parameters and we leave out 3 locations for predictive comparison. %The spatial coordinates are shown in Figure \ref{fgTempMax0}.
Figure \ref{fgTempMax1} (a)-(b) present the boxplots for the maximum temperature (in degrees centigrade) over space and time, respectively. Notice in Figure \ref{fgTempMax1} (c) how large the variability for each station is. %This indicates that the temperature changes a lot over time, sometimes it changes 10 degrees in a couple of days.
In order to model nonstationarities in the mean of the process we considered a mean function that depends on the spatial coordinates and on time. %$s=(s_1,s_2)$ in the UTM coordinate system (divided by $10^5$) and on the time points. Notice that in order to measure distances in kilometers we just need to multiply de coordinates by $100$. The median of the distance between stations is 43 kilometers, the minimum distance is 2.4 kilometers and the maximum distance is 127 kilometers.
In addition, the region considered is quite mountainous with altitudes ranging from 0 to 1188 meters, therefore the altitude ($h$) was included in the mean function as a covariate. The resulting mean function is given by
\begin{equation} \label{meanTemp}
\mu(s,t)=\delta_0+\delta_1 s_1+\delta_2 s_2+\delta_3 h+\delta_4 t+ \delta_5 t^2
\end{equation}
The covariance model considered here is nonseparable as presented in (\ref{CovNonSep}) allowing for interactions between space and time. The model is parameterized as  in Subsection \ref{secTheta} and the chosen covariance in space is the Cauchy type. The parameters estimated are the trend coefficients $(\delta_0,\ldots,\delta_{4})$, the covariance parameters $(\eta_0,\alpha,\beta,a,b,\sigma^2,\tau^2)$ and the mixing parameters $(\nu_1,\nu_2,\nu_h)$. We also generate the auxiliary variables $(\lambda_1,\lambda_2,h)$ in our MCMC algorithm as described in Subsection \ref{SecInf} in order to identify regions in space and time with larger observational variance. The parameters $\eta_1$ and $\eta_2$ are set to $1$. See \citet{Fons08} for a discussion of parameterisation in this class of covariance models.
For the mixing on the nugget effect we only consider lognormally distributed $h_i$, as this seems the most flexible option (see the discussion at the end of Subsection 3.1).

In order to calculate the likelihood function we need to invert a matrix with dimension $2077\times 2077$. Since this is very time demanding we approximate the likelihood by using conditional distributions as described in \citet{Stein04}. In summary, we consider a partition of $Z$ into subvectors $Z_1,... ,Z_{31}$ where $Z_j=(Z(s_1,t_j),\ldots,Z(s_{67},t_j))'$ and we define $Z_{(j)}=(Z_{j-L+1},... ,Z_j)$. Then, taking $\phi=(\theta,\lambda_1,\lambda_2,h)$, we use
\begin{equation}
p(z|\phi)\approx p(z_1|\phi)\prod_{j=2}^{31}p(z_j|z_{(j-1)},\phi).
\end{equation}
Since we have a natural time ordering of $Z$, this means the distribution of $Z_{j}$ will only depend on the observations in space for the previous $L$ time points. In this application we used $L=5$ to make the MCMC feasible.

We estimated the parameters for the following models: the Gaussian model, the non-Gaussian model with $\lambda_1$ only, the non-Gaussian model with lognormal $h$, the non-Gaussian model with lognormal $h$ and $\lambda_1$, the non-Gaussian model with $\lambda_2$ only, the non-Gaussian model with $\lambda_2$ and $h$, the non-Gaussian model with $\lambda_1$ and $\lambda_2$ and finally the non-Gaussian model with lognormal $h$, $\lambda_1$ and $\lambda_2$. Notice that all the models considered here have a nugget effect in space parameterized by $\tau^2>0$ that accounts for measurement errors. The estimated Bayes Factors presented in Table \ref{ComparaTemp} indicate the non-Gaussian models are much more adequate for this dataset than the Gaussian model. The complete model that includes $\lambda_1$, $\lambda_2$ and $h$ is the best one according to both estimators. We expected the models that include $\lambda_1$ and $\lambda_2$ to give good results since the empirical variances seems to change a lot in space and time as shown in Figure \ref{fgTempMax1} and this was indeed the case.

Table \ref{StatMaxDelta1} which presents some posterior summaries for the parameters of interest. The models with $\lambda_1$ give very different results for the smoothness parameter in space $\alpha$, the range in space $a$ and the nugget effect $\tau^2$. The models without $\lambda_1$ suggest smoother processes with much larger values of $\alpha$ than the models that include $\lambda_1$. This may be related to the double use of $M_1(-\gamma_1(s))$ for the correlation in space of the observables and $\ln(\lambda_1)$.  The same does not happen with $\beta$ for the models including $\lambda_2$. Inference on the separability measure $c$ in (\ref{correl}) is relatively unaffected by the model choice and indicates that the data are fairly close to separable. Notice that the process is very rough in time with very small values of $\beta$. For $\beta<1$ the process is not even mean square continuous in time.

The posterior distributions of the parameters of the mixing distributions, which drive the tail behaviour, are depicted in Figure \ref{fgnus}. It is clear that the posteriors are very different from the priors (indicated by dashed lines), and point very strongly towards non-Gaussian behaviour. Remember that the Gaussian model corresponds to the limiting case where all $\nu$ parameters tend towards zero. The posteriors clearly suggest that all three forms of scale mixing in the model are supported by the data, especially for the smooth spatial process. Of course, this is in line with the Bayes factors discussed earlier.

Inference on the coefficients in the mean function (\ref{meanTemp}) show that altitude is an important covariate with a similar effect for all the models (in particular, it indicates a drop in mean level of about 0.8 of a degree centigrade per 100 meters altitude). %On the other hand, only the model with $\lambda_1$ and $\lambda_2$ and the model with $h$, $\lambda_1$ and $\lambda_2$ suggest that the latitude $s_2$ is significant, with the Gaussian model indicating that this covariate is not significant at all. For all the models the longitude $s_1$ is significant. For the most complete model we conclude that there is a increase in the maximum temperature when we move north to south and also when we move east to west.

In the most complete model, the variance of $\tilde Z(s_i,t_j)$ is
\[
\Var[\tilde Z(s_i,t_j)]\equiv \sigma_{ij}^2=\frac{\sigma^2}{\lambda_{2j}}\left[\frac{1-\tau^2}{\lambda_{i1}}+\frac{\tau^2}{h_i}\right],
\]
from which we can deduce the variance structures over stations and time points. In particular, if we marginalise over space by assigning the spatial mixing variables an ``average'' constant value which we can take to be the prior mean, {\it i.e.}~$\lambda_{i1}=h_i=1$ then we can trace the temporal evolution of the variance as $\sigma^2/\lambda_{2j}$. Similarly, if we assume $\lambda_{2j}=1$, we can decompose the spatial variance into a part for the correlated process $\sigma^2(1-\tau^2)/\lambda_{1i}$ and a nugget part $\sigma^2\tau^2/h_{i}$. Figure \ref{fgVarMax12h} gives the posterior boxplots for the variance evolution over time and space for the full model, whereas   Figure \ref{fgVarMax1} does the same for the model with only $h$ and $\lambda_2$ (here the constant contribution to the variance of the spatially smooth process is not depicted).
Notice the clear separation of the effects showing that stations $17$, $18$ and $66$ (both panels (b)) are outliers whereas some regions in time have larger observational variance, especially  (panels (a)). The $3$ stations identified as outliers are all located in the north but quite far from each other, indicating that it is not the case of a region in space with larger variance. The temporal variance pattern in both models is in line with the empirical variances over time shown in Figure \ref{fgTempMax1} (d). The component containing $\lambda_1$ (Figure \ref{fgVarMax12h} (c)) does indicate appreciable differences in variances over spatial regions, in line with the Bayes factor in favour of this model and the posterior inference on $\nu$. The Gaussian model is, of course, not as flexible, assuming just one global variance with posterior 95\% credible interval given by $(8.64,17.03)$.
Figure \ref{fgVarMax2} presents the total posterior variance $\sigma_{ij}^2$ for all time points at stations $1$ (small nugget effect) and $17$ (large nugget effect) using the complete model, while horizontal lines indicate the posterior credible interval for the variance for the Gaussian model. This clearly illustrates the inadequacy of the Gaussian model in adapting to the temporal changes in the variance. For instance, for station 1 and time points $5$ to $10$ it is clear that the Gaussian model overestimates the variance and for time point $31$ it is clear that the variance should be larger than the one estimated by the Gaussian model. This underestimation is even more pronounced for location 17 where the nugget effect is more important, and even time period 22 seems problematic then.

In view of the lack of smoothness of the process in time, we do not present any predictions in time, but we will conduct interpolation to unobserved sites. In our case, we left some stations out of the estimation sample in order to compare the predictions with the actually observed maximum temperature. Table \ref{tbTempMax1} shows %the average width of the 95\% predictive interval and
the estimated value of the LPS as in (\ref{lps}) using the Gaussian and non-Gaussian models for the $3$ stations left out of the estimation step. The non-Gaussian models that include $\lambda_2$ give very good predictions when we compare with the Gaussian model. Specially in the tails, the extreme events are predicted in a much more adequate way. This can be seen in Figure \ref{fgMaxPred2} which plots the predictions (medians and $95\%$ credible intervals) obtained for various models against actual observations at out-of-sample stations for all $31$ time periods. The non-Gaussian models provide interval widths that are similar to the Gaussian on average but the intervals are smaller for many points in time and larger in the tails (e.g. $j=31$) where the posterior prediction interval for the Gaussian model misses the observed value for the process at stations $2^*$ and $3^*$. Station $1^*$ is very close to other stations and even the Gaussian model gives reasonable predictions. This suggests that scale mixing in time is essential in this application in order to produce good interpolations in space, especially for extreme temperatures.


In order to verify how the Gaussian and non-Gaussian models would predict in other regions of the spatial domain we repeated the estimation and prediction steps for 5 different partitions  of the data, selecting at random three testing locations and leaving the remaining 67 for estimation. The results obtained for the Gaussian model and the 2 best non-Gaussian models are presented in Table \ref{tabRepPred}. Notice that the complete model has the smallest log predictive score, but is closely followed by the model with $\lambda_2$ and $h$. %When we compare the $95\%$ credible interval for the predictions the model with only $\lambda_2$ and $h$ is the best one with smaller intervals even when we compare with the Gaussian model.

%\begin{figure}[htb]
%\center
%\begin{tabular}{cc}
%\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/TemperatureNew/SpainMap.pdf} &
%\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/TemperatureNew/CoordPred.pdf} \\
% (a) Spain and France Map. & (b) Basque Country (Zoom).  \\
%\end{tabular}
%\caption{Position of the $67$ stations used for estimation and the $3$ out-of-sample %stations used only for prediction ($1^*$, $2^*$ and $3^*$).}
%\label{fgTempMax0}
%\end{figure}

\begin{figure}[htb]
\center
\begin{tabular}{cc}
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/TemperatureNew/boxplotSpace.pdf} &
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/TemperatureNew/boxplotTime.pdf} \\
 (a) Temperature over space & (b) Temperature over time \\
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/TemperatureNew/empVarSpace.pdf} &
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/TemperatureNew/empVarTime.pdf} \\
 (c) Empirical variance for each location & (d) Empirical variance for each time point  \\
\end{tabular}
\caption{Data summaries.}
\label{fgTempMax1}
\end{figure}

\begin{table}[htb]
\center
\begin{tabular}{c|ccccccc}
            & $h$ & $\lambda_1$ & $\lambda_1$ $\&$ $h$ & $\lambda_2$ & $\lambda_2$ $\&$ $h$ & $\lambda_1$ $\&$ $\lambda_2$ & $\lambda_1$, $h$ $\&$ $\lambda_2$\\
\hline
Shifted Gamma estimator & 172 & 148 & 345 & 138 & 279 & 417 & 547 \\
Raftery estimator $p_4$ & 115 & 116 & 227 & 116 & 203 & 302 & 327 \\
\end{tabular}
\caption{Log Bayes factor in favour of the model in the column versus the Gaussian model using Shifted-Gamma ($\lambda=0.98$) and Newton-Raftery ($d=0.01$) estimators for the marginal likelihood.} \label{ComparaTemp}
\end{table}


\begin{table}[htb]
\center
\begin{tabular}{c|ccccc}
 &  $a$ & $b$ & $c$ & $\alpha$ & $\beta$\\
\hline
Gaussian & (0.99, 2.04) &(3.50, 8.13) & (0.00, 0.15) & (1.31, 1.83) & (0.33, 0.50)\\
$h$ & (0.99, 1.98) & (3.29, 8.22)  & (0.01, 0.13) & (1.39, 1.74) & (0.31, 0.48)\\
$\lambda_1$  & (4.17, 24.08) & (3.49, 7.94)  & (0.01, 0.13) & (0.58, 0.94) & (0.38, 0.55)\\
$\lambda_1$ $\&$ $h$  & (1.44, 3.15) &  (3.30, 7.70) & (0.01, 0.12) & (1.20, 1.55) &  (0.36, 0.53)\\
$\lambda_2$  & (0.95, 2.01) & (4.28, 9.74)  & (0.01, 0.10) & (1.35, 1.84) & (0.32, 0.46) \\
$\lambda_2$ $\&$ $h$ & (1.10, 2.33) & (3.33, 7.88)  & (0.01, 0.13) & (1.29, 1.68) & (0.30, 0.45)\\
$\lambda_1$ $\&$ $\lambda_2$  & (7.28, 21.20) & (0.78, 4.41)  & (0.01, 0.15) & (0.75, 1.00) & (0.29, 0.52)\\
$\lambda_1$, $h$ $\&$ $\lambda_2$ & (6.36, 18.14) & (0.39, 2.58) & (0.01, 0.20) & (0.89, 1.11) & (0.31, 0.68)\\
\end{tabular}
\caption{$95\%$ posterior credible intervals for some parameters of interest.} \label{StatMaxDelta1}
\end{table}

\begin{figure}[htb]
\center
\begin{tabular}{ccc}
\includegraphics[height=5.3cm, angle=360]{FiguresNonGaussian/TemperatureNew/TempPriorPostnu.pdf} & \includegraphics[height=5.3cm, angle=360]{FiguresNonGaussian/TemperatureNew/TempPriorPostnuh.pdf} & \includegraphics[height=5.3cm, angle=360]{FiguresNonGaussian/TemperatureNew/TempPriorPostnu2.pdf}\\
 (a) $\nu$ & (b) $\nu_h$ & (c) $\nu_2$ \\
\end{tabular}
\caption{Prior and posterior densities of the parameters of the mixing distributions for the model with $h$, $\lambda_1$ and $\lambda_2$. Prior densities are given by dashed lines and posterior densities are indicated by drawn lines.}
\label{fgnus}
\end{figure}

\begin{figure}[htb]
\center
\begin{tabular}{cc}
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/TemperatureNew/Var1_ng2h.pdf} & \includegraphics[height=6cm, angle=360]{FiguresNonGaussian/TemperatureNew/Var2_ng2h.pdf} \\
 (a) $\sigma^2/\lambda_{2j},j=1,\ldots,31$ & (b) $\sigma^2\tau^2/h_{i},i=1,\ldots,67$ \\
\end{tabular}
\caption{Posterior boxplots of the variance structure over time and space corresponding to the model  with $h$ and $\lambda_2$, where  $\sigma_{ij}^2=\sigma^2[(1-\tau^2)+\tau^2/h_i]/\lambda_{2j}$. The left panel describes the temporal evolution and the right panel the spatial nugget effect.}
\label{fgVarMax1}
\end{figure}

\begin{figure}[htb]
\center
\begin{tabular}{ccc}
\includegraphics[height=5.3cm, angle=360]{FiguresNonGaussian/TemperatureNew/Var1_ng12h.pdf} & \includegraphics[height=5.3cm, angle=360]{FiguresNonGaussian/TemperatureNew/Var2_ng12h.pdf} & \includegraphics[height=5.3cm, angle=360]{FiguresNonGaussian/TemperatureNew/Var3_ng12h.pdf}\\
 (a) $\sigma^2/\lambda_{2j},j=1,\ldots,31$ & (b) $\sigma^2\tau^2/h_{i},i=1,\ldots,67$ & (c) $\sigma^2(1-\tau^2)/\lambda_{1i},i=1,\ldots,67$ \\
\end{tabular}
\caption{Posterior boxplots of the variance structure over time and space for the model with $h$, $\lambda_1$ and $\lambda_2$, where $\sigma_{ij}^2=\sigma^2[(1-\tau^2)/\lambda_{1i}+\tau^2/h_i]/\lambda_{2j}$. The left panel describes the temporal evolution, the middle panel the smooth spatial process and the right panel the spatial nugget effect.}
\label{fgVarMax12h}
\end{figure}


\begin{figure}[htb]
\center
\begin{tabular}{cc}
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/TemperatureNew/Var3_ng2h.pdf} & \includegraphics[height=6cm, angle=360]{FiguresNonGaussian/TemperatureNew/Var4_ng2h.pdf} \\
 (a) station 1 & (b) station 17 \\
\end{tabular}
\caption{Posterior boxplots of the total variance for two stations at each time point corresponding to the model  with $h$, $\lambda_1$ and $\lambda_2$, where  $\sigma_{ij}^2=\sigma^2[(1-\tau^2)/\lambda_{1i}+\tau^2/h_i]/\lambda_{2j}$. The left panel is for station 1 (small nugget effect) and the right panel for station 17 (large nugget effect). Horizontal lines indicate the 95\% credible interval for the variance of the Gaussian model.}
\label{fgVarMax2}
\end{figure}

\begin{figure}[htb]
\center
\begin{tabular}{ccc}
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/TemperatureNew/Pred1g.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/TemperatureNew/Pred2g.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/TemperatureNew/Pred3g.pdf}\\
(a) Gaussian Model & (b) Gaussian Model & (c) Gaussian Model\\
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/TemperatureNew/pred1ng2h.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/TemperatureNew/pred2ng2h.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/TemperatureNew/pred3ng2h.pdf}\\
(d) Model with $\lambda_2$ $\&$ $h$ & (e) Model with $\lambda_2$ $\&$ $h$ & (f) Model with $\lambda_2$ $\&$ $h$\\
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/TemperatureNew/pred1ng12h.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/TemperatureNew/pred2ng12h.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/TemperatureNew/pred3ng12h.pdf}\\
(g) Model with $h$, $\lambda_1$ $\&$ $\lambda_2$ & (h) Model with $h$, $\lambda_1$ $\&$ $\lambda_2$ & (i) Model with $h$, $\lambda_1$ $\&$ $\lambda_2$\\
\end{tabular}
\caption{Posterior predictive median of the maximum temperature versus the observed maximum temperatures (points) at out-of-sample stations for times $j=1,\ldots,31$ . The dashed line is the $95\%$ credible predictive interval.}
\label{fgMaxPred2}
\end{figure}

%\begin{table}[htb]
%\center
%\begin{tabular}{cccc} \hline
%model & Average width & $\bar{IS}$ & LPS\\ \hline
%Gaussian & 3.78 & 4.35 & 103.81 \\
%$h$ & 3.83 & 4.34 & 102.04 \\
%$\lambda_1$  & 3.74 & 4.36 & 105.09 \\
%$\lambda_1$ $\&$ $h$ & 3.75 & 4.48 & 103.79 \\
%$\lambda_2$ & 3.73 & 3.94 & 87.33 \\
%$\lambda_2$ $\&$ $h$ & 3.73 & 3.87 & 86.57 \\
%$\lambda_1$ $\&$ $\lambda_2$ & 4.51 & 4.65 & 85.89 \\
%$\lambda_1$, $h$ $\&$ $\lambda_2$ & 3.84 & 4.02 & 83.78 \\ \hline
%\end{tabular}
\begin{table}[htb]
\center
\begin{tabular}{cc} \hline
model  & LPS\\ \hline
Gaussian & 103.81 \\
$h$ & 102.04 \\
$\lambda_1$   & 105.09 \\
$\lambda_1$ $\&$ $h$  & 103.79 \\
$\lambda_2$  & 87.33 \\
$\lambda_2$ $\&$ $h$  & 86.57 \\
$\lambda_1$ $\&$ $\lambda_2$  & 85.89 \\
$\lambda_1$, $h$ $\&$ $\lambda_2$ &  83.78 \\ \hline
\end{tabular}
\caption{Log predictive score (LPS) for the predicted maximum temperature at the out-of-sample stations.} \label{tbTempMax1}
\end{table}

%\begin{table}[htb]
%\center
%\begin{tabular}{cccc} \hline
%model & Average width & $\bar{IS}$ & $\bar{LPS}$\\ \hline
%Gaussian & 3.88 & 4.25 & 92.26 \\
%$\lambda_2$ $\&$ $h$ & 3.82 & 4.03 & 81.76 \\
%$\lambda_1$, $h$ $\&$ $\lambda_2$ & 4.20 & 4.29 & 81.43 \\ \hline
%\end{tabular}
\begin{table}[htb]
\center
\begin{tabular}{cc} \hline
model  & $\overline{\rm LPS}$\\ \hline
Gaussian & 92.26 \\
$\lambda_2$ $\&$ $h$  & 81.76 \\
$\lambda_1$, $h$ $\&$ $\lambda_2$  & 81.43 \\ \hline
\end{tabular}
\caption{Average log predictive score ($\overline{\rm LPS}$) for the predicted maximum temperature at the out-of-sample stations for 5 partitions of the data set.} \label{tabRepPred}
\end{table}


\clearpage
\newpage

\section{Conclusions and future work}

We present a non-Gaussian spatiotemporal model that is able to capture departures from Gaussianity in terms of outlier contamination and regions in space or time with larger observational variance. The proposed finite dimensional distributions have heavier tails than the normal distribution and have the normal as a limiting case. The general model includes correlated mixing in the spatiotemporal process (both in time and space) and in the nugget effect. This model is quite flexible, combining nonseparability and non-Gaussian behaviour,  and performed well on simulated data. This was also illustrated in an application to Spanish temperature data, where simultaneous mixing in the nugget effect, in space and in time seems an essential feature. Prediction is straightforward, using the fact that we are using scale mixtures Gaussian processes, and efficient MCMC algorithms for posterior and predictive inference also immediately allow for the identification of outliers and regions in time and space with inflated variances.

As a subject of future research, it might be interesting to explore the effect of using $C^*_1(s)=M_1(-\gamma_1(s))$ and $C^*_2(t)=M_2(-\gamma_2(t))$ in the model for posterior and predictive inference, and to investigate ways of separately modelling correlation structures for the mixing variables and the observables. In the simulated examples we examined here, it did not seem to be an important restriction since the correlation is satisfactorily recovered and this is what matters for prediction.

\appendix

\appsection{Proof of Proposition \ref{Prop1} \label{ApA}}

The Gaussian-Log-Gaussian process $\{\tilde Z(s,t);(s,t)\in D\times T \}$ as defined in (\ref{TildeProcess}) and (\ref{covSpace}) with no nugget effect can be rewritten as
$$\tilde Z(s,t)=\frac{1}{\sqrt{\lambda_1(s)}}Z_1(s;U)Z_2(t;V)=\frac{1}{\sqrt{\lambda_1(s)}}Z(s,t).$$
Then for a fixed time point $t_0$ we have
$$\tilde C(s,0)=\sigma^2C(s,0)f(s),$$
where $f(s)=\exp\{\frac{\nu}{4}[M_1(-\gamma_1(s))-1]+\nu\}$ and $C(s,0)=\sigma^2M_0(-\gamma_1(s))M_1(-\gamma_1(s))$. Therefore, $\tilde C^{(2m)}(s,0)=\sum_{i=0}^{2m}{2m \choose i} C^{(i)}(s,0)f^{(2m-i)}(s)$.

By Fa\'a di Brunos's formula, termwise differentiation of $C(s,0)$ results in
\begin{eqnarray}\label{CovPurSpatDer}
C^{(2m)}(s,0)&=&\sum_A \frac{m!}{k_1!k_2!...k_{2m}!}y^{(k)}(-\gamma_1(s))\prod_{k_j\not = 0}\left\{\frac{-\gamma_1^{(j)}(s)}{j!}\right\}^{k_j}\\ \nonumber
& = & \{y^{(1)}(-\gamma_1(s))[\gamma_1^{(2m)}(s)]+\ldots +y^{(2m)}(-\gamma_1(s))[\gamma_1^{(1)}(s)]^{2m}\}
\end{eqnarray}
where $A=\{k_1,k_2,...,k_{2m}:k_1+2k_2+...+2mk_{2m}=2m\}$, $k=k_1+k_2+...+k_{2m},\;k_i\geq0,\;i=1,2,\ldots,2m$ and $y^{(k)}(x)=\sum_{i=0}^{k}{k \choose i}M_0^{(k-i)}(x)M_1^{(i)}(x).$

In the expression (\ref{CovPurSpatDer}), the highest order derivative of $\gamma_1(s)$ is $2m$ obtained when $k_{2m}=1$ and $k_1=\ldots=k_{2m-1}=0$. Thus, in this case, the behavior of $C^{(2m)}(s,0)$ as $s\to 0$ depends only on the local behavior of $\gamma_1^{(2m)}(s)$ and $M_1^{(1)}(r)$ as $s,r\to 0$. And the highest order derivative of $M_1(-\gamma_1(s))$ is $2m$ obtained when $i=2m$. Thus, the behavior of $C^{(2m)}(s,0)$ as $s\to 0$ depends only on the local behavior of $M_1^{(2m)}(r)$ and $\gamma_1^{(1)}(s)$ as $s,r\to 0$. Notice that $M_0^{(k)}(0)=\E[X_0^k]$ always exists since $X_0\sim {\rm Gamma}(\lambda_0,1)$ in our model.

This proves that $C^{(2m)}(s,0)$ exists and is finite at $0$ if and only if $M_1^{(2m)}(r)$ and $\gamma_1^{(2m)}(s)$ are finite and exist at $0$.
Also, $f^{(2m)}(s)=\frac{\gamma}{4}f(s)\left\{ M_1^{(2m)}(-\gamma_1(s))+g(\nu,s) \right \}$ exists and is finite under the same conditions \citep[see][]{PSteel06}. The result follows.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{ref}
\bibliographystyle{jasa}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\end{document}

\appsection{Large space-time data \label{ApB}}

Notice that in order to calculate the likelihood functions (\ref{likelihoodSpace}) or (\ref{likelihoodTime}) we need to invert $\Sigma$ with dimensions $IJ\times IJ$. Several approaches were proposed to overcome this computational limitation. \citet{Stein04} proposes the use of conditional distributions where we consider a partition of $Z$ into subvectors $Z_1,... ,Z_b$ of possibly different lengths, then
\begin{equation}\label{cond}
p(z|\phi )=p(z_1|\phi )\prod_{j=1}^{b}p(z_j|z_{(j-1)},\phi ),
\end{equation}
where $Z_{(j)}=(Z_1,... ,Z_j)$. In particular, if $S_{(j)}$ is some subvector of $Z_{(j)}$ then
$$p(z|\phi)\approx p(z_1|\phi)\prod_{j=1}^{b}p(z_j|s_{(j-1)},\phi).$$
To use this approximation, we need to order the observations in some manner, and it is suggested to use the rank of the projections of the observations along some axis. Another issue is how to choose the conditioning vectors. It is pointed out in \citet{Stein04} how these choices influence the parameter estimation (using restricted likelihood). \citet{Jones97} use this approach for space-time processes and they suggest defining the nearest neighbors based on a preliminary estimate of the space-time correlation. If there are not too many points in the spatial domain, the definition of these subvectors is straightforward. We can define $z_k=(z(s_1,t_k),\ldots,z(s_I,t_k))'$, $k=1,\ldots,J$ and the subvector $s_{(k)}=(z_{k-L},\ldots,z_{k-1})'$ considers only $L$ points in time.

Another proposal is Tapering (\citet{Nyck06}) where the basic idea is to taper the covariance function to zero beyond a certain range using a positive definite but compactly supported function $K_{\theta}$. The effect of tapering is analyzed using the infill asymptotic theory.
\begin{equation}\label{tap}
K_{tap}(x,y)=K(x,y)K_{\theta}(x,y)
\end{equation}
The Schur product $K_{tap}$ preserves some of the shape of $K$ and it is identically zero outside a fixed range. A possible objection to tapering is that it may not be effective for a spatial covariance with long range correlations. \citet{Fuentes07} also gives an approximation for the likelihood of Gaussian processes but her approach is restricted to the spectral domain. For the simulated data examples provided in the previous sections we didn't use any of these approaches but for the real data application we used the approach based on conditional distributions proposed in \citet{Stein04}.







\subsection{Brazilian air temperature data}

To illustrate the proposed models we analyze the Brazilian air temperature data available in the website {http://satelite.cptec.inpe.br/PCD/}. The data consist of time series of air temperature in degree Celsius obtained every $3$ hours at 37 meteorological stations in Ceara (Brasil) during the period 25/sep/2008 to 05/oct/2008. We use UTM coordinates (divided by $10^5$). Notice that the scale of the spatial coordinates times $100$ is in kilometers. Figure \ref{stationsBr} shows the location of the 37 stations which are scattered irregularly in the state of Ceara in the northeast of Brazil. The median of the distance between stations is 220 kilometers, the minimum distance is 22 kilometers and the maximum distance is 560 kilometers. Figure \ref{fgTempBr1} (a) presents the boxplot of temperature for each station. It looks like the temperature varies significantly across stations, thus we consider the spatial coordinates and altitude ($X_{alt}$) as covariates in the mean function. We also considered the relative humidity in $\%$ ($X_{humid}$) in the mean function. Notice in Figure \ref{fgTempBr1} (b) that there is a marked daily cycle in the time series for the air temperature. This will be considered in the mean function by introducing some sines and cosines. The mean function considered is given by
\begin{eqnarray}\label{meanTempBr}
\mu(s,t) & = & \delta_0+\delta_1 t+\delta_2 s_1+\delta_3 s_2+\delta_4 X_{alt}+\delta_5sin\left (2\pi\frac{10t}{J}\right )+\delta_6cos\left (2\pi\frac{10t}{J}\right )\\ \nonumber
& + & \delta_7sin\left (2\pi\frac{20t}{J}\right )+\delta_8cos\left (2\pi\frac{20t}{J}\right )+\delta_9 X_{humid},
\end{eqnarray}
where $t=1,\ldots,J$, $J=80$.

We calculated the empirical variance for the residual of a classical regression model over space and time as a exploratory data analysis. The empirical variances over time (figure \ref{fgTempBr1} (d)) seems to vary more than over space (figure \ref{fgTempBr1} (c)) with some very high peaks as the one at $t=55$ corresponding to 02/Oct at 9:00 am. Also there are some clusters of large variance in time such as the one from time 28 to 34 corresponding to 29/Sep, from 0:00 to 18:00. There are some large variances across space and some of them are close together (figure \ref{fgTempBr1} (e)), but it is not so clear that they are clusters in space since the distances in space are quite large.

The covariance model considered here is nonseparable as presented in (\ref{CovNonSep}) allowing for interactions in space and time. The model is parameterized as presented in Subsection \ref{secTheta} and the chosen covariance in space is the Cauchy type in time and space-time and exponential in space. The parameters estimated are the trend coefficients $(\delta_0,\ldots,\delta_{9})$, the covariance parameters $(\eta_0,\beta,a,b,\sigma^2,\tau^2)$ and the mixing parameters $(\nu_1,\nu_2,\nu_h)$. We also generate the auxiliary variables $(\lambda_1,\lambda_2,h)$ in our MCMC algorithm as described in Subsection \ref{SecInf} in order to identify regions in space and time with larger observational variance. The parameters $\eta_1$ and $\eta_2$ are set to $1$. See \citet{Fons08} for a discussion of parameterisation in this class of covariance models.

We estimated the parameters for several models by approximating the likelihood as discussed in the application to the Spanish temperature data and here we also used the lag in time equal to $L=5$. The models considered are: the Gaussian model, the non-Gaussian model with $\lambda_1$ only, the non-Gaussian model with lognormal $h$, the non-Gaussian model with lognormal $h$ and $\lambda_1$, the non-Gaussian model with $\lambda_2$ only, the non-Gaussian model with $\lambda_1$ and $\lambda_2$ and finally the non-Gaussian model with lognormal $h$, $\lambda_1$ and $\lambda_2$. Notice that all the models considered here have a nugget effect in space parameterized by $\tau^2>0$ that accounts for measurement errors. The comparison criteria using Bayes Factors presented in  Table \ref{ComparaWindBr} indicates the non-Gaussian models as more adequate for this application when compared to the Gaussian model. The model with only $h$ is better than the model that only includes $\lambda_1$. This suggests that the different variances are mostly due to outliers in space instead of regions with larger variability. But the inclusion of only $\lambda_2$ gives an improvement even more significant for this application. The inclusion of only $\lambda_2$ is even better than including $h$ and $\lambda_2$ suggesting that in this application the heterogeneity is mainly present in time rather than in space. The other models including also $\lambda_1$ and $h$ give a certain improvement but not as large as when $\lambda_2$ is present in the model.

Table \ref{StatTempBrDelta1} presents some summaries for the marginal posterior distributions for some parameters of interest. The posterior distributions for the smoothness parameters obtained by the non-Gaussian models that do not include $\lambda_1$ suggest that the process is much smoother both in space. This might be because these models are able to explain part of the variability in space in a more adequate way. Moreover, the non-Gaussian models present very similar posterior distributions for $\beta$ suggesting smoother processes also in time, except for the models that include $h$ and $\lambda_1$ at the same time. These two models result in very different covariance structures with much larger values of the spatial range parameter. It seems like the inclusion of $h$ and $\lambda_1$ is not a good idea in this example. For all the models it looks like the covariance structure is significantly nonseparable in space and time as the posterior distribution of $c$ does not give much mass for values close to 0 as we can see in Figure \ref{postcTempBr}. The posterior distribution for $c$ in the Gaussian model seems to overestimate the value of $c$.

Table \ref{StatTempBrDelta2} present some posterior summaries for the coefficients in the mean function (\ref{meanTempBr}). The effect of the altitude is very significant and similar for all the models indicating, as expected, that in higher altitudes we find lower air temperatures. It is suggested by all the models that there is a increase in the air temperature when we move north to south and also when we move east to west, that is, as we get more distant from the sea. The humidity was also very significant in all models indicating that the less the humidity the higher the air temperature.

We can separate the variances for each station and each time point in two components that in the most complete model will have component 1 given by $\sigma^2(1-\tau^2)/(\lambda_{1i}\lambda_{2j})$ and component 2 given by $\sigma^2\tau^2/h_{i}$. The second component would give the dimension of the nugget effect in the model. The Gaussian model is not as flexible given just one global variance that has posterior credible interval given by $(2.62,3.41)$. The posterior distribution for components $1$ and $2$ for the model with $\lambda_2$ and $h$ are shown in figure \ref{fgVarTempBr1}.



Notice the clear separation of the effects showing that stations $17$, $18$ and $66$ (Figure \ref{fgVarMax1} (b)) are outliers whereas some regions in time have larger observational variance (Figure \ref{fgVarMax1} (a)). The $3$ stations pointed as outliers are located in the north region but quite far from each other, indicating that it is not the case of a region in space with larger variance. Figure \ref{fgVarMax1} (c-d) present the posterior variance for all time points at stations $1$ (small nugget effect) and $17$ (large nugget effect) and lines in the posterior credible interval for the variance for the Gaussian model. For instance, for station 1 and time points $5$ to $10$ it is clear that the variance is less than the one estimated by the Gaussian model and for time point $31$ it is clear that the variance is larger than the one estimated by the Gaussian model. This will be very important in the interpolation step.

\begin{figure}[htb]
\center
\begin{tabular}{cc}
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/map1.pdf} &
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/CoordReal.pdf} \\
 (a) Brazilian Map. & (b) Ceara (Zoom).  \\
\end{tabular}
\caption{Position of the $33$ stations used for estimation.}
\label{stationsBr}
\end{figure}

\begin{figure}[htb]
\center
\begin{tabular}{cc}
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/WindBr/boxplotSpace.pdf} &
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/WindBr/series.pdf}\\
 (a) Temperature over space. & (b) Temperature over time.  \\
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/WindBr/empVarSpace.pdf} &
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/WindBr/empVarTime.pdf} \\
 (c) Empirical variance for each location. & (d) Empirical variance for each time point.  \\
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/WindBr/empVarSpaceMap.pdf} &
\\
 (e) Empirical variance for each location. &   \\
\end{tabular}
\caption{Data summaries.}
\label{fgTempBr1}
\end{figure}


\begin{table}[htb]
\center
\begin{tabular}{c|ccccccc}
            & $h$ & $\lambda_1$ & $\lambda_1$ $\&$ $h$ & $\lambda_2$ & $\lambda_2$ $\&$ $h$ & $\lambda_1$ $\&$ $\lambda_2$ & $\lambda_1$, $h$ $\&$ $\lambda_2$\\
\hline
Shifted Gamma estimator & 590 & 546 & 703 & 938 & 986 & 1049 & 1268\\
Raftery estimator $p_4$ & 562 & 527 & 634 & - & - & - & - \\
\end{tabular}
\caption{The natural logarithm of the Bayes factor in favor of the model in the column versus Gaussian model using Shifted-Gamma ($\lambda=0.98$) and Newton-Raftery ($d=0.01$) estimators for the predictive density of $z$.} \label{ComparaWindBr}
\end{table}


\begin{table}[htb]
\center
\begin{tabular}{c|ccccc}
 &  $a$ & $b$ & $c$ & $\alpha$ & $\beta$\\
\hline
Gaussian & (0.05,1.26) &(0.78,3.64) & (0.06,0.48) & (0.37,0.77) & (0.46,0.72)\\
$h$ & (1.17,2.61) & (1.23,2.77)  & (0.02,0.30) & (1.02,1.79) & (0.66,0.92)\\
$\lambda_1$  & (0.12,0.55) & (1.35,3.46)  & (0.04,0.38) & (0.44,0.78) & (0.63,0.90)\\
$\lambda_1$ $\&$ $h$  & (3.06,20.97) & (1.21,2.72)  & (0.02,0.22) & (0.50,0.87) & (0.52,0.80)\\
$\lambda_2$  & (0.60,1.83) & (1.78,3.86)  & (0.02,0.27) & (0.72,1.62) & (0.58,0.79)\\
$\lambda_2$ $\&$ $h$ & (1.09,2.44) & (1.34,2.99)  & (0.02,0.30) & (0.97,1.72) & (0.64, 0.90)\\
$\lambda_1$ $\&$ $\lambda_2$  & (0.25,0.64) & (1.39,3.08)  & (0.04,0.29) & (0.57,0.90) & (0.64,0.90)\\
$\lambda_1$, $h$ $\&$ $\lambda_2$ & (26.59,81.83) & (0.74,0.95)  & (0.01,0.18) & ( 0.55,0.78) & (1.53,1.96)\\
\end{tabular}
\caption{$95\%$ posterior intervals for some parameters of interest.} \label{StatTempBrDelta1}
\end{table}

\begin{table}[htb]
\center
\begin{tabular}{c|ccccc}
 & $\delta_1$ (t) & $\delta_2$ (longitude) &  $\delta_3$ (latitude) & $\delta_4$ (altitude) & $\delta_9$ (humidit)\\
\hline
Gaussian & (0.00,0.03) &  (-0.31,-0.03) & (-0.48,-0.07) & (-0.0065,-0.0052) & (-0.13,-0.11)\\
$h$ & (0.00,0.03) & (-0.33,-0.01) & (-0.37,0.03) & (-0.0060,-0.0049) & (-0.14,-0.13) \\
$\lambda_1$  & (0.01,0.03) & (-0.40, -0.10) & (-0.31,0.06) & (-0.0061,-0.0050) & (-0.14,-0.13) \\
$\lambda_1$ $\&$ $h$ & (0.00,0.03) & (-0.53,-0.13) & (-0.42,0.00) & (-0.0063,-0.0051) & (-0.13,-0.12) \\
$\lambda_2$ & (0.00,0.02) & (-0.14,0.10) & (-0.37,-0.07) & (-0.0061,-0.0053) & (-0.15,-0.14) \\
$\lambda_2$ $\&$ $h$ & (0.00,0.02) & (-0.23,0.00) & (-0.33,-0.05) & (-0.0059,-0.0051) & (-0.14,-0.13) \\
$\lambda_1$ $\&$ $\lambda_2$ & (0.01,0.02) & (-0.39,-0.14) & (-0.33,-0.05) & (-0.0062,-0.0054) & (-0.14,-0.13) \\
$\lambda_1$, $h$ $\&$ $\lambda_2$ & (0.01,0.03) & (-0.53,-0.19) & (-0.24,0.22) & (-0.0062,-0.0052) & (-0.14,-0.13) \\
\end{tabular}
\caption{$95\%$ posterior intervals for some parameters of interest.} \label{StatTempBrDelta2}
\end{table}


\begin{figure}[htb]
\center
\begin{tabular}{c}
\includegraphics[height=8cm, angle=360]{FiguresNonGaussian/WindBr/postC.pdf}\\
\end{tabular}
\caption{Marginal posterior distribution for the separability parameter.}
\label{postcTempBr}
\end{figure}


\begin{figure}[htb]
\center
\begin{tabular}{cc}
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/WindBr/Var1_ng2h.pdf} & \includegraphics[height=6cm, angle=360]{FiguresNonGaussian/WindBr/Var2_ng2h.pdf} \\
 (a) Component 1: $\sigma^2(1-\tau^2)/\lambda_2$. & (b) Component 2: $\sigma^2\tau^2/h$. \\
\includegraphics[height=6cm, angle=360]{FiguresNonGaussian/WindBr/Var3_ng2h.pdf} & \includegraphics[height=6cm, angle=360]{FiguresNonGaussian/WindBr/Var4_ng2h.pdf} \\
 (c) $\sigma_{ij}^2$, $i=1$ (large nugget effect). & (d) $\sigma_{ij}^2$, $i=5$ (small nugget effect). \\
\end{tabular}
\caption{Posterior of the two components of the variance for each station and each time point given by $\sigma_{ij}^2=[\sigma^2(1-\tau^2)/\lambda_{2j}]+[\sigma^2\tau^2/h_i]$. Model with $h$ and $\lambda_2$.}
\label{fgVarTempBr1}
\end{figure}





\begin{figure}[htb]
\center
\begin{tabular}{cc}
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/var_ng.pdf} & \includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/Med_ng.pdf} \\
 (a) $\sigma^2_i=\sigma^2\left \{\frac{1-\tau^2}{\lambda_{1i}}+\tau^2\right \}$. & (b) Posterior median of the variance.\\
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/series_5.pdf} & \includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/series_29.pdf} \\
 (c) Station 5 (small variance). & (d) Station 29 (large variance). \\
\end{tabular}
\caption{Posterior of the variance for each location $\sigma^2_i$ for the non-Gaussian model with $\lambda_1$ only.}
\label{fgWindng}
\end{figure}

\begin{figure}[htb]
\center
\begin{tabular}{cc}
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/var_ng2.pdf} &
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/series_t3.pdf}\\
 (a) $\sigma^2_j=\sigma^2\left \{\frac{1-\tau^2}{\lambda_{2j}}+\tau^2\right \}$. &
 (b) $t_3$ (09:00, 25/Sep/08) (small variance).\\
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/series_t71.pdf} & \includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/series_t35.pdf} \\
 (c) $t_{71}$ (09:00, 03/Oct/08) (medium variance). & (d) $t_{35}$ (09:00, 29/Sep/08) (large variance). \\
\end{tabular}
\caption{Posterior of the variance for each time point $\sigma^2_j$ for the non-Gaussian model with $\lambda_2$ only.}
\label{fgWindng2}
\end{figure}


\begin{figure}[htb]
\center
\begin{tabular}{cc}
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/varst1ng12h.pdf} &
\includegraphics[height=7cm, angle=360]{FiguresNonGaussian/WindBr/varst3ng12h.pdf}\\
 (a) $\sigma^2_{ij},i=1$. & (a) $\sigma^2_{ij},i=3$.
\end{tabular}
\caption{Posterior of the variance for each stations and each time point $\sigma^2_{ij}$ for the non-Gaussian model with $\lambda_1$, $\lambda_2$ and $h$.}
\label{fgWindng3}
\end{figure}


\begin{figure}[htb]
\center
\begin{tabular}{ccc}
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred1gvalue.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred2gvalue.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred3gvalue.pdf}\\
(a) Gaussian Model. & (b) Gaussian Model. & (c) Gaussian Model.\\
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred1ng12value.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred2ng12value.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred3ng12value.pdf}\\
(b) Model with $\lambda_1$ and $\lambda_2$. & (c) Model with $\lambda_1$ and $\lambda_2$. & (d) Model with $\lambda_1$ and $\lambda_2$.\\
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred1ng2hvalue.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred2ng2hvalue.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred3ng2hvalue.pdf}\\
(e) Model with $\lambda_2$ and $h$. & (f) Model with $\lambda_2$ and $h$. & (g) Model with $\lambda_2$ and $h$.\\
\end{tabular}
\caption{Predictions at out-of-sample stations for $t=1,\ldots,31$ and observed values (points).}
\label{fgTempBrPred1}
\end{figure}

\begin{figure}[htb]
\center
\begin{tabular}{ccc}
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred1g.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred2g.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred3g.pdf}\\
(a) Gaussian Model. & (b) Gaussian Model. & (c) Gaussian Model.\\
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred1ng12.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred2ng12.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred3ng12.pdf}\\
(b) Model with $\lambda_1$ and $\lambda_2$. & (c) Model with $\lambda_1$ and $\lambda_2$. & (d) Model with $\lambda_1$ and $\lambda_2$.\\
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred1ng2h.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred2ng2h.pdf} &
\includegraphics[height=5.3cm,angle=360]{FiguresNonGaussian/WindBr/pred3ng2h.pdf}\\
(e) Model with $\lambda_2$ and $h$. & (f) Model with $\lambda_2$ and $h$. & (g) Model with $\lambda_2$ and $h$.\\
\end{tabular}
\caption{Posterior median for the maximum temperature versus the observed maximum temperatures (points) at out-of-sample stations for $t=1,\ldots,31$ . The dashed line is the $95\%$ credible interval for the predictions versus the observed maximum temperature.}
\label{fgTempBrPred2}
\end{figure}

\begin{table}[htb]
\center
\begin{tabular}{ccccc} \hline
model & Average width & Median width & Average interval score & LPS*\\ \hline
Gaussian & 6.00 & 6.01 & 6.40 & 338 \\
$h$ & 6.53 & 6.53 & 6.90 & 321\\
$\lambda_1$  & 6.46 & 6.47 & 6.75 & 317 \\
$\lambda_1$ $\&$ $h$ & 6.26 & 6.21 & 6.77 & 320 \\
$\lambda_2$ & 5.48 & 4.42 & 5.81 & 302 \\
$\lambda_2$ $\&$ $h$ & 5.06 & 4.05 & 5.42 & 288 \\
$\lambda_1$ $\&$ $\lambda_2$ & 5.52 & 4.61 & 5.80 & 282 \\
$\lambda_1$, $h$ $\&$ $\lambda_2$ & 5.11 & 4.20 & 5.44 & 294\\ \hline
\end{tabular}
\caption{Comparison of $95\%$ interval forecast and Log predictive score (LPS*) for the predicted air temperature at out-of-sample stations.} \label{ComparaPredTempBr2}
\end{table}


\begin{table}[htb]
\center
\begin{tabular}{cccc} \hline
model & Average width & Average interval score & LPS*\\ \hline
Gaussian & 5.66 & 5.88 & 59.91 \\
$h$ &  &  &  \\
$\lambda_1$  & 5.34 & 5.99 & 59.21 \\
$\lambda_1$ $\&$ $h$ &  &  &  \\
$\lambda_2$ & 5.24 & 6.88 & 61.41 \\
$\lambda_2$ $\&$ $h$ &  &  &  \\
$\lambda_1$ $\&$ $\lambda_2$ & 4.93 & 5.66 & 56.99 \\
$\lambda_1$, $h$ $\&$ $\lambda_2$ & 9.58 &  9.61 &  58.02\\ \hline
\end{tabular}
\caption{Comparison of $95\%$ interval forecast and Log predictive score (LPS*) for the predicted air temperature at time $J+1$.} \label{ComparaPredTempBr}
\end{table}

\clearpage
\newpage
