Flexible shrinkage in high-dimensional Bayesian spatial autoregressive models

Flexible shrinkage in high-dimensional Bayesian spatial autoregressive models

\thispagestyle{empty}

Introduction

In the regional economic literature, spatial econometric model specifications have gained momentum in empirical research as a means to explicitly account for spillover effects between geographically structured units. Increasing availability of data often results in (spatial autoregressive) models where the number of observations is relatively small compared to the number of potential covariates. Standard estimation techniques in such environments therefore typically result in imprecise parameter estimates. In the case of severe overparameterization, where the dimensionality of the regressors exceeds the number of observations, standard estimation approaches may even be infeasible.

Bayesian model averaging techniques to alleviate the problems of overparameterization in spatial autoregressive models have been proposed \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{GEAN:GEAN703} and widely applied, particularly in the empirical study of regional economic growth (see, for example, doi:10.1080/17421770802353758, JAE:JAE2277, doi:10.1080/00343404.2012.678824, GEAN:GEAN12057, or cuaresma2018human). These approaches use weighted averages of parameter estimates based on a multitude of potential combinations of explanatory variables, rather than relying on inference based on a single model specification (for extensive discussions see steel2017bma or koop2003). Bayesian model averaging, however, involves the computation of marginal likelihoods for integrating out the underlying model uncertainty. In contrast to classical linear model frameworks, no closed-form solutions for marginal likelihoods in spatial autoregressive specifications are available \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{GEAN:GEAN703}, resulting in a severe computational burden especially in high-dimensional estimation problems.

Recent advances in the spatial econometric literature focus on extensions and generalizations of standard spatial autoregressive specifications. For example, some extensions aim at explicitly controlling for unobserved heterogeneity by finite mixture or threshold specifications (Piribauer2016; CORNWALL2017148), allowing for heterogeneous parameters across space (GEAN:GEAN12152), heteroskedastic specifications of the innovations (doi:10.1177/016001769702000107), considering continuous spatial effects (Laurini2017), or accounting for uncertainty in the spatial weight matrix specification and the underlying nature of spillover processes (LESAGE2007190; JORS:JORS12188; GEAN:GEAN12057), among several others. However, such extensions to the spatial autoregressive modeling framework further increase the emanating computational burden of Bayesian model averaging, rendering the approach computationally intractable.

For spatial autoregressive model specifications, work by Piribauer2016 and doi:10.1080/17421772.2016.1227468, for example, uses Bayesian stochastic search variable selection priors \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[SSVS,][]{doi:10.1080/01621459.1993.10476353} as an alternative to Bayesian model averaging. The comparative flexibility of this approach allows to easily extend and generalize the basic framework to more complex specifications.\footnote{The computational flexibility of these shrinkage priors comes from the fact that they can easily be implemented in standard Bayesian Markov-chain Monte Carlo (MCMC) algorithms without the need for calculating marginal likelihoods \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{doi:10.1080/01621459.1993.10476353}.} Specifically, the proposed hierarchical modeling approach assumes coefficients of a saturated model under scrutiny to come from a mixture of two Gaussians centered on zero with a spike and a slab component (low and high variance). A binary latent indicator thereby identifies promising subsets of predictors by shrinking unimportant components towards zero. Despite its elegance and appealing theoretical properties, difficulties arise for stochastic search variable selection in large data sets due to the necessity of stochastic search over an enormous space. This implies slow mixing and convergence during estimation \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{doi:10.1080/01621459.2014.960967}. Flexible alternatives for shrinkage in high-dimensional econometric frameworks have therefore been advocated \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{polson2010,griffin2017}.

In this paper we aim at generalizing variants of the Normal-Gamma \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[NG,][]{griffin2010} and the Dirichlet-Laplace \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[DL,][]{doi:10.1080/01621459.2014.960967} shrinkage priors to spatial autoregressive specifications. For the purpose of illustrating prior specific properties in the presence of spatial dependence, we carry out an extensive simulation study using synthetic data sets, considering different scenarios of the number of available covariates and degrees of sparsity of the coefficient vector. The paper moreover considers matrix exponential spatial specifications (MESS), introduced by LESAGE2007190 to model global spillover processes. Our results indicate that conventional stochastic search variable selection priors in the spirit of doi:10.1080/01621459.1993.10476353 work well in relatively low-dimensional settings, where the number of potential covariates do not exceed those of the observations. In high-dimensional modeling frameworks, however, both the Normal-Gamma and Dirichlet-Laplace shrinkage prior exhibit stellar empirical properties. They provide a high degree of adaptiveness of shrinkage, with the Normal-Gamma excelling in terms of precision, while the Dirichlet-Laplace performs particularly well in terms of point estimates.

