2.1 An Illustrative Example of the Different Spatial Processes

Before outlining the alternative spatial model specifications, it is useful to contrast the different spatial processes with respect to their substantive implications for empirical political science research.

Spillover effects occur whenever the behavior (endogenous interactions) or certain characteristics (exogenous interactions) of one unit—may this be a country, a (coalition) government, a political party, or any other entity of interest—affect adjacent units (Darmofal Reference Darmofal2015, 5).Footnote ³ For example, the municipalities’ income tax revenues are directly related to their local economic performance. At the same time, the economic prosperity of its neighbors also exerts a positive impact on a municipality’s income tax revenues as its residents might commute to work. These dependencies in income tax revenues produce spillovers between the municipalities: the effect of a change in one unit’s characteristics propagates to its neighbors (e.g., LeSage and Pace Reference LeSage and Pace2009). Hence, adequately understanding the phenomenon of interest—cross-sectional variation in income tax revenues—necessitates the consideration of these exogenous interaction effects among the municipalities that generate substantively meaningful and theoretically relevant indirect spillover effects.

In contrast, cross-sectional dependencies in the disturbances constitute another spatial process that has different substantive implications and requires an alternative model specification. Several circumstances can cause correlation between the units’ residuals, for example spatial clustering in measurement errors. Alternatively, omitting a spatially dependent explanatory variable that is part of the true data-generating process (DGP) from the regression equation leads to correlated errors (Elhorst Reference Elhorst2014b; Darmofal Reference Darmofal2015). With regards to the example given above, it is reasonable to expect that several unobservable but spatially dependent characteristics, like a region’s general appeal as a place of residence, affect a municipality’s revenue from income taxation as well. In contrast to the theorized exogenous interaction effects, these omitted characteristics merely affect the model residuals and there are no relevant indirect effects present in the process that generates the data. Since exactly modeling the true dependence structure is almost impossible (e.g., Juhl Reference Juhl2020b), omitting relevant variables that are spatially correlated can create linkages among the units’ disturbances.

2.2 Common Factors in the Spatial Durbin Model

While the previous example illustrates the crucial differences in the substantive implications of the underlying process that links the observations to one another, model misspecification may cause severe econometric problems as well. In fact, neglecting cross-sectional interdependencies can induce correlation between the regressors and the residuals, resulting in the canonical endogeneity problem (e.g., Gibbons and Overman Reference Gibbons and Overman2012; Betz, Cook, and Hollenbach Reference Betz, Cook and Hollenbach2019).

In order to illustrate the problem of omitted spillover effects, consider the stylized DGP in which a dependent variable $\boldsymbol {y}$ is entirely determined by two uncorrelated regressors, denoted $\boldsymbol {x}$ and $\boldsymbol {z}$, such that $\boldsymbol {y}=\boldsymbol {x}\beta +\boldsymbol {z}$.Footnote ⁴ Assume that $\boldsymbol {z}$ is unobservable and follows a spatial autoregressive process such that $\boldsymbol {z}=\rho \boldsymbol {Wz}+\boldsymbol {u}_{1}$, where $\rho $ is a scalar parameter, $\boldsymbol {W}$ is an exogenously defined connectivity matrix, and $\boldsymbol {u}_{1}$ is a vector of independent and identically distributed normal disturbances with zero mean and a fixed variance.Footnote ⁵ This scenario leads to the spatial error model (SEM) specification:

(1)$$ \begin{align} \boldsymbol{y} &= \boldsymbol{x}\beta+\underbrace{(\boldsymbol{I}_{n}-\rho\boldsymbol{W})^{-1}\boldsymbol{u}_{1}}_{\boldsymbol{z}}. \end{align} $$

Due to the uncorrelatedness of $\boldsymbol {x}$ and $\boldsymbol {z}$, omitting the spatially autocorrelated variable does not lead to endogeneity concerns. Furthermore, the relationship depicted in Equation (1) implies no indirect spillover effects since the cross-unit interactions are confined to the residuals. Consequently, even a nonspatial OLS model specification would provide unbiased but inefficient parameter estimates (Lacombe and LeSage Reference Lacombe and LeSage2015).

