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.
177,701 characters · 23 sections · 97 citation commands
A Review of Cross-Sectional Matrix Exponential Spatial Models
\def\spacingset#1{ {#1}} \spacingset{1}
{\it Keywords:} Matrix exponential spatial specification, MESS, spatial autoregression, SAR, heteroskedasticity, Bayesian estimation, model selection, impact measures.
\spacingset{1.45}
Spatial econometric models deal with estimation and inference problems that arise from (weak) cross-sectional dependence or correlation in data marked with location stamps. The spatial autoregressive (SAR) specification has been a widely used approach for modeling spatial dependence since its inception in Whittle:1954 and Cliff:1969,Cliff:1973. The matrix exponential spatial specification (MESS) was introduced by Lesage:2007 as an alternative to the SAR specification, primarily due to its computationally appealing properties in likelihood-based estimation schemes. Despite various estimation and inference methods proposed in the econometrics literature for models using either specification, the empirical literature is predominantly populated with papers utilizing the SAR specification.
This highly skewed preference towards the SAR specification by applied researchers is unfortunate in the sense that the MESS attains some attractive properties. First, we must emphasize that the MESS and the SAR imply different rates of decay for cross-sectional dependence. While it is a geometric rate in the case of SAR, the MESS implies an exponential rate of decay for spatial correlation. Consequently, they also imply different reduced forms for the cross-sectional models. Second, contrary to the SAR specification, the MESS does not require any restrictions on the parameter space of the spatial autoregressive parameters as the reduced form of the MESS always exists. In particular, the MESS always yields a positive definite covariance matrix for the outcome variable. Third, in the likelihood based estimation, the SAR specification results in Jacobian determinant terms that can be difficult to compute when the number of cross-sections is large. The MESS likelihood, on the other hand, does not involve such Jacobian terms.
In this paper, our aim is to provide a complete comprehensive review of the econometric literature on the estimation, inference and model selection methods for the cross-sectional MESS-type models.\footnote{We focus on cross-sectional MESS models as there are relatively few papers on panel data MESS models in the literature. See, e.g., Chih:2018, Zhang2019 and yangunbalanced.} More specifically, we aim to present the existing results from the literature in a more accessible way so that they can be utilized easily in empirical applications by applied researchers. Furthermore, we extend the existing literature in some important aspects. Firstly, we propose a new estimation and inference methodology for the cross-sectional MESS-type models with an unknown form of heteroskedasticity. Second, for the model selection problems involving cross-sectional MESS-type models, we consider a new method for computing the marginal likelihoods of the competing models in a Bayesian framework. In an Monte Carlo study, we also assess the finite sample properties of some existing estimators from the literature along with our proposed method. The simulation results show that the suggested M-estimator performs satisfactorily in finite samples.
Various estimation methods for the MESS models have been considered in the literature Lesage:2007, Jin:2015, yangfast, yangunbalanced. Lesage:2007 consider both the maximum likelihood and Bayesian estimation approaches. Jin:2015 formally investigate the large sample properties of the quasi maximum likelihood estimator (QMLE) and the generalized method of moments estimator (GMME). Although both estimators have the standard large sample properties, the GMME can be more efficient than the QMLE when the innovations are non-normal or heteroskedastic. Jin:2015 showed that, unlike the SAR-type models, in the presence of an unknown form of heteroskedasticity, the QMLE of the MESS model can remain consistent if the spatial weights matrices used in the model are commutative. In this paper, we extend on their results by introducing an M-estimation methodology that is robust to heteroskedasticity when the spatial weights matrices do not commute. We also formally establish the large sample properties of the resulting M-estimator.
We provide Bayesian estimation algorithms for the MESS models in the case of both homoskedastic and heteroskedastic error terms Dogan:2023, yangfast, Lesage:2007. In the case of heteroskedasticity, we assume that the error terms follow a scale mixture of normal distributions, where the latent scale variables generate distributions with different variances. The latent variable representation facilitates the estimation through the data augmentation techniques. In both homoskedastic and heteroskedastic models, the conditional posterior distributions of parameters take known forms, except for those of the spatial parameters. The posterior draws for the spatial parameters can be generated either by using the random-walk or the independence-chain Metropolis-Hastings algorithms Lesage:2009,yangfast, LeSage:2007b, Lesage:1997. We also consider the estimation of the MESS with endogenous and Durbin's regressors. Jin.Lee2018 show that the popular nonlinear two stage least squares (N2SLS) estimator, although consistent, may suffer from slow rates of convergence, and may attain nonstandard asymptotic distributions when the true value of a subset of model parameters is zero. We highlight how an adaptive group lasso estimator can provide a solution to these problems.
One important issue for the estimation of MESS-type model relates to the computation of the matrix exponential terms. Although there are various methods suggested in the literature, there is no single method that outperforms the rest in all cases Moler:2003. As such, we visit the computation of the matrix exponential terms and exhibit how the matrix-vector product approach originally suggested by Lesage:2007 can be utilized for quick computation of these terms. A further issue for the MESS models relates to the interpretation of the coefficient estimates for the explanatory variables. In spatial models, the interpretation of the coefficient estimates for the explanatory variables become more complicated due to the cross-sectional interactions. To this end, we review the existing results on the estimation and inference results for the impact measures for the MESS-type models Jin.Lee2018, Lesage:2009, Arbia:2020.
When modeling spatial dependence, researchers may encounter specification problems related to choosing a spatial weights matrix from a pool of candidates, or choosing between nested or non-nested alternative model specifications. Often, modeling is done in an ad hoc manner and there is no guidance from an underlying structural model to address these issues. To this end, we provide a complete review of the literature on testing based, criterion based and marginal likelihood based approaches for model selection problems involving MESS-type models HAN2013250, LIU2019434, Dogan:2023, yangms, Lesage:2009. In this regard, we also visit the Bayesian approaches and consider the modified harmonic mean method of Dey:1994 for the computation of the marginal likelihoods of competing models.
The rest of this paper is organized as follows. Section (ref) reviews several cross-sectional MESS-type models. This section also introduces the main properties of a matrix exponential term. Sections (ref)--(ref) discuss various estimation and inference techniques for the MESS-type models. Section (ref) details the matrix-vector product approach for the efficient computation of matrix exponential terms. Section (ref) presents the impact measures for the MESS-type models and illustrates the inference methods. Section (ref) considers various kinds of model selection approaches involving the MESS-type models. Section (ref) presents results of an Monte Carlo study, focusing on the M-estimation of the MESS-type models. Section (ref) ends our review with some concluding remarks for future research. Some technical results are relegated to an appendix.
We consider the following first order matrix exponential spatial model (for short MESS$(1,1)$)
where $\mathbf{Y}=\left(y_{1},\hdots,y_{n}\right)^{'}$ is the $n\times1$ vector of observations on the dependent variable, $\mathbf{X}$ is the $n\times k$ matrix of non-stochastic exogenous variables with the associated parameter vector $\boldsymbol{\beta}_0$, $\mathbf{U}=\left(u_{1},\hdots,u_{n}\right)^{'}$ is the $n\times1$ vector of regression error terms, and $\mathbf{V}=\left(v_{1},\hdots,v_{n}\right)^{'}$ is the $n\times1$ vector of idiosyncratic error terms. The matrix exponential term $e^{\lambda_0 \mathbf{W}}$ is defined by $e^{\lambda_0 \mathbf{W}}=\sum_{i=0}^\infty\frac{\lambda^{i}_0 \mathbf{W}^{i}}{i!}$, where $\mathbf{W}$ is an $n\times n$ spatial weights matrix with zero diagonal elements and $\lambda_0$ is a scalar spatial parameter. The matrix exponential $e^{\rho_0 \mathbf{M}}$ is defined in a similar way, where $\mathbf{M}$ is another $n\times n$ spatial weights matrix and $\rho_0$ is a scalar spatial parameter.
The MESS(1,1) in (ref) can be considered as the matrix exponential counterpart of the spatial autoregressive model with spatial autoregressive disturbances (SARAR(1,1)),
where $\alpha_0$ and $\tau_0$ are scalar spatial autoregressive parameters. The MESS(1,1) specification is obtained from (ref) by replacing $(\mathbf{I}_n - \alpha_0 \mathbf{W})$ and $(\mathbf{I}_n - \tau_0 \mathbf{M})$ with $e^{\lambda_0 \mathbf{W}}$ and $e^{\rho_0 \mathbf{M}}$, respectively. The matrix exponential terms satisfy the following properties Lesage:2007:
Because of these properties, the spatial models formulated with the matrix exponential terms have some advantages over the spatial models formulated with the spatial autoregressive processes. The first and second properties ensure that the reduced form of matrix exponential models always exists and does not require any restrictions for the spatial parameters. In the context of (ref), the reduced form can be expressed as
The third property implies that $\left|e^{\lambda \mathbf{W}}\right|=e^{\lambda \text{tr}(\mathbf{W})}=1$ because $\mathbf{W}$ has zero diagonal elements. This property ensures that the likelihood function of matrix exponential models is free of any Jacobian terms that need to be computed many times during estimation (see Section (ref) for the details). On the other hand, the likelihood functions of spatial models specified in terms of spatial autoregressive processes is not free of Jacobian terms. For example, the likelihood function of the SARAR(1,1) model includes $|\mathbf{I}_n - \tau \mathbf{M}|$ and $|\mathbf{I}_n - \alpha \mathbf{W}|$, which need to be computed in each iteration during the estimation process.
The MESS(1,1) specification nests two alternative specifications, namely, the MESS(1,0) and MESS(0,1), which can be obtained by setting $\lambda_0=0$ and $\rho_0=0$, respectively. A spatial Durbin extension can be obtained by including the spatial lags of the explanatory variables as regressors:
where $\mathbf{W}\mathbf{X}$ denotes the spatial lag of $\mathbf{X}$ and $\boldsymbol{\delta}_0$ is the corresponding vector of coefficients.
In the MESS(1,1) model, spatial interactions in the outcome variable arise only trough $\mathbf{W}$, and in the disturbance terms only through $\mathbf{M}$. In some cases, spatial dependence may arise from different sources, requiring different matrix exponential terms formulated with different spatial weights matrices. Let $\{\mathbf{W}_i\}_{i=1}^p$ and $\{\mathbf{M}_j\}_{j=1}^q$ be two sequences of spatial weights matrices. Then, a high-order version including the matrix exponential terms formulated with $\{\mathbf{W}_i\}_{i=1}^p$ and $\{\mathbf{M}_j\}_{j=1}^q$ can be specified as
where $\{\lambda_{i0}\}_{i=1}^p$ and $\{\rho_{j0}\}_{j=1}^q$ are sequences of spatial parameters. This model can be called the MESS$(p,q)$ model.
Finally, we specify the distribution of the elements of $\mathbf{V}$. We can consider both homoskedastic and heteroskedastic error terms as specified in the following assumptions.
Both assumptions require that the error terms have more than the fourth moment, which is required by the central limit theorem (CLT) considered by KP:2001, KP:2010 for the linear and quadratic forms of $\mathbf{V}$ (see Lemma 4 in the Appendix).
In this section, we will consider the quasi maximum likelihood estimation of the MESS(1,1) model under Assumption (ref). Let $\boldsymbol{\theta}=(\boldsymbol{\gamma}^{'},\sigma^2)^{'}$, $\boldsymbol{\gamma}=(\boldsymbol{\beta}^{'},\boldsymbol{\zeta}^{'})^{'}$ and $\boldsymbol{\zeta}=(\lambda,\rho)^{'}$. Also let $\boldsymbol{\theta}_0=(\boldsymbol{\gamma}_0^{'},\sigma_0^2)^{'}$ denote the true values of the parameters. Then, the quasi log-likelihood function for the MESS(1,1) is given by
Since $\ln\left|e^{\lambda \mathbf{W}}\right|=\ln(e^{\lambda\text{tr}(\mathbf{W})})=\ln1=0$ and $\ln\left|e^{\rho \mathbf{M}}\right|=\ln(e^{\rho\text{tr}(\mathbf{M})})=\ln1=0$, the two Jacobian terms disappear in the quasi log-likelihood function. Thus, the quasi log-likelihood function simplifies to
We can concentrate out $\sigma^2$ from the quasi log-likelihood function to obtain the concentrated quasi log-likelihood function only involving $\boldsymbol{\gamma}$. From the first order condition with respect to $\sigma^2$, the quasi maximum likelihood estimator of $\sigma^2$ is given by
Substituting (ref) into (ref), we obtain the concentrated quasi log-likelihood function as
Then, the QMLE $\hat{\boldsymbol{\gamma}}$ of $\gamma_0$ is defined as
which is equivalent to
where $Q(\boldsymbol{\gamma})=(e^{\lambda \mathbf{W}} \mathbf{Y}-\mathbf{X}\boldsymbol{\beta})^{'}\ermpe^{\rho \mathbf{M}}(e^{\lambda \mathbf{W}} \mathbf{Y}-\mathbf{X}\boldsymbol{\beta})$. Substituting $\hat{\boldsymbol{\gamma}}$ into (ref), we obtain the QMLE of $\sigma^2$ as $\hat{\sigma}^2=\hat{\sigma}^2(\hat{\boldsymbol{\gamma}})$.
The large sample properties of the QMLE $\hat{\boldsymbol{\gamma}}$ can be established under some regularity conditions. For consistency, the necessary conditions are identifiable uniqueness of $\boldsymbol{\gamma}_0$ and the uniform stochastic convergence of the quasi maximum likelihood function to its population counterpart white:1994. For asymptotic normality of $\hat{\boldsymbol{\gamma}}$, the CLT for linear and quadratic forms can be utilized KP:2001, KP:2010. The low level assumptions guaranteeing the large sample properties of $\hat{\boldsymbol{\gamma}}$ are (i) the existence of moments of disturbance terms up to the fourth moment, (ii) a manageable degree of spatial correlation, (iii) a compact parameter space for $\boldsymbol{\zeta}$, (iv) the non-singularity of certain matrices in large samples, and (v) certain restrictions to guarantee identification of $\boldsymbol{\gamma}_0$ in large samples.\footnote{See Jin:2015 for a complete formal list of these low level assumptions.}
The score functions with respect to the elements of $\boldsymbol{\gamma}$ are given by
where $\mathbf{V}(\boldsymbol{\gamma}) = e^{\rho \mathbf{M}}(e^{\lambda \mathbf{W}} \mathbf{Y}-\mathbf{X}\boldsymbol{\beta})$. Define $\mathbf{B}= \mathrm{Var}\left(\frac{1}{\sqrt{n}}\frac{\partial Q(\boldsymbol{\gamma}_0)}{\partial\boldsymbol{\gamma}}\right)$ and $\mathbf{A}=\mathrm{E}\left(-\frac{1}{n} \frac{\partial^2 Q\left(\boldsymbol{\gamma}_0\right)}{\partial \boldsymbol{\gamma} \partial \boldsymbol{\gamma}^{\prime}}\right)$. To introduce the closed-forms of $\mathbf{A}$ and $\mathbf{B}$, let $\mu_3=\operatorname*{E}(v_i^3)$, $\mu_4=\operatorname*{E}(v_i^4)$, $\mathbb{W}=e^{\rho_0 \mathbf{M}} \We^{-\rho_0 \mathbf{M}}$, $\mathbf{H}^s=\mathbf{H}+\mathbf{H}^{'}$ for any square matrix $\mathbf{H}$ and $\mathrm{vec}_D(\mathbf{H})$ be a vector containing the diagonal elements of $\mathbf{H}$. Then, using (ref) and Lemma 2 in Appendix (ref), we obtain
where the elements are defined as $\mathbf{A}_{\lambda\lambda}=\sigma_0^2 \text{tr}\left(\mathbb{W}^s \mathbb{W}^s\right)+2\left(\mathbb{W} e^{\rho_0 \mathbf{M}} \mathbf{X} \boldsymbol{\beta}_0\right)^{'}\left(\mathbb{W} e^{\rho_0 \mathbf{M}} \mathbf{X} \boldsymbol{\beta}_0\right)$ and $\mathbf{B}_{\rho\rho}=\left(\mu_4-3 \sigma_0^4\right) \mathrm{vec}_D^{\prime}\left(\mathbb{W}^s\right) \mathrm{vec}_D\left(\mathbb{W}^s\right)+4 \mu_3\left(\mathbb{W} e^{\rho_0 \mathbf{M}} \mathbf{X} \boldsymbol{\beta}_0\right)^{\prime} \mathrm{vec}_D\left(\mathbb{W}^s\right)$. The asymptotic distribution of $\hat{\boldsymbol{\gamma}}$ can be derived by applying the mean value theorem to $\frac{\partial Q(\hat{\boldsymbol{\gamma}})}{\partial \boldsymbol{\gamma}}$ around $\boldsymbol{\gamma}_0$. By the mean value theorem, we can write $\sqrt{n}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_0)=-\left(\frac{1}{n}\frac{\partial^2 Q(\bar{\boldsymbol{\gamma}})}{\partial\boldsymbol{\gamma}\partial\boldsymbol{\gamma}^{'}}\right)^{-1}\frac{1}{\sqrt{n}}\frac{\partial Q(\boldsymbol{\gamma}_0)}{\partial\boldsymbol{\gamma}}$, where $\bar{\boldsymbol{\gamma}}$ lies between $\hat{\boldsymbol{\gamma}}$ and $\boldsymbol{\gamma}_0$ elementwise. Then, the desired result follows by showing that $\frac{1}{n}\frac{\partial^2 Q(\bar{\boldsymbol{\gamma}})}{\partial\boldsymbol{\gamma}\partial\boldsymbol{\gamma}^{'}}-\frac{1}{n}\operatorname*{E}\left(\frac{\partial^2 Q(\boldsymbol{\gamma}_0)}{\partial\boldsymbol{\gamma}\partial\boldsymbol{\gamma}^{'}}\right)=o_p(1)$ and the asymptotic normality of $\frac{1}{\sqrt{n}}\frac{\partial Q(\boldsymbol{\gamma}_0)}{\partial\boldsymbol{\gamma}}$ by Lemma 4 in Appendix (ref). Thus, it follows that
Note that there are two cases that yield $\mathbf{B}=\sigma^2_0\mathbf{A}$. The first case arises when $\mathbf{W}$ and $\mathbf{M}$ commute. Under the commutative property, we have $\mathbb{W}=\mathbf{W}$ and $\mathrm{vec}_D(\mathbb{W}^s) = \mathbf{0}$, suggesting that $\mathbf{B}=\sigma^2_0\mathbf{A}$. The second case occurs when the disturbance terms are normally distributed. Under the normality, we have $\mu_4=3\sigma^4_0$ and $\mu_3=0$, yielding $\mathbf{B}=\sigma^2_0\mathbf{A}$. In either case, the result in (ref) reduces to $\sqrt{n}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_0)\xrightarrow{d}N(\mathbf{0},\lim_{n\rightarrow\infty} \sigma_0^2\mathbf{A}^{-1})$.
Finally, for inference, the plug-in estimators of $\mathbf{A}$ and $\mathbf{B}$ can be utilized. To that end, $\sigma_0^2$ can be consistently estimated by evaluating (ref) at $\hat{\boldsymbol{\gamma}}$, and $\mu_3$ and $\mu_4$ can be consistently estimated by their sample analogs using the residuals $\mathbf{V}(\hat{\boldsymbol{\gamma}})$. Thus, the standard error of $\hat{\boldsymbol{\gamma}}$ can be obtained as the square root of the diagonal elements of $\frac{1}{n}\mathbf{A}^{-1}(\hat{\boldsymbol{\gamma}})\mathbf{B}(\hat{\boldsymbol{\gamma}})\mathbf{A}^{-1}(\hat{\boldsymbol{\gamma}})$, where $\mathbf{A}(\hat{\boldsymbol{\gamma}})$ and $\mathbf{B}(\hat{\boldsymbol{\gamma}})$ are the plug-in estimators of $\mathbf{A}$ and $\mathbf{B}$, respectively.
In this subsection, we consider the quasi maximum likelihood estimation of the MESS(1,1) under the assumption of heteroskedastic disturbance terms. Let $\boldsymbol{\Sigma}$ be the variance covariance matrix of the disturbance terms, i.e., $\boldsymbol{\Sigma}=\operatorname*{Diag}(\sigma_1^2,\hdots,\sigma_n^2)$, the diagonal matrix formed by $\sigma_i^2$'s. The score functions of the quasi likelihood function evaluated at $\boldsymbol{\gamma}_0$ are given by
The expectation of the score functions with respect to $\boldsymbol{\beta}$ and $\rho$ at $\boldsymbol{\gamma}_0$ are zero by Lemma 2 in Appendix (ref). However, the expectation of the score function with respect to $\lambda$ at $\boldsymbol{\gamma}_0$ is $\text{tr}(\mathbb{W}\boldsymbol{\Sigma})$. By Lemma 1 in Appendix (ref), the order of this term is $O(n)$ under the assumption that $\mathbf{W}$ and $\mathbf{M}$ are bounded in matrix column sum and row sum norms. Hence, the QMLE $\hat{\boldsymbol{\gamma}}$ may not be consistent. However, when $\mathbf{W}$ and $\mathbf{M}$ commute, we have $\mathbb{W}=\mathbf{W}$, yielding $\text{tr}(\mathbf{W}\boldsymbol{\Sigma})=0$. Hence, when $\mathbf{W}$ and $\mathbf{M}$ commute, the QMLE of MESS(1,1) may remain consistent under the assumption of heteroskedastic disturbance terms.
The consistency and asymptotic normality of the QMLE $\hat{\boldsymbol{\gamma}}$ can be proved similarly to the homoskedastic case. Let $\mathbf{D}=\mathrm{E}\left(-\frac{1}{n} \frac{\partial^2 Q\left(\boldsymbol{\gamma}_0\right)}{\partial \boldsymbol{\gamma} \partial \boldsymbol{\gamma}^{\prime}}\right)$ and $\mathbf{F}= \mathrm{Var}\left(\frac{1}{\sqrt{n}}\frac{\partial Q(\boldsymbol{\gamma}_0)}{\partial\boldsymbol{\gamma}}\right)$. Using Lemma 2 in Appendix (ref), we obtain
where
Then, it can be shown that
For inference, the standard error of $\hat{\boldsymbol{\gamma}}$ can be obtained as the square root of the diagonal elements of $\frac{1}{n}\mathbf{D}^{-1}(\hat{\boldsymbol{\gamma}})\mathbf{F}(\hat{\boldsymbol{\gamma}})\mathbf{D}^{-1}(\hat{\boldsymbol{\gamma}})$, where $\mathbf{D}(\hat{\boldsymbol{\gamma}})$ and $\mathbf{F}(\hat{\boldsymbol{\gamma}})$ are the plug-in estimators of $\mathbf{D}$ and $\mathbf{F}$, respectively. Also, note that $\mathbf{D}$ and $\mathbf{F}$ involve the unknown diagonal matrix $\boldsymbol{\Sigma}$. As in White:1980, the terms involving $\boldsymbol{\Sigma}$ can be consistently estimated by replacing $\boldsymbol{\Sigma}$ with $\hat{\boldsymbol{\Sigma}}=\operatorname*{Diag}\left(v_1^2(\hat{\boldsymbol{\gamma}}),\hdots,v_n^2(\hat{\boldsymbol{\gamma}})\right)$, where $v_i(\hat{\boldsymbol{\gamma}})$ is the $i$th element of $\mathbf{V}(\hat{\boldsymbol{\gamma}})$.
Under heteroskedasticity, when $\mathbf{W}$ and $\mathbf{M}$ do not commute, we can use the M-estimation method to formulate a consistent estimator of $\boldsymbol{\gamma}$ based on the adjusted score functions. Denote $\mathbf{V}(\boldsymbol{\beta},\boldsymbol{\zeta})=e^{\rho \mathbf{M}}(e^{\lambda \mathbf{W}} \mathbf{Y}-\mathbf{X}\boldsymbol{\beta})$ and $\mathbf{V}=\mathbf{V}(\boldsymbol{\beta}_0,\boldsymbol{\zeta}_0)$. Then, the score functions based on (ref) can be determined as\footnote{Under heteroskedasticity, we ignore $\sigma^2$ and aim to construct the adjusted score functions such that $\operatorname*{E}\left(S(\boldsymbol{\beta}_0,\boldsymbol{\zeta}_0)\right)=0$. }
The essential reason why the QMLE is not consistent is $\mathrm{plim}_{n\rightarrow\infty}\frac{1}{n}S(\boldsymbol{\gamma}_0)\ne0$. In the case of the score functions with respect to $\boldsymbol{\beta}$ and $\rho$, we have $\operatorname*{E}(\mathbf{X}^{'}e^{\rho_0 \mathbf{M}'}\mathbf{V})=0$, and $\operatorname*{E}(\mathbf{V}^{'}\mathbf{M}\mathbf{V})=\text{tr}(\boldsymbol{\Sigma}\mathbf{M})=0$ because $\boldsymbol{\Sigma}$ is a diagonal matrix and $\mathbf{M}$ has zero diagonal elements. In the case of the score function with respect to $\lambda$, we have
If $\mathbf{W}$ and $\mathbf{M}$ commute, i.e., $\mathbf{W}\mathbf{M}=\mathbf{M}\mathbf{W}$, then we have $\mathbb{W}=\mathbf{W}$, which yields $\text{tr}(\boldsymbol{\Sigma}\mathbb{W})=\text{tr}(\boldsymbol{\Sigma}\mathbf{W})=0$. Thus, when the commutative property holds, we have $\mathrm{plim}_{n\rightarrow\infty}\frac{1}{n}S(\boldsymbol{\gamma}_0)=0$, suggesting that the QMLE can be consistent under heteroskedasticity. However, if $\mathbf{W}\mathbf{M}\ne\mathbf{M}\mathbf{W}$, then we have $\text{tr}(\boldsymbol{\Sigma}\mathbb{W})=O(n)$ and $\mathrm{plim}_{n\rightarrow\infty}\frac{1}{n}S(\boldsymbol{\gamma}_0)\ne0$ in general, indicating that the QMLE may not be consistent under heteroskedasticity. We will adjust the score function with respect to $\lambda$ so that $\mathrm{plim}_{n\rightarrow\infty}\frac{1}{n}S(\boldsymbol{\gamma}_0)=0$ holds in all cases.
To adjust the score function with respect to $\lambda$, we use the trace property $\text{tr}(\mathbf{D}\mathbf{A})=\text{tr}(\mathbf{D}\operatorname*{Diag}(\mathbf{A}))$, where $\mathbf{D}$ is an $n\times n$ diagonal matrix and $\mathbf{A}$ is a conformable matrix. Using this property, we can express $\operatorname*{E}(\mathbf{Y}^{'}e^{\lambda_0 \mathbf{W}'}\mathbf{W}^{'}e^{\rho_0 \mathbf{M}'}\mathbf{V})$ as
Then, subtracting the last term from the second term in (ref), we obtain
where $\mathbb{W}_D=\mathbb{W}-\operatorname*{Diag}(\mathbb{W})$. Thus, we suggest using the sample counter part of $\operatorname*{E}(\mathbf{Y}^{'}e^{\lambda_0 \mathbf{W}'}e^{\rho_0 \mathbf{M}'}\mathbb{W}_D\mathbf{V})$ as the adjusted score function with respect to $\lambda$. Then, our suggested adjusted score functions take the following form:
where $\mathbb{W}_D(\rho)=\mathbb{W}(\rho)-\operatorname*{Diag}(\mathbb{W}(\rho))$ and $\mathbb{W}(\rho)=e^{\rho \mathbf{M}}\We^{-\rho \mathbf{M}}$. Note that $\operatorname*{E}\left(S^{*}(\boldsymbol{\gamma}_0)\right)=0$ holds by construction. We first derive the estimator of $\boldsymbol{\beta}_0$ for a given $\boldsymbol{\zeta}$ value, which is given by
Then, substituting $\hat{\boldsymbol{\beta}}_M(\boldsymbol{\zeta})$ into the $\lambda$ and $\rho$ elements of (ref), we obtain the concentrated score functions as
where $\hat{\mathbf{V}}(\boldsymbol{\zeta})=\mathbf{V}(\hat{\boldsymbol{\beta}}_M(\boldsymbol{\zeta}),\boldsymbol{\zeta})$. Then, the M-estimator (ME) of $\boldsymbol{\zeta}_0$ is defined by
Substituting $\hat{\boldsymbol{\zeta}}_M$ into (ref), we get the M-estimator for $\boldsymbol{\beta}$ as $\hat{\boldsymbol{\beta}}_M=\hat{\boldsymbol{\beta}}_M(\hat{\boldsymbol{\zeta}}_M)$. To prove the consistency of $\hat{\boldsymbol{\gamma}}_M=(\hat{\boldsymbol{\beta}}_M^{'},\hat{\boldsymbol{\zeta}}_M^{'})^{'}$, we only need to prove the consistency of $\hat{\boldsymbol{\zeta}}_M$ since $\hat{\boldsymbol{\beta}}_M=\hat{\boldsymbol{\beta}}_M(\hat{\boldsymbol{\zeta}}_M)$. To that end, we let $\bar{S}^{*}(\boldsymbol{\beta},\boldsymbol{\zeta})=\operatorname*{E}(S^{*}(\boldsymbol{\beta},\boldsymbol{\zeta}))$ be the population counterpart of the adjusted score functions in (ref). Given $\boldsymbol{\zeta}$, we can write $\bar{\boldsymbol{\beta}}_M(\boldsymbol{\zeta})=(\mathbf{X}^{'}\ermpe^{\rho \mathbf{M}}\mathbf{X})^{-1}\mathbf{X}^{'}e^{\rho \mathbf{M}'}\erme^{\lambda \mathbf{W}}\operatorname*{E}(\mathbf{Y})$, which can be substituted into the $\lambda$ and $\rho$ elements of $\bar{S}^{*}(\boldsymbol{\beta},\boldsymbol{\zeta})$ to obtain
where $\bar{\mathbf{V}}(\boldsymbol{\zeta})=\mathbf{V}(\bar{\boldsymbol{\beta}}_M(\boldsymbol{\zeta}),\boldsymbol{\zeta})$. To investigate the asymptotic properties of $\hat{\boldsymbol{\zeta}}_M$, we maintain the following assumptions.
Assumption (ref) provides the essential properties of the spatial weights matrices. It ensures that the spatial correlation is limited to a manageable degree KP:2001,KP:2010. Assumption (ref) requires that the parameter space of the parameters in the matrix exponential terms is compact. Assumption (ref) and Assumption (ref) imply that the matrix exponential terms are uniformly bounded in both row sum and column sum matrix norms. This can be seen from $\left\Vert e^{\lambda \mathbf{W}}\right\Vert=\left\Vert \sum_{i=0}^{\infty}\lambda^i\mathbf{W}^i/i!\right\Vert\leq\sum_{i=0}^{\infty}|\lambda|^{i}\Vert \mathbf{W}\Vert^{i}/i!=e^{|\lambda|\Vert \mathbf{W}\Vert}$, which is bounded if $|\lambda|$ and $\Vert \mathbf{W}\Vert$ are bounded, where $\Vert\cdot\Vert$ is either the row sum or the column sum matrix norm. Assumption (ref) provides some regularity conditions and corresponds to Assumption 4 of Jin:2015. Assumption (ref) is a high-level assumption and ensures the identification of $\boldsymbol{\zeta}_0$. In Appendix C, we provide two low-level conditions that are sufficient for Assumption (ref).
The uniform convergence $\sup_{\boldsymbol{\zeta}\in\boldsymbol{\Delta}}\frac{1}{n}\left\Vert S^{* c}(\boldsymbol{\zeta})-\bar{S}^{* c}(\boldsymbol{\zeta})\right\Vert\xrightarrow{\enskip p\enskip}0$ and Assumption (ref) ensure the consistency of $\hat{\boldsymbol{\zeta}}_M$.
To derive the asymptotic distribution of $\hat{\boldsymbol{\gamma}}_M$, we apply the mean value theorem to $S^{*}(\hat{\boldsymbol{\gamma}}_M)=0$ at $\boldsymbol{\gamma}_0$, to obtain $\sqrt{n}(\hat{\boldsymbol{\gamma}}_M-\boldsymbol{\gamma}_0)=-\left(\frac{1}{n}\frac{\partial S^{*}(\overline{\boldsymbol{\gamma}})}{\partial \boldsymbol{\gamma}^{'}}\right)^{-1}\frac{1}{\sqrt{n}}S^{*}(\boldsymbol{\gamma}_0)$, where $\overline{\boldsymbol{\gamma}}$ lies between $\boldsymbol{\gamma}_0$ and $\hat{\boldsymbol{\gamma}}_M$ elementwise. By substituting the reduced form $\mathbf{Y}=e^{-\lambda_0 \mathbf{W}}\left(\mathbf{X}\boldsymbol{\beta}_0+e^{-\rho_0 \mathbf{M}}\mathbf{V}\right)$ into $S^{*}(\boldsymbol{\gamma}_0)$, we obtain a linear-quadratic form in $\mathbf{V}$:
where $\mathbb{W}_D=\mathbb{W}_D(\rho_0)$. Thus, the CLT for the linear-quadratic forms of $\mathbf{V}$ in Lemma 4 of the Appendix can be used to establish the asymptotic normality of $\frac{1}{\sqrt{n}}S^{*}(\boldsymbol{\gamma}_0)$. Also, our assumptions ensure that $\frac{1}{n}\frac{\partial S^{*}(\overline{\boldsymbol{\gamma}})}{\partial \boldsymbol{\gamma}^{'}}-\frac{1}{n}\operatorname*{E}\left(\frac{\partial S^{*}(\boldsymbol{\gamma}_0)}{\partial \boldsymbol{\gamma}^{'}}\right)=o_p(1)$. Using these results, we determine the asymptotic distribution of $\hat{\boldsymbol{\gamma}}_M$ in Theorem (ref).
To conduct inference, we need consistent estimators of $\boldsymbol{\Psi}(\boldsymbol{\gamma}_0)$ and $\boldsymbol{\Omega}(\boldsymbol{\gamma}_0)$. For $\boldsymbol{\Psi}(\boldsymbol{\gamma}_0)$, we can use its observed counterpart given by $\boldsymbol{\Psi}(\hat{\boldsymbol{\gamma}}_M)=-\frac{1}{n}\frac{\partial S^*(\boldsymbol{\gamma})}{\partial \boldsymbol{\gamma}^{'}}|_{\boldsymbol{\gamma}=\hat{\boldsymbol{\gamma}}_M}$. The elements of $\boldsymbol{\Psi}(\boldsymbol{\gamma})$ are given by
where $\dot{\mathbb{W}}_D(\rho)=\frac{\partial \mathbb{W}_D(\rho)}{\partial \rho}=\mathbf{M}\mathbb{W}_D(\rho)-\mathbb{W}_D(\rho)\mathbf{M}-\operatorname*{Diag}\left(\mathbf{M}\mathbb{W}_D(\rho)-\mathbb{W}_D(\rho)\mathbf{M}\right)$ and $\mathbf{Y}(\boldsymbol{\zeta})=e^{\rho \mathbf{M}}\We^{\lambda \mathbf{W}}\mathbf{Y}$. In the proof of Theorem (ref), we show that $\boldsymbol{\Psi}(\hat{\boldsymbol{\gamma}}_M)$ is a consistent estimator of $\boldsymbol{\Psi}(\boldsymbol{\gamma}_0)$. Using Lemma 2 in the Appendix, we determined the closed form of $\boldsymbol{\Omega}(\boldsymbol{\gamma}_0)$ as
where $\boldsymbol{\Omega}_{22}=\boldsymbol{\beta}_0^{'}\mathbf{X}^{'}e^{\rho_0 \mathbf{M}'}\mathbb{W}_D\boldsymbol{\Sigma}\mathbb{W}_D^{'}e^{\rho_0 \mathbf{M}}\mathbf{X}\boldsymbol{\beta}_0+\text{tr}(\boldsymbol{\Sigma}\mathbb{W}_D\boldsymbol{\Sigma}\mathbb{W}_D^s)$. Let $\boldsymbol{\Omega}(\hat{\boldsymbol{\gamma}}_M)$ be the plug-in estimator of $\boldsymbol{\Omega}(\boldsymbol{\gamma}_0)$, where we replace $\boldsymbol{\Sigma}$ with $\hat{\boldsymbol{\Sigma}}=\operatorname*{Diag}\left(v_1^2(\hat{\boldsymbol{\gamma}}),\hdots,v_n^2(\hat{\boldsymbol{\gamma}})\right)$ and $v_i(\hat{\boldsymbol{\gamma}}_M)$ is the $i$th element of $\mathbf{V}(\hat{\boldsymbol{\gamma}}_M)$.
Thus, the standard error of $\hat{\boldsymbol{\gamma}}_M$ can be obtained as the square root of the diagonal elements of $\frac{1}{n}\boldsymbol{\Psi}^{-1}(\hat{\boldsymbol{\gamma}}_M)\boldsymbol{\Omega}(\hat{\boldsymbol{\gamma}}_M)\boldsymbol{\Psi}^{-1'}(\hat{\boldsymbol{\gamma}}_M)$.
In this section, we consider the GMM estimation of the MESS(1,1) model under Assumption (ref). Recall again from the definition of MESS(1,1) that $\mathbf{V}(\boldsymbol{\gamma})=e^{\rho \mathbf{M}}(e^{\lambda \mathbf{W}} \mathbf{Y}-\mathbf{X}\boldsymbol{\beta})$, where $\boldsymbol{\gamma}=(\boldsymbol{\zeta}^{'},\boldsymbol{\beta}^{'})^{'}$ and $\boldsymbol{\zeta}=(\lambda,\rho)^{'}$. We consider the following vector of moment functions consisting of $k_p$ quadratic moments and $k_f$ linear moments:
where $\mathbf{P}_{m}$'s are $n\times n$ matrices constants with $\text{tr}(\mathbf{P}_m)=0$ for $m=1,\hdots,k_p$, and $\mathbf{F}$ is the $n\times k_f$ matrix of instrumental variables (IV). Given an arbitrary symmetric weighting matrix $\boldsymbol{\Phi}$ with rank greater than or equal to $k+2$, the GMM objective function is given by $g^{'}(\boldsymbol{\gamma})\boldsymbol{\Phi} g(\boldsymbol{\gamma})$. Then, an initial GMME (IGMME) can be obtained by
Define $\mathbf{G}=\operatorname*{E}\left(\frac{\partial g(\boldsymbol{\gamma}_0)}{\partial\boldsymbol{\gamma}^{'}}\right)$ and $\mathbf{H}=n \mathrm{E}\left(g\left(\boldsymbol{\gamma}_0\right) g^{\prime}\left(\boldsymbol{\gamma}_0\right)\right)$. Let $\mathrm{vec}(\mathbf{A})$ denote the column vector formed by stacking the columns of matrix $\mathbf{A}$ and recall that $\mathrm{vec}_D(\mathbf{A})$ denotes the column vector formed by the diagonal elements of the matrix $\mathbf{A}$. Then, by Lemma 2 in Appendix (ref), we obtain
and
where $\boldsymbol{\omega}=(\mathrm{vec}(\mathbf{P}_1^s),\hdots,\mathrm{vec}(\mathbf{P}_{k_p}^s))$ and $\boldsymbol{\omega}_d=(\mathrm{vec}_D(\mathbf{P}_1^s),\hdots,\mathrm{vec}_D(\mathbf{P}_{k_p}^s))$.
The large sample properties of the IGMME $\hat{\boldsymbol{\gamma}}$ can be established under some regularity conditions. For consistency, the necessary conditions are identification of $\boldsymbol{\gamma}_0$ from the population moments and the uniform stochastic convergence of the generalized method of moments objective function to its population counterpart. For the asymptotic normality of $\hat{\boldsymbol{\gamma}}$, the central limit theorem for linear and quadratic forms can be utilized.\footnote{ The low level assumptions guaranteeing the large sample properties are provided in Jin:2015.} The asymptotic distribution of $\hat{\boldsymbol{\gamma}}$ can be derived by applying the mean value theorem to $\frac{\partial g^{'}(\hat{\boldsymbol{\gamma}})}{\partial \boldsymbol{\gamma}}\boldsymbol{\Phi} g(\hat{\boldsymbol{\gamma}})=0$ at $\boldsymbol{\gamma}_0$ to get $\sqrt{n}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_0)=-\left(\frac{\partial g^{'}(\hat{\boldsymbol{\gamma}})}{\partial\boldsymbol{\gamma}}\boldsymbol{\Phi} \frac{\partial g(\bar{\boldsymbol{\gamma}})}{\partial\boldsymbol{\gamma}^{'}}\right)^{-1}\frac{\partial g^{'}(\hat{\boldsymbol{\gamma}})}{\partial\boldsymbol{\gamma}}\boldsymbol{\Phi}\sqrt{n} g(\boldsymbol{\gamma}_0)$, where $\bar{\boldsymbol{\gamma}}$ lies between $\hat{\boldsymbol{\gamma}}$ and $\boldsymbol{\gamma}_0$ elementwise. Then, the asymptotic distribution of $\sqrt{n}(\hat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_0)$ follows by applying the CLT in Lemma 4 in the Appendix to $\sqrt{n} g(\boldsymbol{\gamma}_0)$ and showing that $\frac{\partial g(\hat{\boldsymbol{\gamma}})}{\partial\boldsymbol{\gamma}^{'}}-\operatorname*{E}\left(\frac{\partial g(\boldsymbol{\gamma}_0)}{\partial\boldsymbol{\gamma}^{'}}\right)=o_p(1)$. Thus, we have
From the expression for the variance-covariance matrix of IGMME, we can see that the precision of the estimator can be improved by replacing the arbitrary weighing matrix $\boldsymbol{\Phi}$ in the objective function with $\mathbf{H}^{-1}$. The resulting GMME is called the optimal GMME Hansen:1982. However, this estimator is not feasible as $\mathbf{H}^{-1}$ is unknown. To make it feasible, a plug-in estimator $\hat{\mathbf{H}}\equiv\mathbf{H}(\hat{\boldsymbol{\gamma}})$ based on the initial GMME $\hat{\boldsymbol{\gamma}}$ can be formulated. Then, the feasible optimal GMME is defined by
Under some conditions, Jin:2015 show that
Jin:2015 determine the best set of moment functions that provide the most efficient GMME for the MESS(1,1) under homoskedasticity. Their idea is to decompose the components of the inverse of the variance-covariance matrix of the optimal GMME, and then use the Cauchy-Schwarz inequality in such a way that an upper bound on the inverse of the variance-covariance matrix that is free of arbitrary pieces of the moment functions (free of $\mathbf{P}_i$'s and $\mathbf{F}$) can be attained. The resulting GMME is termed as the best GMME (BGMME). When the disturbance terms are normally distributed, the BGMME turns out to be asymptotically as efficient as the QMLE. However, when the disturbance terms are not normally distributed, and $\mathbf{W}$ and $\mathbf{M}$ do not commute, the BGMME can be asymptotically more efficient than the QMLE. The best set of moment functions is
where $\mathbf{P}^{*}_1=\mathbb{W}$, $\mathbf{P}^{*}_2=\operatorname*{Diag}(\mathbb{W})$, $\mathbf{P}^{*}_3=\operatorname*{Diag}(e^{\rho_0 \mathbf{M}}\mathbf{W}\mathbf{X}\boldsymbol{\beta}_0)^{(t)}$, $\mathbf{P}^{*}_4=\mathbf{M}$, $\mathbf{P}^{*}_{m+4}=\operatorname*{Diag}(e^{\rho_0 \mathbf{M}}\mathbf{X}_{m})^{(t)}$ for $m=1,...,k^{*}$, and $\mathbf{F}^{*}=(\mathbf{F}^{*}_1,\mathbf{F}^{*}_2,\mathbf{F}^{*}_3,\mathbf{F}^{*}_4)$ with $\mathbf{F}_1^{*}=e^{\rho_0 \mathbf{M}}\mathbf{X}^{*}$, $\mathbf{F}_2^{*}=e^{\rho_0 \mathbf{M}}\mathbf{W}\mathbf{X}\boldsymbol{\beta}_0$, $\mathbf{F}_3^{*}=\bm{l}$, $\mathbf{F}_4^{*}=\mathrm{vec}_D(\mathbb{W})$, where $\mathbf{X}^*$ excludes the intercept term in $\mathbf{X}$ if $\mathbf{M}$ is row-normalized so that $\mathbf{F}_1^{*}$ does not contain the intercept term generated in $e^{\rho_0 \mathbf{M}}\mathbf{X}$, $k^*$ is the number of columns in $\mathbf{X}^{*}$, $\mathbf{A}^{(t)}=\mathbf{A}-\mathbf{I}_n\text{tr}(\mathbf{A})/n$ for any $n\times n$ matrix $\mathbf{A}$ and $\bm{l}$ is an $n\times1$ vector of ones.
In this section, we consider the GMM estimation of MESS(1,1) under Assumption (ref). Recall that $\boldsymbol{\Sigma}$ denotes the variance-covariance matrix of the disturbance terms, i.e., $\boldsymbol{\Sigma}=\operatorname*{Diag}(\sigma_1^2,\hdots,\sigma_n^2)$. Similar to the homoskedastic case, we again employ the following vector of moment functions consisting of $k_p$ quadratic moment functions and $k_f$ linear moment functions:
At $\boldsymbol{\gamma}_0$, we have $\mathrm{E}\left(\mathbf{V}^{'}\mathbf{P}_m\mathbf{V}\right) = \text{tr}\left(\mathbf{P}_m\boldsymbol{\Sigma}\right)=\text{tr}\left(\boldsymbol{\Sigma}\operatorname*{Diag}(\mathbf{P}_i)\right)$, which is equal to zero if the diagonal elements of $\mathbf{P}_m$ are zeros. Hence, in the heteroskedastic case, we require that the diagonal elements of $\mathbf{P}_m$ are zeros, i.e., $\operatorname*{Diag}(\mathbf{P}_m)=\mathbf{0}$ for $m=1,2, \hdots, k_p$. Then, an initial GMME based on an arbitrary symmetric weighting matrix $\boldsymbol{\Phi}$, with rank greater than or equal to $k+2$, can be defined as
Let $\mathbf{G}=\operatorname*{E}\left(\frac{\partial g(\boldsymbol{\gamma}_0)}{\partial\boldsymbol{\gamma}^{'}}\right)$ and $\mathbf{H}=n \mathrm{E}\left(g\left(\boldsymbol{\gamma}_0\right) g^{\prime}\left(\boldsymbol{\gamma}_0\right)\right)$. By Lemma 2 in Appendix (ref), we can show that
and
where $ \boldsymbol{\omega}=\mathrm{vec}(\boldsymbol{\Sigma}^{1/2}\mathbf{P}_1^{s}\boldsymbol{\Sigma}^{1/2},\hdots,\boldsymbol{\Sigma}^{1/2}\mathbf{P}_{k_p}^{s}\boldsymbol{\Sigma}^{1/2})$. It follows again that
Note that $\mathbf{G}$ and $\mathbf{H}$ involve the unknown diagonal matrix $\boldsymbol{\Sigma}$. These terms can be consistently estimated by replacing $\boldsymbol{\Sigma}$ with $\operatorname*{Diag}(v_1^2(\hat{\boldsymbol{\gamma}}),\hdots,v_n^2(\hat{\boldsymbol{\gamma}}))$. Let $\hat{\mathbf{H}}$ be the plug-in estimator of $\mathbf{H}$ based on the initial GMME $\hat{\boldsymbol{\gamma}}$. Then, a feasible optimal robust GMME (RGMME) can be obtained as
It can be shown that
In the heteroskedastic case, the best set of moment functions is not feasible because the moment functions involve the unknown $\boldsymbol{\Sigma}$, which cannot be consistently estimated. In practice, we can formulate the RGMME based on the following vector of moment functions Jin:2015:
where $\mathbf{F}=(\hat{\mathbb{W}}e^{\hat{\tau}\mathbf{M}}\mathbf{X}\hat{\boldsymbol{\beta}},e^{\hat{\tau}\mathbf{M}}\mathbf{X})$ with $\hat{\mathbb{W}}=e^{\hat{\tau}\mathbf{M}} \mathbf{W} e^{-\hat{\tau}\mathbf{M}}$.
Following Lesage:2007, we assume the following independent prior distributions: $\lambda\sim N(\mu_\lambda ,V_{\lambda})$, $\rho\sim N(\mu_\rho,V_{\rho})$, $\boldsymbol{\beta} \sim N(\boldsymbol{\mu}_{\boldsymbol{\beta}}, \mathbf{V}_{\boldsymbol{\beta}})$ and $\sigma^2 \sim IG(a, b)$, where $IG$ denotes the inverse-gamma distribution. Under these prior distributions, the posterior distribution of parameters can be expressed as\footnote{We use $p(\cdot)$ to denote the relevant density functions, and ignore $\mathbf{X}$ in the conditional sets for the sake of simplicity.}
where $p(\boldsymbol{\theta})$ is the joint prior distribution of $\boldsymbol{\theta}$ and $p(\mathbf{Y}|\boldsymbol{\theta})$ is the likelihood function given as
Algorithm (ref) describes a Gibbs sampler that can be used to generate random draws from $p(\boldsymbol{\theta}|\mathbf{Y})$.
In Algorithm (ref), the conditional posterior distributions of $\boldsymbol{\beta}$ and $\sigma^2$ are determined from $p(\boldsymbol{\beta}|\mathbf{Y}, \lambda, \rho, \sigma^2)\propto p(\mathbf{Y}|\boldsymbol{\theta})p(\boldsymbol{\beta})$ and $p(\sigma^2 | \mathbf{Y}, \lambda, \rho, \boldsymbol{\beta})\propto p(\mathbf{Y}|\boldsymbol{\theta})p(\sigma^2)$, respectively. Since we assume conjugate priors for $\boldsymbol{\beta}$ and $\sigma^2$, these conditional posterior distributions take known forms as shown in Algorithm (ref). The Bayesian argument used to determine these conditional posterior distributions is analogous to the one used for a linear regression model. On the other hand, the conditional posterior distributions of spatial parameters are non-standard because the likelihood function is non-linear in terms of these parameters. To sample these parameters, we use the random walk Metropolis-Hastings algorithm suggested by Lesage:2009.
Following Lesage:1997 and Lesage:2009, we assume that the disturbance terms have a scale mixture of normal distributions such that the scale mixture variables generate different distributions with distinct variance terms. Thus, we have $v_i|\eta_i\sim N(0,\eta_i\sigma^2)$, where $\eta_i$'s are independent scale mixture variables with $\eta_i\sim \text{IG}(\nu/2,\nu/2)$ for $i=1,\hdots,n$. Let $\boldsymbol{\eta}=(\eta_1,\hdots,\eta_n)^{'}$ and $\mathbf{H}(\boldsymbol{\eta})=\text{Diag}\left(\eta_1,\hdots,\eta_n\right)$ be the $n\times n$ diagonal matrix with the $i$th diagonal element $\eta_i$. Then, we can derive the conditional likelihood function $p(\mathbf{Y}|\boldsymbol{\theta},\boldsymbol{\eta})$ as
To introduce a Bayesian estimation approach, we adopt the prior distributions assumed in Section 6.1 for $\boldsymbol{\beta}$, $\lambda$, $\rho$ and $\sigma^2$. In the heteroskedastic case, we also need to determine a prior distribution for $\nu$. To that end, we note that the marginal distribution of $v_i$ is a $t$ distribution with mean zero, scale parameter $\sigma^2$ and $\nu$ degrees of freedom, i.e., $v_i\sim t_\nu(0,\sigma^2)$. Thus, we assume the following prior $\nu\sim \text{Uniform}(2,\bar{\nu})$, where $\text{Uniform}(a,b)$ denotes the uniform distribution over the interval $(a,b)$, and $\bar{\nu}$ is a known positive number. This prior distribution ensures that the variance of $v_i$ exists because $\nu>2$. Also, we can set $\bar{\nu}$ to a large positive number so that the $t$ distribution is allowed to approximate the normal distribution well-enough.
The posterior distribution of parameters then takes the following form:
where $p(\mathbf{Y}|\boldsymbol{\theta},\boldsymbol{\eta})$ is the conditional likelihood function stated in (ref) and $p(\boldsymbol{\theta},\boldsymbol{\eta})$ is the joint prior distribution of $\boldsymbol{\theta}$ and $\boldsymbol{\eta}$. Algorithm (ref) describes a Gibbs sampler that can be used to generate random draws from $p(\boldsymbol{\theta},\boldsymbol{\eta}|\mathbf{Y})$.
The conditional posterior distributions of $\boldsymbol{\beta}$, $\boldsymbol{\eta}$, and $\sigma^2$ take known forms as shown in Algorithm (ref). In the case of spatial parameters, we again resort to the random walk Metropolis-Hastings algorithm suggested by Lesage:2009. The conditional posterior distribution of $\nu$ can be determined from $p(\nu|\mathbf{Y},\boldsymbol{\beta},\boldsymbol{\eta}, \lambda,\rho,\sigma^2)=p(\nu|\boldsymbol{\eta})\propto p(\boldsymbol{\eta}|\nu)p(\nu)$. However, this distribution does not take any known form. Since $\nu$ has support over $(2,\bar{\nu})$, we suggest using the Griddy-Gibbs sampler to sample this parameter. Algorithm (ref) describes the Griddy-Gibbs sampler.
The preceding sections consider a regression model with spatial dependence specified by the MESS processes, where no endogenous regressors are included. In this section, we consider a MESS model with endogenous and Durbin regressors. The popular nonlinear two-stage least squares (N2SLS) estimation method in such a setting can have some irregular features.
Consider the following model:
where $\bm{l}$ is an $n\times 1$ vector of ones, $\mathbf{X}_1 $ excludes the intercept term from the exogenous variable matrix $\mathbf{X}$, $\mathbf{Z}$ is an $n\times k_z$ matrix of endogenous regressors, and $\mathbf{X}^*=\mathbf{X}=[\bm{l}, \mathbf{X}_1]$ if $\mathbf{W}$ is not row-normalized to have row sums equal to one, and $\mathbf{X}^{*}=\mathbf{X}_1$ otherwise. The $\bm\beta_{10}$, $\bm\beta_{20}$, $\bm\beta_{30}$ and $\bm\beta_{40}$ are conformable true parameters, and $\mathbf{W}$, $\mathbf{Y}$ and $\mathbf{V}$ have the same meanings as those in (ref). The Durbin regressors $\mathbf{W} \mathbf{X}_1$ are neighbors' characteristics and capture exogenous externalities. When $\mathbf{W}$ is row-normalized, $\mathbf{W} \bm{l}=\bm{l}$ is the intercept term; when $\mathbf{W}$ is not row-normalized, $\mathbf{W} \bm{l}$ is also a Durbin regressor. In particular, if $\mathbf{W}$ is not row-normalized and has binary elements, $\mathbf{W} \bm{l}$ is a vector of out-degrees that measure the overall numbers of links for each spatial unit. Model (ref) includes Durbin regressors explicitly since the MESS structure and the Durbin regressors lead to some irregular features of the N2SLS estimator. Model (ref) has not considered Durbin regressors explicitly but can allow for that, where the theoretical analysis will not be affected although the related expressions for estimators need to be modified accordingly. To focus on the N2SLS estimation, a MESS process for the disturbances is not considered in (ref).\footnote{If there is a MESS process for the disturbances, then as in Jin.Wang2022, the GMM estimation with both linear and quadratic moments can be considered, since instrumental variables alone are not enough to identify parameters for the disturbance process.}
Let $\mathbf{F}$ be an $n\times k_f$ full rank IV matrix for the N2SLS estimation, where $k_f$ is not smaller than the total number of parameters in $\bm\theta=(\lambda, \bm\beta')'$ with $\bm\beta = (\bm\beta_1',\bm\beta_2,\bm\beta_3',\bm\beta_4')'$. For example, $\mathbf{F}$ can be the matrix formed by the independent columns of $[\bm{l},\mathbf{X}_1,\mathbf{W}\bm{l},\mathbf{W}\mathbf{X}_1,\mathbf{W}^2\bm{l},\mathbf{W}^2\mathbf{X}_1,\bar\mathbf{Z}]$, where $\bar\mathbf{Z}$ is the IV matrix for $\mathbf{Z}$.\footnote{If $\mathbf{W}$ is row-normalized, then $\mathbf{W}\bm{l}$ and $\mathbf{W}^2\bm{l}$ are redundant.} Assume that the elements of $\mathbf{V}$ are independent conditional on $\mathbf{F}$ but can have different conditional variances so that $\bm\Sigma=\operatorname*{E}(\mathbf{V} \mathbf{V}'|\mathbf{F})$ is a diagonal matrix of conditional variances. Denote $\mathbf{D}=[\mathbf{X}^*,\mathbf{W}\bm{l},\mathbf{W}\mathbf{X}_1,\mathbf{Z}]$ and $\bm\Pi=\mathbf{F}'\bm\Sigma\mathbf{F}$. The infeasible N2SLS estimation, as if $\bm\Sigma$ were known, has the objective function
The N2SLS estimator $\hat{\bm\theta}$ derived by minimizing $Q(\bm\theta)$ is consistent under regularity conditions.
Let $\bm\delta =(\bm\beta_1',\bm\beta_2)'$ and $\bm\xi=(\bm\beta_3',\bm\beta_4')'$ when $\mathbf{W}$ is row-normalized, and let $\bm\delta=\bm\beta_1$ and $\bm\xi=(\bm\beta_2,\bm\beta_3',\bm\beta_4')'$ when $\mathbf{W}$ is not row-normalized. Then, $\bm\xi$ contains the coefficients for the Durbin and endogenous regressors. When $\bm\xi_0\ne 0$, all components of $\hat{\bm\theta}$ are $\sqrt n$-consistent and $\hat{\bm\theta}$ has the asymptotic distribution
However, some components of $\hat{\bm\theta}$ have a rate of convergence slower than $\sqrt n$ and are not asymptotically normal in the case that $\bm\xi_0= 0$, i.e., the Durbin and endogenous regressors are irrelevant, which is unknown when estimation is considered.
When $\bm\xi_0=0$, we have
where $k^*$ is the number of columns in $\mathbf{X}^*$, $\bm\delta_{20}$ is the last element of $\bm\delta_{0}$ and $\bm\delta_{10}$ contains the remaining elements. Thus, $\frac{1}{\sqrt{n}}\frac{\partial Q(\bm\theta_0)}{\partial\lambda}$ and $\frac{1}{\sqrt{n}}\frac{\partial Q(\bm\theta_0)}{\partial\bm\beta}$ are linearly dependent with probability approaching one (w.p.a.1.). As a result, $\frac{1}{n}\frac{\partial Q(\bm\theta_0)}{\partial\bm\theta} \frac{\partial Q(\bm\theta_0)}{\partial\bm\theta'}$ is singular w.p.a.1. In addition, we can show that $\frac{1}{n}\frac{\partial^2 Q(\bm\theta_0)}{\partial\bm\theta \partial\bm\theta'} $ is also singular for large $n$. Hence, the usual method of deriving the asymptotic distribution of an estimator based on the mean value theorem expansion of the first order condition will not work.
The asymptotic distribution of $\hat{\bm\theta}$ in the case with $\bm\xi_0=0$ can be derived by first raparameterizing the model so that the derivative of the new N2SLS objective function with respect to a new parameter is exactly zero and then investigating a third order Taylor expansion of the first order condition at the true parameter vector. Let $\bar{\bm\Pi}=\operatorname*{E}(\bm\Pi)$, $k_d$ be the number of columns in $\mathbf{D}$, $J$ be a random vector that follows the normal distribution $N(0,\bm\Delta)$, where
and $L= J_2 - \lim_{n\to\infty}[\frac{2}{n}\operatorname*{E}(\mathbf{D}'\mathbf{F})\bar{\bm\Pi}^{-1} \operatorname*{E}(\mathbf{F}'\mathbf{D})]^{-1} \frac{1}{n}\operatorname*{E}(\mathbf{D}'\mathbf{F})\bar{\bm\Pi}^{-1} \operatorname*{E}(\mathbf{F}'\mathbf{W}^2\mathbf{X})\bm\delta_0 J_1$, where $J_1$ is the first element of $J$ and $J_2$ contains the remaining elements of $J$. Then, in the case with $\bm\xi_0=0$, the N2SLS estimator $\hat{\bm\theta}=(\hat\lambda,\hat{\bm\beta}_1',\hat{\bm\beta}_2,\hat{\bm\beta}_3',\hat{\bm\beta}_4')'$ has the asymptotic distribution
where $I(\cdot)$ denotes the indicator function, $J_{2x^*}$ and $L_{x^*}$ are vectors consisting of the first $k^*$ elements of $J_2$ and $L$ respectively, $J_{2z}$ and $L_z$ are vectors consisting of the last $k_z$ elements of $J_2$ and $L$ respectively, and $B$ is a Bernoulli random variable with success probability described in Jin.Lee2018. Thus, only $\hat{\bm\beta_1}$ and $\hat{\bm\beta_4}$ are $\sqrt n$-consistent, and the remaining components of $\hat{\bm\theta}$ have a slow rate $n^{1/4}$ of convergence and follow non-standard asymptotic distributions.
The above N2SLS estimator is an infeasible estimator as $\bm\Pi$ is unknown. A feasible N2SLS estimator can be derived as follows. We may first derive an initial consistent but inefficient N2SLS estimator, e.g., the minimizer $\check{\bm{\theta}}$ of $(e^{\lambda \mathbf{W}}\mathbf{Y} -\mathbf{D}\bm{\beta})'\mathbf{F}(\mathbf{F}'\mathbf{F})^{-1}\mathbf{F}'(e^{\lambda \mathbf{W}}\mathbf{Y} -\mathbf{D}\bm{\beta})$, and then consider the feasible N2SLS estimation with the objective function $\check Q(\bm\theta) =(e^{\lambda \mathbf{W}}\mathbf{Y} -\mathbf{D}\bm{\beta})'\mathbf{F}(\mathbf{F}'\check{\bm\Sigma}\mathbf{F})^{-1}\mathbf{F}'(e^{\lambda \mathbf{W}}\mathbf{Y} -\mathbf{D}\bm{\beta})$, where $\check{\bm\Sigma}=\text{Diag}(\check v_1^2,\cdots,\check v_n^2)$ with $\check v_i$ the $i$th element of $e^{\check\lambda \mathbf{W}}\mathbf{Y} -\mathbf{D}\check{\bm{\beta}}$. The feasible N2SLS estimator $\tilde{\bm\theta}$ has the same asymptotic distribution as the infeasible estimator $\hat{\bm\theta}$.
As $\bm\xi_0=0$ and $\bm\xi_0\ne 0$ lead to different asymptotic distributions of $\hat{\bm\theta}$, Jin.Lee2018 propose several tests for the hypothesis that $\bm\xi_0=0$. Depending on whether $\bm\xi_0=0$ is rejected or not, inference can be based on (ref) or (ref). Consider the case with $\bm{\xi}_0\ne 0$ as an example. By (ref), the variance of $\tilde{\bm\theta}$ can be estimated by $[(-\mathbf{W}\mathbf{D}\tilde{\bm\beta},\mathbf{D})'\mathbf{F} (\mathbf{F}'\tilde{\bm\Sigma}\mathbf{F})^{-1} \mathbf{F}'(-\mathbf{W}\mathbf{D}\tilde{\bm\beta},\mathbf{D})]^{-1}$, where $\tilde{\bm\Sigma}=\text{Diag}(\tilde v_1^2,\cdots,\tilde v_n^2)$ with $\tilde v_i$ the $i$th element of $e^{\tilde\lambda \mathbf{W}}\mathbf{Y} -\mathbf{D}\tilde{\bm{\beta}}$.
An interesting alternative estimation method is the adaptive group LASSO (AGLASSO), which can implement model selection and estimation simultaneously. The resulting estimator has the oracle properties Fan.Li2001, so that the true model can be selected w.p.a.1.\ and the estimator always has the $\sqrt n$-rate of convergence and asymptotic normal distribution. The AGLASSO objective function to be minimized is
where $\alpha_n$ is a tuning parameter that is positive and converges to zero, $\check{\bm\xi}$ is an initial consistent estimator, and $\mu$ is some positive number such as $1$ or $2$. Under regularity conditions, the AGLASSO estimator $\dot{\bm\theta}$ is consistent. In the case that $\bm{\xi}_0=0$, the probability that $\dot{\bm\xi}=0$ goes to one as $n$ goes to infinity, that is, $\dot{\bm\theta}$ has the sparsity property, and for the remaining parameters $\bm\psi=(\lambda,\bm\delta')'$, the AGLASSO estimator has an asymptotic normal distribution as if $\bm\xi_0$ were known:
in the case that $\bm\xi_0\ne 0$, under the condition that $\alpha_n =o(n^{-1/2})$ and other regularity conditions, $\dot{\bm\theta}$ has the same asymptotic normal distribution as that stated in (ref). Similar to the variance estimation of $\tilde{\bm\theta}$, the variance of $\dot{\bm\psi}$ for the case with $\bm{\xi}_0=0$ can be estimated by $[(-\mathbf{W}\mathbf{X}\dot{\bm\delta},\mathbf{X})'\mathbf{F} (\mathbf{F}'\dot{\bm\Sigma}\mathbf{F})^{-1} \mathbf{F}'(-\mathbf{W}\mathbf{X}\dot{\bm\delta},\mathbf{X})]^{-1}$, where $\dot{\bm\Sigma}$ is defined similarly to $\tilde{\bm\Sigma}$.
A practical question for the AGLASSO estimator is the selection of the tuning parameter $\alpha_n$. We can use an information criterion to choose $\alpha_n$. To make the dependence of $\dot{\bm\theta}$ on $\alpha_n$ explicit, denote the minimizer of $\frac{1}{n}\check Q(\bm\theta) + \alpha \|\tilde{\bm\xi}\|^{-\mu}\|\bm\xi\|$ by $\dot{\bm\theta}_\alpha$. Correspondingly, the AGLASSO estimator of $\bm\xi$ is $\dot{\bm\xi}_\alpha$. Consider the following information criterion:
where $\Gamma_n>0$ satisfies $\Gamma_n\to 0$ and $n^{1/2}\Gamma_n\to \infty$ as $n\to\infty$. For example, we may take $\Gamma_n=O(n^{-1/4})$. The tuning parameter chosen by minimizing $h_n(\alpha)$ can achieve model selection consistency.
The Monte Carlo results presented in Jin.Lee2018 show that the N2SLS and AGLASSO estimators have similar performance in the regular case with $\bm{\xi}_0\ne 0$, but the AGLASSO estimator performs significantly better in the irregular case with $\bm{\xi}_0= 0$. Thus, we suggest the use of the AGLASSO estimator.
There are alternative methods in the literature that can be used to compute $e^{\alpha \mathbf{A}}$, where $\alpha$ is a scalar parameter and $\mathbf{A}$ is an $n\times n$ matrix. The computation methods include the Taylor series approximation, Pad{\'e} approximation, ordinary differential equation methods, polynomial methods, matrix decomposition methods, splitting methods and Krylov space methods. Moler:1978, Moler:2003 assess the effectiveness of nineteen methods according to the following attributes: (i) generality, (ii) reliability, (iii) stability, (iv) accuracy, (v) efficiency, (vi) storage requirements, (vii) ease of use, and (viii) simplicity. They conclude that though “none (of the methods in their paper) are completely satisfactory,” a scaling and squaring method with either the rational Pad{\'e} or Taylor approximants can be the most effective one to compute the matrix exponential terms. Popular software such as Python, R, MATLAB and Mathematica provides functions that can be used to compute the matrix exponential of a given matrix. For example, MATLAB (function expm), Mathematica (function MatrixExp) and Python (function scipy.linalg.expm) utilize a scaling and squaring method combined with a Pad{\'e} approximation for the computation of matrix exponential terms.
As pointed out by Moler:1978, Moler:2003, all methods suggested in the literature are “dubious” in the sense that a sole method may not be entirely reliable for all applications. For example, in the context of MESS-type models, the scaling and squaring method combined with the Pad{\'e} approximation as implemented in MATLAB through its expm function can be highly costly in terms of computation time yangfast. Our analysis on the estimation of the MESS(1,1) model indicates that we need to compute terms such as $e^{\lambda \mathbf{W}} e^{\rho \mathbf{M}} \mathbf{Y}$ and $e^{\rho \mathbf{M}} \mathbf{X}$. Lesage:2007 suggest that we should provide approximations to $e^{\lambda \mathbf{W}} e^{\rho \mathbf{M}} \mathbf{Y}$ and $e^{\rho \mathbf{M}} \mathbf{X}$ in terms of the matrix-vector product terms instead of providing approximations to $e^{\lambda \mathbf{W}}$ and $e^{\rho \mathbf{M}}$. This matrix-vector product method can reduce the computation time significantly.
In the following, we show how to apply the matrix-vector product method to $e^{\lambda \mathbf{W}} e^{\rho \mathbf{M}} \mathbf{Y}$ and $e^{\rho \mathbf{M}} \mathbf{X}$. Let $\text{Diag}(a_1,\hdots,a_n)$ be the $n\times n$ diagonal matrix with the $i$th diagonal element $a_i$. We first consider $e^{\lambda \mathbf{W}} e^{\rho \mathbf{M}} \mathbf{Y}$. We can truncate the matrix exponential terms at the $(q+1)$th order and express $e^{\lambda \mathbf{W}} e^{\rho \mathbf{M}} \mathbf{Y}$ as
where
The result in (ref) expresses $\erme^{\lambda \mathbf{W}}\mathbf{Y}$ in terms of $\mathbf{Y}_j$ and $\mathbf{D}_j$ for $j\in\{1,2,3\}$. These terms can be computed once, and then supplied as the inputs of the objective function in an optimization solver.
Let $\mathbf{X}=[\mathbf{X}_1,\mathbf{X}_2,\hdots,\mathbf{X}_k]$, where $\mathbf{X}_i$ is the $i$th column of $\mathbf{X}$. Then, we can express $e^{\rho \mathbf{M}} \mathbf{X}$ as
where
The approximation in (ref) indicates that the computation of $e^{\rho \mathbf{M}}\mathbf{X}$ also requires only the matrix-vector product operations. We can compute $\mathbb{X}$ and $\mathbf{D}_4$ only one time and then pass these terms as the inputs of the objective function in an optimization solver.
In an extensive Monte Carlo simulation study, yangfast compared the computation time required by the matrix-vector product method with the expm function of MATLAB. For the QMLE, they demonstrated that the matrix-vector product method reduced computation time by $98\%$ to $99\%$ compared to the expm function. In the case of GMME, the computation time decreased by $95\%$ to $97\%$. In the context of the Bayesian estimator, the computation time was reduced by at least $99\%$.
For the MESS(1,1) model in (ref), the marginal effect of a change in $\mathbf{X}_k$ on $\operatorname*{E}(\mathbf{Y})$ is given by $e^{-\lambda_0 \mathbf{W}}\boldsymbol{\beta}_{0k}$, where $\boldsymbol{\beta}_{0k}$ is the $k$th element of the true coefficient vector $\boldsymbol{\beta}_0$. Lesage:2009 define three scalar measures for the marginal effect to ease the interpretation and presentation of this marginal effect:
The ADI, AII and ATI are, respectively, the average of the main diagonal elements of $e^{-\lambda_0 \mathbf{W}}\boldsymbol{\beta}_{0k}$, the average of the off-diagonal elements of $e^{-\lambda_0 \mathbf{W}}\boldsymbol{\beta}_{0k}$, and the average of all the elements of $e^{-\lambda_0 \mathbf{W}}\boldsymbol{\beta}_{0k}$.
There are alternative ways that can be used to determine the dispersions of these scalar impact measures Arbia:2020. In the Bayesian estimation approach, a sequence of random draws for each impact measure can be obtained by using the posterior draws. Then, the mean and the standard deviation calculated from each sequence of impact measures can be used for inference.
In the classical estimation approaches, the delta method can be used to determine the asymptotic distributions of impact measures. Applying the mean value theorem to ADI yields
where $\mathbf{A}_{1}=\left(-\frac{1}{n}\text{tr}(e^{-\lambda_0 \mathbf{W}}\mathbf{W} \boldsymbol{\beta}_{0k}),\frac{1}{n}\text{tr}(e^{-\lambda_0 \mathbf{W}})\right)$ and $\mathbf{B}$ is the asymptotic covariance of $\sqrt{n}(\hat{\lambda}-\lambda_0,\hat{\boldsymbol{\beta}}_{k}-\boldsymbol{\beta}_{0k})$. Thus, we can estimate the asymptotic variance of the direct impact as $\frac{1}{n}\hat{\mathbf{A}}_{1}\hat{\mathbf{B}}\hat{\mathbf{A}}^{'}_{1}$, where $\hat{\mathbf{A}}_{1}=\left(-\frac{1}{n}\text{tr}(e^{-\hat{\lambda}\mathbf{W}}\mathbf{W}\hat{\boldsymbol{\beta}}_{k}), \frac{1}{n}\text{tr}(e^{-\hat{\lambda}\mathbf{W}})\right)$, and $\hat{\mathbf{B}}$ is the estimated asymptotic covariance of $\sqrt{n}(\hat{\lambda}-\lambda_0,\hat{\boldsymbol{\beta}}_{k}-\boldsymbol{\beta}_{0k})$. Applying the mean value theorem to ATI$=\frac{1}{n}\hat{\boldsymbol{\beta}}_{k}\bm{l}' e^{-\hat{\lambda}\mathbf{W}}\bm{l}$, we obtain
where $\mathbf{A}_{2}=\left(-\frac{1}{n}\boldsymbol{\beta}_{k}\bm{l}' e^{-\lambda_0\mathbf{W}}\mathbf{W}\bm{l},\frac{1}{n}\bm{l}' e^{-\lambda_0\mathbf{W}}\bm{l} \right)$. Thus, Var$(\frac{1}{n}\hat{\boldsymbol{\beta}}_{k}\bm{l}' e^{-\hat{\lambda}\mathbf{W}}\bm{l} )$ can be estimated by $\frac{1}{n}\hat{\mathbf{A}}_{2}\hat{\mathbf{B}}\hat{\mathbf{A}}^{'}_{2}$, where $\hat{\mathbf{A}}_{2}=\bigl(-\frac{1}{n}\hat{\boldsymbol{\beta}}_{k}\bm{l}' e^{-\hat{\lambda}\mathbf{W}}\mathbf{W}\bm{l},\frac{1}{n}\bm{l}' e^{-\hat{\lambda}\mathbf{W}}\bm{l}\bigr)$. Finally, applying the mean value theorem to the estimator of AII, we obtain
Then, an estimate of $\mathrm{Var}\left(\frac{1}{n}\left(\hat{\boldsymbol{\beta}}_{k}\bm{l}' e^{-\hat{\lambda}\mathbf{W}}\bm{l} -\text{tr}(e^{-\hat{\lambda}\mathbf{W}}\hat{\boldsymbol{\beta}}_{k})\right)\right)$ is given by $\frac{1}{n}(\hat{\mathbf{A}}_{2}-\hat{\mathbf{A}}_{1})\hat{\mathbf{B}}(\hat{\mathbf{A}}_{2}-\hat{\mathbf{A}}_{1})^{'}$.
Various approaches have been proposed in the literature to implement model selection. In this section we review the available approaches.
The classical tests, such as the Wald, Lagrange multiplier (Rao score) and likelihood ratio tests, for inference on spatial parameters can be formulated by using the results on the asymptotic distributions of the estimators Anselin:1988, Anselin:1996,Anselin:2001, Lesage:2009, Elhorst:2014, Dogan:2018. In the literature, to test non-nested hypotheses, the Cox statistic and the J statistic are adapted for mainly spatial autoregressive models Anselin:1984, Anselin:1986, Kelejian:2008, Kelejian:2011, Burridge:2012, Jin:2013. These non-nested testing approaches can also be used for the model selection problem between the spatial autoregressive models and the MESS models.
In the J-test approach, we augment the null model with the predictor from the alternative model and then check whether the predictor can add significantly to the explanatory power of the augmented model Davidson:1981. HAN2013250 consider the J-test for the model selection problem between the SARAR(1,0) and MESS (1,0) models. When the SARAR(1,0) model is the null model, we can formulate the null and the alternative hypotheses as
where $\mathbf{S}^{ex}(\lambda)=e^{\lambda \mathbf{W}}$ and $\boldsymbol{\beta}^{ex}$ is a conformable parameter vector for $\mathbf{X}$ in the alternative model. As in Kelejian:2011, HAN2013250 consider two predictors based on the alternative model. These predictors are $\hat{\mathbf{Y}}_{1}=\mathbf{S}^{ex}(\hat{\lambda})^{-1}\mathbf{X}\hat{\boldsymbol{\beta}}^{ex}$ and $\hat{\mathbf{Y}}_{2}=(\mathbf{I}_n-\mathbf{S}^{ex}(\hat{\lambda}))\mathbf{Y}+\mathbf{X}\hat{\boldsymbol{\beta}}^{ex}$, where $\hat{\lambda}$ and $\hat{\boldsymbol{\beta}}^{ex}$ are the QMLEs of $\lambda$ and $\boldsymbol{\beta}^{ex}$. Note that the first predictor is based on the reduced form of the alternative model while the second predictor is derived from the identity $\mathbf{Y}=(\mathbf{I}_n-\mathbf{S}^{ex}(\lambda))\mathbf{Y}+\mathbf{X}\boldsymbol{\beta}^{ex}+\mathbf{V}$. Then, the null model can be augmented with these predictors to obtain the following testing equation:
for $r_1=1,2$. Denote $\mathbf{V}(\boldsymbol{\eta}_{r_1})=(\mathbf{I}_n-\alpha \mathbf{W})\mathbf{Y}-\mathbf{X}\boldsymbol{\beta}-\hat{\mathbf{Y}}_{r_1}\delta_{r_1}$, where $\boldsymbol{\eta}_{r_1}=(\alpha,\boldsymbol{\beta}^{'},\delta_{r_1})^{'}$. To estimate the augmented model, HAN2013250 consider a GMME based on the following vector of linear and quadratic moment functions:
where $\mathbf{F}$ is a full-column rank matrix of IVs and $\mathbf{P}_m$'s are $n\times n$ matrices of constants with $\text{tr}(\mathbf{P}_m)=0$ for $m=1,\hdots,q$. Following KP:2010, the IV matrix $\mathbf{F}$ can consist of the linearly independent columns of $\left(\mathbf{X}, \mathbf{W}\mathbf{X},\hdots, \mathbf{W}^d\mathbf{X}\right)$, where $d$ is a positive constant. Let $\boldsymbol{\Xi}=\operatorname*{E}[g(\boldsymbol{\eta}_{0r_1})g^{'}(\boldsymbol{\eta}_{0r_1})]$, where $\boldsymbol{\eta}_{0r_1}=(\alpha_0,\boldsymbol{\beta}^{'}_0,0)^{'}$ is the true parameter vector under $H_0$. Then, using Lemma (ref), it can be shown that
where $\boldsymbol{\omega}=[\text{vec}_D(\mathbf{P}_1),\hdots,\text{vec}_D(\mathbf{P}_q)]$. Let $\frac{1}{n}\hat{\boldsymbol{\Xi}}$ be a consistent estimator of $\frac{1}{n}\boldsymbol{\Xi}$. Then, the feasible optimal GMME of $\boldsymbol{\eta}_{0r_1}$ is defined by $\hat{\boldsymbol{\eta}}_{r_1}=\operatorname*{arg\,min}_{\boldsymbol{\eta}_{r_1}} g^{'}(\boldsymbol{\eta}_{r_1})\hat{\boldsymbol{\Xi}}^{-1}g(\boldsymbol{\eta}_{r_1})$. Let $\lambda_{sar}^{*}$ and $\boldsymbol{\beta}^{ex*}_{sar}$ be the pseudo true values of $\lambda_0$ and $\boldsymbol{\beta}^{ex}_0$ under the null model, respectively. Define $\mathbf{S}^{ex*}_{sar}=e^{\lambda_{sar}^{*}\mathbf{W}}$, $\mathbf{S}=\mathbf{I}_n-\alpha_0\mathbf{W}$ and $\mathbf{G}=\mathbf{W}\mathbf{S}^{-1}$. Then, under some assumptions, HAN2013250 show that
for $r_1=1,2$, where
We summarize the estimation of $\boldsymbol{\eta}_{r_1}$ in Algorithm (ref).
The result in (ref) can be used to construct the J statistic in three different ways: (i) the Wald (W) statistic, (ii) the distance difference (DD) statistic, and (iii) the gradient (G) statistic Newey:1987. Let $\mathbf{R}=(\boldsymbol{0}_{1\times(k+1)},1)$ and $\hat{\mathbf{D}}_{r_1}$ be the plug-in estimator of $\mathbf{D}_{r_1}$ based on $\hat{\boldsymbol{\eta}}_{r_1}$ for $r_1=1,2$. Then, the first two statistics are given as
Let $\tilde{\boldsymbol{\eta}}_{r_1}=\operatorname*{arg\,min}_{\{\boldsymbol{\eta}_{r_1}|\delta_{r_1}=0\}} g^{'}(\boldsymbol{\eta}_{r_1})\hat{\boldsymbol{\Xi}}^{-1}g(\boldsymbol{\eta}_{r_1})$ be the restricted optimal GMME. Then, the gradient test statistic is defined by
where $\tilde{\mathbf{D}}_{r_1}$ is the plug-in estimator of $\mathbf{D}_{r_1}$ based on $\tilde{\boldsymbol{\eta}}_{r_1}$ for $r_1=1,2$. Under $H_0$, these statistics have a chi-squared distribution with one degree of freedom. Thus, we will reject $H_0$ at the $5\%$ significance level if the test statistics are larger than $3.84$.
When using the MESS model as the null model, the null and the alternative hypotheses take the following form:
Let $\hat{\alpha}$ and $\hat{\boldsymbol{\beta}}$ be the QML estimates of $\alpha_0$ and $\boldsymbol{\beta}_0$ from the alternative model. Again, we consider two predictors $\hat{\mathbf{Y}}_{1}=(\mathbf{I}_n-\hat{\alpha}\mathbf{W})^{-1}\mathbf{X}\hat{\boldsymbol{\beta}}$ and $\hat{\mathbf{Y}}_{2}=\hat{\alpha}\mathbf{W}\mathbf{Y}+\mathbf{X}\hat{\boldsymbol{\beta}}$. Thus, the augmented model is given by
for $r_2=1,2$. Let $\boldsymbol{\psi}_{r_2}=(\lambda,\boldsymbol{\beta}^{ex'},\delta_{r_2})^{'}$, $\boldsymbol{\psi}_{0r_2}=(\lambda_0,\boldsymbol{\beta}^{ex'},0)^{'}$ for $r_2=1,2$, and $\alpha_{ex}^{*}$ and $\boldsymbol{\beta}_{ex}^{*}$ be the pseudo true values of $\alpha$ and $\boldsymbol{\beta}$ under the null model. HAN2013250 consider the non-linear 2SLS estimator (N2SLSE) for the estimation of the augmented model. Let $g(\boldsymbol{\psi}_{r_2})=\mathbf{F}^{'}\mathbf{V}(\boldsymbol{\psi}_{r_2})$ be the vector of linear moment functions, where $\mathbf{V}(\boldsymbol{\psi}_{r_2})=\mathbf{S}^{ex}(\lambda)\mathbf{Y}-\mathbf{X}\boldsymbol{\beta}^{ex}+\hat{\mathbf{Y}}_{r_2}\delta_{r_2}$ for $r_2=1,2$. Then, the N2SLSE is defined by
Under some assumptions, it can be shown that
where $\mathbf{D}_{1}=\mathbf{F}^{'}\left(\mathbf{W}\mathbf{X}\boldsymbol{\beta}^{ex}_0,\mathbf{X},\mathbf{S}_{ex}^{*-1}\mathbf{X}\boldsymbol{\beta}^{*}_{ex}\right)$ and $\mathbf{D}_{2}=\mathbf{F}^{'}\left(\mathbf{W}\mathbf{X}\boldsymbol{\beta}^{ex}_0,\mathbf{X},\alpha_{ex}^{*}\mathbf{W} \mathbf{S}^{ex-1}\mathbf{X}\boldsymbol{\beta}_{0}^{ex}+\mathbf{X}\boldsymbol{\beta}^{*}_{ex}\right)$ with $\mathbf{S}_{ex}^{*}=\mathbf{I}_n-\alpha_{ex}^{*}\mathbf{W}$ and $\mathbf{S}^{ex}=e^{\lambda_0\mathbf{W}}$. Algorithm (ref) summarizes the estimation of the augmented model in (ref).
Similar to the previous case in which the SARAR(1,0) model was the null model, we can use the result in (ref) to derive the three test statistics. When the disturbance terms are heteroskedastic, robust methods are necessary to derive consistent estimators. However, the process to derive the three test statistics are similar to the homoskedastic case.
Instead of using the critical value $3.84$ from the asymptotic distribution, we can use the bootstrap method to generate the empirical distribution of the test statistics. In this approach, we can report the bootstrapped p-value, which is the percentage of test statistics based on the bootstrapped samples that are greater than the corresponding test statistic obtained from the actual sample, to decide between $H_0$ and $H_1$ Mac:2009. The bootstrap procedure for testing $H_0: \mathbf{Y}=\alpha \mathbf{W}\mathbf{Y}+\mathbf{X}\boldsymbol{\beta}+\mathbf{V}$ against $H_1:\mathbf{S}^{ex}(\lambda)\mathbf{Y}=\mathbf{X}\boldsymbol{\beta}^{ex}+\mathbf{V}$ is described in Algorithm (ref).
In the heteroskedastic case, besides using the heteroskedasticity robust estimation methods, we also need to use a wild bootstrap approach to generate the bootstrapped versions of the test statistics. The details of this approach are summarized in HAN2013250. The extensive simulation results reported in HAN2013250 indicate that all versions of the J-statistic can perform satisfactorily when the sample size is large.
LIU2019434 propose a non-degenerate Vuong-type model selection test for the model selection between the SARAR(1,1) and MESS(1,1) models. The log-likelihood function of the SARAR(1,1) model in (ref) can be expressed as
where $\boldsymbol{\theta}_1=(\boldsymbol{\beta}^{'},\alpha,\tau,\sigma^2)^{'}$ and $z_{i}(\boldsymbol{\theta}_1)=y_{i}-\alpha \mathbf{W}_{i\cdot}\mathbf{Y}-\tau\mathbf{M}_{i\cdot}\mathbf{Y}+\alpha\tau\sum_{k=1}^{n}m_{ik}\mathbf{W}_{k\cdot}\mathbf{Y}-\mathbf{X}_{i}\boldsymbol{\beta}+\tau\mathbf{M}_{i\cdot}\mathbf{X}\boldsymbol{\beta}$, with $\mathbf{X}_{i}$ being the $i$th row of $\mathbf{X}$, $m_{ik}$ being the $(i,k)$th element of $\mathbf{M}$, and $\mathbf{W}_{i\cdot}$ and $\mathbf{M}_{i\cdot}$ being the $i$th row of $\mathbf{W}$ and $\mathbf{M}$, respectively. Then, we can write the log-likelihood function as $\ln L_1(\boldsymbol{\theta}_1)=\sum_{i=1}^nl_{1i}(\boldsymbol{\theta}_1)$, where $l_{1i}(\boldsymbol{\theta}_1)=-\frac{1}{2}\ln(2\pi)-\frac{1}{2}\ln\sigma^2+\frac{1}{n}\ln\left|\mathbf{I}_n-\alpha\mathbf{W}\right|+\frac{1}{n}\ln\left|\mathbf{I}_n-\tau\mathbf{W}\right|-\frac{1}{2\sigma^2}z_{i}(\boldsymbol{\theta}_1)^2$. Similarly, we can express the log-likelihood function of the MESS(1,1) as
where $\boldsymbol{\theta}_2=(\boldsymbol{\beta}^{'},\lambda,\rho,\sigma^2)^{'}$, $l_{2i}(\boldsymbol{\theta}_2)=-\frac{n}{2}\ln(2\pi)-\frac{n}{2}\ln\sigma^2-\frac{1}{2\sigma^2}h_{i}(\boldsymbol{\theta}_2)^2$ and $h_{i}(\boldsymbol{\theta}_2)$ is the $i$th element of $e^{\rho \mathbf{M}}(e^{\lambda \mathbf{W}}\mathbf{Y}-\mathbf{X}\boldsymbol{\beta})$. LIU2019434 first show that the QMLEs of both models, where one of the models or both models are possibly misspecified, are consistent estimators of their pseudo-true values and are asymptotically normal.
LIU2019434 assume that the true data generating process is unknown, and one of the two models or both models might be misspecified. Let $\boldsymbol{\theta}^{*}_1$ and $\boldsymbol{\theta}^{*}_2$ be the pseudo-true parameter vectors in the SARAR(1,1) and MESS(1,1) models, respectively. Then, the null hypothesis and alternative hypotheses are given by
Let $\text{LR}(\hat{\boldsymbol{\theta}}_1,\hat{\boldsymbol{\theta}}_2)=\ln L_{1}(\hat{\boldsymbol{\theta}})-\ln L_{2}(\hat{\boldsymbol{\theta}})$, where $\hat{\boldsymbol{\theta}}_1$ and $\hat{\boldsymbol{\theta}}_2$ are the QMLEs of the two models. Define $\omega^2=\mathrm{Var}(\frac{1}{\sqrt{n}}\text{LR}(\hat{\boldsymbol{\theta}}_1,\hat{\boldsymbol{\theta}}_2))$ and $g_i(\boldsymbol{\theta}_1,\boldsymbol{\theta}_2)=l_{1i}(\boldsymbol{\theta}_1)-l_{2i}(\boldsymbol{\theta}_2)$. Then, following Shi:2017, LIU2019434 consider the following test statistic:
where $\hat{\sigma}$ is a data-dependent scalar, $\hat{\omega}^2$ is an estimator of $\omega^2$ and $U\sim N(0,1)$. Under some regularity assumptions, it is shown that the test statistic converges to the standard normal distribution under the null hypothesis, i.e., $\hat{T}\xrightarrow{d}N(0,1)$ under $H_0$. Under the alternative hypotheses, they show that $\hat{T}\rightarrow+\infty$ under $H_1$, and $\hat{T}\rightarrow-\infty$ under $H_2$. In a Monte Carlo study, LIU2019434 show that the test statistic has good size and power properties.
The predictive accuracy of a model is usually measured through an information criterion, which is typically defined based on the deviance term $-2\ln p(\mathbf{Y}|\boldsymbol{\theta})$ Gelman:2003. The widely used Akaike information criterion (AIC) takes the following form:
where $\hat{\boldsymbol{\theta}}$ is an estimate of $\boldsymbol{\theta}$ and $p$ is the dimension of $\boldsymbol{\theta}$. In a Bayesian context, Spiegelhalter:2002 suggest another criterion called the deviance information criterion (DIC):
where $\bar{D}(\boldsymbol{\theta})$ is called the posterior mean deviance and $p_D$ is a measure of the effective number of parameters in the model. The posterior mean deviance is defined by $\bar{D}(\boldsymbol{\theta})=-2\operatorname*{E}\left(\ln p(\mathbf{Y}|\boldsymbol{\theta})|\mathbf{Y} \right)$, where the expectation is taken with respect to the posterior distribution of $\boldsymbol{\theta}$. This term serves as a Bayesian measure of model fit. The effective number of parameters is defined by $p_D =\bar{D}(\boldsymbol{\theta})-D(\bar{\boldsymbol{\theta}})=-2\operatorname*{E}\left(\ln p(\mathbf{Y}|\boldsymbol{\theta})|\mathbf{Y} \right)+2\ln p(\mathbf{Y}|\bar{\boldsymbol{\theta}})$, where $\bar{\boldsymbol{\theta}}$ is the posterior mean of $\boldsymbol{\theta}$. Thus, the DIC can be written as
Let $\{\boldsymbol{\theta}^r\}_{r=1}^R$ be a sequence of posterior draws. Then, the first term $\operatorname*{E}\left(\ln p(Y|\theta)|Y \right)$ in the $\text{DIC}$ can be computed by $\operatorname*{E}\left(\ln p(\mathbf{Y}|\boldsymbol{\theta})|\mathbf{Y} \right)\approx\frac{1}{R}\sum_{r=1}^R\ln p(\mathbf{Y}|\boldsymbol{\theta}^r)$. The second term $\ln p(\mathbf{Y}|\bar{\boldsymbol{\theta}})$ in the DIC is computed by evaluating the log-likelihood function at the posterior mean $\bar{\boldsymbol{\theta}}$. Using a decision-theoretic perspective, it can be shown that both AIC and DIC choose the model whose predictive distribution is close to the true data generating process Li:2020.
In our heteroskedastic model, there are alternative likelihood functions: (i) the conditional likelihood function denoted by $p(\mathbf{Y}|\boldsymbol{\theta},\boldsymbol{\eta})$, (ii) the complete-data likelihood function denoted by $p(\mathbf{Y},\boldsymbol{\eta}|\boldsymbol{\theta})$, and (iii) the integrated (or observed) likelihood function denoted by $p(\mathbf{Y}|\boldsymbol{\theta})=\int p(\mathbf{Y},\boldsymbol{\eta}|\boldsymbol{\theta})\text{d}\boldsymbol{\eta}$. The log-conditional likelihood function is readily available and given by
where $\mathbf{H}(\boldsymbol{\eta})=\text{Diag}\left(\eta_{1},\hdots,\eta_{n}\right)$ is the $n\times n$ diagonal matrix with the $i$th diagonal element $\eta_{i}$. As shown in Section 6.2, this function facilitates the Bayesian estimation of the heteroskedastic model. Both the conditional likelihood function and the complete-data likelihood function depend on the high-dimensional latent scale mixture variables. Since these high-dimensional variables can not be estimated precisely, the AIC and DIC formulated with the conditional and complete-data likelihood functions may not perform satisfactorily in model selection exercises Chan:2016. Indeed, the latent variable models violate the conditions of the decision-theoretic perspective, indicating that the AIC and DIC cannot be used as a measure of predictive accuracy Li:2020. Hopefully, the log-integrated likelihood function can be obtained analytically by integrating out the scale mixture variables $\boldsymbol{\eta}$ from the complete-data likelihood function, i.e., $p(\mathbf{Y}|\boldsymbol{\theta})=\int p(\mathbf{Y},\boldsymbol{\eta}|\boldsymbol{\theta})\text{d}\boldsymbol{\eta}=\int p(\mathbf{Y}|\boldsymbol{\eta},\boldsymbol{\theta})p(\boldsymbol{\eta}|\boldsymbol{\theta})\text{d}\boldsymbol{\eta}$. This function can be derived as Dogan:2023
where $y_{i}(\boldsymbol{\delta})$ is the $i$th element of $\mathbf{Y}(\boldsymbol{\delta})=e^{\rho\mathbf{M}}\left(e^{\lambda \mathbf{W}}\mathbf{Y}-\mathbf{X}\boldsymbol{\beta}\right)$ with $\boldsymbol{\delta}=(\lambda,\rho,\boldsymbol{\beta}^{'})^{'}$. This function can be used to formulate AIC and DIC in the heteroskedastic model.
Another popular criterion is the Bayesian information criterion, which can be derived from a large sample approximation to the log-marginal likelihood of a candidate model. Let $\{M_k\}_{k=1}^K$ be a sequence of candidate models. Then, the marginal likelihood of the model $M_k$ can be expressed as $p(\mathbf{Y}_k|M_k)=\int_{\boldsymbol{\Theta}_k}p(\mathbf{Y}|\boldsymbol{\theta}_k,M_k)p(\boldsymbol{\theta}_k|M_k)\text{d}\boldsymbol{\theta}_k$, where $\boldsymbol{\theta}_k$ is the $p_k\times1$ vector of parameters in $M_k$. Then, the Laplace approximation to $\ln p(\mathbf{Y}_k|M_k)$ yields the following BIC measure Schwarz:1978:
The Laplace approximation to $\ln p(\mathbf{Y}_k|M_k)$ can also be used to show that Kass:1995
where $\epsilon>0$ is an arbitrary number and $\text{BF}_{kl}=p(\mathbf{Y}|M_k)/p(\mathbf{Y}|M_l)$ is the Bayes factor of $M_k$ against $M_l$. The result in (ref) indicates that the BIC is also a consistent model selection criterion like the Bayes factor. Moreover, both BIC and the Bayes factor can be interpreted as the measures of predictive accuracy because the marginal likelihood function can be interpreted as the predictive density evaluated at $\mathbf{Y}$ Chan:2016.
In the Bayesian setting described in Section 6.2, Dogan:2023 investigate the performance of AIC, DIC and BIC for both nested and non-nested model selection problems through simulations. They consider four popular MESS specifications and aim to see whether the information criteria can select correct model specification and the correct spatial weights matrix from a pool of candidates. Their extensive simulation results show that these criteria perform satisfactorily and can be useful for selecting the correct model in the specification search exercises.
yangms suggest using a Mallows $C_p$ type selection criterion for selecting a spatial weights matrix from a pool of candidates. Let $\mathcal{W}=\left\{\left(\mathbf{W}_s, \mathbf{M}_s\right): s \in \{1, 2, \hdots, S\}\right\}$ be the pool of spatial weights matrices. The quasi log-likelihood function based on the tuple $(W_{s}, M_{s})$ can be expressed as
where $\Vert\cdot\Vert$ denotes the Euclidean norm. For a given $(\hat{\alpha}_s,\hat{\tau}_s)$ value, the first order conditions of (ref) with respect to $\boldsymbol{\beta}$ and $\sigma^2$ yield
Let $\boldsymbol{\mu}=\operatorname*{E}(\mathbf{Y})=e^{-\alpha_0 \mathbf{W}}\mathbf{X}\boldsymbol{\beta}_0$. Substituting (ref) into $\hat{\boldsymbol{\mu}}_s=e^{-\hat{\alpha}_s\mathbf{W}_s}\mathbf{X}\hat{\boldsymbol{\beta}}_s$, we obtain
where $\widetilde{\mathbf{P}}_s=e^{-\hat{\alpha}_s\mathbf{W}_s}e^{-\hat{\tau}_s\mathbf{M}_s}\widehat{\mathbf{P}}_se^{\hat{\tau}_s \mathbf{M}_s}e^{\hat{\alpha}_s\mathbf{W}_s}$ with $\widehat{\mathbf{P}}_s=e^{\hat{\tau}_s \mathbf{M}_s}\mathbf{X}\left(\mathbf{X}^{'}e^{\hat{\tau}_s \mathbf{M}_s^{'}}e^{\hat{\tau}_s \mathbf{M}_s}\mathbf{X}\right)^{-1}\mathbf{X}^{'}e^{\hat{\tau}_s\mathbf{M}_s^{'}}$. Then, yangms consider the following selection criterion function:
where $\boldsymbol{\Omega}=\sigma_0^2 e^{-\lambda_0 \mathbf{W}}e^{-\rho_0 \mathbf{M}}e^{-\rho_0 \mathbf{M}^{'}}e^{-\lambda_0 \mathbf{W}^{'}}$ is the variance of $\mathbf{Y}$, and the closed forms of $\frac{\partial\hat{\lambda}_s}{\partial \mathbf{Y}^{'}}$, $\frac{\partial \widetilde{\mathbf{P}}_s}{\partial \hat{\lambda}}$ and $\frac{\partial \hat{\rho}_s}{\partial \mathbf{Y}^{'}}$ can be found in yangms. Given an estimator of $\boldsymbol{\Omega}$, we can compute $C_s$ for each $s$. Thus, the selected model is defined as $\hat{s}=\operatorname*{arg\,min}_{s\in\{1, \hdots, S\}}C_s$. Under certain assumptions, yangms show that the selection estimator $\hat{\boldsymbol{\mu}}_{\hat{s}}$ is asymptotically optimal in the sense that it is as efficient as the infeasible estimator that uses the best candidate spatial weights matrix. They also show that the selection procedure is selection consistent in the sense that it chooses the true tuple of weight matrices with probability approaching one as $n\rightarrow\infty$.
Instead of selecting the asymptotically optimal model, it is also possible to use a model averaging scheme that compromises across a set of candidate models. Let $\mathbf{w}=(w_1, \hdots, w_S)^{'}$ be a vector of weights, and $\mathcal{N}=\left\{\mathbf{w}\in[0,1]^S: \sum_{s=1}^Sw_s=1\right\}$ be the set of model weights vectors. Let $\widetilde{\mathbf{P}}(\mathbf{w})=\sum_{s=1}^Sw_s\widetilde{\mathbf{P}}_s$ be the weighted average of $\left\{\widetilde{\mathbf{P}}_1, \hdots, \widetilde{\mathbf{P}}_S\right\}$. Then, the model averaging estimator for $\boldsymbol{\mu}$ is given by
Then, yangms consider the following model weights choice criterion function:
The optimal model weights vector is thus given by $\hat{\mathbf{w}}=\operatorname*{arg\,min}_{\mathbf{w}\in\mathcal{N}}\widehat{C}(\mathbf{w})$. Similar to the model selection procedure, the model averaging estimator $\hat{\boldsymbol{\mu}}(\hat{\mathbf{w}})$ is also asymptotically optimal.
The selection and averaging estimators can also be considered for the high order MESS models. In the case of heteroskedastic models, yangms use a heteroskedasticity robust GMM estimator to formulate the selection and model averaging criterion functions. The extensive simulation results in yangms indicate that the model selection and averaging estimators perform satisfactorily.
In the Bayesian approach, the Bayes factor can be used for both nested and non-nested model selection problems. As shown in Section (ref), the Bayes factor for two models is simply the ratio of the corresponding marginal likelihood functions: $\text{BF}_{kl}=p(\mathbf{Y}|M_k)/p(\mathbf{Y}|M_l)$, where $p(\mathbf{Y}|M_j)=\int_{\boldsymbol{\Theta}_j}p(\mathbf{Y}|\boldsymbol{\theta}_j,M_j)p(\boldsymbol{\theta}_j|M_j)\text{d}\boldsymbol{\theta}_j$ for $j\in\{k,l\}$. Thus, the Bayes factor chooses $M_k$ if $p(\mathbf{Y}|M_k)$ is larger than $p(\mathbf{Y}|M_l)$. If the data is generated from $M_k$, then the Bayes factor will consistently choose $M_k$ over $M_l$. To see this, consider the expectation of the log-Bayes factor under $p(\mathbf{Y}|M_k)$:
which is simply the Kullback-Leibler divergence between $p(\mathbf{Y}|M_k)$ and $p(\mathbf{Y}|M_l)$. Thus, the expectation is strictly positive, unless $p(\mathbf{Y}|M_k)=p(\mathbf{Y}|M_l)$ in which case it is zero.
The Bayes factor reduces to the Savage-Dickey density ratio (SDDR) for the nested model selection problems VW:1995. For example, consider the following null and alternative hypotheses: $H_0:\lambda=0$ against $H_1:\lambda\ne0$ or $H_0:\rho=0$ against $H_1:\rho\ne0$. Let $M_R$ and $M_U$ be respectively the restricted and the unrestricted model. Then, the Bayes factor in favor of the unrestricted model is
where $p(\mathbf{Y}|M_j)$ for $j\in\{U, R\}$ is the corresponding marginal likelihood function. Since our prior distributions are independent, the Bayes factor in (ref) reduces to the SDDR given by
where $p(\lambda=0|M_U)$ and $p(\lambda=0|\mathbf{Y},M_U)$ are respectively the prior and the marginal posterior densities of $\lambda$ evaluated at $\lambda=0$. The $\text{BF}_{UR}$ indicates that if $\lambda=0$ is more likely under the prior relative to the marginal posterior, then the $\text{BF}_{UR}$ provides evidence in favor of $H_1$. Under the prior $\lambda\sim N(\mu_\rho,V_\rho)$, we have $p(\lambda=0|M_U)=(2\pi V_\lambda)^{-1/2}\exp(-\mu^2_\lambda/2V_{\lambda})$. Let $\{\boldsymbol{\beta}^r,\lambda^r, \rho^r,\sigma^{2r}\}_{r=1}^R$ be a sequence of posterior draws. Then, one way to estimate the marginal posterior $p(\lambda=0|\mathbf{Y},M_U)$ is to use the following Rao-Blackwell estimator Smith:1990:
This Rao-Blackwell estimator cannot be used in our case because the conditional posterior density of $\lambda$ does not take a standard form. If we assume that the parameter space of $\lambda$ is contained in the interval $(-\tau,\,\tau)$, where $\tau$ is a finite positive constant, then we may resort to a Griddy-Gibbs sampler to estimate $p(\lambda=0|\mathbf{Y},M_U)$. Algorithm (ref) describes how we can use this approach.
The marginal likelihood function of the MESS type models does not take a closed form. There are alternative methods that can be used to estimate or to approximate the marginal likelihood function. In the homoskedastic case, we can analytically integrate out $\boldsymbol{\beta}$ and $\sigma^2$ under the following priors: (i) $\boldsymbol{\beta}|\sigma^2N(\boldsymbol{\mu}_{\beta},\sigma^2\mathbf{V}_{\beta})$ and $\sigma^2\sim IG(a,b)$ or (ii) $p(\boldsymbol{\beta},\sigma^2)\propto1/\sigma$. However, in order to get the marginal likelihood function, we also need to integrate out the spatial parameters, which is not possible analytically. Then, one approach for computing the marginal likelihood function can be based on a numerical integration method Lesage:2009, Han:2013, Hepple:1995. In the case of MESS(1,1), this approach requires a double numerical integration over the parameter space of $\lambda$ and $\rho$. It is clear that this approach may not be feasible for high order MESS models and for models with heteroskedasticity.
Alternatively, since the conditional posterior distributions of the spatial parameters are in non-standard forms, we may resort to the method suggested by Chib:2001 to estimate the marginal likelihood function. This approach is general enough and only requires the MCMC draws of parameters. In the heteroskedastic case, this approach requires the MCMC draws of the high-dimensional scale mixture variables, therefore it may not produce precise estimates.
The modified harmonic mean method of Dey:1994 can also be used to estimate the marginal likelihood function. This method requires a probability density function $g$ whose support lies in the support of the posterior distribution. The method produces an approximation based on $\operatorname*{E}\left(\frac{g(\boldsymbol{\theta})}{p(\mathbf{Y}|\boldsymbol{\theta})p(\boldsymbol{\theta})}\big|\mathbf{Y}\right)$, where the expectation is taken with respect to $p(\boldsymbol{\theta}|\mathbf{Y})$. The expectation gives the following relationship:
Thus, the marginal likelihood function $p(Y)$ can be estimated by the following estimator:
where $\{\boldsymbol{\theta}^r\}_{r=1}^R$ is a sequence of the posterior draws from $p(\boldsymbol{\theta}|\mathbf{Y})$. Under the condition that $g(\boldsymbol{\theta})/\left(p(\mathbf{Y}|\boldsymbol{\theta})p(\boldsymbol{\theta})\right)$ is bounded above over the support of the posterior distribution, it can be shown that this estimator is a simulation consistent estimator when $R$ goes to infinity Geweke:1999. To guarantee this boundedness condition, following Geweke:1999, we can consider a truncated multivariate normal density for $g$. Let $ A=\{\boldsymbol{\theta}\in\mathbb{R}^p:(\boldsymbol{\theta}-\hat{\boldsymbol{\theta}})^{'}\hat{\boldsymbol{\Omega}}^{-1}(\boldsymbol{\theta}-\hat{\boldsymbol{\theta}})<\chi^2_{\alpha,p}\}$ be the truncation set, where $\hat{\boldsymbol{\theta}}$ is the posterior mean of $\boldsymbol{\theta}$, $\hat{\boldsymbol{\Omega}}$ is the posterior covariance of $\boldsymbol{\theta}$, and $\chi^2_{\alpha,p}$ is the $(1-\alpha)$ quantile of the $\chi^2_p$ distribution. Then, $g$ takes the following form:
where $\mathbf{1}_A(\boldsymbol{\theta})$ is the indicator function taking value $1$ if $\boldsymbol{\theta}\in A$, otherwise $0$.
Note that the computation of the modified harmonic mean estimator requires the integrated likelihood function which is available for both homoskedastic and heteroskedastic models. In the context of spatial autoregressive models, Dogan:2023b investigates the finite sample performance of this estimator along with some other popular information criteria for both nested and non-nested model selection problems. His simulation results show that the modified harmonic mean estimator performs satisfactorily, and can be useful for the specification search exercises in spatial econometrics.
In this section, we conduct a Monte Carlo study to investigate the finite sample properties of the estimators considered in Sections (ref) through (ref). To this end, we consider the following data generating process:
where $\mathbf{X}$ contains two explanatory variables with $\boldsymbol{\beta}_0 =(\beta_{10},\beta_{20})'=(1,\,1)^{'}$. The observations for the first explanatory variable are drawn independently from the standard normal distribution, while the observations for the second explanatory variable are drawn independently from Uniform($0, \sqrt{12}$). The spatial parameters $\lambda_0$ and $\rho_0$ can take values from the set $\{(-2, -1), (-2, 1), (0.5, -1), (0.5, 1)\}$. For $\mathbf{V}$, in the homoskedastic scenario, the elements $v_i$ are i.i.d.\ draws from either (i) the standard normal distribution or (ii) the standardized chi-squared distribution with three degrees of freedom, i.e., $(\chi^2_3-3)/\sqrt{6}$. In the heteroskedastic scenario, we set $\mathbf{V}\sim N\left(\mathbf{0},\text{Diag}(\gamma_1,\hdots,\gamma_n)\right)$ with either (i) $\gamma_{i}=2 \vartheta_i / (\sum_{j=1}^n \vartheta_j/n)$, where $\vartheta_i$ is the number of neighbors for unit $i$ using the description of $\mathbf{W}_1$ below, or (ii) $\gamma_{i}=\exp(0.1+0.35X_{2i})$ and $X_{2i}$ is the $i$th element of $X_2$.
For the spatial weights matrix $\mathbf{W}$, we consider the interaction scenario described in Arraiz:2010. To this end, let $n$ entities be distributed across four quadrants of a square grid in such a way that the number of entities in each quadrant can be arranged to allow for sparse or dense quadrants. The location of each entity across the grid is determined by the $xy$-coordinates on the grid. Let $\underline{c}$ and $\overline{c}$ be two integers. The entities in the northeast quadrant of the grid have discrete coordinates satisfying $(\underline{c}+1)\leq x\leq \overline{c}$ and $(\underline{c}+1)\leq y\leq \overline{c}$, with an increment value of $0.5$. For the other quadrants, the location coordinates are integers satisfying $1\leq x\leq \underline{c}$, $1\leq y\leq \overline{c}$, and $1\leq x\leq \overline{c}$, $1\leq y\leq \underline{c}$. The distance $d_{ij}$ between any two entities $i$ and $j$, located respectively at $(x_1,y_1)$ and $(x_2,y_2)$, is measured by the Euclidean distance given by $d_{ij}=\left[(x_1-x_2)^2+(y_1-y_2)^2\right]^{1/2}$. Then, the $(i,j)$th element of $\mathbf{W}$ is set to $1$ if $0\leq d_{ij}\leq 1$, and to $0$ otherwise. We then row normalize $\mathbf{W}$. In this scenario, varying the values for $\underline{c}$ and $\overline{c}$ leads to a different sample size and a different share of units in the northeast quadrant. We consider the following two combinations: $(\underline{c},\overline{c})=(5, 15)$ and $(\underline{c},\overline{c})=(14, 20)$. The first combination produces a sample size of $486$ and locates $75$ percent of the entities in the northeast quadrant ($\mathbf{W}_1$), whereas the second combination generates a sample size of $485$ and locates $25$ percent of the entities in the northeast quadrant ($\mathbf{W}_2$).
For the spatial weights matrix $\mathbf{M}$, we consider a nearest neighbors scheme. To this end, using the Euclidean distances ($d_{ij}$'s) from the construction of $\mathbf{W}$ above, we let entity $i$ be dependent on its 5 nearest neighbors so that the weights corresponding to these neighbors in the $i$'th row of $\mathbf{M}$ are set to $1$ and the rest are set to zero. Then, $\mathbf{M}$ is row normalized. We use the “makeneighborsw” function from the Spatial Econometrics Toolbox to generate $\mathbf{M}$ Lesage:2009.
We evaluate the performance of the following estimators: (i) the QMLE in (ref), (ii) the ME in (ref), (iii) the IGMME in (ref), (iv) the BGMME in (ref), (v) the RGMME in (ref), (vi) the Bayesian estimator (BE) based on Algorithm 1, and (vii) the robust Bayesian estimator (RBE) based on Algorithm 2. For classical estimation methods, we conduct $1000$ repetitions. In the case of Bayesian estimation, we choose the following priors: $\lambda\sim N(0,100)$, $\rho\sim N(0,100)$, $\boldsymbol{\beta}\sim N(\mathbf{0},100\mathbf{I}_2)$ and $\sigma^2\sim\text{IG}(0.01,0.01)$. We set the number of repetitions to $100$, the number of draws to $1500$, and the burn-ins to $500$. For each method, we report bias, root mean squared error (RMSE) and empirical coverage ratio of a 95% confidence interval.
Tables (ref) and (ref) present the simulation results for two homoskedastic cases: (i) $v_i\sim N(0,1)$ and (ii) $v_i\sim(\chi^2_3-3)/\sqrt{6}$. Similarly, Tables (ref) and (ref) report the simulation results for the heteroskedastic cases. Below, we summarize our main findings from these tables.
In this article, we provide an extensive review of cross-sectional MESS models. We mainly focus on a first-order MESS model to discuss specification, estimation, model selection, and interpretation issues. The primary characteristic of a MESS-type model lies in its use of matrix exponential terms to specify spatial dependence in both the dependent variable and the disturbance term. These models possess several distinctive properties:
We have provided a comprehensive description of various estimation methods, namely, the QML approach, the M-estimation approach, the GMM approach, and the Bayesian estimation approach. This detailed overview enables practitioners to easily choose and adapt a method that aligns with their specific needs. Additionally, we address estimation in the presence of endogenous explanatory variables and the Durbin terms.
In future studies, it might be interesting to consider the MESS in a social interactions scenario, and compare its implications with the SAR-type social interaction models. As the QMLE of the MESS can still remain consistent under an unknown form of heteroskedasticity, allowing for such heteroskedasticity in a social interactions model would be a significant contribution. We also think that the literature on nonlinear spatial models, such as the spatial extensions of the limited dependent variable data models, still holds some open questions, and estimation strategies for the the MESS-type limited dependent variable data models must be studied carefully. Finally, although the matrix-vector product approach to compute the matrix exponential terms can reduce the computation time significantly, we think that a faster and more reliable computation approach would be a significant contribution to the literature.