The remainder of this paper is organized as follows. \Autoref{sec:econometrics} presents the Bayesian spatial econometric framework, followed by the introduction of Normal-Gamma and Dirichlet-Laplace shrinkage priors to spatial autoregressive specifications in \Autoref{sec:shrinkage}. A simulation study as a means to comparing the properties of the proposed shrinkage priors is presented in \Autoref{sec:simulation}. \Autoref{sec:application} illustrates the performance of the proposed shrinkage priors using pan-European regional economic growth data. \Autoref{sec:conclusions} concludes.

Econometric framework

We start by considering a model of the form

where $\bm{y}$ is an $N$-dimensional vector of dependent variables and $\bm{S}(\bullet)$ describes a linear transformation dependent on a not yet specified parameter. $\bm{X}$ is an $N\times K$ matrix of explanatory variables (with a vector of ones in the first column), and $\bm{\beta} = (\beta_1, \hdots, \beta_K)'$ is a $K$-dimensional vector of parameters. We assume the error term $\bm{\epsilon}$ to be normally distributed with zero mean and $N\times N$ variance-covariance matrix $\bm{\Omega}$.

Standard econometric models capturing spatial spillover effects are thoroughly discussed in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{lesage-pace_introductionspatial}. Conventional spatial autoregressive models (SAR) typically set $\bm{S}(\xi) = (\bm{I}_N - \xi \bm{W})$, where $\bm{W}$ is an $N\times N$ exogenous right stochastic spatial weights matrix. If observations $i=1,\dots, N$ and $j\neq i$ are considered neighbors, then $W_{ij} \neq 0$, otherwise $W_{ij} = 0$. By convention $W_{ii} = 0$, meaning that a region or other spatial unit is not considered a neighbor to itself \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{lesage-pace_introductionspatial}. $\xi$ is a (scalar) spatial parameter with stability condition $|\xi|<1$. The inverse of $\bm{S}(\xi)$ under the specifying assumptions for $\bm{W}$ and $\xi$ can be expressed as $(\bm{I}_N-\xi\bm{W})^{-1}= \sum_{l=0}^{\infty} \xi^l \bm{W}^l$, implying global geometric decay over space. An alternative to such a geometric decay pattern is given by the matrix exponential spatial specification (MESS), as proposed by LESAGE2007190, where we set

with the (scalar) parameter $\rho$ taking the role of measuring spatial dependence. Consequently, these modeling approaches nest the classical linear regression model in the case where the respective spatial dependence parameter is equal to zero. The main difference between MESS and SAR models is that spatial externalities are modeled by an exponential rather than a geometric decay.\footnote{It is moreover worth noting that a simple model extension frequently used in empirical applications also incorporates an explicit spatial structure in the exogenous variables by including a spatial lag of $\bm{X}$ in (ref). Such a specification is commonly referred to as a spatial Durbin model specification \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{lesage-pace_introductionspatial}.}

MESS is advantageous to the conventional SAR approach especially in a Bayesian framework. This is due to the properties of matrix exponentials presented in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{LESAGE2007190},

which facilitate the likelihood function given in (ref). Albeit both specifications typically produce rather similar estimates and inference (see, for example, LESAGE2007190, GEAN:GEAN12057, or STRAU2017221), matrix exponential spatial specifications have some potential computational advantages, specifically in high-dimensional data environments. A correspondence between spatial coefficients in the conventional SAR and MESS is stated by \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{LESAGE2007190}, where $\xi \approx 1 - \exp(\rho)$. However, it is worth noting that \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{DEBARSY20151} stress that no precise one-to-one correspondence can be established and the two spatial specifications should not be considered perfect substitutes, due to $\rho \in (-\infty,\infty)$ while stability conditions in the conventional SAR case require the parameter space of $\xi$ to be constrained \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{lesage-pace_introductionspatial}. Assuming a normally distributed error term as in (ref), the likelihood of the model is given by

where we define $\bm{\varepsilon} = \left(\bm{S}(\rho)\bm{y} - \bm{X}\bm{\beta}\right)$ for notational convenience.

Flexible shrinkage in high-dimensional Bayesian spatial autoregressive models

One conventional approach for cases where the vector of coefficients is expected to be sparse, but without prior knowledge which variables to exclude, is given by penalized least squares. A prominent example is the least absolute shrinkage and selection operator (LASSO) introduced by \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet[][]{10.2307/2346178}. This paper stresses the correspondence between the conventional LASSO and a Bayesian approach involving independent double-exponential priors on regression coefficients, a notion further elaborated on in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{doi:10.1198/016214508000000337}. Given the choice of a specific mixture distribution $\text{D}$, following \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{griffin2010}, these priors can typically be expressed as

