\section{Introduction}\label{intro}

Regional Climate Models (RCM) describe the medium scale global  atmospheric  and  oceanic  dynamics and serve as downscaling models. RCMs use atmospheric interactions in General Circulation Models (GCM) to develop a higher  resolution  output,  which  is  influenced  by  a smaller scale of topographical features. They represent a tool for assessing climate change impacts,  to develop  seasonal  climate  predictions,  and  are a  powerful  tool  in general for  improving  our  understanding  of regional climate processes \citep{wang04}.  

The  North  American  Regional  Climate  Change  Assessment Program (NARCCAP) is an international program that has as a main goal generating climate scenarios for use in impacts research \cite{Mearns2009}. They have proposed and applied a dynamical downscaling technique that consists of  embedding  RCMs  within  GCMs  to  obtain  higher  resolution evaluations of the climate model over the domain of interest, using different combinations of RCMs and GCMs. As \cite{Wilby1997} state: ``RCMs are computationally demanding and require orders of magnitude more computer time than statistical downscaling to compute equivalent scenarios''. Since this dynamical downscaling technique is computationally  expensive,  a downscaling emulator can serve as a fast statistical approximation to  perform sensitivity analysis  or  for  making  informed  prior judgments relating the GCM drivers for an RCM, is of high value. 

Previous work on statistical emulators from \cite{Castruccio2014}, \cite{overstall2016multivariate}, \cite{Mearns1999}, \cite{OHagan2006} and \cite{Hernanz2022} are concentrated on climate models projections and on computer experiment output. In this work, we describe how to use a spatio-temporal statistical model with varying coefficients (VC), as a downscaling emulator for a RCM. Other experiments to emulate regional models using the NARCCAP output and downscaling techniques have been performed, as described in \cite{Wood2004} and \cite{Laflamme2016}. Nevertheless, to the best of our knowledge there are no attempts to emulate RCM using varying coefficient models. 

Most of the models that use varying coefficients are based on the seminal work by \cite{Gelfand2003a}. A myriad of applications for this model have been developed in epidemiology, ecology, geography, environmental sciences and biology, and some in climate models, such as \cite{liying2019}. These models allow for marginal effects to be non-stationary over space, time or both, and thus offer a higher degree of flexibility than a simple fixed coefficient spatio-temporal model. Due to the computational burden that these models can have with big data, a proposal that uses varying coefficients in space, time or both, needs to be computational efficient. 

Comparisons between methods for analyzing large spatial data, such as the one presented in \cite{Heaton2018}, have compared the speed and accuracy of spatial models but have not included VC methods. In this work, three methods for estimating spatio-temporal downscaling emulators that use VC models are being compared: i) A Multiresolution (MRA) approach such as the one originally described in \cite{katzfuss_multi-resolution_2017}, and expanded in \cite{Katzfuss2020}; INLA, as presented in \citep{Rue2009}; and the approach from \cite{Dambon2021}. The first one is implemented in R for this work, and the other two use R packages that are already published: \textit{INLA} and \textit{varycoef}, respectively.

A Multiresolution (MRA) approach \cite{katzfuss_multi-resolution_2017} can help with the computational burden, and it has been shown to work with massive spatially distributed data \cite{Katzfuss2016} in the context of regular Gaussian processes. To the best of our knowledge, MRA methods have not been applied to estimate VC models with covariates. INLA has limited work on VC models, but \cite{Blangiardo2013} have shown how to work with linear models with random effects. Finally, \cite{Dambon2021} method is the only R package that has been exclusively proposed for spatially VC models to use on medium-size data sets. The objective of this paper is to compare the performance of INLA, varycoef, and MRA for building a statistical downscaling emulator for RCM, and then show that the emulator works properly for NARCCAP data.

The outline of this article is as follows: Section 2 introduces the data from NARCCAP and the applied problem. Section 3 describes the spatio-temporal downscaling emulator, and give an overview on the estimation details for each method. Section 4 describes two simulation studies: the first one investigates the performance of each of the methods when estimating the statistical emulator, and the second one explores the fitting of two different models to emulate data with similar characteristics as the real data application. Section 5 presents the NARCCAP data, the motivation for the problem, and the results from a real data application of the model. Finally, Section 6 presents the discussion and conclusions. 

\section{Data Application}

Our main source of data is NARCCAP \cite{https://doi.org/10.5065/D6RN35ST}. Our data set consists of temperature (ts) measured in Kelvin degrees from the Canadian Regional Climate Model (CRCM) \cite{caya1995description}, seasonal mean pressure vertical velocity measured in Pa/s (OMEGA) and temperature (TREFHT) measured in Kelvin degrees from the Community Climate System Model (CCSM), which is the global model embedded in the CRCM \footnote{The complete data sets can be accessed through: \url{https://www.earthsystemgrid.org/project/NARCCAP.html}}. The two climate models have different resolutions and projections, as presented in Figure \ref{fig1:data}, for the study area of interest.

\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.8]{resolution2.png}
	\caption{Coordinates for the Canadian Regional Climate Model (CRCM) (red points) and the Community Climate System Model (CCSM) in the North American Region (black points).}
	\label{fig1:data}
\end{figure}

The area of study was defined using the intersection between NARCCAP data and the Monsoon area as defined by Higgins et al (2006), given that it is an important climate region that drives weather changes in Central America and the Caribbean, and could be used as an input for weather forecast in the region. 

