Moran's I 2-Stage Lasso: for Models with Spatial Correlation and Endogenous Variables

Moran's $I$ 2-Stage Lasso: for Models with Spatial Correlation and Endogenous Variables.

Introduction

The main aim of structural economic modeling is to explain the evolution of endogenous variables of interest, given fundamental processes such as productivity, taste, and policy. It has long been known that Ordinary Least Squares (OLS) estimation of the coefficients of such endogenous variables is invalidated by endogeneity bias and that instrumental variables (IV) offer a way around this problem wright28. This paper considers the case where the researcher is similarly interested in estimating the parameters on endogenous variables, but where in addition both the structural equation being estimated and the endogenous variables themselves spatial processes based on a given spatial weights matrix (SWM).\footnote{A spatial weights matrix is an $n\times n$ matrix that describes the pair-wise relationship between each of the $n$ cross-sectional units.} Crucially, while SWM is assumed to be known, we do assume that the exact functional forms of the spatial processes are unknown, and possibly include higher-order powers of the SWM. Because the researcher is only interested in estimating the direct effect of the right-hand-side variable(s), the corresponding spatial parameters are thus considered nuisance parameters.

This setup is arguably a realistic situation in applied research: testing for cross-sectional/ spatial dependence is relatively easy, for example using a Moran's $I$ test moran50, but determining the exact form of the spatial process is much more challenging, and might not form the focus of the research. Similarly, spatial dependence and endogeneity are common in many economic models. Some examples include modelling the relationship between economic growth and energy consumption or pollution, employment and migration, and the effect of policing on crime. Many papers in the econometrics literature have shown how to incorporate endogenous variables into a given spatial model.\footnote{Some recent examples include H18_JBES, J16_ET, XL13_ER, FLG08.} The Generalised Method of Moment (GMM) based estimation techniques such as Generalised Spatial Two-Stage Least Squares (GS2SLS) are commonly used by applied researchers when estimating a spatial model with an endogenous variable. However, to use any of the proposed GMM-based estimation techniques, the researcher must specify (1) a spatial economic model and (2) define the spatial structure, i.e., the SWM. A misspecified model will yield inconsistent estimates, and this problem is more acute if the SWM is also misspecified lsK_sens14.

Given that the spatial process is assumed to be of a lesser interest to the researcher than the direct economic impact of the endogenous variable, i.e. the spatial parameters are considered nuisance parameters, we propose relying on the Eigenvector Spatial Filtering (ESF) approach developed by grif2000, grif03. This has the key advantage of being agnostic to the underlying functional form of the spatial process. Instead of explicitly modelling the underlying spatial process, ESF uses a subset of eigenvectors from the SWM as controls in a linear regression framework to control of the spatial dependence, removing the need to specify and estimate a spatial process.

Leaving aside the issue of endogeneity for a moment, the main downside of ESF is that estimation using the full set of eigenvectors is infeasible using OLS. Given $k$ covariates, the addition of the $n$ eigenvectors produced by the spectral decomposition of the SWM necessarily produces a rank-deficient Gram matrix with $n+k$ parameters and $n$ observations. This problem can be mitigated by making a sparsity assumption, i.e. assuming that only a subset of the eigenvectors are relevant and will have non-zero coefficients. This generates a separate problem, however, which is the selection of the relevant subset of eigenvectors. To solve this selection problem, we propose using a Lasso-based procedure that uses information contained in the Moran's $I$ statistic to determine a point estimate for the Lasso tuning parameters. The proposed estimator, called Moran's $I$ two-stage Lasso (Mi-2SL), is a three-step procedure: the first and second stages of a general two-stage least squares (2SLS) specification are separately estimated by using this Moran's $I$ based Lasso, in order to extract the relevant eigenvectors. The union of the two sets of selected eigenvectors is then used to provide supplementary covariates in a standard 2SLS regression. This 2SLS specification deals with the endogenous variables, with the additional eigenvectors selected via Moran's $I$ based Lasso dealing with the (weak) cross-sectional dependence.\footnote{We will use the terms cross-sections dependence and spatial dependence interchangeably.}

Several studies have already used two-stage Lasso procedures in a spatial setting. For example, peng19 estimates a spatial autoregressive model (SAR) by a two-stage Lasso procedure to allow heterogeneous peer effects and the identification of the influential individuals in a network. As both stages are high-dimensional, they are both estimated by Lasso. Ahren15 estimate the effect of conflict risk on economic growth using bcch12 two-stage procedure, where Lasso estimates the high-dimensional first stage and the second is a low-dimensional panel SAR model. Additionally, arhbha15, LS16_SWM, LS20_SWM all use two-stage Lasso-based procedures to estimate/select a SWM. We are the first, however, to consider a two-stage Lasso procedure for ESF.

The specific contribution we bring is to derive theoretical results on consistent and asymptotically normal parameter estimation. Proving consistency and asymptotically normality is tricky: as the eigenvectors are derived from the SWM, which itself encodes the pair-wise dependence between the observations, one cannot rely on the standard assumption of row-wise independence. To get around this problem we rely instead on the LT_KMS20 notion of $\psi$-dependence and corresponding limit theorems to derive our results. These theoretical results are supported by a set of Monte Carlo simulations, where the estimator is tested against competing methodologies for varying degrees of correlation between the first and second-stage errors as well as varying levels of spatial dependency in the covariates. The analysis shows that Mi-2SL performs well relative to competitors in small samples, and out-performs them in terms of bias and mean squared errors in the presence of spatially correlated covariates.