for $r=1,\dots,K$. In this setting, the marginal distribution for $\beta_r$ has heavier than normal tails but places substantial mass on zero. \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{griffin2010} indicate that the standard discrete spike-and-slab prior (where a variant is presented in (ref)) may be represented in this form, dependent on the mixture distribution being chosen accordingly.

In addition, various other cases such as the double-exponential mentioned above (leading to the Bayesian LASSO, which is closely tied to the Normal-Gamma prior presented later on), arise for different mixing distributions \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{doi:10.1198/016214508000000337}. Note that in this basic framework, the specific shape of $\text{D}$ depends on the deterministic choice of prior hyperparameters. More flexibility and adaptive shrinkage can be achieved by introducing additional hierarchical layers of priors at this stage. Approaches in this spirit are termed global-local shrinkage priors \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{polson2010}, and seem to perform well in high-dimensional settings \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[for instance in time-series analysis, see][]{bittosfs,kastner,doi:10.1080/07350015.2016.1256217,feldkircherkastnerhuber}. In the following, we discuss two alternatives of continuous global-local shrinkage, the Normal-Gamma and the Dirichlet-Laplace shrinkage prior.

Normal--Gamma shrinkage prior

A variant of the Normal-Gamma global-local shrinkage prior, as proposed by griffin2010, is given by a scale mixture of Gaussians,

where $\psi_r$ is an idiosyncratic scaling parameter following a Gamma distribution that involves parameter specific shrinkage and can be collected in a vector $\bm{\psi} = (\psi_1,\hdots,\psi_K)'$. Overall shrinkage is governed by the global parameter $\lambda$ that also follows a Gamma distribution. According to \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{polson2010}, this setup reflects the necessity of noise reduction by shrinking all coefficient means to zero based on the global parameter, while allowing for signals to override this effect using local scaling parameters. Hyperparameters must be set by the researcher, where standard approaches in the literature include $d_0 = d_1 = 0.01$, and $\theta = 0.1$ as default. This implies heavy overall shrinkage of the parameters stemming from the global component but provides enough flexibility to detect individual non-zero coefficients if necessary.\footnote{In addition, it is worth mentioning that $\theta$ may be also integrated out, as for instance described in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{doi:10.1080/07350015.2016.1256217}. This could easily be implemented within the given structure. However, for the sake of simplicity, we refrain from doing so.} Setting $\theta = 1$ would present the case leading to the Bayesian LASSO \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{doi:10.1198/016214508000000337}. Deriving the posterior distribution of $\psi_r$, one finds that it follows a generalized inverse Gaussian distribution \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{griffin2010},

with $\text{GIG}(\zeta,\chi,\varrho)$ being parameterized such that its density $f(x) \propto x^{\zeta-1}\exp\{-(\chi/x + \varrho x)/2\}$, while the conditional posterior distribution of the global parameter is a Gamma distribution with

This setup defines the prior variance-covariance matrix on $\bm{\beta}$, denoted by $\ubar{\bm{\Sigma}}$, that is updated during estimation. We set the prior coefficient vector denoted by $\ubar{\bm{\beta}} = \bm{0}$. In the case of the NG prior, the variance-covariance matrix is given by $\text{diag}(\ubar{\bm{\Sigma}}) = \bm{\psi}$. The conditional posterior for coefficients $\bm{\beta}$ is of typical form \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{koop2003},

Dirichlet--Laplace prior

Even though working well empirically, \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{doi:10.1080/01621459.2014.960967} note that many aspects of global-local shrinkage priors based on scale mixtures of Gaussians are poorly understood theoretically and stress the necessity of simultaneous consideration of marginal properties of shrinkage priors and the joint distribution of obtained parameters.

As an alternative, we continue with the DL prior suggested by doi:10.1080/01621459.2014.960967. As opposed to the three hyperparameters to be set by the researcher in the NG case, a single hyperparameter suffices for establishing the DL setup. Similar to the NG prior, however, it is composed hierarchically of global and local shrinkage parameters, and may be structured as follows:

The local parameters are denoted by $\varphi_r$ and assigned an exponential prior distribution. In contrast to the single global parameter $\lambda$ in the case of the NG prior, the DL approach uses a vector of scales $(\phi_1 \tau, \hdots, \phi_K \tau)$ to provide more flexibility regarding idiosyncratic coefficient shrinkage, where $\bm{\phi} = (\phi_1,\hdots,\phi_K)$ is defined to lie in the ($K-1$)-dimensional simplex $\mathcal{S}^{K-1} = \{\mathbf{x} = (x_1, \hdots, x_K)':\ x_r\geq 0,\ \sum_{r=1}^{K} x_r = 1\}$ and is assigned a Dirichlet prior distribution with hyperparameter $a$ controlling the tightness of the prior. This quantity may again be integrated out as shown in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{doi:10.1080/01621459.2014.960967}. Based on theoretical results obtained by \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{doi:10.1080/01621459.2014.960967}, standard deterministic choices include $a = 1/K$ (the default setting used for the simulation study), but different choices such as $a = 1/2$ are justifiable theoretically. Notice that this allows the DL prior setup to be dependent on the dimensionality of the underlying model on theoretical grounds, different to both the SSVS and NG prior.

In the case of the DL prior, \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{doi:10.1080/01621459.2014.960967} show that $\varphi_r$ can be sampled efficiently involving independent inverse Gaussian distributions. This is done by obtaining $\tilde{\varphi}_r|\phi,\bm{\beta}$ from a reparameterization of the generalized inverse Gaussian distribution with $\mu_r = \phi_r\tau/|\beta_r|$,

where we subsequently set $\varphi_r = 1/\tilde{\varphi}_r$ to obtain draws from the conditional posterior distribution of $\varphi_r$. The global shrinkage component $\tau$ is again sampled from a generalized inverse Gaussian distribution

In order to sample $\bm{\phi}|\bm{\beta}$ we rely on Theorem 2.1 in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{doi:10.1080/01621459.2014.960967} and draw auxiliary variables $T_1,\hdots,T_K$ independently, with $T_r\sim\text{GIG}(a-1, 2|\beta_r|, 1)$. Afterwards, we set $\phi_r= T_r/\sum_{j=1}^{K} T_j$. \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{doi:10.1080/01621459.2014.960967} indicate this step as an important feature of their setup, as it significantly accelerates mixing and convergence. Consequently, the structure of the prior variance-covariance matrix $\ubar{\bm{\Sigma}}$ is updated, where $\text{diag}(\ubar{\bm{\Sigma}}) = (\varphi_1 \phi_1^2 \tau^2, \hdots, \varphi_K \phi_K^2 \tau^2)'$. The coefficient vector is then sampled analogously to the NG prior, based on (ref).

Priors for remaining parameters

Without loss of generality, we assume a homoskedastic error variance $\bm{\Omega} = \sigma^2 \bm{I}_N$, where $\bm{I}_N$ is an $N\times N$-dimensional identity matrix.\footnote{The generic notation in (ref), however, allows for various specifications of the variance-covariance matrix $\bm{\Omega}$, for instance reflecting heteroskedastic error terms \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, for instance,][in a spatial context]{doi:10.1177/016001769702000107}.} To complete the prior setup, we have to elicit prior distributions for $\sigma^2$ and the spatial dependence parameter $\rho$. Here we use standard specifications in the literature \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{lesage-pace_introductionspatial}. Specifically, we impose an inverse Gamma prior on the variance of the error term, $\sigma^2 \sim \text{G}^{-1}(\ubar{a},\ubar{b})$ with scalar hyperparameters $\ubar{a}$ and $\ubar{b}$ that we set equal to $0.01$, rendering them rather uninformative. Conditioning on all other parameters of the model and the data we obtain a conditional posterior density for the error variances,

The quantities for $\bm{\beta}$ and $\sigma^2$ are standard and can be sampled efficiently using Gibbs sampling (koop2003). For $\rho$, we follow LESAGE2007190 and lesage-pace_introductionspatial, by eliciting a normal prior $\rho \sim \text{N}(0,c)$, where $c$ may be chosen dependent on the prior belief of the researcher regarding the strength of spatial association in the data. For the simulation study, we again use a rather uninformative specification and choose $c = 10$. Since the posterior distribution of $\rho$ conditional on all other quantities of the model is in general not of well-known form,

it is sampled by a Metropolis-within-Gibbs step. We follow the standard approach by proposing a new value from a Normal distribution, $\rho^{\ast}\sim\text{N}(\rho_{t-1},\varsigma)$, where the subscript $t-1$ denotes the value of the parameter from the previous iteration of the sampling algorithm and $\varsigma$ is a tuning parameter. The acceptance probability of the proposal is calculated using

If the proposal is accepted, we set $\rho_{t} = \rho^{\ast}$. Otherwise, the proposal is rejected, and we retain the value from the previous draw. The process of generating proposals for $\rho$ is tuned during half of the burn-in phase of the algorithm by incrementally increasing or decreasing the variance of the proposal distribution $\varsigma$ to yield an acceptance rate for the parameter between $0.2$ and $0.4$. An overview of the algorithm employed can be found in (ref).