Now consider a slightly modified scenario in which the included and the omitted regressors are no longer independent from one another but correlated. To induce correlation between the variables, suppose that the random variable $\boldsymbol {u}_{1}$ in Equation (1) is replaced by $\boldsymbol {u}_{2}$ which is an additive linear function of $\boldsymbol {x}$ and a stochastic disturbance term $\boldsymbol {v}$ such that $\boldsymbol {u}_{2}=\boldsymbol {x}\gamma +\boldsymbol {v}$. In addition to the spatial autocorrelation, the unobserved covariate $\boldsymbol {z}$ is now also correlated with $\boldsymbol {x}$ and the scalar $\gamma \in (0,1]$ as well as the dispersion of $\boldsymbol {v}$ ($\sigma _{v}^{2}$) jointly determine the strength of the correlation. In this slightly modified scenario, the true DGP becomes

(2)$$ \begin{align} \boldsymbol{y} & = \boldsymbol{x}\beta+\underbrace{(\boldsymbol{I}_{n}-\rho\boldsymbol{W})^{-1}(\boldsymbol{x}\gamma+\boldsymbol{v})}_{\boldsymbol{z}}. \end{align} $$

In contrast to the DGP shown in Equation (1), the correlation between the included regressor $\boldsymbol {x}$ and the spatially clustered variable $\boldsymbol {z}$ causes an endogeneity problem if $\boldsymbol {z}$ is omitted from the regression’s systematic part. Importantly, the effect of a change in regressor $\boldsymbol {x}$ on the outcome $\boldsymbol {y}$ in Equation (2) is more complex and both a nonspatial OLS model as well as a SEM model do not yield unbiased estimates since they ignore indirect effects produced by the spatial patterning of the correlated omitted variable. Consequently, the effect of a change in $x_{i}$ is not confined to unit i’s outcome $y_{i}$ but rather pertains to all nonisolated units in the entire system through indirect spillover and instantaneous feedback effects (e.g., LeSage and Pace Reference LeSage and Pace2009; Betz et al. Reference Betz, Cook and Hollenbach2019).

In order to address the endogeneity problem and to identify meaningful spillover effects, the unrestricted SDM model plays an important role. In fact, it serves as a general nesting model in many specification search strategies since it comprises several simpler spatial regression models frequently employed in empirical studies (e.g., Mur and Angulo Reference Mur and Angulo2009; Elhorst Reference Elhorst2010; Angulo and Mur Reference Angulo and Mur2011).Footnote ⁶ By allowing researchers to test different parameter restrictions, the SDM model facilitates the specification of an econometric model that reflects the spatial process generating the data most appropriately. For the hypothetical scenario with one regressor $\boldsymbol {x}$ and the (possibly correlated) unobserved variable $\boldsymbol {z}$ discussed here, the SDM model to be estimated takes the following form:

(3)$$ \begin{align} \boldsymbol{y} & = \rho\boldsymbol{Wy} + \boldsymbol{x}\beta + \boldsymbol{Wx}\theta + \boldsymbol{v}. \end{align} $$

While it is easy to verify that the SDM model reduces to the popular spatial autoregressive (SAR) model, also known as spatial lag model (e.g., Elhorst Reference Elhorst2014b, 5), that features global spillover effects if $\theta = 0$ and to the spatial lag of X (SLX) model featuring local spillovers when $\rho = 0$, it also subsumes the SEM model that rules out any substantive indirect effects by assumption.Footnote ⁷ To illuminate the relationship between these models, it is useful to restate the SEM DGP displayed in Equation (1) by multiplying both sides of the equation by $(\boldsymbol {I}_{n}-\rho \boldsymbol {W})$ and rearranging terms which results in the following structural form (Burridge Reference Burridge1981; Anselin Reference Anselin2003):

(4)$$ \begin{align} (\boldsymbol{I}_{n}-\rho\boldsymbol{W})\boldsymbol{y} & = (\boldsymbol{I}_{n}-\rho\boldsymbol{W})\boldsymbol{x}\beta + \boldsymbol{u}_{1}\nonumber\\ \boldsymbol{y} - \rho\boldsymbol{Wy} & =(\beta\boldsymbol{I}_{n}-\rho\beta\boldsymbol{W})\boldsymbol{x} + \boldsymbol{u}_{1} \nonumber\\ \boldsymbol{y} & = \rho\boldsymbol{Wy} + \boldsymbol{x}\beta + \boldsymbol{Wx}(-\rho\beta) + \boldsymbol{u}_{1}. \end{align} $$