Finally, as a motivating application, we apply our methodology to ck16, who analyse the impact of Mexican worker mobility on local labour market outcomes of natives in the US, using a standard IV strategy to correct for endogeneity. Despite having an explicit spatial dimension in their data, their analysis does not allow for spatial dependence in their specification. A standardised Moran's $I$ test on the first and second-stage residuals indicates significant spatial correlation for most demographic groups, with a higher spatial correlation level in the first stage than the second. This forms an idea use-case for Mi-2SL, as the functional form of the spatial process is uncertain, and it is not the main focus of the research question. We re-estimate their model using Mi-2SL to account for the unknown spatial structure and find that while Mi-2SL does not change the overall conclusion of ck16, it substantially improves the strength of the Bartik instrument in the first stage, thus improving the precision of the second stage estimates.

The rest of the paper is structured as follows. Section (ref) presents the underlying structural model, the notation, and the proposed the Mi-2SL procedure. In section (ref) we drive the theoretical properties of Mi-2SLS under perfect selection. Section (ref) provides Monte Carlo studies to evaluate the finite sample properties of the proposed estimator and in Section (ref), we apply the proposed procedure to ck16. Finally, Section (ref) offers our concluding remarks.

Structural model and estimation procedure

Underlying structural model

Consider the following structural equation where the endogenous $n\times 1$ vector $\bm{y}$ which depends on an $n\times k_1$ matrix of exogenous regressors $\bm{X}_1$, an $n\times 1$ endogenous vector $\bm{x}_2$ and follows some spatial process:

where $f(\bm{W},\bm{y},\bm{X}_1)$ is a linear combination of spatial lags of $\bm{y}$ and $\bm{X}_1$ obtained with $\bm{W}$, a $n \times n$ symmetric weights matrix, and $\bm{\varepsilon}$ is an $n\times 1$ vector of innovations. $f(\bm{W},\bm{y},\bm{X}_1)$ is allowed to contain higher-order spatial lags $\bm{W}^i\bm{y}$ and $\bm{W}^i\bm{X}_1$ with $i>1$. An example of a common special case of this process is:

where $\rho_{i,0}$'s and $\bm{\psi}_0$ are unknown parameters that represent the degree of spatial correlation in the endogenous variable $\bm{y}$ and the predetermined exogenous variables $\bm{X}_1$ with moment conditions $\operatorname{\mathbb{E}}[\bm{X}_1'\bm{\varepsilon}]=0$ and $\operatorname{\mathbb{E}}[(\bm{W}\bm{X}_1,\bm{X}_1)'\bm{\varepsilon}]=0$. The exact spatial process is unknown, in the sense that some of these spatial parameters, including $p$, are allowed to be zero-valued.\footnote{The data generating process of $\bm{y}$ could also include spatial autoregressive disturbances; however this is excluded from the model for simplicity.}

The regressor $\bm{x}_2$ in (ref) is endogenous, in the sense that $\operatorname{\mathbb{E}}(\bm{x}_{2}'\bm{\varepsilon} )\neq 0$, and $\beta_{2,0}$ is the parameter of interest to the researcher. The extension to the case where $\bm{x}_2$ is a matrix is straightforward and omitted for simplicity. We assume that $\bm{x}_2$ also follows some unknown spatial process:

where $\bm{Z}_2$ is a $n\times q$ matrix of instrument variables with $q\geq 1$ and moment conditions $\operatorname{\mathbb{E}}(\bm{Z}_2'\bm{\varepsilon})= 0$. $\bm{u}_2$ is a vector of disturbances with $\operatorname{\mathbb{E}}[(\bm{X}_1,\bm{Z}_2,\bm{W}\bm{X}_1,\bm{W}\bm{Z}_2)'\bm{u}_2]=0$ and $\operatorname{\mathbb{E}}[(\bm{X}_1,\bm{Z}_2,\bm{W}\bm{X}_1,\bm{W}\bm{Z}_2)'\bm{\varepsilon}]=0$. Again, the spatial process $g(\bm{W},\bm{x}_2, \bm{X}_1,\bm{Z}_2)$ is some linear combination of spatial lags of $\bm{x}_2$, $\bm{X}_1$ and $\bm{Z}_2$, obtained with $\bm{W}$. An example of such a process is:

Let $N_n=N=\{1, \ldots, n\}$ be the set of cross-sectional unit indices with $n\in\operatorname{\mathbb{N}}$ denoting the number of observations. For reasons of generality, we allow the elements of $\bm{\varepsilon}=\bm{\varepsilon}_n$, $\bm{y}=\bm{y}_n$, $\bm{W}=\bm{W}_{n}$, $\bm{Z}_2=\bm{Z}_{2,n}$, $\bm{u}_2=\bm{u}_{2,n}$, $\bm{X}_1=\bm{X}_{1,n}$ and $\bm{x}_2=\bm{x}_{2,n}$ to be dependent on $n$, that is to form triangular arrays. However, to simplify the notation, the $n$ index is omitted.

Equation (ref) contains two sources of endogeneity, first $\bm{x}_2$ because $\operatorname{\mathbb{E}}(\bm{u}_2'\bm{\varepsilon} )\neq 0$, which implies $\operatorname{\mathbb{E}}(\bm{x}_2'\bm{\varepsilon} )\neq 0$. Second, $\bm{y}$ itself is endogenous as it appears on both sides of (ref), via $\bm{W}^i\bm{y} \; \forall i$. Both sources of endogeneity cause the OLS estimate of $\bm{\beta}_0=(\bm{\beta}_{1,0}, \beta_{2,0})'$ to be inconsistent ($\hat{\bm{\beta}}_{ols}\not\to_p \bm{\beta}_0$).

Substituting (ref) into (ref) gives the reduced forms for $\bm{y}$:

where $\bm{\pi}_{1,0}=\beta_{2,0}\bm{\zeta}_{1,0}$, $\bm{\pi}_{2,0}=\beta_{2,0}\bm{\zeta}_{2,0}$, $\bm{\pi}_{3,0}=\beta_{2,0}\bm{\zeta}_{3,0}$, $\bm{\pi}_{4,0}=\beta_{2,0}\bm{\zeta}_{4,0}$, $\bm{d}=\bm{S}_2^{-1}\bm{u}_2\beta_{2,0} + \bm{\varepsilon}$ and both $\bm{S}_1\equiv(\bm{I}-\sum^p_{i=1}\rho_{i,0}\bm{W}^i)$, and $\bm{S}_2\equiv(\bm{I}-\sum^l_{i=1}\bm{W}^i\zeta_{i,3,0})$ are non-singular.

Moran's $I$ 2-Stage Lasso

The existence of valid instruments $\bm{Z}_2$ for the endogenous $\bm{x}_2$ implies that we can deal with the problem of endogeneity, leaving the key challenge of controlling for the unknown underlying spatial processes in (ref) and (ref). Even if the exact underlying spatial process were known, estimation of (ref) would be feasible, albeit non-trivial. One method would be to first estimate (ref) by GS2SLS, first developed by kelpru98 and extended by DEP19_HOG2SGS to allow for higher-order spatial lags, which would use higher order spatial lags of the exogenous variables in (ref) as instruments for $\bm{W}^i\bm{x}_2 \; \forall i$. The resulting fitted values can then be used to estimate (ref). GS2SLS has the advantage that it can be easily extended to include other right-hand-side endogenous variables. However, the procedure requires that the researcher specify which spatial parameters to estimate, and given this extra layer of estimation, the standard GS2SLS standard errors would be invalid.

Given the additional assumed uncertainty regarding the true functional form of the spatial process in the model, we propose using eigenvectors $\bm{E}_{n}=\bm{E}$ from a spectral decomposition of $\bm{W}$ to represent $f(\bm{W},\bm{y},\bm{X}_1)$ and $g(\bm{W},\bm{x}_2, \bm{X}_1,\bm{Z}_2)$ i.e. $f(\bm{W},\bm{y},\bm{X}_1)=\bm{E}\bm{\gamma}_{y,0}$ and $g(\bm{W},\bm{x}_2, \bm{X}_1,\bm{Z}_2)=\bm{E}\bm{\gamma}_{x,0}$ where $\bm{\gamma}_{y,0}$ and $\bm{\gamma}_{x,0}$ are vectors of unknown constants. This methodology has the key advantage that it is agnostic to the exact form of $f(\bm{W},\bm{y},\bm{X}_1)$ and $g(\bm{W},\bm{x}_2, \bm{X}_1,\bm{Z}_2)$, including the presence of higher-order lags, stemming from the spectral property that the eigenvectors from $\bm{W}$ and $\bm{W}^i$ $\forall i \in \mathbb{Z}^+$ are the same. Using this linear representation, one could in principle estimate the following system instead of (ref) and (ref):

where $\bm{G}=[\bm{X}_1,\bm{x}_2,\bm{E}]$, $\bm{\Upsilon}_0=[\bm{\beta}_{1,0},\beta_{2,0}, \bm{\gamma}_{y,0}]'$, $\bm{Z}=[\bm{X}_1,\bm{Z}_2, \bm{E}]$ and $\bm{\zeta}_0=[\bm{\zeta}_{1,0}, \bm{\zeta}_{2,0}, \bm{\gamma}_{x,0}]'$ with $\operatorname{\mathbb{E}}[\bm{G}'\bm{\varepsilon}]=0$ and $\operatorname{\mathbb{E}}[\bm{Z}'\bm{u}_2]=0$.

The practical obstacle is that (ref) and (ref) are both high-dimensional linear regressions, as in each equation the number of parameters is greater than the number of observations. This means both the (re-scaled) Gram matrices $\bm{G}'\bm{G}/n$ and $\bm{Z}'\bm{Z}/n$ are necessarily rank deficient. Thus, neither (ref) nor (ref) cannot be estimated by OLS nor (ref) by 2SLS. grif2000 argues, however, that in most cases only a subset of eigenvectors are relevant to the data generating process (DGP) of $\bm{y}$ and $\bm{x}_2$, i.e. the parameter vectors $\bm{\gamma}_{y,0}$ and $\bm{\gamma}_{x,0}$ are sparse. The intuition behind this sparsity assumption is each of the $n$ eigenvectors can be viewed as an orthogonal spatial pattern, and only a specific subset of these patterns are relevant to the DGP of $\bm{y}$ and $\bm{x}_2$ grif03. Thus, the estimation problem turns into a selection problem.

We propose addressing this selection problem with an extension of the Moran's $I$ based Lasso first proposed in c23. This procedure considers a single structural equation where all the covariates are exogenous, i.e., (ref) with $\beta_2=0$. It only penalises the $\bm{\gamma}_y$ coefficients on the eigenvectors $\bm{E}$ and set the Lasso tuning parameter to $z^{-2} \; \forall \; z\neq 0$ where $z$ is the standardised Moran's $I$ ($z$) of the residual $\hat{\bm{h}}=\bm{M_Xy}$, with $\bm{M_X}=\bm{I}-\bm{X}_1(\bm{X}_1'\bm{X}_1)^{-1}\bm{X}_1'$.

with

Given that the aim of ESF is to directly control for spatial correlation patterns in the regression, the intuition behind calibrating the tuning parameter this way is that when the level of spatial correlation in the residuals is low, only a small set of eigenvectors is necessary, thus a high level of regularization (large tuning parameter) is required. In contrast, when the level of spatial correlation is high, a larger set of eigenvectors will be necessary, thus a low level of regularization (small tuning parameter) is required. As $z$ gives a large value when the overall correlation is high and small values when the overall correlation is low, they propose using the inverse square of the standardised Moran's $I$ as the tuning parameter.\footnote{A positive tuning parameter is required for the Lasso solution to be unique. Thus, the squared value of $z$ is used.}

Our proposed Moran's $I$ 2-stage Lasso (Mi-2SL) procedure, outlined in Algorithm (ref), can handle both endogenous covariates and cross-sectional dependence. The procedure is straightforward: first a spectral decomposition of the SWM is performed to get the candidate set of eigenvectors. The standardised Moran's $I$ on the naïve first stage residuals (ignoring the spatial correlation) provides the tuning parameter for a Lasso (or post-Lasso) estimation of (ref) to get $\hat{\bm{x}}_2$ as well as the selected eigenvectors $\hat{\bm{E}}_x$. Subsequently, $\hat{\bm{x}}_2$ is used instead of $\bm{x}_2$ to calculate standardised Moran's $I$ for the naïve second stage residuals (ignoring the spatial correlation), which serves as the tuning parameter for a lasso estimation of (ref), providing a second set of selected eigenvectors $\hat{\bm{E}}_y$. As a final step, $\beta_2$ is estimated by standard 2SLS using the union of $\hat{\bm{E}}_x$ and $\hat{\bm{E}}_y$ as additional controls.

Theoretical results

Assumptions

We will now derive some theoretical properties of the proposed Mi-2SL procedure. This requires two sets of assumptions, the first of which applies to the underlying data generating processes (ref) and (ref).

Assumption (ref).1 is standard in the spatial econometrics literature kelpru98,kp99, lee04. Note, assumption (ref).1 (a) is required for the spectral decomposition and Assumption (ref).1 (d) ensures that the elements of the eigenvectors have the same dependence coefficient as the elements of the SWM. Assumption (ref).2 and (ref).3 are standard assumptions in the instrument variables literature.

The second set of assumptions relates to the ESF approximation itself. Assumption (ref).1 says there exists a set of linearly dependent eigenvectors and corresponding unknown constants that will approximate the functions $f(\bm{W,y,X}_1)$ and $g(\bm{W},\bm{x}_2, \bm{X}_1,\bm{Z}_2)$. Assumption (ref).2 assumes this approximation is weakly sparse and Assumption (ref).3 assumes the approximation errors go to zero at a sufficient speed.

Under these assumptions, the high-dimensionality ESF system ((ref)) and ((ref)) can be expressed as the following low dimensional reduced form system of equations:\footnote{While these assumptions cannot be verified in practice or even in simulations, they are common feature in the ESF literature, as well as related methodology such as factor or principal component analysis}

where $\bm{G}_{\Omega}=[\bm{X}_1,\bm{x}_2,\bm{E}_{\Omega}]$, $\bm{Z}_{\Omega}=[\bm{X}_1,\bm{Z}_2, \bm{E}_{\Omega}]$ , $\bm{\Upsilon}_{\Omega}=[\bm{\beta}_{1,0}', \beta_{2,0}, \bm{\gamma}_{\Omega}']'$, $\bm{\zeta}_{\Omega}=[\bm{\zeta}_{1,0}', \bm{\zeta}_{2,0}', \bm{\zeta}_{3,\Omega}']'$, $\bm{\varepsilon}_\pi=\bm{\varepsilon}+\pi_y$ and $\bm{u}_\pi=\bm{u}_2+\pi_x$ .

Even assuming that the subset of relevant vectors $\Omega$ is known, establishing that (ref)-(ref) can be estimated by 2SLS is non-trivial, for two reasons. First, we have the two additional approximation errors $\pi_y$ and $\pi_x$ in the first and second stage errors and second, the standard weak law of large numbers (LLN) and central limit theorem for triangular arrays used for spatial models requires assuming the row-wise independence. This is not realistic here as $\bm{G}_{\Omega}$ and $\bm{Z}_{\Omega}$ contain elements of $\bm{E}$, constructed from a linear transformation of $\bm{W}$, a matrix which itself encapsulates the spatial dependence of observations. Establishing the theoretical properties of the procedure therefore requires formalising this dependence and applying the appropriate limit theorems.

To do so, we use the notion of $\psi$-dependence first proposed by dl99 for time-series data and adapted by LT_KMS20 to allow for cross-sectional dependence. This allows us to use the limit theorems proposed by LT_KMS20. Roughly speaking, $\psi$-dependence measures the strength of dependence between two sets of random variables by the covariance of non-linear functions of the random variables.

Let $\{w_{ij},1\leq i\leq n,n\geq 1\}$, $j=1,\ldots,n$, be a triangular array of random variables, where $w_{ij}=w_{ij,n}$ denotes the $i,j$th element of matrix $\bm{W}$ which is derived from a spatial structure as follows. For any $a\in \operatorname{\mathbb{N}}$, we endow $\operatorname{\mathbb{R}}^{a}$ with distance: \[ \mathtt{d}_a(\bm{q},\bm{h}) = \sum_{l=1}^a|q_l-h_l| \] where $\bm{q}=(q_1,\ldots, q_a)$ and $\bm{h}=(h_1,\ldots,h_a)$ are points in $\operatorname{\mathbb{R}}^{a}$. Let $\pazocal{L}_{a}$ denote the family of real valued, bounded Lipschitz functions, with $\operatorname*{Lip}(f)$ the Lipschitz constant of $f$ and $||f||_\infty=\sup_x|f(x)|$ its sup-norm. \[ \pazocal{L}_{a}=\{f:\operatorname{\mathbb{R}}^{a}\to \operatorname{\mathbb{R}}:||f||_\infty <\infty ;\; \operatorname*{Lip}(f)<\infty\} \]

Now consider two sets of cross-sectional units (of size $a$ and $b$ $\in\operatorname{\mathbb{N}}$) with a distance between each other of at least $r>0$. Let $\pazocal{P}_{a,b;r}$ denote the collections of all pairs \[ \pazocal{P}_{a,b;r} = \{(A,B):A,B\subset N, |A|=a,\; |B|=b,\; d_{A,B}\geq r\} \] where $d_{A,B}=\min_{i\in A}\min_{j\in B}d_{ij}$.\footnote{Note that $\pazocal{P}_{a,b;r}$, $d_{A,B}$ and $d_{ij}$ are also implicitly indexed by $n$, but we again omit the index to simplify the notation} $\forall$ sets $A$ of positive integers, define $w_A=\{w_{ij}:i\in A\}$.

We take $\{\operatorname{\pazocal{C}}_n=\operatorname{\pazocal{C}}\}$ be a sequence of given $\sigma$-fields, such that for each $n\geq 1$, the spatial weights matrix $\bm{W}_n=\bm{W}$ is $\operatorname{\pazocal{C}}$-measurable. Definition (ref) gives the exact definition of conditional $\psi$ dependence we use.

We will now explicitly specify the latent spatial formation process. We consider binary connectivity based on physical distance plus some stochastic elements. Specifically, the connection for each pair of spatial units $i$ and $j$ ($i\neq j$) is randomly realised if and only if:

where the $\phi_{ij}$'s and $\eta_{ij}$'s are random variables such that $\phi_{ij,n}=\phi_{ij}=\phi_{ji}$, $\eta_{ij}=\eta_{ji}$ and $\{\eta_{ij} \; : \; i < j \}$ is $i.i.d.$ and independent of $\phi = (\phi_{ij})_{i<j}$. The random variable $\phi_{ij}$ which determines the formation probabilities is assumed to be a function of observable characteristics $\bm{l}_{ij,n}=\bm{l}_{ij}$(e.g., the physical distance between the spatial units) and unit specific unobservable characteristics $\bm{t}_{i,n}=\bm{t}_i$ (i.e. $\phi_{ij}=f(\bm{t}_i, \bm{t}_j, \bm{l}_{ij})$ where $f(\cdot )$ is some function). Thus, the $\sigma$-field $\operatorname{\pazocal{C}}$ is generated by $\bm{t}_i$, $\bm{t}_j$ and $\bm{l}_{ij}$ for all $i$ and $j$.

Let us introduce following additional notations. Let $\tilde{\bm{E}}$ be either equal to $\bm{M}_H \bm{E}$ or $\bm{M}_{\hat{\bm{X}}} \bm{E}$, and $\bm{C}_{\Omega k\Omega k}=n^{-1}\tilde{\bm{E}}_{\Omega k}'\tilde{\bm{E}}_{\Omega k}$, $\bm{C}_{\Omega k\grave{\Omega k}}=n^{-1}\tilde{\bm{E}}_{\Omega k}'\tilde{\bm{E}}_{\grave{\Omega k}}$, $\bm{C}_{\grave{\Omega k}\Omega k}=n^{-1}\tilde{\bm{E}}_{\grave{\Omega k}}'\tilde{\bm{E}}_{\Omega k}$ and $\bm{C}_{\grave{\Omega k}\grave{\Omega k}}=n^{-1}\tilde{\bm{E}}_{\grave{\Omega k}}'\tilde{\bm{E}}_{\grave{\Omega k}}$ where $\tilde{\bm{E}}_{\Omega k}$ is an $n\times s_k$ matrix with columns corresponding to the active set $\Omega k$. $\grave{\Omega k}$ is the complement set and the $n\times q_k$ matrix $\tilde{\bm{E}}_{\grave{\Omega k}}$ is defined accordingly with $q_{nk}=q_k=s_k-n$. Now the (re-scaled) Gram matrix $\bm{C}_n=\bm{C}=n^{-1}\tilde{\bm{E}}'\tilde{\bm{E}}$ can be expressed in block-wise form as:

Similarly we define $\bm{\gamma}=[\bm{\gamma}_{\Omega k},\bm{\gamma}_{\grave{\Omega k}}]'=[\gamma_1,\ldots,\gamma_{s_k},\gamma_{s_k+1},\ldots,\gamma_{n}]'$, with $k=x$ or $y$

Assumption (ref) is similar to Assumption 4 in c23. These are conditions on the eigenvectors and eigenvalues to assure consistent selection.

Consistent Eigenvector Selection

We now derive conditions under which Algorithm (ref) selects the relevant eigenvectors in steps 3 and 4. c23 discusses the conditions for consistent selection in a SAR model, which involve some restrictions on eigenvalues and the level of sparsity $s_1$ and $s_2.$

Proof: The proof of the Lemma (ref) follows immediately form the application of Theorem 2 from c23.

Estimation consistency

We will now derive a consistency proof for estimating $\bm{\Upsilon}_{\Omega}$ by 2SLS, assuming $\Omega$ is known. In scalar notation (ref) can be rewritten as:

for $i=1,\ldots, n$.

We now state the additional assumptions for consistent estimation of $\beta_{2,0}$ by 2SLS.

Assumption (ref).1 is from LT_KMS20 and the function $\psi_{a,b}$ satisfies Assumption (ref).1 if: \[ \sup_{n\geq 1}\max_{i\in N}\operatorname{\mathbb{E}}[|w_{ij}|^q|\operatorname{\pazocal{C}}_n] \enspace < \enspace \infty \quad \operatorname*{a.s.} \] for some $q>4$ and $\forall j$. Assumption (ref).2 states that all variables have conditional finite second moments, so all are $\operatorname{\pazocal{C}}$ measurable. Assumption (ref).3 is also from LT_KMS20 and puts a restriction on the denseness of the spatial structure and the rate of decay of dependence with regards to the distance between the spatial units. In the mixing literature, it is common to assume the mixing coefficients can be summed $n^{-1}\sum_{r=1}^\infty \mu_{r}=O_p(1)$ as $n\to\infty$. A sufficient condition for Assumption (ref).3, in this case, is if the average number of neighbours at distance $r$ grows slower than the sample size $n$, i.e. $\sup_{r\geq 1} \delta^d_r = o_p(n)$. Intuitively, this assumption requires that the number of spatial connections at distance $r$ not grow too fast as $r$ increases. However, as the precise condition (ref) includes the dependence coefficient $\mu_{r}$, this assumption can be relaxed if $\mu_{r}$ itself decreases at an appropriate rate relative to $r$. This assumption seems reasonable as the literature on estimating SWMs often assumes a sparse spatial structure arhbha15,LS16_SWM,LS20_SWM. An example of where Assumption (ref).3 could fail is if one unit is connected to all other units, such as in the star network. This is because the distance between any two units is never larger than 2.\footnote{$\delta^d_1 = 2(n-1)/n$, $\delta^d_2 = (n-2)(n-1)/n$ and $\delta^d_r = 0$ for $r\geq 3$} Assumption (ref).4 requires the instruments (including $\bm{E}_\Omega$) be uncorrelated with the structural error, conditional on $\operatorname{\pazocal{C}}$.

Lemma (ref) below establishes that assumption (ref).1, which requires that $w_{ij}$ are $\psi$-dependent triangular arrays, carries over to the eigenvector elements $e_{ik}$. Given the eigendecomposition $\bm{W} = \bm{E} \bm{\Lambda} \bm{E}^T$, these are generated by a linear combination of $w_{ij}$, $\lambda_k$ and $e_{jk}$ as follows:

for all $\lambda_k \neq 0$

Proof: This is provided in appendix (ref).

Lemma (ref) shows that as long as the largest eigenvalue is bounded and the eigenvectors are mutually orthogonal (both of these requirements are satisfied by Assumption (ref).1) the eigenvector elements will have the same dependence coefficients $\{\mu_r\}$ as the elements of the SWM.

Proof: This is provided in appendix (ref).

Theorem (ref) shows that under an appropriate mixing condition, some additional regularity conditions and if $\Omega$ is known, we could estimate $\bm{\Upsilon}_{\Omega}$ consistently by 2SLS. The proof of Theorem (ref) uses the weak LLN for triangular arrays, which gives convergence in probability, and the strong LLN for cross-sectionally dependent random variables of LT_KMS20 which gives almost sure convergence, thus, overall gives convergence in probability. An almost sure convergence result could be obtained similarly by using the strong LLN for triangular arrays instead of the weak LLN for triangular arrays.

Asymptotic Distribution

In order to derive the asymptotic distribution of the 2SLS estimator for a known support $\Omega$ of the relevant eigenvector set, we need some additional assumptions:

Assumption (ref).1 states that all variables have at least conditional fourth finite moment, so are all $\operatorname{\pazocal{C}}$ measurable, which is in line with many spatial and 2SLS models. Assumption (ref).2 is from LT_KMS20 and limits the extent of the spatial dependence of the random variables through restrictions on the spatial structure. When the spatial structure is given $c_{r,m;k}$ can be computed, it is composed of two parts $\Delta_{r,m;k\alpha}$ and $\delta^d_{r,\alpha/(1-\alpha)}$, which capture the denseness of the spatial structure through the average size of neighbourhoods and the average shell size of the neighbourhood. Note that after $r$ goes beyond a certain level $\Delta_{r,m;k}$ tends to decrease fast, as the set $N_{j,r-1}$ becomes large quickly. For (ref) to be satisfied $\mu_r$ (the spatial dependence) needs to decay fast enough as $r$ becomes large, this is because it will become increasingly difficult to find a slowly increasing sequence $m$ to satisfy the condition.

Proof: This is provided in appendix (ref).

Theorem (ref) shows that if $\Omega$ is known, then under an appropriate mixing condition, restriction on the denseness of the spatial structure, and some additional regularity conditions, the 2SLS estimate of $\bm{\Upsilon}_{\Omega}$ and thus, $\beta_{2,0}$ will be asymptotically normal, with a convergence rate of $n^{-1/2}$.

Simulation

In this section, we provide simulation evidence to assess the finite sample performance of the Mi-2SL estimator and compare its performance to some commonly used estimator for spatial models. We generate the following system of equations (ref) - (ref) where the structural equation includes a SAR(1) with spatial lags of the exogenous/endogenous variables, and the endogenous variable follows a SAR(2) with spatial lags of the exogenous variable/instrument:

with $\bm{z}_2\sim N(0,\bm{I})$, $\bm{x}_1\sim N(0,\bm{I})$ and $u_i$, $v_i$ (the $i$th elements of $\bm{u}$ and $\bm{v}$) are given by: \[ (u_i,v_i)\sim N\left(0,

\right) \]

We set the non-spatial parameters to $\zeta_1=\zeta_2=\beta_1=\beta_2=1$ and $\sigma^2_{v,u}=0.9$, and the spatial parameters are combinations of the following values: $\rho\in\{0,0.4,0.8\}$, $\zeta_{3,1}\in\{0.4,0.8\}$, $\zeta_{3,1}\in\{0,0.4\}$ and $\omega\in\{0,0.4,0.8\}$.

The SWM $\bm{W}$ is generated using a watts1998collective small world network model. Small world networks are a popular way of modelling cross-sectional dependency in social networks, and have been used in many economic applications, particularly the economics of innovation diffusion and industrial clusters jackson2005economics, cassi2008opportunity, maggioni2011networks, ter2011co, gulati2012rise, bagley2019small. The number of neighbours is set to 10 and the rewiring probabilities to $p\in\{0.4,0.8\}$. This allows to see the difference been a higher level of clustering ($p=0.4$) and a lower level of clustering ($p=0.8$). Each SWM is normalised by the largest row sum and the eigenvectors are from the normalised SWM. Sample sizes considered are $n\in\{100,250,500\}$, and we run 1000 Monte Carlo replications for each experiment.

The estimators and specifications compared are:

{5pt}

{6pt}

Note, 2SLS-SAR is based on the procedure proposed by kelpru98 (also commonly refereed to as Generalised Spatial two Stage Least Squares). This is included in the comparison set as it is a common spatial model used by applied researchers, and provides a more challenging benchmark for Mi-2SL, as it is a specification where a genuine attempt is made at controlling for both the endogeneity of $\bm{x}_2$ and the presence of a spatial process.

Tables (ref), (ref), and (ref) presents the bias, mean squared error (MSE) and average asymptotic standard error (AASE) of $\beta_2$ for $\omega=0.4$ and sample sizes of 100, 250 and 500 respectively.\footnote{Tables (ref), and (ref) are provided in appendix (ref). Additional extended results for $\omega=\{0,0.8\}$ and the $\rho=0$ case can be found in the supplementary material} These tables exhibit the standard bias-variance trade off between naïve OLS and naïve IV. OLS has the smallest AASE but the largest bias, whereas IV eliminates a substantial part of the bias but has the largest AASE. As both ignore the presence of a spatial process, the 2SLS-SAR is able to decrease both the bias and AASE compared to IV.

The Mi-2SL estimators have the second smallest AASE (smaller than simpIV and 2SLS-SAR) overall, and when $\rho$ is large Mi-2SL has the same or slightly larger AASE as OLS while having a small bias similar to the SAR(1). Generally when the rewiring probability is small ($p=0.4$) Mi-2SL has a larger absolute bias compared to when the rewiring probability is large, regardless of the sample size. In terms of eigenvector selection behaviour, the number of selected eigenvectors increases with $\rho$ and as the sample size increases. When the sample size is small ($n=100$) more eigenvectors are selected when $p=0.8$, whereas for larger sample sizes more eigenvectors are selected when $p=0.4$. More eigenvector are selected when the first stage fitted values come from Lasso estimate (Mi-2SLl), this is because more eigenvectors are selected in the second stage.

In summary Mi-2SL performs well compared to OLS, naïve IV and the SAR(1) estimated by 2SLS. It has a smaller AASE than Classical IV and the SAR(1). In particular, when the level of spatial correlation of the dependent variable in the structural equation is high, Mi-2SL AASE is similar to that of OLS. In terms of bias Mi-2SL generally performs better than both OLS and IV, and similarly to the SAR(1).

Application on impact of migration on labour markets

This section revisits the empirical application of ck16. Using an IV strategy to control for the endogeneity of labour market decisions, their main finding is `that low-skilled Mexican-born immigrants’ location choices respond strongly to changes in local labour demand, which helps equalize spatial differences in employment outcomes for low-skilled native workers. ck16 starts from the observation that over the Great Recession low-educated Mexican-born male immigrants were more mobile than their native counterparts. Given this observation, their aim was to test if location choice of migrants was being driven by local labour market conditions, leveraging the geographic variation in employment changes during the Great Recession as a natural experiment. Their argument rests on fact that changes in labour market conditions during the Great Recession can be approximately measured by changes in employment, as traditionally sticky-downwards wages were essentially fixed during that period. Thus, they look at the effect of changes in employment on population changes for 20 different demographic groups, split by gender (males and females), education (`high school or less' and `some college or more'), and location of birth (native-born, foreign-born, Mexican-born, and other foreign-born). The unit of observation is a metropolitan statistical area (MSA), 95 of which are included in their IV analysis. Their empirical specification is:

where $i$ indexes the MSA, $\Delta pop_i$ is the proportional change in working-age population from 2006–2010, $\Delta emp_i$ is the proportional change in employment from 2006–2010, $mex_i$ is the share of Mexicans-born population in 2000, $policy_i$ and $287g_i$ are both immigration policy controls, and $\Delta bartik_i$ is the `Bartik instrument' b91, which predicts changes in local labour demand by assuming that in each industry national employment changes are proportionately allocated across cities, based on each cities initial industry composition of employment. For the reader's convenience, Table (ref) replicates their main IV results, Table 4 in ck16. We have also added the first stage (full) F-statistic, so this can be compared to the partial F-statistic and give further insight into the actual impact of the Bartik instrument in their estimates. Table (ref) shows the full F-statistic is always smaller than the partial F-statistic, implying that in their specification, the Bartik, which is supposed to give the identification, is not helping the first stage.

A potential issue with the estimation of (ref) is the potential existence of spatial dependency between the MSAs used in the analysis. Figure (ref) shows the 95 MSAs included in the ck16 IV analysis, revealing clear spatial heterogeneity. In order to account for potential spatial correlation we construct a SWM using a binary distance-based cut-off, where $w_{ij}=1$ if the distance between the metropolitan areas is less than $A$ kilometres and zero otherwise. We consider cut-off distances ($A$) of 500km, 600km, and 700km. Note 500km is the smallest distance that ensures every metropolitan area has at least one neighbour.

Table (ref) shows the standardised Moran's $i$ for the first and second stages of ck16 IV regressions obtained for the three SWMs considered. This exercise shows that the standardised Moran's $i$ of the first stage is always significant at the one percent level, and the second stage is also significant at the ten percent level in most configurations. In almost all cases, the first-stage has a substantially higher level of spatial correlation than the second-stage. For low-educated Mexican-born male migrants, the standardised Moran's $i$ is significant at the five percent level in both stages for all three SWMs, with a test statistic three times larger in the first than second stage. Given the presence of spatial dependence in the data there is a legitimate question as to how the IV estimates in (ref) might be affected. As explained in the introduction, this setting provides a realistic use-case for Mi-SL: the research question focuses on the impact of an endogenous covariate (the movements of Mexican-born lower skill workers) in a context where the presence of spatial dependence between the observational units (MSAs) potentially invalidates IV estimation. The spatial process present in the data can be detected in a straightforward manner, however accurately specifying it would go beyond the scope of the research question, and thus the researcher might well prefer to simply control for it, as a nuisance parameter.

Table (ref) shows the estimates obtained using Mi-2SLl (the first stage fitted values from the Lasso estimates, step 3 in Algorithm (ref)) with the 500km cut-off SWM. No eigenvectors are selected in any of the second stages, due to the lower levels of spatial correlation in each of the second stages, so the fitted values from Lasso and post-Lasso yield the same results. Tables of results obtained with the larger the larger cut-off SWMs, which were run as a robustness check, are provided in appendix (ref). The Mi-2SL results do not change the qualitative conclusion of ck16 that low-educated Mexican-born migrants respond positively to changes in employment. However, for low-educated Mexican-born males we find the magnitude of the coefficient increases by approximately a standard error. More generally, a key general impact of the included eigenvectors is a substantial improvement in the first stage F-statistic and partial F-statistic. For example, for low-educated Mexican-born males, the first-stage F-statistic and partial F-statistic increase from 8.28 and 11.94 to 101.39.35 and 55.84 respectively. This improvement in the first-stage estimates leads to an increase in the precision of the predicted values and this of the second-stage estimates, which can be seen by the reduction in the estimated standard errors on employment change from 0.468 to 0.359. The fact that the partial F-statistic is now smaller that the full F-statistic also implies that the Bartik is now having a stronger positive effect, at least in the case of low-educated Mexican-born migrants.

Conclusion

In conclusion, we have proposed a new two-stage lasso-based procedure, called Moran's $I$ 2 stage Lasso (Mi-2SL), to estimate classical regression parameters of endogenous variables in the presence of spatial correlation of an unknown functional form. Under the assumption that the relevant set of eigenvectors is known, that an appropriate mixing condition holds, that some restriction exists on the spatial structure, and some assuming some additional regularity conditions, we show that the Mi-2SLparameter estimates are consistent and asymptotic normal.

Our simulations results establish that the Mi-2SL estimators offer good performance in small samples against a range of competing estimator in the presence of spatial correlation, both in terms of bias and asymptotic variance. In particular, performance is equivalent to IV estimation when spatial correlation is absent or its impact is small: in such cases the Lasso procedure simply fails to select any eigenvectors and the resulting estimator boils down to a simple 2SLS. Should a researcher need to estimate an IV specification but then detect the presence of spatial dependence with a given SWM, our recommendation is therefore to instead run Mi-2SL as a protection against the adverse effect of that dependence on the estimates. At worst, if the spatial dependence is weak, the two estimators will produce the same estimates and thus Mi-2SL does no harm. At best, Mi-2SL will effectively control for the spatial dependence. Our empirical application, where we replicate the IV results of ck16, demonstrates the benefits of using Mi-2SL in the presence of clear spatial dependence, by improving the first-stage partial F-statistic and full F-statistic, and reducing the second stage standard errors.

Several avenues of further research involve investigating the robustness of Mi-2SL to various misspecifications. The first is the fact that the set of relevant eigenvectors $\Omega$ need to be estimated, and therefore the robustness of consistency and asymptotic normality in the presence of mistakes in eigenvector selection should be investigated. A related direction is robustness to the specification of the SWM. The results obtained here use the true SWM from the data-generating process, however in empirical settings, the true SWM is unobserved and it is likely that the empirical SWM will be misspecified in some way. Clearly, if the empirical SWM is correlated enough to the truth the Moran's $I$ test will have power to detect the correlation. However the methodology would benefit from a greater theoretical understanding of how performance will degrade as the misspecification of the empirical SWM increases.

\onehalfspacing

\doublespacing