Simulation study

In this section, we evaluate the empirical properties and merits of the proposed shrinkage priors. Estimates are obtained using the Normal-Gamma (NG) and Dirichlet-Laplace (DL) shrinkage priors sketched above. As a benchmark specification, we compare the results of both setups with a stochastic search variable selection (SSVS) prior put forward by doi:10.1080/01621459.1993.10476353 and applied to spatial autoregressive models by Piribauer2016 and doi:10.1080/17421772.2016.1227468. The SSVS prior as a means to introducing shrinkage on the slope coefficients mimics Bayesian model-averaging frameworks. Details on the SSVS setup along with the concurrent Markov-chain Monte Carlo (MCMC) sampling algorithm are given in Appendix (ref). To further illustrate and underline the necessity of applying shrinkage in cases where the coefficient vector is sparse, we also provide results for the case where a rather uninformative prior variance-covariance matrix $\ubar{\bm{\Sigma}} = 1000 \bm{I}_N$ is used. This approach represents the basic setup given in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{LESAGE2007190}, and is labeled None, referring to the fact that no shrinkage is applied. For each specification employed, we use $2,000$ iterations, discarding the initial $1,000$ draws as burn-in. Inference is then based on the $T=1,000$ posterior draws for all parameters, where point estimates are calculated using the median of the respective sample.

Data Generating Process

Simulating a synthetic data set for different data generating processes requires decisions on the number of observations $N$ and explanatory variables $K$. Moreover, we have to set parameter values for the coefficients $\tilde{\bm{\beta}}$, the variance of the error terms $\tilde{\sigma}^2$ and the spatial autoregressive parameter $\tilde{\rho}$. We consider varying degrees of sparsity with respect to the coefficient vector. The final data generating process, reflecting the basic modeling approach as in (ref), has the form

where the first column in $\bm{Z}$ contains an intercept. Robustness for different combinations of parameters and generated variables is achieved via repeating each of the basic simulations sketched below for $100$ times and relying on stochastically setting most of the quantities involved. In particular, we set the required parameters as follows:

Diagnostics

We briefly consider diagnostics regarding obtained draws from the posterior distributions by evaluating trace plots and estimated densities. This allows us to present both the obtained marginal posterior densities for the respective parameters, and to discuss mixing and convergence properties. Diagnostic plots of non-zero coefficients typically mirror the well-known picture resulting from standard applications without introducing shrinkage. For illustrative purposes, we thus pick examples with true coefficients being equal to zero. A low-dimensional problem is given by choosing $q = 10$ and $K=50$, while the case of $K=150$ serves as high-dimensional scenario where applying shrinkage is required to find meaningful results on an acceptable level of precision.

The first example, for $K = 50$ exogenous covariates, is presented in (ref). The upper panel shows trace plots for all prior specifications, while the lower depicts the marginal posterior density of the parameter. The dashed red line indicates the true parameter value, the blue line is the median of the obtained posterior distribution. Even in this case, where ten out of the 50 slope coefficients in the parameter vector $\bm{\beta}$ are different from zero, we find that the standard specification without shrinkage (None) of the parameter space results in a comparatively high variance of parameter estimates, observable in (ref). Considering the three shrinkage approaches, we find that the DL and NG priors concentrate more point mass on zero by design, while the SSVS prior shows considerable variation in the interval around zero.

A more interesting case is given in (ref). In this setup we include $150$ covariates, where only ten exhibit non-zero values by construction. As is evident considering None, the variance around parameter estimates is huge when the number of parameters to estimate exceeds the number of observations. Since the semi-automatic setup of the SSVS prior set forth in (ref) relies on these quantities to scale prior hyperparameters, this is influential regarding resulting parameter estimates for the SSVS approach. Notice that the obtained scaling factor renders the variance of the spike component comparatively large due to the specific setup involved, and SSVS does not provide enough shrinkage, as shown in (ref).

This issue is not observable in the case of the DL and NG priors, providing evidence for their excellent adaptive shrinkage properties in high-dimensional settings. For the NG prior, it is worth mentioning that the point mass placed on zero is comparable to the case for $K = 50$. By the fact that the single hyperparameter of the DL prior is set to $1/K$ in the default case, higher dimensional problems result in even stronger shrinkage towards zero for true zero coefficients. Evidence of the increased tightness of the prior can be observed with regard to the density scales in (ref). Moreover, note that we also observe that the shrinkage priors are able to recover coefficients close to zero equally well.

