Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
57,534 characters · 22 sections · 0 citation commands
Fast Two-Stage Variational Bayesian Approach to Estimating Panel Spatial Autoregressive Models with Unrestricted Spatial Weights Matrices
\affil[1]{University of Leicester} \affil[2]{University of Leicester, Bank of Greece and University of Pretoria} \affil[3]{Bank of Greece and Hoover Institution, Stanford University }
\doublespacing
The spatial autoregressive (SAR) models, first proposed by Cliff and Ord (1973), have been widely used in the literature to investigate the spatial dependence in cross-sectional units (e.g., Anselin, 1988; Baltagi et al., 2003, 2013, Lee and Yu, 2010). In practice, the spatial weights matrices are usually set a priori based on geographical distances or economic characters (e.g. Cliff and Ord, 1973, Anselin, 1988, Case, 1991). This is not surprising as a spatial weights matrix potentially involves $N^{2}-N$ interrelationships between $N$ spatial units, which makes it difficult to estimate, especially when $N$ is large.
In recent years, a number of variable selection and parameter shrinkage methods have been developed to estimate the spatial weights matrices of panel SAR models. Among them, methods resorting to least absolute shrinkage and selection operator (Lasso) of Tibshirani’s (1996) and its variants have gained a lot of attention. For example, Basak et al. (2018) propose to estimate a triangular weights matrix under the assumption of recursive ordering. Ahrens and Bhattacharjee (2015) develop a two step Lasso estimator to identify the weights matrix. Lam and Souza (2019) estimate the weights matrix using adaptive Lasso with sparse adjustment in mind. Most of the studies, however, usually impose sometime unrealistic restrictions on the model's coefficients or covariances. Krock et al. (2021) develop a graphical Lasso approach to estimating the unrestricted covariances. Their method, however, does not deal with the impacts of any possible exogenous variables. Moreover, to our knowledge, few of those studies focus on the short panels where $N$ is large while $T$ is small. Only recently, Krisztin and Piribauer (2023) and Piribauer et al. (2023) use a hierarchical prior set-up to identify sparsity when $N$ exceeds $T$ by a large margin, assuming the same spatial parameter for all the cross-sectional units.
This paper contributes to the SAR literature by developing a fast two-stage variational Bayesian (VB) approach to estimating panel SAR models with unknown spatial weights matrices. We do not impose any restrictions on spatial weights matrix or the covariance functions, hence our approach lets the data speak. The prior we used for Bayesian regularization is the Dirichlet–Laplace (D-L) prior of Bhattacharya et al. (2015). With D-L prior, the entire posterior distribution concentrates at the optimal rate. This nice feature remains unchanged when the number of parameters to be estimated is much larger than the number of observations, providing strong theoretical justifications for the two-stage VB's effectiveness in uncovering the spatial dependencies in a short panel.
Our second contribution is to develop VB methods to estimate large covariance matrices with a global-local shrinkage prior. We are among the first in the literature to develop VB estimator for large covariance matrices with unknown sparsity patterns. Our VB methods using D-L prior can be easily extended to allow for other popular priors such as the graphic Lasso of Wang (2012), the half-Cauchy prior of Makalic and Schimidt (2016) and the graphical horseshoe prior of Li et al. (2019). This is not trivial as VB is a more computationally efficient alternative of Markov Chain Monte Carlo (MCMC), and our approach can be used in estimating other popular models involving large covariance matrices such as large Bayesian vector autoregressions (BVARs).\footnote{Matlab code is provided.}
We have conducted a wide range of simulation studies using a traditional panel SAR model and a panel SAR model that takes account of the simultaneous relationships between cross-sectional groups.\footnote{This research used the ALICE High Performance Computing Facility at the University of Leicester.} Monte Carlo experiments show that two-stage VB is accurate and computationally efficient when $T\gg N$, which usually is more pertinent to macroeconomic and financial data. When $N\gg T$, which tends to be more relevant to microeconomic data, two-stage VB estimates tend to have slightly larger biases and empirical standard deviations. Tighter priors can help lessen that problem.
In empirical applications, we use two-stage VB to estimate a two-equation simultaneous spatial model using panel data comprising sovereign bond ratings and spreads of ten eurozone countries. The research is motivated by the importance of a more in-depth understanding in how cross-country contagion works in the southern euro area countries and their northern counterparts, especially in the European debt crises (e.g. Gibson et al., 2021, Hall et al., 2022). Using rolling windows of 2 years, we compare and contrast how a rate or spread change in the south and the north affect each individual member country. Our results provide ample evidence for the marked differences between the two country groups, both in how a shock to the south and a shock to the north impact a particular country in a very different way and in how the shocks to a particular country group affect a southern and a northern country very differently. Our results also highlight the big impacts of the global financial crises and how the south was severely affected by the European debt crises while the north was much less so.
The rest of the paper is organised as follows. Section 2 extends the traditional panel SAR models to an unrestricted panel SAR. Section 3 develops the two-stage VB. Section 4 conducts Monte Carlo studies. Section 5 applies two-stage VB to the sovereign bond ratings and spreads data of ten euro area countries. Section 6 concludes. Sources for data used in the empirical example are listed in the Appendix. Detailed VB derivation formulas and more extensive Monte Carlo results are relegated to Online Supplements.\footnote{Online Supplements can be found at https://github.com/DBayesian/GKT2022}
In this section, we start from a traditional standard panel SAR and then relax the restrictions imposed upon it in steps, with the aim of giving a flavour of the differences between the traditional model and the unrestricted panel SAR model that we set to estimate using two-stage VB.
Let $Y$, $X$ and $V$ denote the $T\times N$ matrix of endogenous variables, $T\times (Nm)$ matrix of exogenous variables, and $T\times N$ matrix of disturbances, respectively. A traditional panel SAR model takes the following form:
where $y_{t}=(Y_{t1},Y_{t2},...,Y_{tN})'$ is the $N\times 1$ vector of observations of the dependent variables, $W_{n}$ is an $N\times N$ known spatial weights matrix with zero diagonal entries, $X_{t}=\left(
\right) $ is the $N\times m$ matrix of exogenous variables, with $X_{t,i}$ denoting the $1\times m$ row vector of exogenous variables associated with dependent variable $y_{ti}$, $\beta$ is an $m\times 1$ vector of parameters, $\lambda$ is a scalar parameter, and $u_{t}=(V_{t1},V_{t2},...,V_{tN})'$ is the $N\times 1$ i.i.d error terms with mean zero and diagonal covariance matrix $ \left(
\right) $.
Model ((ref)) imposes the following unrealistic restrictions on the data generating process: 1) $W_{N}$ is predetermined, in a fashion that is not related to the variations in the data; 2) $\lambda$ and $\beta$ remain the same across equations associated with different dependent variables; and 3) the covariance matrix of $u_{t}$ is diagonal with the same diagonal entries.
Relaxing those restrictions, model ((ref)) can be written as:
where $\widetilde{W}_{N}$ is an $N\times N$ unknown spatial weight matrix with zero diagonal entries, $\tilde{\lambda}$ is a $N\times 1$ parameter vector, $\odot$ is the Hadamard product, $\widetilde{X}_{t}=\left(
\right) $, $\tilde{\beta}$ is a $Nm\times 1$ parameter vector, and $u_{t}$ is i.i.d with mean zero and diagonal covariance matrix $\tilde{\Sigma}$ with diagonal entries that can be different from each other. In this model, the dimensions of $X_{t,i}$ and $X_{t,j}$ for $i\neq j$ can differ from each other. Let the dimension of $X_{t,i}$ to be $1\times m_{i}$. The dimension of the parameter vector $\tilde{\beta}$ is thus $(\sum_{i=1}^{n}m_{i})\times 1$.
Note that model ((ref)) is quite flexible. For example, with appropriate restrictions, it can be easily transformed back into the traditional form described in ((ref)) or a panel SAR containing the simultaneous cross-sectional spatial relationship as described in Yang and Lee (2017) and Liu and Saraiva (2019).
Since our main concerns in panel SAR models are the spillover effects, there is therefore little research interest in separately identifying $\tilde{\lambda}$ and $\widetilde{W}_{N}$. What we care about is the product $(\tilde{\lambda} \otimes l_{N}) \odot \widetilde{W}_{N}$, where $l_{N}$ is a $N\times 1$ column of ones and $\otimes$ is the Kronecker product, as $(\tilde{\lambda} \otimes l_{N}) \odot \widetilde{W}_{N}$ is the $N\times N$ parameter matrix which captures the spillover effects between spatial units.
Let $\Lambda= (\tilde{\lambda} \otimes l_{N}) \odot \widetilde{W}_{N}$. Model ((ref)) can be written as:
where $I_{N}-\Lambda$ is nonsingular and the characteristic roots of $I_{N}-\Lambda$ lie within the unit circle.
The attractiveness of model ((ref)) is that it turns an unrestricted panel SAR model into a system of simultaneous equations (SEM). As shown in Zellner and Thell (1962) and Fox (1979), the $i^{th}$ equation in model ((ref)) is just identified if $ \sum_{i=1}^{N}m_{i}=N-1+m_{i}$ and over-identified if $\sum_{i=1}^{N}m_{i}>N-1+m_{i}$. Under these circumstances, a myriad of estimation methods, such as two-stage least squares (2SLS), three-stage least squares (3SLS), maximum likelihoods methods and simultaneous generalized method (GMM), can be used to uncover the structural parameters $\Lambda$ and $\tilde{\beta}$. Moreover, standard tests can be developed to test the restrictions on $\tilde{\lambda}$, $\widetilde{W}_{N}$ and $\tilde{\beta}$, if those restrictions are of the researchers' concerns.
This paper proposes to estimate $ \Lambda $ and $\tilde{\beta}$ in two stages as it is computationally simple. To estimate the parameters associated with the $i^{th}$ individual dependent variable, in the first stage, we estimate
where $ Y_{/ i}$ is the $T\times (N-1)$ matrix of dependent variables except for the $i^{th}$ dependent variable, and $E_{i}$ is a $T\times (N-1)$ matrix of error terms whose precision matrix might not be diagonal.
Making use of the estimated $\Upsilon_{i}$, in the second stage, we estimate
where $ y_{i}$ is the $T\times 1$ vector of the $i^{th}$ dependent variable, $\widehat{Y}_{/i}=X\Upsilon_{i}$, $\Lambda _{i\bullet}$ is the $1\times (n-1) $ vector of the $i^{th}$ row of $\Lambda$ with $\Lambda_{ii}$ dropped, and $\tilde{\beta}_{i}$ is the corresponding coefficients in $\tilde{\beta}$.
Note that in a panel SAR model, we can have $N\gg T$ and $(N-1+m_{i})\gg T$, which makes it difficult or even impossible to uncover $B_{i}$ and $\Lambda _{i\bullet}$ using traditional estimation techniques.
As detailed in Ormerod and Wand (2010) and Blei et al. (2017), the essence of VB is to use appropriate densities from a mean field variational family to approximate the posterior densities through minimizing the Kullback-Leibler divergence, which is equivalent to maximising the evidence lower bound (ELBO). As a more efficient alternative to MCMC, VB has been increasingly used in sophisticated models involving large data where MCMC is too computationally expensive or even untenable (e.g. Gefang et al., 2020, 2022, Loaiza-Maya et al., 2022).
In the two-stage VB, we identify the parameters in model ((ref)) equation by equation. In each stage, we update the parameters using the approximate $q$ densities by iterations. The convergence of the algorithm can be measured by the changes in ELBO across iterations is less than a convergence criteria. When the number of parameters is large, however, calculating ELBO can be time consuming. It is therefore more convenient to check if convergence has occurred by examining if the VB estimates of parameters stop changing across iterations.
In the first stage, we estimate model ((ref)) in order to construct the predicted value of $ Y_{/i}$.
Let $\gamma=vec(\Upsilon_{i})$. We set hierarchical D-L prior for the $j^{th}$, for $j=1,...,np$, element of $\gamma$ as follows:
where $DE(\bullet)$ denotes Double Exponential or Lapalace distribution, $Dir(\bullet)$ denotes Dirichlet distribution, $G(\bullet)$ denotes Gamma distribution, $n=N-1$, and $p=\sum_{i=1}^{N}m_{i}$.
Next, we set Exponential priors and D-L priors for the elements of $\Omega$, the precision matrix of $E_{i}$, as follows:
where, with a slight abuse of notations, we use $\omega_{ii}$ and $\omega_{ij}$ to denote the diagonal and off-diagonal elements of $\Omega$.
Following Wang's (2012) Block Gibbs sampler to update the relevant parameters and hyperparameters, we use the last column and row of $\Omega$ as an example on how to update $\Omega$.
Let $S=E_{i}'E_{i}$ and $H$ be the $n\times n$ matrix with $0$ diagonal elements and the off diagonal element at $i^{th}$ row and $j^{th}$ column be $\psi_{\omega,ij}\phi_{\omega,ij}^{2}\tau_{\omega}^{2}$. Partition $\Omega$, $S$ and $H$ as follows:
where $-n$ denotes the set of all indices except for $n$.
Relegating technical details to Online Appendix A, we outline the VB approximation densities as follows:
where $$\overline{V}=(V^{-1}+\overline{\Omega}\otimes(X'X))^{-1},$$ $$\overline{\gamma}=\overline{V}(\overline{\Omega}\otimes X')vec(Y_{/i}), $$ and $$V^{-1}=diag(\frac{1}{\overline{\psi_{1}}\hspace{0.1cm}\overline{\phi_{1}^{2}}\hspace{0.1cm}\overline{\tau^{2}}},...,\frac{1}{\overline{\psi_{np}}\hspace{0.1cm}\overline{\phi_{np}^{2}}\hspace{0.1cm}\overline{\tau^{2}}})$$
Let $\chi=\sum_{j=1}^{np}2(\overline{b}_{j}^{2}+\overline{V}_{jj})^{1/2}\frac{1}{\overline{\phi_{j}}}$, we have $$\bar{\tau}=\frac{\sqrt{\chi}K_{npa-np+1}(\sqrt{\chi})}{K_{npa-np}(\sqrt{\chi})}$$ $$\overline{\tau^{2}}=\bar{\tau}^{2}+\chi[\frac{K_{npa-np+2}(\sqrt{\chi})}{K_{npa-np}(\sqrt{\chi})}-(\frac{K_{npa-np+1}(\sqrt{\chi})}{K_{npa-np}(\sqrt{\chi})})^{2}]$$ where $K_{\ast}[\bullet]$ is the modified Bessel functions of the second kind.
where $iG(\bullet)$ denotes s Inverse Gaussian distribution.
Let $\rho=\sqrt{\frac{\overline{\phi_{j}^{2}}\hspace{0.2em} \overline{\tau^{2}}}{\overline{\gamma}_{j}^{2}+\overline{V}_{jj}}}$, $$\overline{\psi_{j}^{-1}}=\rho$$ and $$\overline{\psi}_{j}=1+1/\rho$$
where $giG(\bullet)$ denotes the generalized inverse Gaussian distribution.
Let $\varpi=2 \sqrt{\overline{\gamma}_{j}^{2}+\overline{V}_{jj}}$, we have $$ \overline{ \xi_{j}}=\frac{\sqrt{\varpi}K_{a} (\sqrt{\varpi})} {K_{a-1}(\sqrt{\varpi})}, $$ and $$var(\xi_{j})=\varpi \{\frac{K_{a+1}(\sqrt{\varpi})} {K_{a-1}(\sqrt{\varpi})} -[\frac{K_{a}(\sqrt{\varpi})} {K_{a-1}(\sqrt{\varpi})}]^{2}\},$$ where $var(\bullet)$ denotes the variance.
Scaling $\xi_{i}$, we have $$\overline{\phi_{j}}=\frac{\overline{ \xi_{j}}}{\sum_{j=1}^{np}\overline{ \xi_{j}}},$$ and $$\overline{\phi_{j}^{2}}=\overline{\phi_{j}}^{2}+\frac{var(\xi_{j})}{(\sum_{j=1}^{np}\overline{ \xi_{j}})^{2}}$$
Thus, the optimal $q$ density of $\phi_{i,j}$ takes the following form:
Let $b_{1}=\omega_{n,n}- \omega _{-n,n}' \Omega_{-n,-n}^{-1} \omega_{-n,n}$ .
where $$\overline{s}_{n,n}=\frac{1}{2}(s_{n,n}+tr(X'X\overline{V_{n}})+\underline{s}),$$ and $$\overline{V_{n}}=V_{(n-1)\times p+1:n\times p, (n-1) \times p+1:n \times p}.$$
Hence $$\overline{b_{1}}=\frac{\frac{T}{2}}{\overline{s}_{n,n}}$$
Here we use $b_{2}$ to denote $\omega _{-n,n}$. Let $\overline{s}_{-n,n}=s_{-n,n}+\tilde{s}_{-n,n}$, where $\tilde{s}_{-n,n}$ is a $(n-1)\times 1$ vector with the $j^{th}$ element being $tr(X'XA_{j})$ and $A_{j}=V_{(j-1)\times p+1:j \times p,(j-1)\times p+1:j \times p }$.
where $$\overline{C}=(2\overline{s}_{n,n}\Omega_{-n.-n}^{-1}+\overline{H}^{*-1})^{-1},$$ and $$\overline{b_{2}}=(-\overline{C}\overline{s}_{-n,n}).$$
Note that $\overline{H}^{*}=diag(\overline{h}_{-n,n})$, and $j^{th}$ element of $\overline{h}_{-n,n}$ is $\overline{\psi_{\omega,jn}}\overline{\phi_{\omega,jn}^{2}}\overline{\tau_{\omega}^{2}}$.
Let $\chi_{\omega}=\sum\limits_{j<k}2(\overline{\omega_{jk}}^{2}+\overline{V}_{jj})^{1/2}(\overline{\phi_{\omega,jk}})^{-1}$, we have $$\overline{\tau_{\omega}}=\frac{\sqrt{\chi_{\omega}}K_{\frac{n^{2}-n}{2}(a_{\omega}-1)+1}(\sqrt{\chi_{\omega}})}{K_{\frac{n^{2}-n}{2}(a_{\omega}-1)}(\sqrt{\chi_{\omega}})}$$
$$\overline{\tau_{\omega}^{2}}=\overline{\tau_{\omega}}^{2}+\chi_{\omega}[\frac{K_{\frac{n^{2}-n}{2}(a_{\omega}-1)+2}(\sqrt{\chi})}{K_{\frac{n^{2}-n}{2}(a_{\omega}-1)}(\sqrt{\chi_{\omega}})} -(\frac{K_{\frac{n^{2}-n}{2}(a_{\omega}-1)+1}(\sqrt{\chi_{\omega}})}{K_{\frac{n^{2}-n}{2}(a_{\omega}-1)}(\sqrt{\chi_{\omega}})})^{2}]$$
Let $\rho_{\omega}=\sqrt{\frac{\overline{\phi_{\omega,jn}^{2}}\hspace{0.2em} \overline{\tau_{\omega}^{2}}}{\overline{\omega_{jn}}^{2}+\overline{C}_{jj}}}$, $$\overline{\psi_{\omega,jn}^{-1}}=\rho_{\omega}$$ and $$\overline{\psi_{w,jn}}=1+1/\rho_{\omega}$$
Let $\varpi_{\omega}=2 \sqrt{\overline{\omega_{jn}}^{2}+\overline{C}_{jj}}$, we have $$ \overline{ \xi_{\omega,jn}}=\frac{\sqrt{\varpi_{\omega}}K_{a_{\omega}} (\sqrt{\varpi_{\omega}})} {K_{a_{\omega}-1}(\sqrt{\varpi_{\omega}})}, $$ and $$var(\xi_{\omega,jn})=\varpi_{\omega} \{\frac{K_{a_{\omega}+1}(\sqrt{\varpi_{\omega}})} {K_{a_{\omega}-1}(\sqrt{\varpi_{\omega}})} -[\frac{K_{a_{\omega}}(\sqrt{\varpi_{\omega}})} {K_{a_{\omega}-1}(\sqrt{\varpi_{\omega}})}]^{2}\},$$ where $var(\bullet)$ denotes the variance.
Scaling $\xi_{\omega,jn}$, we have $$\overline{\phi_{\omega,jn}}=\frac{\overline{ \xi_{\omega,jn}}}{\sum_{j<k}\overline{ \xi_{\omega,jk}}},$$
and $$\overline{\phi_{\omega,jn}^{2}}=\overline{\phi_{\omega,jn}}^{2}+\frac{var(\xi_{\omega,jn})}{(\sum\limits_{j<k}\overline{ \xi_{\omega,jk}})^{2}}$$
Thus, the optimal $q$ density of $\phi_{\omega,ij}$ takes the following form:
We explain the technical details of second-stage VB in Online Appendix B. Below we briefly describe the priors of the parameters and hyperparameters and then provide their optimal $q$ densities.
Let $Z=(\widehat{Y}_{ /i} \hspace{0.2cm}X_{\bullet, i})'$ and $\theta=[(\Lambda _{i\bullet})' \hspace{0.2cm}\tilde{\beta}_{i}]$.
We elicit hierarchical D-L prior for $\theta$ as follows:
where $k=N-1+m_{i}$.
Next we set a Gamma prior for $\sigma^{-2}$:
The VB optimal densities can be found as follows:
where $$\overline{\tilde{V}}=(\frac{\frac{T}{2}+\nu}{\overline{\tilde{S}}}Z'Z+\tilde{V}^{-1})^{-1}$$ $$\overline{\theta}=(\frac{\frac{T}{2}+\nu}{\overline{\tilde{S}}})\overline{\tilde{V}}Z'y_{i}$$ $$\tilde{V}^{-1}=diag(\overline{\tilde{\psi}_{1}^{-1}}\hspace{0.1cm} \overline{\tilde{\phi}_{1}^{-2}}\overline{\tilde{\tau}^{-2}},...,(\overline{\tilde{\psi}_{k}^{-1}} \hspace{0.1cm} \overline{\tilde{\phi}_{k}^{-2}}\overline{\tilde{\tau}^{-2}})$$
where $$\overline{S}=\frac{1}{2}[||y_{i}-Z\overline{\theta}||^{^{2}}+tr(Z'Z\overline{\tilde{V}})]+\tilde{S}$$
Let $\tilde{\chi}=\sum_{j=1}^{k}2(\overline{\theta}_{j}^{2}+\overline{\tilde{V}}_{jj})^{1/2}\frac{1}{\overline{\tilde{\phi}_{j}}}$, we have $$\overline{\tilde{\tau}}=\frac{\sqrt{\tilde{\chi}}K_{k\tilde{a}-k+1}(\sqrt{\tilde{\chi}})}{K_{k\tilde{a}-k}(\sqrt{\tilde{\chi}})}$$ and $$\overline{\tilde{\tau}^{2}}=\overline{\tilde{\tau}}^{2}+\chi[\frac{K_{k\tilde{a}-k+2}(\sqrt{\tilde{\chi}})}{K_{k\tilde{a}-k}(\sqrt{\tilde{\chi}})}-(\frac{K_{k\tilde{a}-k+1}(\sqrt{\tilde{\chi}})}{K_{k\tilde{a}-k}(\sqrt{\tilde{\chi}})})^{2}]$$
Let $\tilde{\rho}=\sqrt{\frac{\overline{\tilde{\phi}_{j}^{2}}\hspace{0.2em} \overline{\tilde{\tau}^{2}}}{\overline{\theta}_{j}^{2}+\overline{\tilde{V}}^{jj}}}$, $$\overline{\frac{1}{\tilde{\psi}_{j}}}=\tilde{\rho}$$ and $$\overline{\tilde{\psi}}_{j}=1+1/\tilde{\rho}$$
Let $\tilde{\varpi}=2\sqrt{\overline{\theta}_{j}^{2}+(\overline{\tilde{V}}_{jj})^{2}}$, we have $$ \overline{ \tilde{\xi}}_{j}=\frac{\sqrt{\tilde{\varpi}}K_{\tilde{a}} (\sqrt{\tilde{\varpi}})} {K_{\tilde{a}-1}(\sqrt{\tilde{\varpi}})}, $$ and $$var(\tilde{\xi}_{j})=\tilde{\varpi} \{\frac{K_{\tilde{a}+1}(\sqrt{\tilde{\varpi}})} {K_{\tilde{a}-1}(\sqrt{\tilde{\varpi}})} -[\frac{K_{\tilde{a}}(\sqrt{\tilde{\varpi}})} {K_{\tilde{a}-1}(\sqrt{\tilde{\varpi}})}]^{2}\}.$$
Scaling $\tilde{\xi}$, we have $$\overline{\tilde{\phi}}_{j}=\frac{\overline{ \tilde{\xi}}_{j}}{\sum^{k}\overline{ \tilde{\xi}}_{j}},$$ and $$\overline{\tilde{\phi}_{j}^{2}}=\overline{\tilde{\phi}}_{j}^{2}+\frac{var(\tilde{\xi}_{j})}{(\sum^{k}\overline{ \tilde{\xi}}_{j})^{2}}$$
Thus, the optimal $q$ density of $\tilde{\phi}_{j}$ takes the following form:
In the Monte Carlo studies, we look into two traditional panel SAR models of various sample sizes. The first model is:
where to specify $W_{N}$, we let each cross-sectional unit be connected with the unit ahead of it and the unit behind, and then normalize $W$ by rows. When conducting Monte Carlo, we generate $x_{t}$ and $u_{t}$ independently from $N(0,1)$ and $0.1N(0,1)$, respectively.
The second model is:
where $y_{t,1}$ is a vector of half of the $N$ dependent variables observed at time $t$, and $y_{t,2}$ is the other half. When setting $W_{N_{1}}$ and $W_{N_{2}}$, we assume a variable in $y_{t,1}$ is only spatially related with the unit ahead of it and the unit behind, likewise a variable in $y_{t,2}$. Both $W_{N_{1}}$ and $W_{N_{2}}$ are normalized by rows. In Monte Carlo, we generate each element of $x_{t,1}$ and $x_{t,2}$ independently from $N(0,1)$, then each element of $u_{t,1}$ and $u_{t,2}$ independently from $0.1N(0,1)$.\footnote{We have also experimented on $W_{N}$, $W_{N_{1}}$ and $W_{N_{2}}$ of other forms. In addition, we have looked into models with different coefficients, $u_{t,1}\sim 0.1N(0,1)$ and $u_{t,2}\sim 0.3N(0,1)$. The Monte Carlo results provide further evidence that two-stage VB method works well.}
For both ((ref)) and ((ref)), the sample sizes considered are $N=30\hspace{0.1cm}$ and $\hspace{0.1cm} T=20$, $N=30\hspace{0.1cm}$ and $\hspace{0.1cm} T=80$, $N=50\hspace{0.1cm}$ and $\hspace{0.1cm} T=30$, $N=50\hspace{0.1cm}$ and $\hspace{0.1cm} T=100$, $N=100\hspace{0.1cm}$ and $\hspace{0.1cm} T=50$, and $N=100\hspace{0.1cm}$ and $\hspace{0.1cm} T=200$. For each case, we conduct 1000 Monte Carlo replications and use the changes in parameters instead of that of ELBO to check whether two-stage VB has converged.
Results of Monte Carlo simulations, which are relegated to Online Supplement to save space, provide strong evidence that the two-stage VB is able to recover the true parameters in the data generating process, especially when $T\gg N$. When $T\ll N$, two-stage VB estimates have larger biases and larger empirical standard deviations, which can be reduced by setting tighter priors. More important, there is clear evidence that the true spatial connections can be identified regardless of the length of the panel, long or short.
To give a flavour, tables 1-2 report the Monte Carlo results of parameter $\Lambda$ and $\tilde{\beta}$ for model ((ref)) with $N=30$ and $T=80$ as well as $N=30$ and $T=20$. Since a tablet with 30 columns is too big to fit into a page, we only report the mean and standard deviations, the latter in parenthesis, of the empirical distributions of the $5\times 5$ sub-matrices in the four corners of $\Lambda$, which are associated with $y_{t,1},...,y_{t,5}$ and $y_{t,26},...,y_{t,30}$, and the first and last five elements in $\tilde{\beta}$.
We use high-performance computing (HPC) services for estimation. All computations are done using 1 compute node and 4 processor core. On average, each Monte Carlo replication for model ((ref)) where $N=30$ and $T=80$ takes about 20 seconds, while for the same model, when $N=30$ and $T=20$ it takes about 60 seconds. For model ((ref)), it again takes 20 seconds when $N=30$ and $T=80$. However, estimating model ((ref)) takes 3 minutes when $N=30$ and $T=20$. The same pattern can be observed when estimating models where $N=50$ and $N=100$, that is: two-stage VB takes much longer to converge when $N\gg T$, and the larger the number of nonzero true parameters, the longer the estimation takes. To give a flavour of how fast two-stage VB converges, we would like to mention that, when allowing for parallel computing, each Monte Carlo replication of model ((ref)) where $N=100$ and $T=50$ takes about 20 minutes, and that number reduces to 5 minutes for the same model where $N=100$ and $T=200$.
The relationship between sovereign bond ratings and spreads in eurozone countries is of great interest to researchers and policy makers alike, as explained in papers such as Gibson et al. (2021) and Hall et al. (2022). One of the important questions posed is whether rate changes of the southern euro area countries and their northern counterparts have different impacts on a member country, especially during the European debt crisis, also known as euro area crises, that started in late 2009 and lingered on well into 2015. Equally, it is important to know whether rate changes of a country group affect all the other countries differently. To answer those questions, we apply the two-stage VB approach on a sample containing five southern euro area countries -- Spain, Greece, Ireland, Italy and Portugal -- and five northern euro area countries -- Austria, Belgium, France and Netherlands.
We use monthly data of Gibson et al. (2021). Sovereign rating is denoted by Rate, which is the combined sovereign ratings given by Standard and Poor's, Fitch and Moody's then transformed into numerical scale with `triple A' having value 1 and `selected default' having the value 22. Spread is defined as the difference between the yield on 10-year government bond of a country and that of Germany. Hence a rise in the rate or spread implies a worsening situation. Following Gibson et al. (2017), we use the following variables that capture the economic and political fundamentals as control variables: GDPgrowth, which is the real GDP growth rate; News, which is fiscal news constructed using European Commision forecasts; $\frac{Debt}{GDP}$, which is the ratio of government debt to GDP; $\frac{CA}{GDP}$, which is the ratio of current account balance to GDP; $\frac{P}{P^{*}}$, which is the ratio of a country's harmonised consumer price to that of Germany; and Pol, which is the index of political uncertainty reflecting the climate for foreign investors and political uncertainty. The monthly data runs from January 2000 to April 2019. For brevity, we report the data sources in the Appendix.
Taking account of the feedback loop between sovereign bond ratings and sovereign spreads, our dynamic model takes the following form:
where $W_{i}^{rate}$ and $W_{i}^{spread}$ are the $i^{th} $ row of the $10\times 10$ spatial weights matrices $W^{rate}$ and $W^{spread}$, respectively. Note that $W^{rate}$ and $W^{spread}$ are both with zero diagonals and the rest of the elements to be estimated.
Model ((ref)) can be rewritten as
where $Rate_{t}=(Rate_{1t}, ...,Rate_{10t})'$, $Spread_{t}=(Spread_{1t}, ...,Spread_{10t})'$, $ \varepsilon_{t}^{rate}=( \varepsilon_{1t}^{rate},..., \varepsilon_{10t}^{rate})'$, $\varepsilon_{t}^{spread}=(\varepsilon_{1t}^{spread},...,\varepsilon_{10t}^{spread})'$, $\mathbf{c}_{12}=diag(\mathbf{c}_{12,1},..., \mathbf{c}_{12,10})$, $\mathbf{c}_{22}=diag(\mathbf{c}_{22,1},..., \mathbf{c}_{22,10})$, $\mathbf{c}_{17}=diag(\mathbf{c}_{17,1},..., \mathbf{c}_{17,10})$, $\mathbf{c}_{27}=diag(\mathbf{c}_{27,1},..., \mathbf{c}_{27,10})$, $\mathbf{x}$ is the vector containing all the exogenous variables, and $A$ is the matrix containing their respective parameters. In the spirit of Debarsy et al. (2012), we use equation ((ref)) to calculate the impulse response functions and evaluate how rate and spread changes in one country spill over to the other countries.
We use a rolling window of 24 months to trace how the spillovers fluctuate over time. For the $i^{th}$ country, the average rating spillovers from the south is computed by taking the mean of the cumulative impacts of a 1 notch value increase in ratings of southern euro area countries; and the average rating spillovers from the north is computed by the average of the cumulative impacts of a 1 notch value increase in northern euro area countries. In both cases, the impact responses to a shock of own country are excluded. In the same fashion, we calculate the average spread spillovers from the south and the north. We set the amount of shock to the spread to be 1 basis point.
Figure 1 presents the average spillovers from the southern euro area countries and those from the north, where spillovers are measured by the cumulative impacts of 5 years.\footnote{The impulse response functions usually converge to zero after 2 to 3 months} Two salient patterns leap out. First, spatial spillovers between countries only become apparent after the financial crises in 2007. Before that, the spillovers between countries are barely minimum. Second, shocks to a country's spread tend to always increase other countries' spreads, but the impacts of shocks to sovereign ratings are mixed.
Let us start with the impacts of rating changes in a northern euro area country. We can observe that shocks to rates of the southern countries tends to improve the sovereign ratings of France, Finland and Netherland (causing the values of their rates to decrease), but worsen the ratings of Austria and Belgium (causing their rates to increase), especially during the period of euro area crises. For a northern country, shocks to rates of other northern countries, by contrast, tend to slightly worsen its sovereign ratings.
For a southern country, apart from Greece and Ireland, shocks to rates of the northern countries tend to have little cumulative impacts on its ratings. For Greece and Ireland, especially in the euro area crises, positive shocks to northern countries' ratings will cause their ratings to deteriorate as well. For all the southern countries, shocks that worsen other southern countries' ratings will worsen their own ratings too, in a magnitude that is much higher than the impacts of the same amount of shocks to northern countries' ratings.
For all countries, an increase in the spreads of northern countries is always followed by an increase in their own spreads, making it more costly for the country to borrow, same as the impacts of an increase in the spreads of southern countries. But the latter's impacts are much larger than the former's in all cases. In euro area crises, the spreads of France, Finland and Spain increased a lot after shocks to southern countries' spreads.
Since model ((ref)) is dynamic, the spatial spillovers measured by impulse responses are complicated by time dependence. To better understand the spatial structure, we next plot the average spatial weights associated with each country in Figure 2. For the $i^{th}$ country, the average rating spatial weights from the south is computed by taking the mean of the non-zero elements in $W_{i}^{rate}$ that are associated with the southern euro area countries; and the average rating spatial weights from the north is computed by taking the mean of the non-zero elements in $W_{i}^{rate}$ that are associated with the northern euro area countries. In the same way, we calculate the average spread spatial weights from the south and the north. Plots in Figure 2 provides further evidence for the marked differences in how a country is spatially linked with the north and the south. In general, northern countries are less spatially influenced by other countries than the southern countries, warranting we taking a closer look at how a southern country's rate and spread are spatially related to those of other countries.
In terms of the spatial relationships between one country's ratings and the ratings of other countries, Greece and Ireland stand out during the crises period, with other countries' ratings spatially influencing these two countries' ratings negatively, highlighting the limitations of traditional spatial weights matrix where all the elements are non-negative. Apart from Greece and Ireland, a country's sovereign ratings tend to be positively associated with other countries's ratings. Interestingly, ratings of Spain and Portugal are more closely associated with those of other southern countries than those of the northern countries from 2007 to 2010, but become more influenced by northern countries' ratings in the euro area crises after 2010. Italy's sovereign ratings, by contrast, are more spatially influenced by ratings of the northern countries than those of the south throughout time.
For spreads, negative weights are rare, indicating an increase in one country's spread tend to be associated with spread increases in other countries too. Apart from Spain, whose spreads are more spatially influenced by the northern countries for a brief period during euro area crises, spreads of all southern countries are more closely spatially linked with other southern countries than with the northern ones.
In applied work, if the spatial weights matrix set a priori were to be far from its true value, the empirical analysis using SAR models would be misleading. In this paper, we have developed a two-stage VB approach to estimating panel SAR models with unknown spatial weights matrices so as to let the data speak. Our two-stage VB with D-L priors method can be easily extended by using other popular priors such as the Lasso, Horseshoe and spike and a slab priors. Furthermore, the success of two-stage VB shows the potential of combining VB with more sophisticated methods, such as three-stage least squares, full information likelihood and GMM to estimate panel SAR models, especially those involving $N\gg T$. Monte Carlo experiments show that our two-stage VB is rather fast, and it can recover the spatial impacts well for both the long and short panels.
As an empirical example, we apply the two-stage VB to the sovereign bond ratings and sovereign spreads data of 10 eurozone countries to uncover the impacts of spatial spillovers of the southern countries, which were more severely hit by the eurozone debt crises, and those of the northern countries. Without pre-imposing any spatial weights matrices that might be unrealistic, we are able to shed new lights on the spillover behaviours of the south and the north. To our best knowledge, our findings are among the first in the literature to delineate how an individual country is affected by the spillovers from other countries in the eurozone.