Equation (4) elucidates that the SEM process imposes a nonlinear common factor restriction on the parameter associated with $\boldsymbol {Wx}$ which can be assessed by using estimates from the unrestricted SDM model depicted in Equation (3).Footnote ⁸ More precisely, if the estimates from the unrestricted SDM model satisfy the common factor restriction $\theta =-\rho \beta $, the model can be simplified to a SEM model specification because there are no discernible substantive spillover effects in the data. Given that $E(\hat {\beta })-\beta =0$ in Equation (4), the common factor restriction holds for SEM processes.

However, if $\boldsymbol {x}$ and $\boldsymbol {z}$ are correlated in the true DGP, an endogeneity problem occurs since the variable $\boldsymbol {z}$ is unobserved and the common factor restriction no longer holds. Restating the DGP in Equation (2) that features this correlation in a similar fashion yields:

(5)$$ \begin{align} (\boldsymbol{I}_{n}-\rho\boldsymbol{W})\boldsymbol{y} & = (\boldsymbol{I}_{n}-\rho\boldsymbol{W})\boldsymbol{x}\beta + \boldsymbol{x}\gamma + \boldsymbol{v} \nonumber\\ \boldsymbol{y} & = \rho\boldsymbol{Wy} + \boldsymbol{x}(\beta+\gamma) + \boldsymbol{Wx}(-\rho\beta) + \boldsymbol{v}. \end{align} $$

Given the discussion above, Equation (5) illustrates that, in the presence of an omitted variable that is (i) spatially clustered and (ii) correlated with an included regressor, the true DGP resembles the SDM specification. Although the SDM model shown in Equation (3) provides consistent estimates for $E(\hat {\rho })=\rho $ and $E(\hat {\theta })=-\rho \beta $, the estimate of $\beta $ is asymptotically biased since $E(\hat {\beta }) = \beta + \gamma $. Based on these model estimates, the common factor restriction does not hold because of the endogeneity bias present in $E(\hat {\beta })$. Hence, a violation of the common factor restriction is indicative of the existence of indirect spillover effects that need to be included in the systematic part of the regression model.

As this discussion suggests, the spatially lagged exogenous variables in the SDM model specification can also be understood as instruments for omitted variables that are correlated with included regressors (e.g., Elhorst Reference Elhorst2014b, 18). At the same time, as Gibbons and Overman (Reference Gibbons and Overman2012) emphasize, this strategy provides only weak identification of the model parameters and crucially depends on the exogeneity and the assumed perfect knowledge of $\boldsymbol {W}$. Consequently, the SDM model does not provide a general solution to the omitted variables problem that would allow researchers to identify causal effects. Just like in any other observational study, doing so requires the application of appropriate research designs (see also Betz et al. Reference Betz, Cook and Hollenbach2019; Rüttenauer Reference Rüttenauer2019). Moreover, while the SDM specification features global spillover effects, local spillovers would require a different model specification like the spatial Durbin error model (e.g., Halleck Vega and Elhorst Reference Halleck Vega and Elhorst2015).

3.1 Analytical Derivation and Asymptotic Distribution of the Wald Statistic

Consider a situation in which a test needs to be constructed in order to evaluate a single nonlinear restriction $H_{0}: g(\boldsymbol {\lambda }) = 0$, where $\boldsymbol {\lambda }$ is a parameter vector and $g(\cdot )$ is some function that is continuously differentiable in a neighborhood of $\boldsymbol {\lambda }$. For this general case, the Wald statistic is defined by

(6)$$ \begin{align} w=g(\hat{\boldsymbol{\lambda}})\big[\widehat{V(g(\hat{\boldsymbol{\lambda}}))}\big]^{-1}g(\hat{\boldsymbol{\lambda}}), \end{align} $$

where $\widehat {V(g(\boldsymbol {\hat {\lambda }}))}$ is the estimated variance of $g(\boldsymbol {\hat {\lambda }})$. Under the null hypothesis, w asymptotically follows a $\chi ^{2}$ distribution with the number of degrees of freedom equal to the number of restrictions imposed.