Following this brief description of diagnostic plots, we again use the scenarios above as examples and consider all of the resulting parameter estimates jointly. Plots of the true regression coefficients against their posterior medians are depicted in (ref) and (ref). The black line indicates the 45 degree line, reflecting the objective of perfectly estimated parameters. In the low-dimensional example given in (ref), we find that point estimates based on the median of obtained posterior distributions for the parameters are mostly correct also in the case where no shrinkage is applied. However, note that these estimates are rather imprecise, which renders them suboptimal in terms of significance for interpretation, and also for applications involving predictions. The three shrinkage approaches closely mirror each other, even though it is worth noting that the NG and DL prior allow for even more precise estimates in terms of the variance.

A completely different case emerges for $150$ covariates. (ref) showcases severe problems regarding the standard prior specification without shrinkage (None). Even though the SSVS prior performs slightly better than this approach, it is evident that the default setup regarding the scaling of the hyperparameters is suboptimal in high-dimensions. In particular, and as we will discuss below, the SSVS prior actually tracks non-zero coefficients quite well, but imprecisely. Issues arise mainly from its disability to capture zero coefficients adequately in the given amount of iterations of the MCMC algorithm, pointing towards mixing problems and suboptimal convergence as stated in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{doi:10.1080/01621459.2014.960967}. In this respect, both the NG and DL prior are superior. The extraordinary empirical properties of these priors are evident in (ref), where the majority of the coefficients is recovered almost perfectly. Note that the DL prior provides slightly tighter shrinkage of true zero coefficients.

Evaluation of performance

Turning to an overview of all estimated parameters of the models, with respect to a total of eight different scenarios given by combinations of $K = \{50,100,150,200\}$ and $q = \{10, 20\}$, we present root mean squared error (RMSE) measures.\footnote{We adopt the standard root mean squared error measure for the posterior median $\hat{\bm{\beta}}$, $\text{RMSE}(\hat{\bm{\beta}}) = (\sum_{k=1}^K(\hat{\beta}_k-\tilde{\beta}_k)^2/K)^{\frac{1}{2}}$ and also calculate the corresponding quantity with respect to each draw from the MCMC algorithm, $\text{RMSE}_\text{dr}(\hat{\bm{\beta}}) = \sum_{k=1}^K(\sum_{t=1}^{T}(\hat{\beta}_{k,t}-\tilde{\beta}_k)^2/T)^{\frac{1}{2}}/K$. This serves as a means to illustrate differing degrees of precision of the estimators.} The specific choice of the number of included covariates and number of non-zero coefficients implies that we have degrees of sparsity ranging from 5 to 40 percent for the parameter vector. Average results for $100$ iterations are presented in (ref) and (ref).

The results for the standard RMSE measure closely resemble the notions obtained in the discussion above. In particular, we find that the SSVS prior performs slightly better than other shrinkage approaches in low dimensional settings where $K < N$, both in terms of estimation errors for parameters in $\bm{\beta}$, $\sigma^2$, the spatial dependence parameter $\rho$ and also time elapsed for cycling through the MCMC algorithm. However, both the NG and DL prior show promising values regarding RMSEs and can be considered as feasible alternatives. For the case of $K = N$, we find a different picture. The approach without shrinkage (None) performs rather poorly, evidenced by large errors in terms of all parameters. In particular, based on poor estimates for the coefficients in $\bm{\beta}$, this results in the inability to correctly estimate $\sigma^2$. We find particularly good performance measures for the DL prior with the NG prior being close second. This finding is particularly pronounced in the denser case with respect to the coefficient vector ($q=20$), where the NG prior outperforms the DL prior in terms of precision regarding the error variances and the spatial dependence parameter.

Turning to the higher dimensional cases where $K$ exceeds $N$ it is worth noting that SSVS runs into similar problems than the standard approach, where poor parameter estimates result in huge values for the error variances. Since this problem only occurred in a minor subset of the $100$ iterations per scenario, we chose to exclude these erroneous simulations. In this setup, the approach without shrinkage again appeared to perform poorly, as already observed in the case with $K = 100$. Second, due to the notion that SSVS typically shows suboptimal mixing and convergence properties in high dimensional settings -- stemming from the necessity of stochastic search over an enormous parameter space -- it appears that the default setting of $2,000$ iterations for the MCMC algorithm is not sufficient to obtain reasonable posterior distributions for the parameters. Finally, we find that the DL prior performs best regarding coefficients and error variances -- compared to the NG prior, which is only superior in terms of estimating the spatial dependence parameter $\rho$. Note that estimates in the case of more dense coefficient vectors ($q=20$) are typically slightly worse, due to the necessity of estimating more none-zero parameters with the same number of observations. A similar picture is present in the case of $K \gg N$. Mirroring the results already obtained in the previous scenario, the DL prior appears to be the fastest of the shrinkage priors. Interestingly, the superiority of the DL prior regarding estimates of the error variances gets even larger when compared to the worse estimates resulting from the NG prior.

