Testing Endogeneity of Spatial Weights Matrices in Spatial Dynamic Panel Data Models

Testing Endogeneity of Spatial Weights Matrices in Spatial Dynamic Panel Data Models

\onehalfspacing

Introduction

\onehalfspacing

In a connected world, it is natural that an outcome of interest arises as a complex of interactions among spatial units in some effective area. In this spirit, Spatial econometrics analyzes the effects of spatial interactions among spatial units on their final economic outcome.

An ordinary factor considered to affect such interactions is the physical distances that represent the adjacency with its neighborhood or farness where the interactions are expected to decay at some specific rate as the distances get far. Thus, an intuitive use in spatial weights matrices ($W$) has been found in the predetermined geography, wherefore exogenous $W$ have been conventional in spatial econometrics (Moran, 1950; Cliff and Ord, 1972; Anselin, 1988, 2001; Ertur and Koch, 2007; Kelejian and Prucha, 2010; Elhorst, 2010, 2014; Lee and Yu, 2010; Dogan and Taspinar, 2013).

Empirical literature, however, assert that economic activities such as international trade also carry out knowledge spillover effects (Frankel and Rose, 1998; Baxter and Kouparitas, 2005; Ditzen, 2018). Inspired by this, a special $W$ of the bilateral trade flow is used (Ertur and Koch, 2011; Ho, Wang, and Yu, 2013). For more examples where economic factors form $W$, see Conley and Ligon (2002), Conley and Topa (2002), Parent and LeSage (2008), and Skevas, Skevas, and Cabrera (2021).

$W$ can be endogenous, however, if the space is intrinsically physical or economic (Pinkse and Slade, 2010). This is because economic variables are broadly and closely connected one another and thus the random shock in one variable might be correlated to that in the outcome of interest. If they are significantly dependent each other, then the exogenous assumption on $W$ may not be valid any more and inferences from the ordinary spatial autoregressive (SAR) estimators and test statistics will be invalid.

To tackle this problem, Qu and Lee (2015) modeled the source of endogeneity of $W$ and proposed the maximum likelihood estimator (MLE) that controls endogeneity of $W$, which is known to be consistent and asymptotically normal distributed in cross section (Qu and Lee, 2015). In a spatial dynamic panel data (SDPD) model, the score function is not centered around zero. Also, the ML estimator has a bias with an order of $O(max\{n^{-1},T^{-1}\})$, even with large $n$ and large $T$, due to the two-way fixed effects (Yu, Jong, and Lee, 2008; Lee and Yu, 2010; Qu, Lee, and Yu, 2017). This can be rectified by the bias correction method (Yu, Jong, and Lee, 2008; Lee and Yu, 2010; Qu, Lee, and Yu, 2017; Bera, Dogan, Taspinar, and Leiluo, 2019). I thus firstly introduce the bias-corrected score function centered around zero. The initial condition issue with finite $T$ is not of concern here with large $T$ setting.

I further adjust score functions to avoid the over-rejection of the null hypothesis under presences of local misspecification (Davidson and Mackinon, 1987; Saikkonen, 1989) in contemporaneous dependence over space, dependence over time, or spatial time dependence. This adjustment is in agreement with the existing Robust Rao's Score (RS) test literature: see Bera and Yoon (1993); Dogan, Taspinar, and Bera (2018); Bera, Dogan, and Taspinar (2018, 2019); Bera, Dogan, Taspinar, and Leiluo (2019); Bera, Bilias, Yoon, Taspinar, and Dogan (2020). The RS test is robust in the sense that its asymptotic distribution is a central chi-square distribution regardless of local parametric misspecifications. A Monte Carlo simulation shows nice finite sample properties in size and power.

Another advantage of the RS test is that it is computationally efficient since it only requires the restricted MLE under the null where the parameters above are assumed to be zero, reducing the spatial dynamic panel data models to the simple fixed-effects model. As illustrated in Section 6, the elapsed time (in seconds) for the Robust RS test is less than that for Conditional Lagrange Multiplier test (Qu and Lee, 2015; Cheng and Fei Lee, 2017).

Literature on testing and model specifications have been mostly considered for cross-sectional spatial models (For example, see Anselin, 1988, 2001; Kelejian and Robinson, 1992; Anselin et al., 1996; Baltagi and Li, 2001; Yang, 2010; Bera et al., 2018, 2019; Dogan et al., 2018), whereas only few studies are found for the spatial static panel models (Baltagi et al., 2003, 2007, 2009; Baltagi and Liu, 2008; Debarsy and Ertur, 2010; Baltagi and Yang, 2013) and for the SDPD models (Yang, 2016; Taspinar et al., 2017; Bera et al., 2019). The majority of these studies lie on testing the presence of spatial dependence. To my best knowledge, testing endogeneity of $W$ has been developed only in Bera et al. (2018) in the cross-sectional spatial models and no such test has been developed yet in the SDPD models. A valid score test to test the endogeneity of $W$ is developed in this paper resolving two challenges of uncentered score function and presence of local misspecifications in the SDPD models.