The only complication involved here is that obtaining w necessitates knowledge about the sampling distribution of a nonlinear function. While it is straightforward to compute the value of $g(\boldsymbol {\hat {\lambda }})$ at the parameter estimates, deriving the Wald statistic in Equation (6) additionally requires information about its variability which depends on the estimator $\boldsymbol {\hat {\lambda }}$ and the restriction $g(\cdot )$. However, given the restriction’s nonlinearity, exact distributional results become inapplicable (Greene Reference Greene2012). Instead, the delta method provides an approximation of the restriction’s asymptotic distribution.Footnote ⁹ It is based on a first-order Taylor series expansion of $g(\boldsymbol {\hat {\lambda }})$ around the true parameter vector $\boldsymbol {\lambda }$. Assuming that $\boldsymbol {\hat {\lambda }}$ is a consistent estimator with a limiting distribution defined by ${\sqrt {n}(\boldsymbol {\hat {\lambda }} - \boldsymbol {\lambda }) \overset {d}{\rightarrow } \mathcal {N}(\boldsymbol {0}, \boldsymbol {\Sigma })}$ and that the standard regularity conditions hold (see e.g., Newey and McFadden Reference Newey, McFadden, Engle and McFadden1994), the delta method implies that

(7)$$ \begin{align} \sqrt{n}(g(\boldsymbol{\hat{\lambda}}) - g(\boldsymbol{\lambda})) & \overset{d}{\rightarrow} \mathcal{N}\big(0,\boldsymbol{G}(\boldsymbol{\lambda})\boldsymbol{\Sigma}\boldsymbol{G}{(\boldsymbol{\lambda})'}\big), \end{align} $$

where $\boldsymbol {G}(\boldsymbol {\lambda }) = \partial g(\boldsymbol {\lambda })/ \partial {\boldsymbol {\lambda }'}$ is a row vector of partial derivatives. It follows that the asymptotic distribution of the restriction under $H_{0}$ is $g(\boldsymbol {\hat {\lambda }}) \overset {a}{\sim } \mathcal {N}\big (g(\boldsymbol {\lambda }), \boldsymbol {G}(\boldsymbol {\lambda })n^{-1}\boldsymbol {\Sigma }\boldsymbol {G}{(\boldsymbol {\lambda })'}\big )$. By using the consistent estimates obtained from an unrestricted model, the restriction’s sampling variability derived from the delta method is given by

(8)$$ \begin{align} \widehat{V(g(\boldsymbol{\hat{\lambda}}))} = \boldsymbol{G}(\boldsymbol{\hat{\lambda}})\widehat{\boldsymbol{\Sigma}}\boldsymbol{G}{(\boldsymbol{\hat{\lambda}})'}, \end{align} $$

with $\boldsymbol {G}(\boldsymbol {\hat {\lambda }})$ evaluated at $\boldsymbol {\hat {\lambda }}$ and $\widehat {\boldsymbol {\Sigma }} = n^{-1}\boldsymbol {\Sigma }$ being a consistent estimator of the symmetric, positive definite asymptotic variance–covariance matrix. By substituting Equation (8) into Equation (6), the Wald test statistic can be calculated. Asymptotically, since $\text {plim}_{n\to \infty }\boldsymbol {\hat {\lambda }}=\boldsymbol {\lambda }$, the function $g(\boldsymbol {\hat {\lambda }})$ converges in distribution to $g(\boldsymbol {\lambda })$ with a mean given by $\text {plim}_{n\to \infty }g(\boldsymbol {\hat {\lambda }})=g(\boldsymbol {\lambda })$. At the same time, the necessity to estimate the nonlinear restriction’s sampling variability can cause a mismatch between the Wald statistic’s asymptotic $\chi ^{2}$ distribution and its finite sample distribution which has considerable consequences for hypothesis testing (e.g., Lafontaine and White Reference Lafontaine and White1986; Phillips and Park Reference Phillips and Park1988).

3.2 The Wald Test of Common Factors in Spatial Models

The analytical results derived above are directly applicable to the empirical assessment of the common factor restriction in spatial regression models because the null hypothesis can be expressed as a nonlinear function of the estimates obtained from an unrestricted SDM model. In order to calculate the test statistic, it is necessary to determine the functional representation of the null hypothesis of common factors. Yet, there are numerous algebraically equivalent alternative parameterizations that satisfy $g(\boldsymbol {\hat {\lambda }})=0$, where $\boldsymbol {\hat {\lambda }}=[\hat {\rho },\hat {\beta },\hat {\theta }]$ are the estimates obtained from the SDM model in Equation (3).

Table 1, for example, lists the four alternative expressions of the null hypothesis considered by Gregory and Veall (Reference Gregory and Veall1986) in a time-series context. While all of the alternative statements declare the same restriction, they produce distinct test statistics and p values in finite samples because the approximation of the restriction’s sampling variability is based on the partial derivatives of the parameter estimates. Depending on the exact representation of the null hypothesis, the right part of Table 1 shows that the vector of partial derivatives obtained from these nonlinear functions differ.

Table 1 Algebraically identical formulations of the common factor hypothesis.

As a result, alternative expressions of the common factor hypothesis use different estimators for the nonlinear restriction’s sampling variability which yields distinct test statistics and causes them to converge to the asymptotic $\chi ^{2}$ distribution at individual rates. In large samples, this circumstance is unproblematic as the accuracy of the Taylor series approximation increases in sample size while the contribution of the restriction’s estimated variability to the test statistic becomes negligible. In finite samples, however, the differences between the alternatives can be substantial (e.g., Gregory and Veall Reference Gregory and Veall1985). At worst, alternative and algebraically identical functional representations of the parameter restriction can indicate opposing conclusions regarding its validity despite the fact that the same model estimates are used to calculate the test statistic.

Importantly, while many statistical software packages used to estimate spatial regression models, like Stata or R packages, report results from a Wald test by default, they do not test for common factors. Instead, the null hypothesis these packages evaluate is cross-sectional independence, that is, $\hat {\rho }=0$. Since this is a linear restriction, the Wald statistic’s noninvariance problem does not arise and the tests these packages perform are only indicative of the presence of nonrandom spatial clustering. They do not permit any inferences regarding the spatial process at work. Hence, scrutinizing the common factor restriction requires researchers to amend the Wald test’s null hypothesis.

3.3 Modifications of the Wald Test

Since the application of the Taylor series expansion results in distinct Wald statistics for algebraically equivalent formulations of the null hypothesis, the asymptotic $\chi ^{2}$ distribution might constitute an inappropriate approximation of the statistic’s finite sample distribution for some parameterizations. Therefore, modifications proposed in the literature that attempt to address the Wald test’s noninvariance problem primarily focus on adjusting the statistic’s reference distribution.

Phillips and Park (Reference Phillips and Park1988), for example, show that an Edgeworth expansion of the Wald statistic provides additional information on the statistic’s distribution which can be used to obtain corrected critical values and modified test statistics for each functional representation of the null hypothesis (e.g., de Paula and Cribari-Neto Reference de Paula Ferrari and Cribari-Neto1993; King and Goh Reference King, Goh, Ullah, Wan and Chaturvedi2002). Besides these corrections, simulation techniques allow researchers to generate the empirical distribution under the null hypothesis for each specification of the common factor restriction and base inferences on these reference distributions (e.g., Lafontaine and White Reference Lafontaine and White1986; Goh and King Reference Goh and King1996). In particular, bootstrap methods provide a way to estimate critical values and use them as an alternative to the (corrected) asymptotic critical values which can be unreliable in finite samples (Godfrey and Veall Reference Godfrey and Veall1998). Using the following procedure, it is straightforward to derive bootstrap critical values to test for common factors in spatial regression models:

1. 1. Use the estimates from an unrestricted SDM model and calculate the observed Wald statistic w according to Equation (6) for a given restriction.

2. 2. Estimate the restricted SEM model. With these estimates and the DGP shown in Equation (1), generate 100 bootstrap samples by resampling with replacement from the residual vector to obtain bootstrap disturbances.

3. 3. Repeat step 1 for each bootstrap sample and store the Wald statistic in vector $\boldsymbol {\tilde {w}}$.

4. 4. Sort $\boldsymbol {\tilde {w}}$ in ascending order. The value with rank $(1-\alpha )\times 100+1$ is the estimated bootstrap critical value, $\chi ^{2}_{boot}$, corresponding to a predefined $\alpha $-level (e.g., $\alpha =0.05$).

By comparing w calculated in step 1 to the corresponding bootstrap critical value $\chi ^{2}_{boot}$ from step4, its statistical significance can be assessed. This procedure can be repeated for any functional representation of the common factor restriction in order to obtain individual bootstrap critical values for each restriction and a given region of the parameter space. Thereby, researchers can base inferences on the empirical distribution under the null hypothesis instead of relying on the asymptotic $\chi ^{2}$ distribution. This is especially important since the performance of the Wald test not only depends on the specific expression of the common factor hypothesis but also on the particular region in the parameter space. In fact, previous research shows that there is no single formulation of the restriction that consistently outperforms all alternatives (e.g., Gregory and Veall Reference Gregory and Veall1986; Lafontaine and White Reference Lafontaine and White1986; Phillips and Park Reference Phillips and Park1988).

While Goh and King (Reference Goh and King1996) demonstrate that both asymptotic modifications—the corrected critical values and the improved test statistics—might even deteriorate the Wald test’s power and size properties, they conclude that the bootstrap approach constitutes a useful improvement for applied research (see also Lafontaine and White Reference Lafontaine and White1986; Godfrey and Veall Reference Godfrey and Veall1998; King and Goh Reference King, Goh, Ullah, Wan and Chaturvedi2002). Of course, neither the Edgeworth corrections nor simulation techniques completely resolve the noninvariance problem inherent to the Wald test. Doing so requires the application of alternative tests such as the asymptotically equivalent LR test that is invariant to such reparameterizations (e.g., Mur and Angulo Reference Mur and Angulo2006). Yet, by providing corrections for the Wald test’s empirical size, these modifications reduce the possibility of intentionally manipulating the result by amending the functional expression of the null hypothesis (King and Goh Reference King, Goh, Ullah, Wan and Chaturvedi2002, 260).

4.1 Experimental Setup

In order to investigate the finite sample performance of the Wald test, I conduct Monte Carlo experiments in which I vary the sample size, the strength of the interdependence, and the severity of the omitted variables bias through the degree of correlation between the included and the omitted regressor. Using the spatial process depicted in Equation (2), I generate $1,000$ samples of the outcome vector $\boldsymbol {y}$ for each of the parameter configurations. In the simulations, $\beta = 2$ and $\sigma _{v}^{2}=1$ are held constant and $\boldsymbol {x}$ is drawn from a standard normal distribution. The parameter space of $\gamma $ ranges from $0$ to $1$ in steps of $0.2$ while $\rho $ takes on values between $0$ and $0.8$ in steps of $0.2$.Footnote ¹⁰ This setup includes a nonspatial DGP without omitted variables bias ($\rho =0$ and $\gamma =0$), nonspatial DGPs with omitted variables bias ($\rho =0$ and $\gamma>0$), SEM DGPs ($\rho>0$ and $\gamma =0$), and SDM DGPs ($\rho>0$ and $\gamma>0$). $\boldsymbol {W}$ is a row-stochastic contiguity matrix based on the queen criterion of adjacency. In contrast to the rook connectivity scheme which links spatial units ordered on a lattice to their direct horizontal and vertical neighbors, the queen criterion additionally connects the units to their diagonal neighbors (e.g., Cliff and Keith Ord Reference Cliff and Keith Ord1981).Footnote ¹¹ The sample sizes specified here contain $49$, $100$, $225$, and $400$ observations distributed on regular grids ($7 \times 7$, $10 \times 10$, $ 15 \times 15$, and $20 \times 20$) to realistically reflect small to medium sized samples frequently encountered in political science.Footnote ¹²

Since the consequences of model misspecification for the estimation of unbiased effect estimates have been studied elsewhere (e.g., LeSage and Pace Reference LeSage and Pace2009; Pace and LeSage Reference Pace, LeSage, Páez, Le Gallo, Buliung and Dall’erba2010; Lacombe and LeSage Reference Lacombe and LeSage2015; Rüttenauer Reference Rüttenauer2019), this Monte Carlo analysis focuses on the ability of the Wald test to identify the true spatial model and differentiate between substantive and residual dependence across a range of alternative DGPs.Footnote ¹³ To this end, I investigate the performance of the Wald test using the four alternative null hypotheses of common factors summarized in Table 1.Footnote ¹⁴

4.2 Performance of the Original Wald Test

Table 2 reports the rejection rates of the four expressions of the null hypotheses of common factors at an $\alpha $-level of $0.05$ across the simulations when the true DGP is that of the SEM model ($\gamma =0$). Since there are no omitted spillovers in this scenario, the common factor restriction holds and the four variants of the Wald test are expected to reject the true null hypothesis in about $5\%$ of the simulation trials with a $95\%$ confidence interval of $[3.65\%; 6.35\%]$.

Table 2 Share of false positives (type I error rates) using asymptotic critical values.

Based on the $\chi ^2$ distribution with $df=1$ and $\alpha = 0.05$, the asymptotic critical value used for all variants of the Wald test and across the different levels of spatial autocorrelation is $\chi ^2_{asym} = 3.841$. The theoretically expected rejection rate across the 1,000 simulation iterations is $5\%$ with a $95\%$ (binomial proportion) confidence interval of $[3.65\%; 6.35\%]$.

Although Section 3 analytically shows that the alternative Wald tests are asymptotically equivalent, their type I error rates differ notably in finite samples. Especially $H_{0}(II)$ but also $H_{0}(IV)$ deviate considerably from the expected error rate. Across all sample sizes, $H_{0}(II)$ is too conservative when $\rho $ is small. Since $\hat {\rho }$ appears in the restriction’s denominator (see Table 1), the restriction has no derivative at zero which violates the assumed continuity of derivatives. However, the Wald test based on $H_{0}(II)$ remains valid as its asymptotic distribution is obtained under the null hypothesis which precludes the problematic value (Gregory and Veall Reference Gregory and Veall1985).Footnote ¹⁵

At the same time, incorrectly rejecting a true null hypothesis might be less problematic in this case since the SDM model derives unbiased impact estimates even if only the residuals are spatially clustered and no substantive spillovers exist (e.g., Elhorst Reference Elhorst2010). The only drawback is that the appropriate SEM specification would be more efficient which might affect inferences regarding the statistical significance of a regressor’s impact. Consequently, it is crucial for any test of the common factor hypothesis to have satisfactory power properties in order to reduce concerns about biased effect estimates.

Against this background, Figure 1 compares the performance of the alternative variants of the Wald test for different levels of correlation between the spatially dependent omitted variable and the included regressor by reporting their power. As the correlation increases, the tests should be more likely to reject the null hypothesis. In order to account for the effects of the sample size and the strength of the interdependence on the performance of the tests, Figure 1 is comprised of 16 panels. In each panel, the horizontal axis depicts the different values of $\gamma $ and the vertical axis shows the observed share of rejections across the simulation trials.

Figure 1 Share of null hypothesis rejections at a nominal significance level of $5\%$.

A brief inspection of Figure 1 already confirms that alternative parameterizations of the null hypothesis—although algebraically equivalent—yield strikingly different results in finite samples. Even with a decently sized sample, there are pronounced differences in the rejection rates of the four Wald tests. While $H_{0}(I)$, which is considered to be the common way to express the restriction (Gregory and Veall Reference Gregory and Veall1986, 204), and $H_{0}(III)$ perform comparatively well in these simulations, the specifications based on $H_{0}(II)$ and $H_{0}(IV)$ have inferior power properties. The remarkably low rejection rates of these specifications of the Wald test increase the likelihood that researchers incorrectly infer the absence of meaningful spillover effects.

Moreover, the behavior of $H_{0}(IV)$ differs greatly from the expectation as its rejection rate initially decreases in almost all parameter settings as $\gamma $ increases. This phenomenon—known as nonmonotonicity in the power function—makes the rejection of the null hypothesis even less likely as the difference between the true DGP and the restriction increases (King and Goh Reference King, Goh, Ullah, Wan and Chaturvedi2002, 256–58).Footnote ¹⁶ In practice, these tests would suggest that the data was generated by a DGP with spatially correlated errors even if there are sizable spillover effects. Researchers would incorrectly conclude that a SEM model or even a nonspatial OLS model appropriately represents the unobservable DGP. Given that these model specifications produce biased impact estimates if a SDM process generated the data, the low rejection rates are highly problematic for substantive inferences.

Although the different variants of the Wald test use the same data and identical parameter estimates, this simulation study shows that, depending on the functional representation of the null hypothesis, they can come to contradictory conclusions regarding the validity of the common factor restriction.Footnote ¹⁷ In fact, Breusch and Schmidt (Reference Breusch and Schmidt1988) analytically show that it is possible to obtain any desired Wald statistic by appropriately specifying the restriction which opens up the possibility to intentionally manipulate the test result (see also King and Goh Reference King, Goh, Ullah, Wan and Chaturvedi2002). Therefore, any search strategy utilizing the Wald test, like the basic general-to-specific approach or the multistep procedure suggested by Elhorst (Reference Elhorst2014a), is subject to this malfunctioning. Since there is no theoretically justified functional representation of the common factor hypothesis and given the strikingly large share of inconsistent inferences across a range of parameter settings, the evidence presented here strongly caution against the use of the standard Wald test based on an asymptotic reference distribution.

4.3 Performance of the Modified Wald Test Based on Bootstrap Critical Values

While the simulations performed here illustrate that the standard Wald test based on asymptotic critical values is unreliable for the identification of the unobservable spatial process, this section investigates whether the application of simulated reference distributions improves the test’s performance. To this end, I use the bootstrap procedure outlined in Section 3.3 and compare the observed Wald statistics based on the different formulations of the common factor restriction to their estimated critical values.

The results reported in Table 3 show that the estimated critical values from the bootstrap approach, $\chi ^{2}_{boot}$, displayed in parentheses not only differ from their asymptotic counterpart $\chi ^{2}_{asym} = 3.841$ on which the original Wald test is based. They also reveal sizable discrepancies between the alternative parameterizations of the null hypothesis. While the estimated critical values for $H_{0}(I)$ and $H_{0}(III)$ are always higher than $\chi ^{2}_{asym}$, the simulated null distributions of $H_{0}(II)$ suggest much smaller critical values for this expression of the null hypothesis in most scenarios.

Table 3 Share of false positives (type I error rates) using bootstrap critical values.

For the different levels of spatial autocorrelation, the median bootstrap critical values $\chi ^{2}_{boot}$ for each variant of the Wald test at the nominal significance level of $5\%$ are displayed in parentheses. Again, the theoretically expected rejection rate across the simulation trials is $5\%\; [3.65\%; 6.35\%]$.

Since the functional expression of the nonlinear common factor restriction determines the Wald statistic’s rate of convergence to the asymptotic $\chi ^{2}$ distribution, estimating critical values for each alternative parameterization improves the empirical size of the Wald test in finite samples. Compared to the original tests based on the asymptotically derived critical value, the observed rejection rates of each of the four variants of the Wald test is closer to the nominal significance level of $5\%$ across all sample sizes. Whereas the observed rejection rate of $H_{0}(I)$ ranges from $6.9\%$ to $11\%$ across the different values of $\rho $ for a sample size of $n=49$ when relying on the asymptotic $\chi ^{2}$ distribution (see Table 2), its corresponding range is narrowed to 5.1%–5.8% when using bootstrap critical values. Similarly, the bootstrap critical values even improve the size of $H_{0}(IV)$ which performed poorly under the asymptotic reference distribution. For $n=49$, basing inferences on the simulated null distribution narrows the range of rejection rates from 2.6%–11.5% to 4.1%–6.8% across the different levels of spatial autocorrelation.

In conclusion, the Monte Carlo evidence presented here demonstrate that using the simulated null distribution as a reference distribution and basing inferences on estimated rather than asymptotically derived critical values ameliorates the problems posed by the Wald test’s lack of invariance to alternative parameterizations of the common factor hypothesis.Footnote ¹⁸ Since the bootstrap critical values account for differences in the convergence rates of the Wald statistics, this modification constitutes a superior alternative that facilitates the empirical assessment of the common factor restriction in spatial regression models. Alternatively, the LR test constitutes another option that is invariant to such reparameterizations (Godfrey and Veall Reference Godfrey and Veall1998).Footnote ¹⁹ Hence, irrespective of the empirical model search strategy employed, researchers should utilize the modified Wald test based on the simulated null distribution or the LR test in order to empirically evaluate the appropriateness of the spatial model employed.