The main results discussed above also hold in the case of the adapted RMSE focusing on density predictions in (ref). Recall that we use this measure for the purpose of gaining insight into the precision of the obtained parameter estimates. Interestingly, even though being inferior to the DL prior in terms of point estimates, we find that the NG prior produces more precise estimates around the true values of the parameters in the coefficient vector $\bm{\beta}$. However, this finding does not carry over to estimates of the error variances $\sigma^2$, where the DL prior outperforms the NG approach. The spatial dependence parameter is estimated precisely in all cases regarding sparsity and number of covariates within the NG and DL shrinkage prior framework, and in scenarios where the SSVS prior and the standard approach still perform well.

Empirical illustration

In this section we aim at illustrating the performance of the proposed model specification using real data on pan-European regional economic growth and its empirical determinants. Specifically, we consider a cross-regional spatial Durbin model specification (see, for example, lesage-pace_introductionspatial) which also allows for spatially lagged explanatory variables as potential covariates. The specification used can be written as

where $\boldsymbol{y}$ denotes an $N$-dimensional column vector of regional economic growth rates of per capita gross value added and the $N\times 1$ error term $\bm{u}$ is defined as before as iid normal. The $N\times K$ matrix $\boldsymbol{X}$ contains the set of potential predictors. $\boldsymbol{W}$ is an $N\times N$ row-stochastic spatial weight matrix as defined before. $\boldsymbol{S}(\rho)$ denotes a matrix exponential spatial filter with corresponding scalar parameter $\rho$ as defined before. $\boldsymbol{WX}$ is the spatial lag of the explanatory variables and explicitly incorporates spatially lagged information of $\boldsymbol{X}$ from neighboring regions, resulting in the spatial Durbin model specification mentioned above. In this empirical illustration we follow the typical structure of (spatial) growth regressions by assuming that information in matrix $\boldsymbol{X}$ is measured at the beginning of the sample period (which is the year $2000$).

Regions, spatial weights and data

For the empirical illustration we use data on regional economic growth on a sample of $273$ European NUTS-2 regions of 28 European countries. The dependent variable in the regression framework is the average annual growth rate of per capita gross value added in the period $2001-2010$. A detailed list of the regions used in the application is given in (ref) in (ref). Table (ref) provides detailed information on the set of potential covariates in the matrix exponential spatial growth specification.

The set of predictors is in line with recent empirical applications on regional economic growth. Specifically, the matrix of explanatory variables contains information on the initial level of economic growth rates, human capital endowments, proxies for knowledge capital stocks, regional population structure, infrastructure, the region-specific industry mix, and other socio-economic quantities. It is worth noting that the spatial lag of the (non-constant) explanatory variables in the spatial Durbin framework sketched above results in doubling the set of potential covariates depicted in (ref). For the specification of the spatial weight matrix $\boldsymbol{W}$ we used a $10$-nearest row-stochastic specification.\footnote{Robustness checks using alternative numbers of nearest neighbors appeared to exert negligible effects on the results.}

Empirical results

In this subsection we present the results of the stochastic search variable selection (SSVS) prior, along with Normal-Gamma (NG), and Dirichlet-Laplace (DL) shrinkage setups including alternative prior hyperparameter specifications. Specifically, for the Normal-Gamma shrinkage prior a natural candidate is the Bayesian LASSO \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{doi:10.1198/016214508000000337}, achieved by setting the prior hyperparameter $\theta=1$. Alternatively, we also consider $\theta = 0.1$, which refers to the generalized version characterized by heavier overall shrinkage. For the Dirichlet-Laplace (DL) shrinkage prior, we consider $a = 1/2$, which presents a usual benchmark in empirical research, while $a = 1/K$ is the default prior specification based on theoretical considerations \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{doi:10.1080/01621459.2014.960967}.

Table (ref) depicts estimation results for the competing prior setups under scrutiny. Note that only covariates where the corresponding slope coefficient is statistically significant from zero based on the lower $10$ percent and upper $90$ posterior interval in at least one of the candidate specifications are reported. For each prior setup, (ref) reports the posterior means (labeled Mean) and corresponding posterior standard deviations (labeled SD) for the parameters under scrutiny. Both quantities are based on $3,000$ retained draws of the MCMC algorithms described in (ref).\footnote{Note that for all alternative specifications we used a number of $12,000$ iterations with the first $3,000$ serving as burn-ins. To reduce the autocorrelation of the draws for $\rho$ we have used thinning by considering only every third draw (see, for example, koop2003), resulting in a total of $3,000$ draws for posterior inference.}