The rest of the paper is organized as follows. In Section 2, I review the model specification and adopt the assumptions in Qu, Lee, and Yu (2017). In Section 3, I review the ML estimation which maximizes the concentrated log-likelihood. In Section 4, I develop the robust Rao's score test for testing endogeneity of $W$. In Section 5, I derive the explicit forms of the test statistic. In Section 6, I conduct a Monte Carlo simulation to explore its finite sample properties and an empirical illustration is provided using Penn World Table version 6.1. In Section 7, I conclude. The proofs of Propositions are provided in the Appendix.

Model specification

Following Jenish and Prucha (2009, 2012), consider a generalized spatial processes on a possibly unevenly spaced lattice. Let $\vert\vert Z \vert\vert_p=[E|Z|^p]^{1/p}$ when the absolute $\text{$p^{th}$}$ moment exists.

I adopt the following assumptions as in Qu, Lee, and Yu (2017) on the spatial setting of observations.

Assumption 1 allows the asymptotic inference on increasing domain under the space-time NED. Without loss of generality, assume $\rho_{ct}=1$ for all $t=1,\dots,T$, which means their least physical distance is one unit apart given each period $t.$

Let $\{(\epsilon_{l(i,t)}, v_{l(i,t)}): l(i,t) \in D_L, i \in \mathbb{N}, t \in \mathbb{T}\}$ be a random field of error terms defined on a probability space $(\Omega, \mathcal{F}, P)$, where $D_L \subset D$ is a finite set and $D$ satisfies Assumption 1. To simplify the notation, let $(\epsilon_{l(i,t)}, v_{l(i,t)})$ be denoted by $(\epsilon_{it}, v_{it})$, where $\epsilon_{it}$ is the $p \times 1$ column vector for $t=1,\dots,T$ formulated from the $i^{th}$ row of $n \times p$ matrix $\varepsilon_{nt}$, and $v_{it}$ is the $i^{th}$ element of $n \times 1$ vector $V_{nt}$, for all $t = 1, \cdots, T.$ A SDPD model with individual and time fixed effects for $n$ cross-sectional units can be specified as

where $Y_{nt}=(y_{1t},\dots,y_{nt})'$ is a $n \times 1$ vector of observations on the dependent variable, $W_{nt}=(w_{ij})_{t}$ is a $n \times n$ spatial weights matrix with zero diagonal elements, $\lambda_0$ denotes the autoregressive parameter, $X_{1nt}$ is a $n \times k_1$ matrix of deterministic explanatory variables that are bounded in absolute value, and $\beta$ is a $k_1 \times 1$ vector of parameters. $c_{n10}$ is a $n \times 1$ column vector of individual fixed effects and $\alpha_{t10}$ is the $t^{th}$ element of $T \times 1$ time fixed effects vector $\alpha_{T10}$, respectively, and $V_{nt}=(v_{1t}, \cdots, v_{nt})'$ is a $n \times 1$ vector of error terms with zero mean and variance $\sigma_0^2$. I also generalize $W_{nt}$ to be time varying and endogenous by modeling a construction of $W_{nt}$ as a bounded function of $Z_{nt}$ such that

where $h_{ij,t}(\cdot)$ is a nonnegative and uniformly bounded function and $Z_{nt}=(z_{1t}, \cdots, z_{nt})'$ is a $n \times p$ matrix of dependent variables with $z_{it}=(z_{1it}, \cdots, z_{p it})'$ being a $p \times 1$ vector and $\kappa_0$ is its associated parameter. $X_{2nt}$ is a $n \times k_2$ matrix of deterministic explanatory variables whose elements are bounded in absolute value, $\Gamma_0$ is a $k_2 \times p$ matrix of associated coefficients, and $\varepsilon_{nt}=(\epsilon_{1t}, \cdots, \epsilon_{nt})'$ is a $n \times p$ matrix of errors, where $\epsilon_{it}=(\epsilon_{it,1}, \cdots, \epsilon_{it,p})'$ is a $p \times 1$ vector. $c_{n20}$ is again a $n \times 1$ column vector of individual fixed effects and $\alpha_{t20}$ is the $t^{th}$ element of $T \times 1$ time fixed effects vector $\alpha_{T20}$, respectively. $1_n$ is a $n \times 1$ vector of ones. The initial values in $Y_{n0}$ are assumed to be observable.

For the SDPD model in (1), it might be convenient to adopt the subscript $n$ to simplify notation. Due to non-linearity in $W^{'}_{nt}s$, one may work on the random fields for $\{y_{it}\}$ and $\{(v_{it},\epsilon'_{it})\}$. The settings are similar as in Jenish and Prucha (2012) so that one may apply the spatial-time LLN. At time $t=0,$ the $n$ units are located in the Euclidean space $\mathbb{R}^{d+1}$. At time $t=1$, the $n$ units shift vertically upward to an affine plane parallel to $\mathbb{R}^d$ planes. As time passes, it keeps shifting upward each time by one time unit. For the variables in past periods (history), one may shift the plane $\mathbb{R}^d$ at $t=0$ downward to an affine plane at $t=-1$, and so on. A sample of individuals with $n$ units over $T$ periods are located in the $\mathbb{R}^{d+1}$ space. Using the maximum metric (Jenish and Prucha, 2012)