There are 2482 points in for the regional-model resolution, and 270 points for the global-based one. Monthly data from 1990 to 1999 (inclusive) was used to construct and test the model. The main objective of the emulator is to use the CCSM output to infer the values of the CRCM output. For that propose, it is important to recognize the non-stationary nature of the relationship between the CCSM output and the regional model response, as well as giving the model enough flexibility to describe a non-stationary mean in space and time, together with an acceptable downscaling performance. 

A descriptive analysis to gather evidence about this non-stationary relationship was done with the data, and it was found that the mean difference between both outputs in log-scale was not constant over space (see \ref{fig2:coefficients}), and that defining the mean as constant in space and time could limit the ability of the model to predict the response accurately. A spatio-temporal downscaling emulator that can incorporate spatio-temporal VC is described in the next section.

\begin{figure}[htp!]
	\centering
	\includegraphics[scale=0.8]{promedios.png}
	\caption{Mean difference between the monthly temperature log-outputs among the regional and global models. The difference was computed over the period 1990-1999 inclusive.}
	\label{fig2:coefficients}
\end{figure}

\section{Statistical Methods}\label{sec2}

\subsection{Spatio-temporal Downscaling Emulator}

Consider a climate-related variable $C_t(s)$ indexed by time $t$ and spatial location $s\in \mathcal S$. The spatial set $\mathcal S$ is assumed to be regularly-spaced with a fixed resolution. Assume there exists a set of $q$ covariates linearly related with the dependent variable $C_t(s)$ in the following way:
\begin{align}\label{Global_eq}
	C_t(s)=\alpha+\beta'X_t(s)+\epsilon_t(s),
\end{align} 
where $\alpha$ and $\beta$ are random parameters and $\epsilon_t(s)$ is a white noise in both space and time. The main idea is to model the variable $C_t(\cdot)$ over a set of finer resolution on the same spatial domain $\mathcal S$. Let $\mathcal W$ be such set, and assume that the resolution satisfies that each location $s\in \mathcal S$ can be placed as a center of a rectangle containing several (or none) locations $w_1,\ldots,w_K$ on the finer set $\mathcal W$. Figure \ref{diagrama} shows how the location of both sets can be visualized when $K=4$.
\begin{figure}[htp!]
	\centering
	\includegraphics*[scale=0.35]{diagrama}
	\caption{Illustration of points in $\mathcal{S}$ (blue) and $\mathcal{W}$ (black).}
	\label{diagrama}
\end{figure}

Consider an arbitrary point $w\in \mathcal{W}$ which is located within the rectangle whose center is $s$. We assume that \eqref{Global_eq} holds on this scale, but we adjust the local relationship among the response and predictor as follows:
\begin{align}\label{Regional_eq}
	C_t(w)=[\alpha+\alpha_t^r(w)]+[\beta+\beta_t^r(w)]'X_t(s)+\epsilon_t(s)+\gamma_t(w),
\end{align}
where 
\begin{itemize}
	\item $\alpha_{\cdot}^r(\cdot)\sim N(\beta_0,\Sigma_0(\theta_0))$,
	\item $\beta_{\cdot}^r(\cdot)\sim N(\beta_1,\Sigma_1(\theta_1))$ and
	\item $\gamma_t(\cdot)\stackrel{i.i.d}{\sim} N(0,\tau^2)$
\end{itemize}
and $\Sigma_0$ and $\Sigma_1$ are separable spatio-temporal covariance matrices parameterized with the vectors $\theta_0$ and $\theta_1$ respectively and the parameters $\beta_0$ and $\beta_1$ take into account any resulting bias due to the resolution change. The random terms $\alpha_t^t(w)$, $\beta_t^r(w)$ and $\gamma_t(w)$ account for the increase in uncertainty due to the downscaling process. Based on the above, we are interested in the following model:
\begin{align}
	Y_t(w):=C_t(w)-C_t(s)=\alpha_t^r(w)+\beta_t^r(w)'X_t(s)+\gamma_t(w),
	\label{modelo}
\end{align}
where $C_t(s)$ is treated as a known covariate in this model and hence we are interested in the modeling of the difference $Y_t(w)$ in order to get predictions on the finer variable $C_t(w)$. Note that $\alpha_{\cdot}^r(\cdot)$ and $\beta_{\cdot}^r(\cdot)$ represent random spatio-temporal correction terms on $\alpha$ and $\beta$ in \eqref{Global_eq} due to the use of the predictor $X_t(s)$ over the coarser set $\mathcal S$ and \eqref{Regional_eq} the change in scale between the sets $\mathcal S$ and $\mathcal W$.  

Under a Bayesian approach, we are interested in the conditional distribution of the dependent vector $\mathbf Y:=(Y_t(w))_{t,w}$ given the parameters $\Phi:=(\beta_0,\beta_1,\theta_0,\theta_1, \tau^2)$:
\begin{align*}
	\mathbf Y|\Phi \sim N(\beta_0+\mathbf X^T\beta_1,\Sigma_Y),
\end{align*}
where 
\begin{align*}
	\Sigma_Y=\mathbf X^T \left(\Sigma_1(\theta_1)+\Sigma_2(\theta_2)\right) \mathbf X+\tau^2 I_n
\end{align*}
and $\mathbf X$ is a block-diagonal matrix whose $i$-th sub-matrix is $X^T(s_i)$, $i=1,\ldots,n$ and $n$ is the number of locations in $\mathcal S$. In order to evaluate the likelihood of $\mathbf Y$ we need to evaluate the inverse and determinant of $\Sigma_Y$ which can be time-consuming for a combination of a large $n$ and/or a large number of times $t$. Several methods have been proposed to bypass calculating the inverse and determinant of $\Sigma_Y$ for large data using approximations. In the next section, these methods are explained in light of our model. 