Highly significant estimates for the spatial parameter $\rho$ result in all specifications considered, ranging from $-0.811$ to $-0.972$. Using the correspondence between the spatial parameter ($\xi$) in conventional spatial autoregressive frameworks and its matrix exponential spatial counterpart ($\rho$) stated by \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{LESAGE2007190}, we find an implied spatial autoregressive parameter $\xi$ between $0.56$ and $0.62$. This degree of spatial dependence resembles the findings in recent studies on regional spatial economic growth in Europe \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, for example,][]{doi:10.1080/00343404.2012.678824,doi:10.1080/17421770802353758}.

Turning attention to the (non-constant) potential growth determinants, (ref) shows four slope coefficients that are statistically significant in all prior setups. These variables are the initial level of income, the old-age dependency ratio, as well as both the share of low-educated employment and its corresponding spatial lag. Overall, (ref) shows rather similar results for both the magnitudes of the estimated effects as well as their significance. However, a notable exception is the Bayesian LASSO setup NG($\theta=1$), which highlights a markedly higher amount of significant slope parameters as compared to the competing specifications.

As expected, the initial income variable appears to negatively affect regional economic growth rates in all specifications, providing evidence for conditional convergence among the regions in the sample. However, the SSVS setting also shows significant posterior mean of the spatially lagged initial income variable, pointing towards the existence of positive growth spillovers emanating from the level of income of neighboring regions. This means that regions seem to benefit from the spatial proximity of rich regions in terms of accelerated growth rates. The old-age dependency ratio (measured in terms of the ratio of the number of people aged 65 and over to the working age population) appears to exhibit a negative impact on regional economic growth rates. Except for the Bayesian LASSO specification (NG ($\theta=1$)), the estimated coefficients are rather similar.

Interestingly, the results presented in (ref) corroborate the findings of previous studies (see, for example, GEAN:GEAN12057, or doi:10.1080/00343404.2012.678824) by detecting the share of low educated working age population (Lower education workers) as being more robustly correlated to regional economic growth as compared to a measure of tertiary education attainment (higher education workers). As expected, lower education workers appear to exhibit a negative impact on regional growth in all specifications considered. However, it is worth noting that the spatial lag of this variable appears to exhibit a positive impact, indicating that positive effects on income growth to neighboring regions. Work by olejnik2008, for example, argue that an increase in the lower educated labor force might result in positive growth spillovers, by assuming that such an increase might be primarily due to migration of workers between regions.

Concluding remarks

Dealing with model uncertainty in spatial autoregressive model specifications has been subject to numerous studies in general, especially in the empirical economic growth literature. However, spatial econometric applications typically rely on Bayesian model-averaging techniques, which suffer from severe drawbacks both in terms of computational time and possible extensions to more flexible model specifications. In spatial autoregressive models, the computational burden emanates from the calculation of marginal likelihoods, where no closed form solutions are available. Recent contributions to the literature as a means to alleviating the computational difficulties include a variant of the conventional stochastic search variable selection prior discussed in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{doi:10.1080/17421772.2016.1227468}. However, shortcomings of this approach include slow mixing and convergence properties in the presence of a large number of available covariates and difficulties in choosing prior hyperparameters.

This paper aims at generalizing two absolutely continuous hierarchical shrinkage priors -- the Normal-Gamma and the Dirichlet-Laplace shrinkage prior -- to the matrix exponential spatial specification in order to alleviate the above-mentioned drawbacks of both Bayesian model averaging and standard stochastic search variable selection priors. The proposed frameworks allow for flexible and adaptive, but also computationally efficient stochastic variable selection, where extensions to basic spatial model specification can be easily implemented in Markov chain Monte Carlo sampling algorithms. For illustrative purposes, and to evaluate prior-specific properties in the presence of spatial dependence, the paper carries out an extensive simulation study using synthetic data sets. An empirical illustration is given by a study on economic growth of European regions.

Our results indicate that standard stochastic search variable selection priors work particularly well in relatively low-dimensional settings. However, severe problems occur in high-dimensional settings, where the number of potential covariates relative to the number of available observations becomes large. The proposed global-local shrinkage priors provide the required adaptiveness of shrinkage in a flexible and computationally efficient way. The Normal-Gamma shrinkage prior excels in terms of precision of estimates, while the Dirichlet-Laplace shrinkage prior slightly outperforms the former in terms of point estimations of parameters. Both proposed shrinkage priors can be considered valuable tools as a means to incorporate model uncertainty in spatial autoregressive frameworks, regardless of dimensionality of the problem at hand.

\singlespacing \addcontentsline{toc}{section}{References}

\onehalfspacing