each observation indexed by $(i,t)$ will be located at $\ell(i,t)\in\mathbb{R}^{d+1}$ where $\ell(i,0)$ is the physical location of an individual at time 0. Any two individuals located in the $\mathbb{R}^d$ plane at $t=0$ are assumed at least one unit apart. For each spatial unit $i,$ the location $\ell(i,t)$ is always one unit apart from $\ell(i,t-1)$ with respect to the time dimension. There are $L$ locations in $\mathbb{R}^{d+1}$. For the spatial weight matrix, $w_{ij,nt} \neq 0$ if an individual at spatial unit $i$ at time $t$ directly links to $j$ at time $t$; $w_{ij,n,t-1}\neq 0$ if $i$ links to $j$ at time $t-1$. Any individual at time $t$ does not directly link to any future loads from anyone including itself $t+\tau$ for $\tau \geq 1,$ and not the past $t-\tau$ with $\tau \geq 2$ but indirect links are allowed at $t-\tau$ for $\tau \geq 2$. Note that $i$ indirectly links to itself in past periods. This is a mapping of $(i,t)$ to a location $\ell$ in $D_L$, i.e., $\ell=\ell(i,t) \in D_L$ with $L=nT.$ Since $y_\ell$ corresponds to a $y_{it}$, one may define the spatial NED process with a base $\sigma$-field generated by $\{v_{it},\epsilon_{it}'\}$ as follows:

The $Y=\{y_{it}\}$ is a NED if $\vert\vert y_{\ell,L}-E(y_{\ell,L}\vert\mathcal{F}_{\ell,L}(s))\vert\vert_p \leq d_{\ell,L}\Psi(s)$ for some $\Psi(s),$ which goes to zero as $s \rightarrow \infty.$

Note that since $W_{nt}$ is a function of $\varepsilon_{nt}$, $W_{nt}$ is endogenous in the estimation equation if $\sigma_{v\epsilon 0} \neq 0$, which leads to a biased estimator due to the endogeneity. From the joint distribution of $(v_{it},\epsilon_{it}')',$ one has $E(v_{it}|\epsilon_{it})=\epsilon_{it}'\delta_0$ where $\delta_0=\Sigma_{\epsilon 0}^{-1}\sigma_{v\epsilon 0}$ and $Var(v_{it}|\epsilon_{it})=\sigma_{\xi 0}^2$ where $\sigma_{\xi 0}^2=\sigma_{v0}^2-\sigma_{v \epsilon 0}'\Sigma_{\epsilon 0}^{-1}\sigma_{v \epsilon 0}.$ Let $\xi_{nt}=V_{nt}-\epsilon_{nt}\delta_0.$ Since the expectation of $\xi_{nt}$ conditional on $\varepsilon_{nt}$ is zero, $\xi_{nt}$ is uncorrelated with $\varepsilon_{nt}.$ Hence $\eqref{eq1}$ can be represented as

where $E(\xi_{nt}|\varepsilon_{nt})=0$ and $Var(\xi_{nt}|\varepsilon_{nt})=\sigma_{\xi 0}^2 I_n,$ and the elements of $\xi_{nt}$ are iid across $i$ and $t$. Note that $(Z_{nt}-Z_{n,t-1}\kappa_0-X_{2nt}\Gamma_0-c_{n20}-1_n\alpha_{t20}')$ are control variables in $Y_{nt}$ to control the endogeneity of $W_{nt}$.

I additionally assume the followings.

Assumption 3 requires large $n$ and large $T$ case. As $T$ is large under assumption 3, the initial condition problem would not be an issue here. Assumption 4 excludes an issue of multicollinearity. Also $X_{1nt}$ and $X_{2nt}$ are allowed to overlap and have the common variables because identification is not of interest for our purpose. Assumption 5 imposes features of $W_{nt}$ such that: (i) allows time-varying $W_{nt}$ while the physical distance is fixed over time; For technical purpose, (ii) states that two units are not considered spatially connected if their exogenous distance exceeds a certain threshold, even though their economic/social factors are. This is a popular setting in empirical studies; (iii) is a condition on $h(\cdot)$ so that $z_{it}$ and $z_{jt}$ determine the time NED property for $w_{ij,t}$. Assumption 6 guarantees the stability of the dynamic process by controlling the magnitude of spatial interaction of $W_{nt}$ matrix, which is defined in Section 3. The parameter space is characterized under Assumption 7. \\

ML estimation

Following Qu, Lee, and Yu (2017), one may put the model $\eqref{eq1}$ into a big matrix form. Denote

and $X_{1L}, \epsilon_L, \xi_L, Z_L$ and $X_{2L}$ are defined similarly. Also, denote

with $W_{1L}=

,$ $W_{2L}=

,$ \\ $W_{3L}=

,$ \\

so that

where $\eta=(\lambda, \gamma, \rho)'$ and $\eta_0$ is its true value. Therefore, the model in $\eqref{eq1}$ and in $\eqref{eq2}$ can be written in the matrix form as:

where $c_{1L} = 1_T \otimes c_{1n}$, $c_{2L}=1_T \otimes c_{2n}$. From Assumption 2, $V_L=\varepsilon_L \delta_0 + \xi_L$ where $\xi_L \sim N(0,\sigma_{\xi_0}^2 I_n)$ conditional on $\varepsilon_L$. Hence,

where $ \sbox{\myboxA}{$\m@thc_{1L0}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.75\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi=c_{1L0}-c_{2L0}\delta_0$ and $ \sbox{\myboxA}{$\m@th\alpha_{1L0}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.75\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi=\alpha_{1L0}-\alpha_{2L0}\delta_0.$ \\

For asymptotic analysis, it is useful to have the likelihood function presented in matrix form involving $W_L({\eta}).$ Let $\theta=(\lambda,\phi_1',\delta',\phi_2',\alpha',\sigma_{\xi}^2)',$ where $\phi_1=(\gamma,\rho,\beta')'$, $\phi_2=vec(\Phi_2)$ with $\Phi_2=(\kappa',\Gamma')',$ and $\alpha$ is the $J \times 1$ column vector of distinct elements in $\Sigma_\epsilon.$ Corresponding to $\phi_1$ and $\phi_2$, I define $R_{nt}=[Y_{n,t-1}, W_{n,t-1}Y_{n,t-1}, X_{1nt}]$, $K_{nt}=[Z_{n,t-1},X_{2nt}]$ and their matrix form $R_L$ and $K_L$ with row dimension $L=nT.$ I denote $S_L(\eta)=I_L-W_L(\eta)$, $S_L=I-W_L(\eta_0)$ and correspondingly $G_{jL}(\eta)=W_{jL}S_L^{-1}(\eta)$ and $G_{jL}=W_{jL}S_{L}^{-1}$ for $j=1, 2, 3$, where $I_L$ is the $L \times L$ identity matrix. I also denote $Q_{1L}(\theta_0,c_{1L0},\alpha_{1L0})=G_{1L}(X_{1L}\beta_0+\ell_0(\gamma_0,\rho_0)+\varepsilon_L\delta_0+c_{1L0}+\alpha_{1L0})$. One may derive the log likelihood function from Assumption 2 as it only imposes assumptions only on the first and second moments but third and fourth ones. From Assumption 2 and equation $\eqref{eq4}$, the log likelihood function in the matrix form is

As in Lee and Yu (2010), I use the two orthogonal projectors for taking time deviation and cross sectional deviation from their means, $J_n = I_n-\frac{1}{n}1_n 1_n'$ and $J_T = I_T - \frac{1}{T}1_T 1_T'$. Let $J_L=J_T \otimes J_n$. The concentrated log-likelihood function is then

The robust RS test statistics needs the asymptotic distribution of score functions derived as follows (Qu, Lee, and Yu, 2017):

Due to the two way fixed effects, the score function of $lnL_L^c(\theta)$ can be decomposed of the unbiased score function and bias term (Qu, Lee, and Yu, 2017) as follows:

where

and

The first term in $\frac{\partial lnL_L^c(\theta_0)}{\partial \theta}$, $\frac{\partial \ln L_{1,L}^c(\theta_0)}{\partial \theta}$, has zero mean and is asymptotically distributed with further assumption of the third and fourth moment of $\xi_{it}$ conditional on $\epsilon_{it}$.

Proof. See Appendix 1. \\

Furthermore, the bias term $\Delta_L$ is composed of two types of biases such that

where $a_{1,\theta_0}$ and $a_{2,\theta_0}$ are of $O(1)$ given as

and

The bias $a_{1,\theta_0}$ is from individual effects and $a_{2,\theta_0}$ is due to time effects.

Proof. See Appendix 1.

Proof. See Appendix 1. \\

Now let $\hat{\theta}_{nT}=\underset{\theta \in \Theta}{\text{argmax}} \hspace{1mm}\ln L_L^c(\theta)$ be the MLE. Then the asymptotic distribution of $\hat{\theta}_{nT}$ is (Qu, Lee, and Yu, 2017)

where $b_{1,\theta_0,nT}=\ell_{nT,\theta_0}^{-1}a_{1,\theta_0}$ and $b_{2,\theta_0,nT}=\ell_{nT,\theta_0}^{-1}a_{2,\theta_0}$ are the bias terms with their orders $O(1)$. Note that $\hat{\theta}_{nT}$ has the bias $\frac{1}{T}b_{1,\theta_0,nT}$ and $\frac{1}{n}b_{2,\theta_0,nT}$ with $\frac{n}{T}\rightarrow k$, which implies the bias term is of order $O\left(\text{max}\left(\frac{1}{n},\frac{1}{T}\right)\right).$\footnote{Indeed, a bias term is of order $O\left(\text{max}\left(\frac{1}{n},\frac{1}{T}\right)\right)$ in any cases if $\frac{n}{T}\rightarrow k<\infty, \frac{n}{T}\rightarrow 0,$ or $\frac{n}{T}\rightarrow \infty$.} A bias-corrected estimator can be defined as

where $\hat{B}_{1,nT}=[\ell_{nT,\theta}^{-1} \cdot a_{1,\theta}] \vert _{\theta=\hat{\theta}_{nT}}$ and $\hat{B}_{2,nT}=[\ell_{nT,\theta}^{-1} \cdot a_{2,\theta}]\vert _{\theta=\hat{\theta}_{nT}}.$ With further assumption, $\hat{\theta}_{nT}^1$ is properly centered.

Proof. See Appendix 1.

The Robust Rao's Score Tests for Endogeneity of Spatial Weights Matrices

Consider the following partition of the parameter vector $\theta = (\delta', \eta', \omega')'$, where $\omega=(\beta',\phi_2',\alpha',\sigma_{\xi}^2)'$. In this partition, $(\delta',\eta')'$ represent the parameter vector of interest, while $\omega$ is pure nuisance parameters. From equation $\eqref{eq4}$, I may design a test for endogeneity of $W_{nt}$ by testing if $\delta=0.$ The null hypothesis then can be stated as $H_0: \delta_0 = 0$ versus the alternative hypothesis as $H_1: \delta_0 \neq 0$. Let $I(\theta_0)=\lim\limits_{n,T \rightarrow \infty} \ell_{nT,\theta_0}$ and its estimator $I(\theta)=-\frac{1}{nT}\frac{\partial^2 lnL_L^c(\theta)}{\partial \theta \partial \theta'}.$ Let $L_{\Psi}(\theta)=\frac{1}{nT}\frac{\partial lnL_L^c(\theta)}{\partial \Psi}$ and $L_{\Psi\Psi}(\theta)=\frac{1}{nT}\frac{\partial^2 lnL_L^c(\theta)}{\partial \Psi \partial \Psi'},$ where $\Psi \in\{\delta,\eta,\omega\}.$ I consider the following partition of $I(\theta)$:

and the partition of the bias term $\Delta_L(\theta)$:

such that

as shown in Remark 1. I denote $I\equiv I(\theta_0)$ and let $\tilde{\theta}=(0',0',\tilde{\omega}')'$ be the restricted MLE when the joint null $H_0^{\delta,\eta}: \delta_0=0, \hspace{1mm}\eta_0=0$ holds. I first consider the case of $H_0^{\delta}: \delta_0=0$ when $H_0^{\eta}: \eta_0=0$ holds. Since $L_{\delta}(\theta)$ is not centered around zero as in Proposition 1, I introduce the bias-corrected score function. As derived in the Appendix 3 (Proposition 4 proof), the first-order Taylor expansion of $L_\delta(\tilde{\theta})$ gives the following equation:

which implies

where $I_{\delta \cdot \omega}:=I_{\delta\delta}-I_{\delta\omega}I_{\omega \omega}^{-1}I_{\omega\delta}$. The bias-corrected score function $C_\delta(\tilde{\theta})$ is therefore of the form

where the asymptotic distribution of $C_\delta(\tilde{\theta})$ is centered around zero. The standard RS test statistic is

Now I investigate the asymptotic distribution of $RS_\delta(\tilde{\theta})$ under the sequences of local alternatives of $H_a^\delta: \delta_0=\zeta/\sqrt{nT}$ and $H_a^\eta: \lambda_0=\nu/\sqrt{nT}$, where $\zeta$ and $\nu$ are bounded constant vectors. From the first-order Taylor expansion of $L_\delta(\tilde{\theta})$,

By Proposition 1,

Then, $\eqref{eq11}$ and $\eqref{eq12}$ imply that

where $I_{\delta\eta \cdot \omega}=I_{\delta \eta}-I_{\delta \omega}I_{\omega \omega}^{-1}I_{\omega \eta}.$ That $\sqrt{nT}C_\delta(\tilde{\theta})$ is asymptotically not centered around zero yields that $RS_\delta$ follows $\chi^2$ distribution with non-zero noncentrality parameter. This is a problem to be resolved because the test statistic leads to over-rejection of the null hypothesis (Davidson and Mackinnon, 1987; Saikkonen, 1989). A robust version of RS test can be constructed by adjusting $C_\delta(\tilde{\theta})$ so that the RS test statistic is centered around zero: see Bera and Yoon (1993); Dogan, Taspinar, and Bera (2018); Bera, Dogan, and Taspinar (2018, 2019); Bera, Dogan, Taspinar, Leiluo (2019); Bera, Bilias, Yoon, Taspinar, and Dogan (2020). The asymptotic behavior of the standard and robust RS test statistics under $H_0^\delta$ and $H_a^\delta$ are provided in the following Proposition.

Proof. See Appendix 3. \\

This result indicates that $RS_\delta^*(\tilde{\theta})$ is a robust test since it fixes the over-rejection rate of the null hypothesis. That is, it gives asymptotically correct size with its asymptotic null distribution being centered chi-square distribution under the sequence of alterantives $H_a^\eta: \eta_0=\nu/\sqrt{nT}$. Also note that under $H_a^\delta$ and $H_0^\eta$, the result shows $RS_\delta^*(\tilde{\theta}) \xrightarrow{d} \chi^2_{p}(\varphi_2)$ and $RS_\delta(\tilde{\theta}) \xrightarrow{d} \chi^2_{p}(\varphi_1)$ with $\varphi_1-\varphi_2 \geq 0$, indicating $RS_\delta^*(\tilde{\theta})$ has less asymptotic power than $RS_\delta(\tilde{\theta})$ when there is no local misspecification, i.e., when $\eta_0=0.$ That is, one pays the premium to have less power if no presence of local misspecification is found.

A different approach on testing endogeneity of $W$ could be implemented using the Conditional Lagrange Multiplier (CLM) test following the framework of Qu and Lee (2015) and Cheng and Fei Lee (2017). Let $\hat{\theta}=\underset{\theta: \; \delta=0}{\text{argmax}}\ln L_L^c(\theta)$ be the restricted ML estimator under $H_0^\delta$. Then a valid CLM is formulated as

where $\Psi=(\lambda,\phi_1^{'},\phi_2^{'},\alpha',\sigma_\xi^2)'$. Note that $LM(\hat{\theta})$ is asymptotically central chi-squared under $H_0^\delta$ and has the same form of the non-robust $RS_\delta(\theta)$ test at $\hat{\theta}$. Notably, the robust RS test is computationally efficient in the sense that it does not require for $\eta=(\lambda,\gamma,\rho)'$ to be estimated, while CLM requires the restricted ML estimators obtained under $H_0: \delta_0=0$, i.e., $(\lambda, \phi_1', \phi_2', \alpha', \sigma_\xi^{2'})$ need to be estimated. A comparison over elapsed times between two tests is reported in Section 6.

The Test Statistics

Now I explicitly derive the robust RS test statistics for the employed model (Qu, Lee, and Yu, 2017) using Proposition 4. One may have the following hypotheses:

1. $H_0^\delta: \delta_0=0$ and $H_0^\eta: \eta_0=0.$ \\ 2. $H_a^\delta: \delta_0=\zeta/\sqrt{nT}$ and $H_a^\eta: \eta_0=\nu/\sqrt{nT}$.

The first joint null hypothesis tests if $W_{nt}$ is endogenous when there is no local misspecification. I then consider the asymptotic distribution of the test statistic under the alternative hypothesis of the endogeneity parameter $\delta$ and local presence of misspecification in $\eta$.

Recall that the restricted MLE is denoted by $\tilde{\theta}=(0',0',\omega')'$ under the joint null hypothesis, where $\omega=(\beta',\phi_2',\alpha',\sigma_\xi^2)'$. It is highlighted that the robust RS test has computational advantage in the sense that it only requires $\tilde{\theta}$ and one does not need to estimate other parameters under the alternative hypothesis. I firstly consider the asymptotic distribution of the test statistic under the joint null hypothesis $H_0^\delta$ and $H_0^\eta$. The concentrated log-likelihood function at $\tilde{\theta}$, $\ln L_L^c(\tilde{\theta})$, turns down to

which is decomposed into two unrelated components as $\ln L_L^c(\tilde{\theta})=\ln L_L^{C1}(\tilde{\theta})+\ln L_L^{C2}(\tilde{\theta})$, where

Using the bias-corrected ML estimator in (ref), one obtains the restricted ML estimators and the residuals, $\xi_L(\tilde{\theta})=Y_L-X_{1L}\tilde{\beta}$ and $\varepsilon_L(\tilde{\theta})=Z_L-K_L\widetilde{\Phi_2}$. I now provide explicit expressions for $RS_\delta(\tilde{\theta})$ and $RS_\delta^*(\tilde{\theta})$. The test statistic requires the score functions with respect to $\delta$ and $\eta$. The score functions evaluated at $\tilde{\theta}$ are given as

From the information matrix in equation (8) as well as the second order score functions using Corollary 1, one can find the consistent estimator for the information matrix as

where

Remark that the estimator for the information matrix at $\tilde{\theta}$, $I(\tilde{\theta})$, forms a block diagonal matrix with respect to $(\phi_2',\alpha)'$. Thus one may regard $\omega=(\beta',\sigma_\xi^2)$. Hence the estimators for the information matrices necessary for computing the test statistic under the joint null $H_0^\delta$ and $H_0^\eta$ are

where $(A)_{(i,j:k)}$ represents the entries in a matrix $A$ located in $i^{th}$ row, $j^{th}$ to $k^{th}$ columns.

The bias terms evaluated at $\tilde{\theta}$ are given as

where

and

with the order of $\theta=(\lambda,\phi_1,\delta,\phi_2,\sigma_\xi^2,\alpha)$. Denoting $\Delta_1=\sqrt{\frac{n}{T}}a_{1,\theta_0}$ and $\Delta_2=\sqrt{\frac{T}{n}}a_{2,\theta_0}$, the bias terms necessary for computing the test statistics under the joint null are

Given the above quantities, one can obtain the bias-corrected score functions of $C_\delta(\tilde{\theta})$ and $C_\eta(\tilde{\theta})$.\footnote{The explicit forms are provided in Appendix 2.} Then the standard RS test statistic is

while the robust RS test statistic in Proposition 4 is provided as

where $K(\tilde{\theta})=I_{\delta\eta \cdot \omega}(\tilde{\theta})I_{\eta \cdot \omega}^{-1}(\tilde{\theta})$ and $C_\delta^*(\tilde{\theta})=C_\delta(\tilde{\theta})-K(\tilde{\theta})C_\eta(\tilde{\theta})$, respectively. One may observe that two corrections are made in the score function and its variance with $K(\tilde{\theta})$ in order to manage to rectify over-rejection of the null hypothesis due to the local presence of misspecification in another testing parameter. Note that if $K(\tilde{\theta})=0$ or there is no local misspecification, then the robust RS test is simply the standard RS test.

Now the quantities to compute the CLM test (Qu and Lee, 2015; Cheng and Fei Lee, 2017) are provided. Let $\theta=(\delta',\Psi')'$ where $\Psi=(\lambda,\phi_1^{'},\phi_2^{'},\sigma_\xi^{2},\alpha^{'})'$ and $\hat{\theta}=\underset{\theta:\;\delta=0}{\text{argmax}}\ln L_L^c(\theta)$ be the restricted ML estimator under $H_0^\delta$. Remark that the estimator for the information matrix at $\hat{\theta}$, $I(\hat{\theta})$ forms a block diagonal matrix with respect to $(\phi_2^{'},\alpha')'$, where one thus may regard $\Psi=(\lambda,\phi_1^{'},\sigma_\xi^2)'$. The concentrated log-likelihood under $H_0^\delta$ is then

The centered score function at $\hat{\theta}$ is

with

Monte Carlo study and Empirical illustration

Monte Carlo simulation

I run a Monte Carlo simulation to explore the finite properties of the test statistic. It is conducted for 1,000 times and the sample size for individuals and time periods are set by four cases where $n$ is relatively large than $T$ or vice versa, followed by (approximately) doubled sample size, respectively: (i) $n=100, T=10$; (ii) $n=196, T=20$; (iii) $n=9, T=100$; and (iv) $n=16, T=200$. The type I error is set as 0.05 and the dimension of $Z$ is set as $p=1$. The parameters for the main and auxiliary equations follow the setup in Qu, Lee, and Yu (2017): $\beta_0=1, \kappa_0=0.2, \Gamma_0=0.3, \alpha_0=1$, respectively. Data are generated by (1) and (2). The magnitude of local misspecification in $\eta=(\lambda,\gamma,\rho)$ increases by 0.05 from 0 to 0.3 and $\delta_0$ is set from 0 to 0.2, which also increases by 0.05. The initial values of $Y_0, Z_0$ and the deterministic explanatory variables of $X_{1L}, X_{2L}$ are generated from independent standard normal distributions. To generate the two-way fixed effects, data generated from the multivariate normal distribution of nonzero correlation, 0.5 for both fixed effects, with $X_{iL}$ are averaged over time for each spatial unit ($C_{iL0}$) or over spatial unit for each time ($\alpha_{iL0}$), $i=1,2$. The joint pdf of the disturbance terms, $(v_{it},\epsilon_{it})$, follows bivariate normal distribution of $N\left(

,

\right).$ In summary, we have

I also follow the setup for the spatial weight matrices $W_{nt}$ as in Qu, Lee, and Yu (2017), which is generated by Hadamard product of the physical contiguity and the economic distance: $W_{nt}=W_n^d \circ W_{nt}^e$, i.e. $w_{ij,nt}=w_{ij}^dw_{ij,nt}^e$ and is row normalized afterwards. I explore two kinds of contiguities of Queen & Rook matrices for $W_n^d$ which characterize different adjacency: Queen allows contiguities over edges & corners, whereas Rook only allows edges. The economic distance $W_{nt}^e$ is generated by $w_{ij,nt}^e=1/|z_{it}-z_{jt}|$ if $i \neq j$ or zero otherwise. Table 1 summarizes our notations on the test statistics. The notation of A:B:C represents that the parameter starts from A to C, increased by B. I will estimate the size and power of the test statistics under the local misspecification in $\eta_0=(\lambda_0,\gamma_0,\rho_0)$ when $n$ is larger than $T$: $(n,T)=(100,10)$ or the reverse: $(n,T)=(9,100)$, and finally when they get doubled in each case: $(n,T)=(196,20)$ and $(n,T)=(16,200)$. Here, $n$ is set as a square number so that the Queen & Rook matrices could be constructed accordingly.

The size of the test statistics is summarized in Figure 1 and Table 2. First consider the case when $n$ is larger than $T$. For example, when $(n,T)=(100,10)$, the size of the standard RS generally increases as the magnitude of local misspecification in $\eta_0=(\lambda_0,\gamma_0,\rho_0)$ increases, both in Queen and Rook $W_{nt}^{d'}$s. In particular, the largest increase in size is found in the presence of local misspecification in the contemporaneous spatial dependence ($\lambda_{0}$), whereas those in the spatial time dependence ($\rho_{0}$), and dependence over time ($\gamma_{0}$) slightly increase the size. Results even get clearer when both $n$ and $T$ get doubled. When $(n,T)=(196,20)$, the size explodes under the presence of local misspecification in $\lambda_{0}$. Meanwhile, one may find that the robust RS test stay settled around 0.05 and its performance improves as $n$ and $T$ increase.

Now consider the other case when $T$ is larger than $n$. For example, consider $(n,T)=(9,100)$. Results show that increase in size is remarkable for the presence of local misspecification in $\lambda_{0}$ for the standard RS test as well as decrease in size for the robust RS test for the presence of local misspecification in $\gamma_{0}$, both in Queen and Rook $W_{nt}^{d'}$s. However, the size of the Robust RS test get settled as $n$ and $T$ get doubled, $(n,T)=(16,200)$, while the standard RS test shows even larger increase in size under the local misspecification of $\lambda_{0}$. This result again supports a nice performance of the robust RS test as $n$ and $T$ increase.

The power of the test statistics is summarized in Figure 2 & 3 and in Table 4 & 5. When there is no local misspecification, i.e., $\eta_{0}=(0,0,0)$, the power of the standard RS test is higher than that of the robust RS test, indicating the premium one pays of losing little power when adjusting the standard RS test so that one rectifies the over-rejection of the null hypothesis. For $(n,T)=(100,10)$ or $(n,T)=(9,100)$, the power sharply increases to 1 as $\delta_0$ increases to 0.2. As shown in Figure 2 & 3 and Table 4, the power of the robust RS is as good as the standard RS without local misspecification for all $\delta_0 \in \{0.05,\;0.1,\;0.15,\;0.2\}$, i.e., one does not lose power in using the robust RS. Similar results are found when $n$ and $T$ get doubled, i.e., $(n,T)=(196,20)$ or $(n,T)=(16,200)$ in Table 5.

Now the numerical values of the robust RS test- and the conditional LM test statistics are presented in Table 6 & 7, showing similar results as $n$ and $T$ increase. The robust RS test is expected to have less cost in computation since it only requires the ML estimation of the simple fixed-effects model for panel data, whereas the conditional LM test requires that of the SDPD model. As expected, the elapsed time (in seconds) for the Robust RS test is less than that for Conditional LM test (Table 3).

Empirical illustration

Using the Penn World Tables (PWT version 6.1), I introduce how to test the endogeneity of the spatial weights matrices ($W$) before regular estimation. As in Ertur and Koch (2007), I measure the variables in Solow-Swan growth model of the logarithm of savings (lns), the logarithm of the growth of the working-age population (ages 15 to 64) summed up with the growth rate and the depreciation rate in capital (ln(n+g+$\delta$)) in year 1960-1995. I suppose that $g+\delta=0.05$ as in MRW (1992) and Romer (1989). As a general sense, I consider SDPD growth model to capture the contemporaneous dependence over space, dependence over time, or spatial time dependence.

Now I suspect that $W_{nt}$ is a function of the share of gross consumption in GDP (kc) and real Gross Domestic Income (GDI) for terms of trade changes (rgdptt). As mentioned in Qu, Lee, and Yu (2017), identification is not an issue here and thus one is allowed to have $X_{1nt}$ and $X_{2nt}$ share the common variables. The equations are therefore

where $t=1,\dots,T$, $y_{nt}$ is the real income per worker at time $t$, $X_{1nt}=(lns_{nt},ln(n_{nt}+g+\delta),W_{nt}lns_{nt},W_{nt}ln(n_{nt}+g+\delta))'$, $\beta=(\beta_1,\beta_2,\beta_3,\beta_4)'$, and $Z_{nt}\in \{kc_{nt},rgdptt_{nt},(kc_{nt},rgdptt_{nk})\}$, $\kappa$ is the associated parameter, $X_{2nt}=(lns_{nt},ln(n_{nt}+g+\delta))'$, and $\Gamma=(\Gamma_1,\Gamma_2)$. The equations above can be augmented as

where $\xi_{nt}\sim N(0,\sigma_\xi^2I_{n})$. Recall that $(Z_{nt}-Z_{n,t-1}\kappa-X_{2nt}\Gamma-c_{n2}-1_n\alpha'_{t2})$ is our control variables for the endogenous $W_{nt}$.

Now I use the test statistic, the Robust RS test, to determine if $W_{nt}$ is exogenous, i.e., $H_0^\delta: \delta_0=0$. The results are presented in Table 8. One may find that Robust RS test statistic is far smaller than Standard RS, indicating the standard RS generally leads to over-rejection of the null hypothesis. It also implies the potential problem of $W$ being endogenous when economic distances are effective in constructing $W$, urging the importance of testing the endogeneity of $W$ as a basic work for analysis. \\

Conclusion

Even though conventional uses for the spatial weights matrices ($W$) have been found in the predetermined geography, one may allow $W$ to include economic distances, following the accumulating evidence in economics literature. However, this may lead to the violation of the exogenous assumption for the ordinary spatial autoregressive (SAR) estimators as well as for the related test statistics. For this purpose, I propose Robust Rao's Score (RS) test to determine endogeneity of spatial weights matrices ($W$) in spatial dynamic panel data (SDPD) models.

The robust Rao's Score (RS) test is robust in the sense that it is asymptotically central chi-squared under the null regardless of the presence of local misspecifications in the contemporaneous dependence over space, dependence over time, and spatial time dependence, keeping the type I error fixed. It is also computationally efficient in the sense that it only requires the restricted ML estimators under the null where the parameters above are assumed to be zero, reducing the spatial dynamic panel data models to the simple fixed-effects model.

A Monte Carlo simulation supports the analytics and shows nice finite sample properties as $n$ and $T$ increase. Subsequently, an empirical illustration using Penn World Table version 6.1 shows how large the robust RS test adjusts toward the local misspecifications in parameters compared to the standard RS test. Also, it reaffirms the importance of testing the endogeneity of $W$ as a basic work for analysis due to the potential problem of $W$ being endogenous when economic distances are effective in constructing $W$.

Reference

\scriptsize