\subsection{Approximation Methods}

Depending on both the temporal and spatial dimensions, the inverse of the resulting covariance matrix in its Gaussian likelihood can be time-consuming to compute, specially when any (or both) dimensions are large. The likelihood calculation is required to estimate the parameters of \eqref{modelo}, either by a frequentist or Bayesian approach. For the purpose of this work, we are interested in the Bayesian inference of the model and we aim to remedy the computational issues by means of the following approximations: MRA, INLA and varycoef. 

\subsubsection{Multi-Resolution Approximation}
The model in \eqref{modelo} is a spatially varying coefficient model in the sense of \cite{Gelfand2003a} when the time $t$ is considered fixed. Under this setting, a first type of approximation that can be applied is by using a predictive process, as described in \cite{banerjee_gaussian_2008}. In general terms the predicted process is obtained by the ``kriged'' projection of any of the random spatial effects over a set of knots on $\mathcal S$. 

This approximation has issues capturing local/small scale dependence properly, leading to biased estimations and errors in predictions \citep{Sang2011}. That can be improved by considering an approximation of the difference between the random effect and the predictive process (called residual process) $\epsilon_1:=\omega(s)-\omega_1(s)$ with a new random field whose covariance matrix is given by:
\begin{align*}
	\Sigma_{\epsilon_1}=(\Sigma_\omega-\Sigma_{\omega_1})\circ K,
\end{align*}
where $\circ$ stands for the Shur product of matrices (entry-wise product), K is a modulating matrix and $\Sigma_\omega,\Sigma_{\omega_1}$ are the covariance matrices for $\omega(s)$ and $\omega_1(s)$ respectively (see \citealp{Sang2011} for more details) . A natural way to choose the modulating matrix is to consider a partition in $K$ disjoint sub-regions of $\mathcal{S}$. Then the matrix $K$ is defined as $K(s,s')=1$ if $s,s'$ belong to the same sub-region and 0 otherwise. 
The inverse and determinant of $\Sigma_Y$ can be easily computed using the Sherman-Woodbury-Morrison formula (see \citealp{Harville1997})
The above method is called full-scale approximation (FSA) (see \citealp{Sang2011})

The above procedure can be repeated over $K$ disjoint regions of $\mathcal{S}$ and the number of knots on each region is chosen to guarantee that the resulting predictive process has smaller dimension.
If we compute the predictive process over $\epsilon_1$ and we decompose any of the disjoint sets in $K$ disjoint subsets, we can apply FSA using a certain number of knots over the resulting subset. This process can be repeated $M$ times to get a multi-resolution approximation (MRA) (see \citealp{katzfuss_multi-resolution_2017} and \citealp{Katzfuss2020}) of the covariance matrix of $Y$ up to level $M$:
\begin{align*}
	\Sigma_Y=\sum_{i=1}^M\mathbf X^T \Sigma_{w_i} \mathbf X+\mathbf X^T \Sigma_{\epsilon_M} \mathbf X+\tau ^2 I_n,
\end{align*} 
and we are able to compute the inverse and determinant of $\Sigma_Y$ by applying the Sherman-Woodbury-Morrison formula recursively. To this point, there is a difference with respect to the work in \cite{katzfuss_multi-resolution_2017} since the model in \eqref{modelo} is a linear model with spatially varying slope and intercept, whereas  \cite{katzfuss_multi-resolution_2017} uses the MRA only for a linear model with only spatially varying intercept. 

\cite{Appel2020} extends the MRA approach to spatio-temporal data under different covariance structures, using an appropriate two-dimensional partitioning. Up to now, this article shows the first attempt of MRA for spatio-temporal varying linear models as \eqref{modelo}. 


\subsubsection{Dambon Approach}

As we mentioned before, the model in \eqref{modelo} is a spatially-varying model when the time $t$ is considered as fixed. A second approach that we use in this case is the likelihood approximation according to \cite{Dambon2021}. This approach has two important assumptions, first it assumes mutual prior independence of the random effects $\alpha^r(\cdot)$ and $\beta^r(\cdot)$ and furthermore it assumes that each random effect has tapering on its covariance matrix. The first assumption allows the covariance matrix of $Y$ in model \eqref{modelo} to be written as:

\begin{align*}
	\Sigma_Y=\sum_{i=1}^p \Sigma^{j}  \odot \mathds{X}^{j}+\tau^2 I_n,
\end{align*} 
where $\odot$ is the Hadamard product and $\mathds{X}^{j}:= x^{(j)}\left(x^{(j)}\right)^T$ uses the $j$th covariate $x^{(j)}$. In order to ease the computation burden, \cite{Dambon2021} incorporate the tapered covariance matrix, $\Sigma_{tap}^{j}:=\Sigma^{j} \odot C_{\rho^*}$, where $\rho^*$ is the tapering range. In this case, it is a sparse matrix with $[\Sigma_{tap}^{j}]_{kl}=0$ for $||s_k-s_l||\geq \rho^*$. As a result, the covariance matrix can be simplified as: 

\begin{align*}
	\Sigma_{Y,tap}=\sum_{i=1}^p \Sigma_{tap}^{j}  \odot \mathds{X}^{j}+\tau^2 I_n.
\end{align*} 

The tapering, together with the sparsity of the matrix $\textbf{X}$ let the execution time to be reduced. This approach is fully implemented with maximum likelihood estimation on the R package varycoef, available on CRAN. We adapted the approach by using MCMC estimation. 

\subsubsection{Integrated Nested Laplace Approximation}

The computational challenge of MCMC inference in model \eqref{modelo} is due to the frequent evaluation of the likelihood during the Metropolis-Hastings steps. A possible strategy is to rely on Integrated Nested Laplace Approximation (INLA) to accelerate the MCMC procedure. The INLA approach is applicable to a general specification for which the mean $\eta_i$ of the observations $y_i$ follows a linear structure:
\begin{align}\label{eq:meanINLA}
	\eta_i = \alpha +\sum_{m=1}^M\beta_mx_{mi}+\sum_{l=1}^Lf_l(z_{li}),
\end{align}
where $\alpha$ represents an intercept, the coefficients
$\mathbf{\beta} = (\beta_1,\ldots,\beta_M)$ relate $M$ covariates
$(x_1,\ldots,x_M)$ to $\eta_i$, and $f = \{f_1(\cdot),\ldots,f_L(\cdot)\}$ is a collection of
random effects defined on a set of $L$ covariates $(z_1,\ldots,z_L)$ (see
\citealp{Rue2009} and \citealp{Blangiardo2013}). 

In model \eqref{modelo}, if we consider $y_i=\eta_i=Y_t(w)$ and the varying intercept and slope are the latent variables with fixed effects as their respective means $\beta_0$ and $\beta_1$, then our models fall into the general specification of INLA. The main objective of our Bayesian estimation is to compute the
marginal posterior distribution of each parameter in $\Phi$. To attain computational advantages, INLA assumes that the prior of
vector $\Phi$ is a multivariate normal random vector with a precision matrix that depends on a certain set of hyperparameters. INLA further 
approximates the conditional distribution of $\Phi$ by using a simplified Laplace Approximation together with optimization procedures (see \citealp{Rue2009}). 

\section{Simulation Study}\label{sec3}

Two independent simulations were performed. In order to compare the three approximation methods, a first simulation study was done with the objective of evaluating the goodness of fit of the three approximation methods under study, for a spatial downscaling emulator and for three different spatial resolutions. Once the best approximation method was defined, a second simulation study uses that best approximation method to compare, under a predictive perspective, two different models for the spatio-temporal emulator, using simulated data that resembles the real data application. 

\subsection{Simulation I: Comparing approximation methods}

We study the fitting performance of three approximation methods for a spatially VC model: MRA, varycoef and INLA. For MRA, the method was implemented using R and adapted for a SV model with MCMC estimation. For the Dambon Approach, we used the package \textit{varycoef} with tapering of 0.1, available on CRAN, and adapted its likelihood evaluations to use it with a MCMC. Finally, we used the INLA package available on CRAN.

Our main focus is to describe how well we can estimate a simulated response according to \eqref{Regional_eq}, and we compare goodness-of-fit performance for different structures of spatial correlation when estimating the response variable and the parameters involved, according to elapsed time and the final model estimation in \eqref{modelo}. For that, we use the following statistics: mean elapsed time to estimate the model \eqref{modelo} (in minutes), Interval Score (IS) \citep{Gneiting2007a} and MSE for each estimated model; where the last two metrics are defined as follows:

\begin{align*}
	MSE = \frac 1 n \sum_{t,w}(Y_t(w)-\hat Y_t(w))^2,
\end{align*}

\begin{align*}
	IS_{\alpha}=\frac 1 n \sum_{t,w}\left[(U_t(w)-L_t(w))+\frac{2}{1-\alpha}(L_t(w)-Y_t(w))\cdot 1_{Y_t(w)<L_t(w)}\right. \\
	  \left.+\frac{2}{1-\alpha}(Y_t(w)-U_t(w))\cdot 1_{Y_t>U_t}\right],
\end{align*}
where $n$ is the number of observations, $\hat Y_t(w)$ is the estimated response, $L_t(w)$ is the $(1-\alpha)$\%-predictive lower bound and $U_t(w)$ is the $(1-\alpha)$\%-predictive upper bound at location $w$ and time $t$ over the period of study. 

In the spatial region $\left[0,1\right]^2$, we set three different resolutions for the regularly-spaced set $\mathcal S$ and the finer set $\mathcal W$. If $(g \times g)$ and $(r \times r)$ denote the resolution size for $\mathcal S$ and $\mathcal W$ respectively, we will denote ($g \times r$) as a way to simplify the notation for a single scenario. To be precise, we will use $25 \times 10$, $40 \times 20$ and $55 \times 25$ as three different resolution scenarios for this simulation and for each scenario we generate 10 replications.

At these locations, we simulate two models according to \eqref{Global_eq} and \eqref{Regional_eq} with $q=1$, and $\alpha=1$, $\beta_0=2$, $\epsilon_t(s)\stackrel{i.i.d}{\sim} N(0,\sigma^2)$ with $\sigma^2=2$ and $\tau^2=1$. In this simulation, we consider a spatially varying intercept with a spatial covariance structure following either a Matern$(\phi=3,\sigma=1,\nu=1.5)$ or an Exponential$(\phi=3,\sigma=1)$. 

\subsubsection{Simulation I Results}

Data was analyzed using a PC prior in all cases according to \cite{Franco2019}. Table~\ref{tab:simulation_prediction} presents the goodness of fit metrics by model and resolution, including the elapsed time in minutes for each of the approximation methods.
\begin{table}
		\caption{\label{tab:simulation_prediction}Goodness of fit metrics by approximation method, covariance structure of the data, and resolution. Elapsed time has minutes as units.}
		\centering
		{\scriptsize
		\begin{tabular}{ll|rrp{1cm}|rrp{1cm}|rrp{1cm}}
			\toprule
			\multicolumn{2}{c}{} & \multicolumn{9}{c}{\textbf{Resolution}} \\
			\cmidrule{3-11}
			\multicolumn{2}{c}{ } & \multicolumn{3}{|c|}{$25\times 10$} & \multicolumn{3}{c|}{$40\times 20$} & \multicolumn{3}{c}{$55\times 25$} \\
			\midrule
			Covariance & Model & MSE & $IS_{.95}$ & Elapsed time & MSE & $IS_{.95}$ & Elapsed time & MSE & $IS_{.95}$ & Elapsed time \\
			\midrule
			\multirow{3}{*}{Exponential}  & MRA      & 20.827	  & 11.158  & 37.53   & 22.041&11.242  &496.20  & - & -& -  \\
			& varycoef &12.588   &10.276   &5.80    &20.674  &10.450  &63.44   &28.984&10.593 &450.32 \\
			& INLA     &5e-08	 &0.0837   &8.70 	   &2e-08  &0.0852   &10.44    &1e-08 &0.0847 &17.10  \\
			\midrule
			\multirow{3}{*}{Matern} & MRA      & 14.205 & 11.181  &54.32   &22.011  & 11.242 & 481.89 & - & -&- \\
			& varycoef &12.660	   &10.224   & 5.94  & 20.620	 &10.414  &62.16  &28.963&10.586&430.62  \\ 
			& INLA     &6e-08	&0.0820   & 8.63   &1e-08	  &0.0853 &11.11  &1e-08	&0.0844& 15.81 \\
			\bottomrule
		\end{tabular}
		}
	\end{table}

It is clear from these results that MRA is not an efficient approximation choice, given that it takes more than 6 hours to run a model with a relatively small number of points ($40 \times 20$), and that we were not able to obtain a competitive simulation, in terms of elapsed time, with more points than that. Now, comparing only INLA and varycoef, INLA is a very superior option in both time and accuracy, as shown in Table~\ref{tab:simulation_prediction}.

The parameter estimation did not present much of a difference for the approximation methods, in terms of bias or variability. Given these results, INLA is going to be the approximation of choice in the rest of the paper.

\subsection{Simulation II: Comparing spatio-temporal models}

Once the approximation method is defined, the objective of the second simulation is to compare two alternatives for a spatio-temporal SV model given the same set of spatio-temporal data generated in Simulation I. In this exercise, we consider a spatial-temporal varying intercept model with a spatial covariance structure following a Matern$(\phi=5,\sigma=0.0003,\nu=1)$, and two different temporal structures: random noise and AR(1) processes. For simplicity, the spatio-temporal covariance structure is assumed to be separable. Consequently, the models can be written as follows:

\begin{description}
	\item[Model 1:] Spatio-temporal varying intercept following a random noise process, where the spatial covariance structure follows a Matern.
	\begin{align*}
		\begin{cases}
			Y_t(w)&=\alpha_t(w)+\gamma_t(w)\\
			\alpha_t(w)&= \epsilon_t(w)
		\end{cases}
	\end{align*} 
	\item[Model 2:] Spatio-temporal varying intercept following an $AR(1)$ process with $\rho=0.8$, where the spatial covariance structure follows a Matern.
	\begin{align*}
		\begin{cases}
			Y_t(w)&=\alpha_t(w)+\gamma_t(w)\\
			\alpha_t(w)&= \rho \alpha_{t-1}(w)+\epsilon_t(w)
		\end{cases}
	\end{align*}	
\end{description}

Data was generated in a spatial region $\mathcal S=\left[0,20\right]^2$, we set three different resolutions for the regularly-spaced $(g \times g)$ global and $(r \times r)$ regional locations, and $T$ periods of time for each location. Table \ref{tab:simulation_setting} shows the three spatio-temporal resolution settings. For each simulation setting, we generate 30 replications, and for each simulated data, we divided it into two sets: a training data set containing the first $5/6$ periods of time, and a testing data set, containing the last $1/6$ periods of time.

At these locations, we simulate two models according to \eqref{Global_eq} and \eqref{Regional_eq} with $q=0$, and $\alpha=5.707$, $\beta_0=-0.001$, $\epsilon_t(s)\stackrel{i.i.d}{\sim} N(0,\sigma^2)$ with $\sigma^2=0.0003$ and $\tau^2=1/700000$. We choose these parameter values based on the calculated values for the climate variable of interest in Section \ref{sec4}. Furthermore, we transform ${C_t(\cdot)}$ using an exponential function, to ensure a symmetric distribution for the response variable. 

\subsubsection{Simulation II Results}

\begin{table}
	\caption{\label{tab:simulation_setting}Grid setting and elapsed time (in minutes) by model and global $\times$ regional $\times$ time resolution}
	\centering
	\begin{tabular}[t]{lccc}
		\toprule
		\multicolumn{1}{c}{ } & \multicolumn{3}{c}{\textbf{Resolution}} \\
		\cline{2-4}
		& 1 & 2 & 3\\
		\midrule
		Grid setting $(g\times r\times T)$ & $(25\times 10\times12)$ & $(40\times20\times12)$ & $(25\times10\times72)$  \\
		Model 1 & 0.51 & 0.81 & 1.84\\
		Model 2 & 2.79 & 3.74 & 20.74\\
		\bottomrule
	\end{tabular}
\end{table}

Table \ref{tab:simulation_setting} presents the elapsed time for different resolutions and models. The execution time for the Model 1 increases gently from Resolution 1 to 3, but the Model 2, with an extra parameter $\rho$ of the autoregressive structure, it increases significantly between resolutions 2 and 3.



Table \ref{tab:simulation_predict} shows the predictive metrics for the training and testing data set for the two models by each resolution scenario. For each resolution and model, similar metrics of the testing and training data set are observed. This indicates that the prediction are reliable and we are not overfitting the models.


\begin{table}
	\caption{\label{tab:simulation_predict}Prediction metrics by model and resolution}
	\centering
	\begin{tabular}[t]{>{}ll|rr|rr|rr}
		\toprule
		\multicolumn{2}{c}{ } & \multicolumn{6}{c}{\textbf{Resolution}} \\
		\cline{3-8}
		\multicolumn{2}{c}{ } & \multicolumn{2}{|c|}{1} & \multicolumn{2}{c|}{2} & \multicolumn{2}{c}{3} \\
		\cline{3-4} \cline{5-6} \cline{7-8}
		Model & Data & MSE & $IS_{.95}$ & MSE & $IS_{.95}$ & MSE & $IS_{.95}$\\
		\midrule
		& test & 27.58 & 139.39 & 27.35 & 139.89 & 27.21 & 139.62\\
		\multirow[t]{-2}{*}{\raggedright\arraybackslash \textbf{1}} & 
		train & 26.95 & 138.09 & 27.01 & 139.51 & 26.99 & 139.30\\
		\cline{1-8}
		& test & 27.74 & 148.16 & 27.44 & 147.42 & 27.57 & 148.04\\
		\multirow[t]{-2}{*}{\raggedright\arraybackslash \textbf{2}} & train & 27.16 & 146.20 & 27.24 & 147.00 & 27.01 & 146.28\\
		\bottomrule
	\end{tabular}
\end{table}

Finally, we compute the relative MSE (rMSE) and relative $IS_{\alpha}$ (r$IS_{\alpha}$) for each posterior estimate by dividing its corresponding metric by its true absolute parameter, for ease of comparison. In this way we can compare the ability of each simulated model to identify the true parameters on each resolution setting. Table \ref{tab:simulation_coef} presents the results for each model and resolution. We can observe that for both models, $\beta_0$ and $\phi$ can be estimated very well. $\rho$ estimates can also be reliably estimated for the Model 2. However, the parameters related to the variability: the standard deviation of the random coefficient $\sigma$ and the random noise $\tau$ are not estimated accurately. Although these estimates presents high bias and variability, we argue that from a predictive point of view, these models produce satisfactory results.



\begin{table}
		\caption{\label{tab:simulation_coef}Parameter estimation metrics (relative MSE and $IS_\alpha$) by model and resolution.}
	\centering
	\begin{tabular}[t]{cc|rr|rr|rr}
		\toprule
		\multicolumn{2}{c|}{ } & \multicolumn{6}{c}{\textbf{Resolution}} \\
		\cline{3-8}
		\multicolumn{2}{c}{ } & \multicolumn{2}{c|}{1} & \multicolumn{2}{c|}{2} & \multicolumn{2}{c}{3} \\
		\cline{3-4} \cline{5-6} \cline{7-8}
		Model & Parameter & rMSE & rIS & rMSE & rIS & rMSE & rIS\\
		\midrule
		& $\beta_0$ & 0.0001 & 5.9386 & 0.0000 & 5.8912 & 0.0000 & 5.8728\\
		& $\phi$ & 2.6287 & 12.7381 & 1.8426 & 7.8701 & 0.4986 & 13.5038\\
		& $\sigma$ & 0.0003 & 39.9919 & 0.0003 & 39.9938 & 0.0003 & 39.9951\\
		\multirow[t]{-4}{*}{\raggedright\arraybackslash \textbf{1}} &  $\tau$ & 0.0625 & 8101.7592 & 0.0628 & 8244.6291 & 0.0625 & 8252.1995\\
		\midrule
		& $\beta_0$ & 0.0003 & 4.0196 & 0.0001 & 3.9623 & 0.0001 & 3.9407\\
		& $\phi$ & 0.6798 & 14.3594 & 0.6296 & 13.3838 & 2.0720 & 11.4603\\
		& $\rho$ & 0.0211 & 3.3914 & 0.0131 & 2.7524 & 0.0083 & 0.9989\\
		& $\sigma$ & 0.0003 & 40.0000 & 0.0003 & 40.0000 & 0.0003 & 40.0003\\
		\multirow[t]{-5}{*}{\raggedright\arraybackslash \textbf{2}} & $\tau$ & 0.0629 & 8102.5600 & 0.0633 & 8233.7478 & 0.0621 & 8220.7104\\
		\bottomrule
	\end{tabular}
\end{table}

\section{Analysis of NARCCAP Data}\label{sec4}

A formal application of the spatio-temporal downscaling emulator was performed using the NARCCAP data set described in Section 2.  We chose the regional resolution as the main spatial domain of this analysis. The response variable is the temperature from the CRCM Regional Model, measured in Kelvin degrees, for this variable we choose to keep the model units and also to transform the response variable using a log scale, to attain a normal behavior approximately. Therefore, we put attention on a varying coefficient model according to \eqref{modelo} with:
\begin{itemize}
	\item $C_t(\cdot) = \log T_t(\cdot)$, where $T_t$ is the temperature and the domains $\mathcal S$ and $\mathcal W$ correspond to the global and regional domains respectively. 
	\item $X_t(\cdot) = OMEGA_t(\cdot)$ Vertical rising rate of air parcels. Provides a measure of large scale rising and sinking motions in the atmosphere. It was obtained from the CCSM model. 
\end{itemize}

The physical reasoning behind the model in \eqref{modelo} with the above variables can be justified by the fact that negative values of OMEGA indicate rising and a tendency for convection (and cooler temperatures), and positive values indicate sinking motions associated with warmer temperatures and high pressure (see \cite{Hostetler2011}). 

In order to compare the predictive performance of the model and choose an adequate emulator for our data, we use the following three alternatives:
\begin{description}
	\item[Model 0:] Non varying intercept model.
	\begin{align*}
		Y_t(w)=\alpha+\gamma_t(w)
	\end{align*}
	\item[Model 1:] Spatio-temporal varying intercept following a random noise process, where the spatial covariance structure follows a Matern$\left(\nu=1 \right)$.
	\begin{align*}
		\begin{cases}
			Y_t(w)&=\alpha_t(w)+\gamma_t(w)\\
			\alpha_t(w)&= \epsilon_t(w)
		\end{cases}
	\end{align*} 
	\item[Model 2:] Spatio-temporal varying intercept following an $AR(1)$ process, where the spatial covariance structure follows a Matern$(\nu= 1)$.
	\begin{align*}
		\begin{cases}
			Y_t(w)&=\alpha_t(w)+\gamma_t(w)\\
			\alpha_t(w)&= \rho \alpha_{t-1}(w)+\epsilon_t(w)
		\end{cases}
	\end{align*}	
	
	
	\item[Model 3:] Spatio-temporal varying intercept following an $AR(1)$ process with a fixed-parameter covariate:
	\begin{align*}
		\begin{cases}
			Y_t(w)&=\alpha_t(w)+\beta \cdot OMEGA_t(w)+\gamma_t(w)\\
			\alpha_t(w)&= \rho \alpha_{t-1}(w)+\epsilon_t(w)
		\end{cases}
	\end{align*}	 
\end{description}

The models are fitted using INLA package, in particular by the stochastic partial differential equation (SPDE) approach of \cite{Lindgren2011}. In order to use this method, we train the above models with historical data from 1990 to 1998 (inclusive) at a monthly basis and we use 1999 data as a testing set. The global-model spatial domain is used as knots for the mesh construction in the SPDE approach to reduce the dimensionality of the linear model due to the regional resolution of the RCM. The priors of the non-varying mean coefficients are assumed to be gaussian with the default parameterization of INLA. The prior of the range ($\phi$) and standard deviation ($\sigma$) in each of the Matern structures is assumed to be a Penalized Complexity (PC) prior (see \citealp{Fuglstad2019}) satisfying: $P[\rho < 700]=0.5$ and $P[\sigma > 0.32] = 0.01$ based on descriptive analysis on the temperature from RCM. A PC prior was also used for the correlation coefficient in the AR1 model, using as hyperparameters 0 and 0.9.

In order to compare the predictive performance of the three models, we apply the squared-root MSE and Interval Score between the observed RCM temperatures and the predicted temperatures over the one-year testing period. Table \ref{tabla_N_uno} contains a comparison of the two predictive metrics for the three models
\begin{table}
	\caption{\label{tabla_N_uno}Comparison of models according to the predictive metrics (1990-1998: training period, 1999: testing period) and elapsed time in minutes.}
	\centering
	\begin{tabular}{l|ll|ll|ll|ll}
		\toprule
		& \multicolumn{2}{c|}{\textbf{Model 0}} & \multicolumn{2}{c|}{\textbf{Model 1}} & \multicolumn{2}{c|}{\textbf{Model 2}} & \multicolumn{2}{c}{\textbf{Model 3}} \\ \midrule
		Year & MSE          & $IS_{.95}$           & MSE          & $IS_{.95}$           & MSE          & $IS_{.95}$           & MSE          & $IS_{.95}$           \\
		\midrule
		1990         & 12.52        & 68.86        & 1.41         & 12.67        & 1.42         & 13.21        & 1.42         & 13.19        \\
		1991         & 12.82        & 66.94        & 1.48         & 13.45        & 1.48         & 14.00        & 1.48         & 14.01        \\
		1992         & 10.24        & 59.12        & 1.43         & 12.97        & 1.43         & 13.51        & 1.43         & 13.50        \\
		1993         & 10.41        & 58.88        & 1.47         & 13.29        & 1.47         & 13.83        & 1.47         & 13.83        \\
		1994         & 10.63        & 59.99        & 1.51         & 13.69        & 1.51         & 14.26        & 1.51         & 14.24        \\
		1995         & 10.64        & 63.22        & 1.33         & 12.22        & 1.32         & 12.73        & 1.32         & 12.68        \\
		1996         & 9.72         & 57.17        & 1.58         & 13.92        & 1.59         & 14.47        & 1.59         & 14.48        \\
		1997         & 12.33        & 67.15        & 1.37         & 12.62        & 1.37         & 13.11        & 1.37         & 13.12        \\
		1998         & 9.20         & 52.94        & 1.41         & 12.87        & 1.41         & 13.39        & 1.41         & 13.38        \\
		\midrule
		1999         & 34.59        & 131.03       & 34.22        & 35.79        & 34.11        & 33.34        & 34.47        & 33.93        \\
		\midrule
		Elapsed Time & \multicolumn{2}{c|}{2.87}    & \multicolumn{2}{c|}{9.69}    & \multicolumn{2}{c|}{36.55}   & \multicolumn{2}{c}{36.57}   \\ \bottomrule
	\end{tabular}
\end{table}
and Table \ref{tabla_N_dos} contains the resulting parameters with their 95\%-predictive intervals for the best model according to the metrics (Model 2). Figure \ref{grafico_temps} shows the predicted values (mean and median), observed temperatures and interquantile ranges for three chosen months during the testing period. 	
\begin{figure}
	\caption{\label{grafico_temps}Mean temperatures, median temperatures and observed temperatures (left plot) and interquantile ranges (right plot) for three chosen months during the testing period (Model 2).}
	\centering
	\begin{tabular}{l}
		\includegraphics[scale=0.5]{graficoTemps.png} \\ \includegraphics[scale=0.5]{graficoRango_c.png}
	\end{tabular}
\end{figure}		 

Note that there is significant difference among the $IS$ obtained over the testing period for Model 0 and Model 1, which is an indicator that the spatial-temporal structure over the intercept is a valid assumption. If we make the temporal structure as AR(1), there is small gain in terms of predictive performance, but the addition of OMEGA as covariate is not improving the IS metric at all comparing with Model 2. In terms of computational time, the Model 2 offers an interesting alternative and the additional time obtained by including OMEGA as covariate is not substantial. The downscaling capacity of the model can be visualized in Figure \ref{grafico_temps}. It is interesting to note that the general patterns and trends of the temperature was predicted successfully in most of the months, but when the observed temperature raises during the summer months, the emulator was not able to capture extreme events during the year. Parameter estimates and prediction intervals are shown for Model 2 in \ref{tabla_N_dos}, where we were able to do quite precise inference in terms of standard errors for all the parameters of interest, but due to the identifiability issues found in simulation II, it would not be surprising to have large bias in the variance parameters.     


\begin{table}
	\caption{\label{tabla_N_dos}Parameter estimates and prediction interval for the best scenario (Model 2).}
	\centering
	\begin{tabular}{r|rrr}
		\toprule
		& \textbf{Estimate}  & \textbf{Upper 95\%} & \textbf{Lower 95\%}\\
		\midrule
		$\beta_0$ & -0.0023   & -0.0034    & -0.0011\\
		$\sigma$ & 0.3673    & 0.3396     & 0.3896 \\
		$\phi$   & 302.22    & 289.09     & 312.12 \\
		$\rho$   & 0.9424    & 0.9402     & 0.9451 \\
		$\tau$   & 0.0000195 & 0.0000193  & 0.0000197\\
		\bottomrule
	\end{tabular}
\end{table}

\section{Conclusions}\label{sec5}

A spatio-temporal downscaling emulator for regional models was proposed using spatial varying coefficients. Among the advantages of the said model, we are able to estimate non-stationary marginal effects, which means that the downscaling can vary over space. Furthermore, the model has flexibility to estimate the mean of any variable in space and time, as shown on the real data application, and finally, that the model has good prediction results, as shown on both the simulation results, and the analysis of NARCCAP data. This good results in terms of prediction happens in spite of the possible identification issues observed in the same sections. 

Throughout the simulations, INLA was by far the best approximation method for both the spatial and spatial temporal versions of the proposed model. Moreover, INLA was the fastest method for all the cases, and the approximation with best accuracy to estimate the different parameters from the model and the posterior distribution of the response variable. Computational time could be improved using other programming languages, but the comparison shows us what is possible to do only using R.

For future model extensions, it is important to highlight that when dealing with spatio-temporal models having a starting value for $\tau$ or $\sigma$ is not an easy task, and thus a lot of different values need to be tested in order to avoid identification issues with the model. 

The downscaling emulator showed potential for predicting accurately output from RCM using the output from a GCM. Specifically, temperature from the CRCM can be predicted using mean pressure vertical velocity output from CCSM with an out of sample MSE of approximately 34 units for all models proposed, and with a $95\%$ IS that is significantly lower for Models 2 and 3, compared to the other options, as shown in Table 5. Furthermore, when adding the elapsed time to run the models, and given that Model 1 and Model 2 give very similar predictive performance results, then Model 1 is preferred, since it can be run in less than one third of the time that takes to complete Model 2. 

For future work, there are several directions that can be taken. First, more combinations of RCMs and GCMs can be tested, using the data from NARCCAP. In addition, more variables can be used from the GCM output, in order to improve an emulator with the structure of Model 3. Lastly, another programming language such as Python or Julia can be used with the same methods, in order to improve the computational time. 


\section{Computational Resources and Code}

All the analysis were run using an HP Proliant DL380 Gen9 server, Intel(R) Xeon(R) CPU E5-2630 v3 @ 2.40GHz, with 16 logical cores, and 128GB RAM. The R scripts to read and wrangle the NARCCAP data, as well to run all of the models described in this paper are available at (will be updated for review) and were run using RStudio Server 2022.02.2 Build 485. Packages and R version are specified in the repository file called session.info.

\bibliographystyle{rss}
