Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
131,634 characters · 21 sections · 60 citation commands
Accounting for Unobservable Heterogeneity in Cross Section Using Spatial First Differences
\doparttoc \faketableofcontents
\onehalfspacing
\thispagestyle{empty}
\setcounter{page}{1}
We consider the problem of estimating causal effects in cross-sectional regressions when important covariates, which influence outcomes and are thought to be correlated with the treatment, cannot be observed. It is well understood that the omission of these variables may lead to substantial bias in standard regression approaches to inference. Here we propose a new cross-sectional research design that is capable of recovering such causal effects even in the presence of omitted variables. We demonstrate the performance of this approach in simulation and in two real data sets by intentionally withholding important covariates during estimation, thereby mimicking contexts with omitted variables. The first application demonstrates the ability of the procedure to recover a well established relationship (returns to schooling) while the second recovers new relationships (geographical determinants of agricultural productivities). The core insight of our approach is that unobserved heterogeneity in many cross-sectional contexts is captured by trends in outcomes across space, which can be understood as a non-parametric component of partially linear semiparametric models robinson1988root. Recognizing this, we suggest that omitted variables bias due to this heterogeneity can be eliminated from estimates using a simple and general differencing approach yatchew1997elementary in situations where the spatial position of observations can be located.
When units of observation are organized and densely packed across physical space---such as gridded data or county-level data---we propose an estimator that only compares observations to their immediately adjacent neighbors and simultaneously compares all observations to a neighbor. This approach assumes that immediately adjacent observational units are comparable to one another but does not assume that distant units are comparable, as is assumed in standard cross-sectional approaches. By restricting comparisons to adjacent neighbors in our procedure, the influence of all omitted variables that are common to neighboring units are differenced out. This fundamentally transforms the identifying assumption of cross-sectional analysis into one that matches modern quasi-experimental research designs, such as regression discontinuity designs, in terms of its mathematical structure and plausibility. Conceptually, this approach is similar to using first differences over time in a panel regression to purge data of unobserved factors specific to a panel unit, however in our case the unobserved factors are shared by two observations that are adjacent in space rather than adjacent in time. In fact, our approach is essentially identical, mathematically, to the well known first differences (FD) estimator where the key alteration is to exchange the time index of observations to an index describing the position of observations in space. For this reason, we call our research design “spatial first differences" (SFD).
In the standard cross-sectional multiple regression research design, we often study situations where an outcome of interest $y$ is influenced by a vector of $K$ observable variables $\mathbf{x}$ (“treatments”) and might possibly be influenced by $M$ unobservable variables $\mathbf{c}$ as well:
where $\epsilon$ is an i.i.d. disturbance term with mean zero. $N$ observational units are indexed by $i$. The $K$ parameters of interest are estimates of the causal effects in the vector $\mathbf{\beta}$. It is well known that if $\mathbf{c}_i$ is omitted from the cross-sectional ordinary least-squares (OLS) regression in levels (denoted by subscript $L$)
then the “omitted variables bias" in the vector of parameter estimates is $$ E[\hat{\beta}_{L}-\beta] = E[(\mathbf{x}'\mathbf{x})^{-1}(\mathbf{x}'\mathbf{c}\alpha)] $$ which may be large if the covariance between $\mathbf{x}_i$ and $\mathbf{c}_i$ is large and/or if any of the $M$ elements in $\alpha$ are large. However, because $\mathbf{c}$ is unobserved, it is not generally possible to know whether either of these conditions apply. Due to this fact, the specter of omitted variables bias now looms large over cross-sectional regression analyses and $\hat{\beta}_{L}$ is often assumed to be biased unless corroborated using an alternative research design. Thus, in many fields, $\hat{\beta}_L$ is no longer used as a basis for causal inference leamer1983let,holland1986statistics,clarke2005phantom, angrist2010credibility, regardless of how many covariates are included in the regression model.
The weakness of the standard cross-sectional research design results from how it addresses the fundamental challenge of causal inference, i.e. the estimation of plausible counterfactual outcomes holland1986statistics. For a change of $\mathbf{x}$ from $\mathbf{x}_j$ to $\mathbf{x}_i$, the average treatment effect of interest from Eq. ((ref)) is
where $E[y_i|\mathbf{x}_j]$ is the expected potential outcome for observational unit $i$ if it were treated with the $\mathbf{x}$'s of observation $j$. However, since this term is never observed in the real world, a researcher estimating Eq. ((ref)) assumes
which states that the expected potential outcome for $i$ and the outcome for $j$ would be the same if both units were treated with $\mathbf{x}_j$, which in reality was only received by $j$ and not $i$. This Conditional Independence Assumption is a relatively strong assumption in many contexts because it assumes all observational units in a cross section of data are comparable. Substituting Eq. ((ref)) into Eq. ((ref)) delivers the standard cross-sectional research design in levels, which provides an unbiased estimate of treatment effects if Eq. ((ref)) is true. However, in the presence of unobserved heterogeneity, such as the variables described by $\mathbf{c}$ in Eq. ((ref)), then the assumption in Eq. ((ref)) will not be true since units $i$ and $j$ will not longer be comparable when conditioned only on their $\mathbf{x}$'s.
We propose that the treatment effect in Eq. ((ref)) can sometimes be credibly identified in the presence of unobserved $\mathbf{c}$'s by reformulating the estimation procedure to only compare small differences in $\mathbf{x}$ and $y$ between adjacent observational units. This approach exploits a conditional independance assumption that is dramatically weaker than Eq. ((ref)) because observational units are only compared to their immediately adjacent neighbors. If $i$ is an index that matches the rank-order of observations across space in an arbitrary coordinate system, such that observations $i$ and $i-1$ are immediately adjacent to one another, then the SFD research design replaces the assumption of Eq. ((ref)) with the substantially weaker assumption
which states that the expected potential outcome for two immediate neighbors $i$ and $i-1$ are equal if they were to receive the same treatment $\mathbf{x}_{i-1}$. Eq. ((ref)) is a strictly weaker assumption than Eq. ((ref)) because the latter holds for all pairs of observations in the sample, whereas Eq. ((ref)) states that the same conditions hold only for the subset of pairs where the observations are adjacent.\footnote{It may be tempting to suggest that Eq. ((ref)) necessarily implies Eq. ((ref)) because of transitivity, but this is not the case. To see why, note that $E[y_i|\mathbf{x}_{i-1}]=E[y_{i-1}|\mathbf{x}_{i-1}]$ and $E[y_{i+1}|\mathbf{x}_{i}]=E[y_{i}|\mathbf{x}_{i}]$ do not share common terms, since the former is conditioned on $\mathbf{x}_{i-1}$ and the latter on $\mathbf{x}_{i}$.} Because Eq. ((ref)) imposes that units are only conditionally independent with respect to their local neighbors, we denote it the Local Conditional Independence Assumption.
Conditions under which Eq. ((ref)) is plausible are discussed below, but it is worth noting at the outset that this assumption is conceptually similar to (i) the assumption that immediately sequential observations within a time-series are comparable
the assumption exploited in event-study research designs and many FD time-series models; (ii) the assumption that sequential observations within a panel unit are comparable $$ E[y_{it}|\mathbf{x}_{i,t-1}]=E[y_{i,t-1}|\mathbf{x}_{i,t-1}] \;\; \forall \;\; \{t,t-1 \mid i\}, $$ an assumption exploited in differences-in-differences panel research designs\footnote{These are a subset of the assumptions required for many differences-in-differences research designs, since these approaches also often assume common trends across panel units.} (e.g. panel fixed-effects estimators); and (iii) the assumption that observations just above and just below a treatment discontinuity are comparable
the assumption exploited in regression discontinuity research designs. In fact, the SFD research design is, mathematically speaking, almost identical to the FD approach in times-series and panel analysis wooldridge2010econometric, except the one-for-one transposition of time and space indices---a similarity that allows researchers to easily implement SFD by “tricking” software packages into using time series operators on cross-sectional data sets by substituting the spatial indices for time indices. We also note that, if it is helpful, one may also think of the SFD research design “{as if}” the researcher is simultaneously running a large number of regression discontinuity research designs in space, in the sense of black1999better, with one “discontinuity” for every pair of adjacent observations in the SFD setup. Thus, overall, we argue that the Local Conditional Independence Assumption required by the SFD research design is at least as valid as corresponding assumptions in other widely accepted identification strategies, when each is applied to the appropriate context.
To illustrate how SFD brings the number of identifying assumptions needed for cross-sectional analyses into parity with the research designs described above, Figure (ref) graphically depicts the comparisons exploited to identify causal effects in different research designs. Each grid depicts a different research design, and each observation in a data set appears on both a row and column of that grid. A square is shaded grey if the two observations corresponding to that row and column are assumed to be comparable when using the associated research design (pairs are only shaded once). Panel a illustrates how only observations that are adjacent in time are compared to one another in FD time-series models (Eq. (ref)). Panel b displays how only observations with running variable values just below and just above the cutoff value $x^*$ are assumed to be comparable in a regression discontinuity design (Eq. (ref)). Panel c shows the large number of $(N-1)\frac{N}{2}$ comparisons made when using the standard “levels” approach to cross-sectional research designs (Eq. (ref)), where every observation is compared to every other observation. Panel d demonstrates how SFD reduces the number of comparisons to a strict subset of the comparisons in the levels model, since observations are only compared to those immediately adjacent in space (Eq. (ref)). The necessary $N-1$ assumptions in SFD regarding the comparability of neighbors (panel d) resembles the $T-1$ assumptions in FD time-series models (panel a); or, alternatively, $N-1$ different regression discontinuity analyses (panel b) executed across space. In contrast, the strong Conditional Independence Assumption necessary for the cross-sectional levels approach (panel c) requires exactly $\frac{N}{2}$ times as many pair-wise assumptions (compared to SFD) in order to identify $\hat{\beta}$.
We exploit the Local Conditional Independence Assumption by comparing each observation to a spatially adjacent neighbor. Critically, we difference each pair of neighboring observations to purge the data of unobserved factors that are common to the pair. We thus construct the SFD approach by writing Eq. ((ref)) for spatially adjacent observations $i$ and $i-1$, and then differencing the equations
where the $\Delta$ operator is analogous to the difference operator in time series analysis. Here we use the convention of denoting differences by the index of the higher-valued index of the pair of observations ($i$ rather than $i-1$). Because $\Delta\mathbf{c}_i$ cannot be observed, it does not appear in the SFD regression model
which can be estimated using any regression solution, such as maximum likelihood or ridge regression. We focus here on OLS, since it is straightforward to solve and its properties are widely understood. The OLS estimator for an SFD research design is
where as usual a vector of ones is included in the $\Delta \mathbf{x}$ matrix to guarantee variables are centered and errors are mean zero, although this is not explicit in our notation for parsimony.\footnote{This “constant” term is a nuisance parameter, discussed below.}
If the local conditional independence assumption Eq. ((ref)) holds, then it implies
where $0_{K,M}$ is the $K\times M$ null matrix. Eq. ((ref)) is the core identifying assumption of SFD when it is solved by OLS, stating that the covariance between changes in $\mathbf{x}$ and changes in $\mathbf{c}$ between immediately adjacent neighbors is not systematically correlated within the population.\footnote{See the Appendix for an alternative interpretation of Eq. ((ref)) based on isotropy.} The plausibility of identification via SFD depends on the validity of Eq. ((ref)), which is actually weaker than the local conditional independence assumption; although we emphasize the latter above because it is a more easily conceptualized heuristic for considering the validity of SFD in applied settings and, if true, it implies Eq ((ref)) holds.
Note that the orthogonality condition in Eq. ((ref)) describes the covariance between the local derivatives of $\mathbf{x}$ and $\mathbf{c}$ with respect to space. Thus violation of Eq. ((ref)) only occurs if the second derivatives of these variables with respect to space are sufficiently synchronized that they systematically induce changes in these derivatives in the same locations. The assumption in Eq. ((ref)) therefore differs fundamentally from the identifying assumption in the levels model (Eq. (ref)):
which describes covariance between the levels of $\mathbf{x}$ and $\mathbf{c}$ across the entire sample. In general, the failure of Eq. ((ref)) implies nothing about the validity of Eq. ((ref)), thus a levels-based cross-sectional approach may be invalid at the same time that the analogous SFD cross-sectional approach is valid. We discuss this point further below.
If the identifying assumption in Eq ((ref)) is true, then we have
so $\hat{\beta}_{SFD}$ will be unbiased. Intuitively, if $\mathbf{c}$ is common between neighbors, its influence on $y$ will be differenced out and $\Delta \mathbf{c}$ will be very near or equal zero. Should there exists a component of $\mathbf{c}$ that is not common between neighbors, $\hat{\beta}_{SFD}$ will still be unbiased so long as the non-zero component of $\Delta \mathbf{c}$ is uncorrelated with changes in $\mathbf{x}$ between neighbors ($\Delta \mathbf{x}$). Thus $\hat{\beta}_{SFD}$ is generally robust to unobservable to heterogeneity in factors that are spatially correlated ($\mathbf{c}_i\approx\mathbf{c}_{i-1}$) and factors that are i.i.d with respect to spatial position (Eq. (ref)). As we demonstrate in later sections, the SFD approach effectively purges estimates of the influence of unobserved factors $\mathbf{c}$ across a variety of applied contexts.
As described here, $\hat{\beta}_{SFD}$ falls within a class of difference estimators explored by yatchew1997elementary, yatchew1999differencing. To our knowledge, Yatchew did not discuss applying those results to an explicitly spatial context to identify causal effects---nonetheless, Yatchew's results apply here since our context is a specific case of that more general problem. yatchew1997elementary demonstrated that under mild conditions, estimation of $\hat{\beta}_{SFD}$ via OLS is consistent, yielding an asymptotic distribution
as the number of observations increases to infinity and the spatial distance between observations vanishes. Here $\sigma^2_{\epsilon}$ is the population variance of $\epsilon$ and $\Omega_{x}={E}[Cov(\mathbf{x}|\ell_i)]$ where $\ell_i$ is the spatial position of observation $i$. yatchew1997elementary also demonstrated that these variances can be estimated consistently:
As can be seen in Eq. ((ref)), $\hat{\beta}_{SFD}$ does not achieve the Cramer-Rao lower bound\footnote{yatchew1997elementary suggests an optimal combination of higher-order differencing estimates to achieve this bound asymptotically, a result that should in theory apply to the SFD context, but whose practical exploration we leave to future work.} and instead has an efficiency of 66.7% relative to this bound. In many applied context where ample data is available, we think this sacrifice of efficiency may be reasonable in order to obtain an unbiased cross-sectional estimate.
Importantly, in finite samples, the usual OLS estimator $\hat{s}^2_{\epsilon}$ in Eq. ((ref)) is not appropriate because ${\Delta {\mathbf{\epsilon}}}$ will be autocorrelated between adjacent units (eg. $\hat{\Delta {\epsilon}}_i$ and ${\hat{\Delta {\epsilon}}}_{i+1}$ both contain ${\epsilon}_i$). Thus, in practice, we recommend that the variance of $\hat{\beta}_{SFD}$ be estimated using the autocorrelation-robust approaches described by newey87 and conley1999gmm, allowing for autocorrelation in disturbances between nearby observations (in one and two-dimensions, respectively) after differencing.\footnote{Differencing requires the use of a kernel that spans at least one adjacent unit in each direction when estimating the covariance matrix for $\hat{\Delta \epsilon}$, but larger kernels may be appropriate if $\hat{\Delta {\epsilon}}$ is spatially correlated across larger scales, as in the maize example we consider below.} In the Appendix, we show that this approach generally provides the most conservative inferences relative to alternatives in our empirical application.
In the simplest case, the physical space in which observations are located has only one dimension, such as households located along a road. Panel a in Figure (ref) depicts this setup, where the $i$th observation in the differenced data set contains the change in both the treatment ($\Delta \mathbf{x}_i$) and the outcome ($\Delta y_i$) between immediately adjacent neighbors in positions $i$ and $i-1$. In this situation, the setup is directly analogous to FD in time series analysis, where the only change is that the position of an observation in time is replaced by the position of an observation along the one-dimensional space. When data is arranged in this way, it is straightforward to estimate SFD by applying basic time series functions that are standard in most statistical packages. For example, a researcher might estimate the effect of years_of_education on wages among individuals living at position house_number along a single road. Implementing Eq. ((ref)) via OLS in the statistical package Stata would then only require the two commands \footnote{In the statistical package R, the same procedure is implemented by the two commands: \\
$\quad \quad$ dplyr::arrange(data, house_number) \\ $\quad \quad$ lm(diff(wages) $\sim$ diff(years_of_education), data) \\
and in Python, one could implement this procedure using pandas after importing statsmodels.formula.api as sm: \\
$\quad \quad$ data.sort_values(by=[`house_number']) \\ $\quad \quad$ diff_data = data.diff() \\ $\quad \quad$ model = sm.ols(formula= "wages $\sim$ years_of_education", data=diff_data).fit()} \\
tsset house_number\\ regress D.wages D.years_of_education\\
where the first command “tricks” the software by telling it that the data is a “time series” where the time variable is house_number. The second command then exploits the difference operator D which computes first differences in both variables along the road and estimates the SFD regression. Whether the software is informed that “time” moves forward as one travels up or down the house_number variable is irrelevant, the SFD estimate will be the same.
Implementing SFD in two-dimensional space is similarly simple if data are “gridded” on a regular lattice, such as pixels describing topographical ruggedness or night lights. Gridded data sets of this sort are rapidly growing in availability donaldson2016, making this a particularly useful case to consider. Panel b and c in Figure (ref) depict two ways that SFD can be implemented using such gridded data: differences can be computed between neighbors defined in the East-West sense (panel b) or in the North-South sense (panel c). Neighbors that are adjacent along the dimension that is not differenced (e.g. North-South neighbors in panel b) are simply not compared. In this case, SFD is implemented in a manner analogous to FD applied to panel data, where the row (panel b) or column (panel c) of each sequence of differenced observations (indexed by $k$ in Figure (ref)) are analogous to the panel units in the FD model. A researcher interested in the effect of ruggedness on night_lights in a gridded data set where pixels are indexed by latitude and longitude could implement the East-West SFD model in Stata using the two commands\\
xtset latitude longitude\\ regress D.night_lights D.ruggedness\\
where the first command tells the software to treat latitude as if it were the panel variable and longitude as if it were the time variable in a normal panel dataset. The North-South SFD model could be similarly estimated, but switching which dimension is declared analogous to time\\
xtset longitude latitude\\
prior to estimating Eq. ((ref)).
The ability to estimate $\hat{\beta}_{SFD}$ twice along two orthogonal dimensions of a gridded data set---exploiting entirely different variation in the independent variables---provides a natural and appealing check on the robustness and validity of the two estimates since spatial patterns in omitted variables along one dimension might be different than along the other dimension.
Implementing SFD in two dimensional data when the data are gridded is straightforward, although it is somewhat more difficult to implement on data sets with irregular spatial structure. Below we demonstrate one approach that produces similar “panel-like” data structures (similar to Figure (ref)b) for the cross section of US counties. Although we first attempt to develop the readers intuition for why the procedure works, provide practical guidance, and consider the performance of SFD under simpler one-dimensional scenarios.
A central benefit of the SFD approach is that it eliminates bias due to all spatially correlated unobserved variables, which in many cross-sectional contexts represents most or all of the important omitted factors $\mathbf{c}$. For example, in a cross-sectional regression of earnings on years of schooling, if households that have high levels of education tend to live in areas with more Whites and Whites tend to earn more than other races, then race will be a spatially correlated omitted variable if it is not included in the model. SFD eliminates spatially correlated unobserved heterogeneity at two levels: the procedure filters out the influence of all factors that vary at low spatial frequencies (any factor that affects observations that are not immediately adjacent) and it differences out all common influences that idiosyncratically affect any two observations that are immediately adjacent to one another.
We think a useful way to see the benefit of SFD's high-pass filtering is to consider how $\mathbf{x}$ and $\mathbf{c}$ vary as an observer traverses the physical space in the sample, leveraging intuition and language from thinking about cross sections spanning space as analogous to time-series spanning time. Let us define some arbitrary initial position as $i=1$ (analogous to $t=0$ in time-series) and observe how $\mathbf{x}$ and $\mathbf{c}$ evolve as we move sequentially from adjacent neighbor to neighbor away from $i=1$ (analogous to moving forward in time). Then we see that $\mathbf{x}$ and $\mathbf{c}$ evolve in a “unit-root-like” manner across space because each variable is equal to the sum of its “spatial history”---i.e. all evolutions of the variables that have occurred since the initial position---and the change in the variable that occurred between the immediately previous position and the current position. Call $\tilde{\mathbf{x}}_i$ the spatial history of $\mathbf{x}_i$ where
and $\mathbf{\tilde{x}}_i$ represents the accumulation of changes since the arbitrarily defined starting point $$ \tilde{\mathbf{x}}_i = \sum_{s=1}^{i-1}\Delta \mathbf{x}_s. $$ Define the spatial histories $\tilde{\mathbf{c}}_i$ and $\tilde{y}_i$ analogously. These terms are the cumulative effect of all changes in each variable as a path from position $1$ to $i-1$ is traced out through space, analogous to a line-integral of first differences through space. Panel a of Figure (ref) illustrates this decomposition.
Taking the standard cross-sectional regression shown in Eq. ((ref)), which we hereafter refer to as the “levels” model, we know the OLS estimate for $\beta$ is $$ \hat{\beta}_{L} = \beta + \underbrace{(\mathbf{x}'\mathbf{x})^{-1}\mathbf{x}'(\mathbf{c}\alpha+\epsilon)}_{\mathrm{bias\;in\;levels\;model}} $$ where the key term that generates omitted variables bias is $\mathbf{x}'\mathbf{c}\alpha$. The bias originating from this term can be decomposed into contributions from spatial histories and spatial first differences
where the total bias depends on the size of elements in the matrices $\mathbf{W}_1$, $\mathbf{W}_2$ and $\mathbf{W}_3$, each of which is $K\times M$. $\mathbf{W}_1$ is the sample covariance between the spatial histories of $\mathbf{x}$ and all omitted factors, $\mathbf{W}_3$ is the sample covariance between their spatial first differences, and $\mathbf{W}_2$ is the sum of the cross-covariances.
In most contexts, the most important source of bias is $\mathbf{W}_1$, which is the sample covariance between the spatial histories of observable and unobservable factors. This term may be quite large, since any realizations in $\Delta\mathbf{x}_j$ and $\Delta\mathbf{c}_l$ that occur within the spatial history of observation $i$ (i.e. $j<i$ and $l<i$) and are correlated will induce correlation in $y_i$ and $\mathbf{x}_i$. Because each realization of $\Delta\mathbf{x}$ and $\Delta\mathbf{c}$ affect all observations that occur in their “spatial future,” these effects accumulate, causing correlations between $\mathbf{x}$ and $\mathbf{c}$ to sometimes grow large in finite samples, even if there is no causal or otherwise mechanical relationship between the two variables. Such correlation is clear in Panel a of Figure (ref), where the accumulation of i.i.d. realizations (Panel b) generate large correlations between $x$ and $y$ across space that have no causal meaning. In a large number of cross-sectional contexts, the bulk of correlation between $\mathbf{x}$ and $\mathbf{c}$ is captured by their spatial histories $\mathbf{\tilde{x}}$ and $\mathbf{\tilde{c}}$. In contrast, when the SFD estimator is used, spatial histories and their associated bias are eliminated by differencing (since $\mathbf{\tilde{x}}_{i}=\mathbf{x}_{i-1}$), thus all omitted variables biases attributable to $\mathbf{W}_1$ are purged. Biases due to $\mathbf{W}_2$ are also purged in the SFD approach, although this is usually only a modest improvement since this term is likely small in most contexts. Thus the magnitude of of the omitted variables bias that is eliminated from the cross-sectional regression by differencing out spatial histories is
which we show below may be substantial. If the condition in Eq. ((ref)) is satisfied, which is true if the local conditional independance assumption is valid, then $E[\mathbf{W}_3]=0$ and $\hat{\beta}_{SFD}$ is identified even if $\hat{\beta}_{L}$ is not. Note that the identifying assumption in Eq. ((ref)) does not constrain the magnitude of $\mathbf{W}_1$, so the validity of the SFD estimator provides no support for the validity of the analogous levels estimator.
The low-frequency spatial correlations between $\mathbf{x}$ and $\mathbf{c}$ is a major source of omitted variables bias in many cross-sectional settings, however spatial correlations between $\mathbf{x}$ and $\mathbf{c}$ with high spatial frequencies (analogous to local “shocks”) may also be problematic in the traditional levels model. For example, if a single unobserved hospital is particularly good at providing healthcare, then average health for individuals in adjacent neighborhoods may be idiosyncratically high, possibly confounding regressions of health on $\mathbf{x}$. Because the SFD estimator only exploits changes in $\mathbf{x}$ and $y$ that occur between neighboring observations, any localized change in $\mathbf{c}$ that affects both locations $i$ and $i-1$ (e.g. the presence of a hospital) is differenced out when $\Delta y_i$ is constructed.
Thus, in the SFD model, there remains no spatial correlations in unobservables left to bias the estimated parameters. Spatial correlations that affect observations more than one unit away from one another (e.g. $i$ and $i-2$) are purged because the spatial histories of observations are eliminated from the data (observation $i$ does not “know” about the existence of observation $i-2$) and spatial correlations that affect adjacent observations (e.g. $i$ and $i-1$) are differenced out.
SFD is a distinct research design that differs from, but sometimes links to, other approaches that also consider spatial relationships. Conceptually, the assumptions of SFD generalize the assumptions used in spatial regression discontinuity (RD) black1999better so that obtaining arbitrary boundaries is not a requirement and variation throughout a sample can be exploited. SFD uses spatial relationships strictly to organize observations for the purpose of identification and does not rely on the tools used to handle spatial dependence or spatial autocorrelation developed elsewhere anselin1988spatial. However, we recommend estimating SFD standard errors using the approach described by newey87 or conley1999gmm, in one and two-dimensional contexts, respectively, to account for the possibility that $\hat{\Delta \epsilon}$ is spatially autocorrelated. Additionally, in the presence of spatial spillovers, it is important that spatial lags are included in the regression model (as additional elements in the vector $\mathbf{x}$) before spatial first differences are computed. At a high level, SFD can be thought of as a simple and general approach to identifying partially linear semiparametric models robinson1988root with unknown omitted variables in contexts where observations are dense and organized across space. We describe the connections of SFD to these different models in greater detail below.
\paragraph{Spatial regression discontinuity research designs} Spatial RD designs, following holmes1998effect, black1999better and others, exploit arbitrary borders in space (e.g. state borders or school districts) and are now widely used in applied research for causal inference. Spatial RD designs rely on the assumption that observations on immediately opposing sides of an arbitrary spatial boundary are conditionally comparable to one another in factors that are not determined by the boundary, an assumption that is the same as the Local Conditional Independence Assumption (Eq. (ref)) required by SFD, albeit restricted to apply only to observations near the exploited border. In spatial RD designs, this assumption is often stated as the requirement “that all confounding factors vary smoothly `enough' in space across the boundary.” This assumption is important because it implies that changes in treatment $\mathbf{x}$ across the boundary are not correlated with changes in confounders $\mathbf{c}$, a statement which written mathematically is precisely the orthogonality condition in Eq. ((ref)).
Thus the underlying assumptions of SFD can be interpreted as a generalization of the assumptions in spatial RD, since in SFD the comparability of immediately adjacent neighbors---only explicitly assumed to apply near a the boundary in spatial RD---is assumed to hold everywhere. Intuitively, if one accepts the spatial RD approach, then extending the assumed local comparability of adjacent observations to units far from the boundary may be natural, since it is often the case that adjacent units far from a boundary (e.g. adjacent counties inside the same state) are at least as comparable to one another as adjacent units on opposite sides of a boundary (e.g. adjacent counties in different states). Thus, SFD can be thought of “as if” it generalizes spatial RD, where small “discontinuities” between all pairs of neighbors throughout the sample (not only at borders) are exploited; although the correspondence is more true conceptually than it is exact mathematically, since spatial RD is implemented in a variety of creative ways that do not all map exactly to SFD.\footnote{The generalization of spatial RD using SFD can be seen most easily in one-dimensional space. The univariate spatial RD approach can be rewritten as the SFD estimate, where the treatment $\mathbf{x}=1$ if observations are on one side of a border and $\mathbf{x}=0$ on the other side of the border. However, not all implementations of spatial RD are so simple. For example, in some cases spatial RD is implemented by discarding observations far from the boundary, while in other cases these observations are kept and the outcome is allowed to depend flexibly on the running variable (space) far from the boundary, perhaps using splines or nonparametric approaches.}
\paragraph{Models of spatial dependence} SFD is fundamentally distinct from the class of methods that collectively compose the field of “spatial econometrics” in which regression models are specifically structured to account for relationships that manifest over space across different observational units\footnote{anselin1988spatial explicitly defines spatial econometrics as “the collection of techniques that deal with the peculiarities caused by space in the statistical analysis of regional science models”.} anselin1988spatial, lesage2009introduction. Core methods developed in spatial econometrics account for spatial dependence---where the outcome in one location (${y}_i$) influences the outcome in a nearby location (${y}_{i+1}$), and visa versa---as well as spatial spillovers---where the treatment at one location ($\mathbf{x}_i$) influences the outcome in a nearby location (${y}_{i+1}$)---and numerous tools to measure and model various patterns of spatial autocorrelation of disturbances explicitly LeSage2010. In contrast, SFD is simply a research design that exploits the spatial structure of data instrumentally in order to identify the average causal effect of a unit's observable treatments ($\mathbf{x}_i$) on that same unit's outcome (${y}_i$). Unlike spatial econometric models, Eq. ((ref)) does not contain any explicit relationships that depend on space, since neither neighbors' treatments nor neighbors' outcomes are on the right-hand side. Space is only used to identify $\beta$ by informing how this equation is transformed via Eq. ((ref)). Notably, however, in principle one could write down a model from spatial econometrics and then apply SFD to that model as an identification strategy in the appropriate context.
\paragraph{Models of spatial spillovers} In many contexts, a treatment at location $i$ may alter an outcome at some nearby location $j$, where the magnitude of the effect decays with distance. For example, construction of a police station may affect crime locally and in nearby locations. Such spatial spillovers are generally accounted for using “spatial lag models" that recode the neighbor's treatment as a new type of treatment variable which is then introduced as an additional independent variable in the original regression. For example, in the model $y_i={x}_i\beta + {x}_{i-1}\gamma+\epsilon_i$, ${\gamma}$ is the spillover from location $i-1$ to $i$. In a more general model, lags in multiple spatial directions and over different distances may be accounted for. Descriptions of spatial lag models and SFD may both draw on language or intuition from time-series\footnote{Descriptions of spatial lag models sometimes invoke similarities to distributed lag time-series models.}, but are employed for unrelated reasons. Spatial lag models may account for a particular data-generating process in which spillovers are present, while SFD is an identification strategy that may be employed regardless of whether spillover are present or not.
One might imagine that the presence of spatial spillovers necessarily invalidates the the SFD research design because contamination of nearby control groups (i.e. violation of the Stable Unit Treatment Value Assumption) is an obstacle to identification in some spatial settings, such as spatial RD designs, but this is not a problem for SFD. Unlike spatial RD designs, there may be substantial variation in neighbor's treatments within a sample so that spatial spillover effects can themselves be identified separately from own treatment effects. To see this, note that $i$'s neighbor's treatment (e.g. $x_{i-1}$), which generates a spillover effect to $i$, differs from $i$'s neighbor's {neighbor's treatment} (e.g. $x_{i-2}$), which generates a spillover effect onto $i$'s neighbor. These small differences in adjacent values for the neighbor's treatment variable are then exploited in the SFD design to identify the spillover effect. This intuition generalizes for multiple directions and any distance. Thus, the SFD research design can be used in combination with a spatial lag model without complication if spillovers are thought to be present, all that is required is that the standard spatial lag specification (in levels) is differenced across space.\footnote{ The general spatial lag model is $$y_i={x}_i\beta + {L}_i(\mathbf{X})\gamma+\epsilon_i$$ where $\mathbf{X}$ is the matrix of all treatments for all locations, ${L}_i(.)$ is a spatial lag operator relative to the position of $i$ and $\gamma$ is a vector of estimated spillover effects. Note that ${L}_i(.)$ could encode spatial lags in multiple directions and distances and spillovers need not be isotropic. The SFD model in this general case is then $$\Delta y_i=\Delta {x}_i\beta + \Delta {L}_i(\mathbf{X})\gamma+\Delta \epsilon_i.$$ } For example, in the simple model of spillovers above, the SFD estimating equation is simply $\Delta y_i=\Delta {x}_i\hat{\beta}_{SFD} + \Delta {x}_{i-1}\hat{\gamma}_{SFD}+ \hat{\Delta\epsilon}_i$ which will recover unbiased estimates of ${\beta}$ and ${\gamma}$. In general, the precise form of spillovers need not be known ex ante to achieve unbiasedness, since one can always include “more than enough” spatial lags and irrelevant lags will be uncorrelated with the outcome. With sufficient spatial lags included, the orthogonality condition in Eq. ((ref)) will be satisfied.
Importantly, however, a mistaken omission of spatial lags in the SFD research design when spatial spillovers are present may generate bias. For example, suppose the data generating process is $y_i={x}_i\beta + {x}_{i-1}\gamma+\epsilon_i$ where $\gamma \neq 0$. If the SFD model is not the model described above but instead omits the spatial lag term, estimated as $\Delta y_i=\Delta{x}_i\hat{\beta}_{SFD}^\dagger +\hat{\Delta\epsilon}_i$, then the estimate $\hat{\beta}_{SFD}^\dagger$ will be biased. This occurs because ${x}_{i-1}$ behaves as a source of unobserved heterogeneity (${x}_{i-1}=\mathbf{c}_i$ in Eq. (ref)) and the orthogonality condition Eq. ((ref)) fails for this particular form of heterogeneity. Unbiasedness of $\hat{\beta}_{SFD}^\dagger$ would require the orthogonality condition $\mathrm{E}[\Delta {x}_{i}\Delta {x}_{i-1}]=2\mathrm{E}[{x}_{i}{x}_{i-1}]-\mathrm{E}[{x}_{i}{x}_{i-2}]-\mathrm{E}[{x}_{i-1}^2]=0$, which will almost never be true; even if there is no spatial autocorrelation in ${x}$, the last term is a variance that will not equal zero. Similar issues arise with spatial spillovers that have other spatial structures. Conceptually, this issue arises because first differencing embeds information regarding $\mathbf{x}_{i-1}$ into the independent variable $\Delta\mathbf{x}_i$, with estimated coefficient $\hat\beta_{SFD}$, while at the same time information from $\mathbf{x}_{i-1}$ persists in the first differenced error term (because of the spillover) if $\Delta\mathbf{x}_{i-1}$ is not also made a regressor. This contrasts with the standard levels model (Eq. (ref)) because $\mathbf{x}_{i-1}$ is not incorporated into the regressor with coefficient $\hat\beta_L$, so omission of $\mathbf{x}_{i-1}$ from the regression does not necessarily induce this mechanical bias. Notably, though, if there is any spatial auto-correlation in $\mathbf{x}$---which is almost always true in spatially dense data---and spatial spillovers are present, then estimating the levels regression without spatial lags will also be biased but for a different reason.\footnote{For example, if the data generating process is $ y_i={x}_i{\beta} + {x}_{i-1}\gamma+\epsilon_i$ as above, then estimating a levels models that omits the spatial lag $ y_i={x}_i\hat{\beta}_{L}^\dagger +\hat{\epsilon}_i$ requires the orthogonality condition $\mathrm{E}[{x}_{i}{x}_{i-1}]=0$, which fails if $x$ is spatially auto-correlated.} Thus, inclusion of spatial lags is important for preventing omitted variables bias when spatial spillovers are thought to be important, regardless of whether a cross-sectional model is estimated in levels or SFD.
\paragraph{Spatial autocorrelation robust standard errors} Recently in applied work, there has been increasing attention to the role of spatial autocorrelation in disturbances when estimating uncertainty in cross-sectional models, particularly in the approach developed by conley1999gmm. In Conley's procedure, the structure of spatial autocorrelation is accounted for by computing $\hat\epsilon_i\hat\epsilon_j$ for nearby observations when estimating a covariance matrix, similar to the analogous time-series procedure developed by newey87. However, in the SFD research design, all of these cross-covariances in levels are immaterial because common information between units has been differenced out, as discussed above. Nonetheless, it is possible that the off-diagonal terms in $E[\hat{\Delta \epsilon}\hat{\Delta \epsilon}']$ could be nonzero if there are higher-order aspects of the data-generating process that generate spatial correlations in the gradients of disturbances. In such a scenario, it is appropriate to apply the procedure in conley1999gmm to the spatially first-differenced data, which in one-dimensional space is identical to the procedure in newey87. To avoid possible over-rejection of the null, we recommend such an approach when one is unsure about the spatial autocorrelation of $\hat{\Delta \epsilon}$ and we use this approach in empirical examples below.
\paragraph{Semiparametric regression models} The SFD design can be understood as a specific case of partially linear semiparametric regression where the vector of unobservable variables $\mathbf{c}$ is unknown and the dependance of $y$ on $\mathbf{c}$ is governed by an unknown function $g(\mathbf{c})$. In its usual formulation, the semiparametric model that would replace Eq. ((ref)) is written $y = \mathbf{x}\beta + g(\mathbf{c})+\epsilon$ robinson1988root, carroll1997generalized. We point out that if unobserved covariates $\mathbf{c}$ are functions of space (with observation $i$ at position $\ell_i$, for consistency with the sections above), then a general solution is to rewrite $g(\mathbf{c})=g(\mathbf{c}(\ell_i))=g^*(\ell_i)$ and account for unobservables by estimating a partially linear model that is nonparametric over positions:
Thus, the unobservable cross-sectional heterogeneity described by $\mathbf{c}_i\alpha$ in Eq. ((ref)) is captured by the nonparametric component $g^*(\ell)$ in the semi-parametric model. SFD leverages the idea that physical space can be used as a metric in which to organize and index observations in order to remove any confounding influence of $g^*(\ell)$. The SFD solution to estimating Eq. ((ref)) is to first-difference across space, so that $g^*(\ell)$ is differenced out. The idea of estimating the linear component of partially linear models through differencing was first proposed by yatchew1997elementary based on similar intuition, although, to our knowledge, previous literature has not proposed indexing observations based on physical location for this purpose.
An alternative approach to estimating Eq. ((ref)) would be the procedure proposed by robinson1988root. Implementing Robinson's approach would involve estimating smooth non-parametric “trends” in $y$ and $\mathbf{x}$ across positions $\ell_i$ using kernel estimators, and then regressing the resulting residuals of $y$ on the residuals of $\mathbf{x}$ in order to recover $\beta$. For example, the “spatial fixed effects" approach proposed by conley2010 is a creative implementation of Robinson's approach applied in space, where the kernel used to estimate $g^*(.)$ has uniform mass for all observations within a cutoff radius.\footnote{We demonstrate this implementation in Appendix (ref) and compare the results to SFD for our one-dimensional empirical example.} Under such an approach, to achieve identification of $\beta$ one would need to assume that all units $j$ near enough to $i$ to inform these kernel estimates at $\ell_i$ are comparable to $i$. This assumption may serve well in some contexts, but may be difficult to defend if the elements of $\frac{\partial \mathbf{c}}{\partial \ell}$ (and thus possibly $\frac{\partial g^*}{\partial \ell}$) are hypothesized to be highly variable across space, to exhibit unknown discontinuities, or to be anisotropic in the neighborhood of $\ell_i$. SFD circumvents this assumption, which may be useful if, for example, crucial dummy variables that would otherwise capture level-shifts in $g^*(\ell)$ were omitted from the model (in the SFD approach, effects of omitted dummy variables will nonetheless be absorbed even if they are not explicitly accounted for). Differencing can be thought of as restricting the bandwidth of the kernel estimator in the first stage of Robinson's procedure to its very smallest possible value, such that it contains only a single neighboring observation. In this sense, one could think of SFD “as if” it applies Robinson's procedure by estimating $g^*(\ell)$ using a completely nonparametric “spatial trend” that includes a single dummy variable for every single pair of neighboring observations. However, critically, implementing such a procedure using dummy variables is not identified for $N$ adjacent observations, since it would require estimating at least $N$ parameters ($N-1$ dummy variables and $\hat{\beta}$). Thus, SFD is the only feasible approach that also allows for this level of flexibility in the possible structure of $g^*(\ell)$. Although, notably, the cost of this flexibility is a loss of efficiency. As shown by yatchew1997elementary, Robinson's approach converges faster.
The design of the SFD estimator emerges naturally from recognizing that adjacent neighbors in a sample may be comparable, although exact comparability across every single pair of neighbors is not actually essential for SFD to provide identification. If the Local Conditional Independence assumption (Eq. (ref)) is true, then the identifying assumption that $\Delta \mathbf{x}$ and $\Delta \mathbf{c}$ are orthogonal (Eq. (ref)) will hold. However, Eq. ((ref)) may remain valid under weaker conditions in actual data. Here we discuss some practical issues to consider when determining whether SFD may provide causal identification in different contexts.
\paragraph{The spatial density of observation} The assumption that adjacent observations in a data set are comparable relies on the notion that observations that are adjacent in a data set are actually “nearby” one another in the real world. If adjacent observations in a data set are extremely far from one another in actual space, then it may not be reasonable to assume they are comparable. For example, in a sample of countries, China and Russia are adjacent, but from the standpoint of many economic questions, they may not be directly comparable. One reason they do not seem comparable, intuitively, is that there are many portions of Russia that are extremely distant from many portions of China, despite the existence of a common boundary between the two units. Thus, when assessing whether SFD is an appropriate approach for a given sample, it is important to consider whether adjacent observational units are nearby one another in their entirety, especially with respect to the range of spatial coverage across the sample. For this reason, it is likely that local conditional independance will be most nearly true, and SFD well identified, if the spatial extent of individual observational units is limited and their spatial density is high. We cannot provide precise guidance on how dense is “dense enough” in practice, rather, it is up to the judgement of the econometrician to determine whether the spatial density of observations is sufficiently high that Eq. ((ref)) (or at least Eq. (ref)) is likely to be satisfied.
\paragraph{Considering potential common causes of both regressors and omitted variables}
A natural question to consider is whether the SFD estimator will be identified in cases where some unobserved variable $z$ influences both the regressors $\mathbf{x}$ and the unobserved variable $\mathbf{c}$: $$\mathbf{x}_i = \mathbf{x}(z_i),\;\;\; \mathbf{c}_i = \mathbf{c}(z_i).$$ In this case, as the sample is traversed in physical space, $z$ will evolve as a function of location, driving changes in both $\mathbf{x}$ and $\mathbf{c}$ along the way. For example, higher elevations ($z$) across counties might cause temperatures (${x}$) to fall and the air to be thinner (${c}$), both of which might in turn affect crop yields ($y$). In a levels model, the existence of such an external factor $z$ might generate correlation in $\mathbf{x}$ and $\mathbf{c}$, thereby inducing bias in $\hat{\beta}_{L}$. However, $\hat{\beta}_{SFD}$ is substantially more robust to this scenario, in the sense that such a bias is much less likely even if an unknown common cause exists, because a violation of the identifying assumption occurs only if a fairly restrictive condition on the functions $\mathbf{x}(z)$ and $\mathbf{c}(z)$ is met. Specifically, a common cause $z$ generates bias if the curvature of the functions $\mathbf{x}(z)$ and $\mathbf{c}(z)$ mirror one another throughout the support of $z_i$ in a sample. For example, if $\mathbf{c}(z)$ is concave in $z$ and $\mathbf{x}(z)$ increases linearly in $z$ (Figure (ref)a) or exhibits higher order variations as a function of $z$ (Figure (ref)b), then it is likely SFD will be unbiased---even though $\mathbf{x}$ and $\mathbf{c}$ are correlated in levels across $z$.
To see why this is true in general, note that for sufficiently small changes in $z$ between neighbors we obtain the first-order Taylor approximations $$ \Delta \mathbf{x}_i = \frac{\partial \mathbf{x}(z_{i-1})}{\partial z} \Delta z_i, \;\;\; \Delta \mathbf{c}_i = \frac{\partial \mathbf{c}(z_{i-1})}{\partial z} \Delta z_i. $$ We use these approximations to rewrite the local conditional independence condition in Eq. ((ref)) as
where, for clarity, we explicitly write out the subtraction of the sample average of $\Delta\mathbf{x}$, denoted $\overline{\Delta\mathbf{x}}$, which is partialled out by the vector of ones in the SFD regression. If $\Delta z_i$ is held fixed\footnote{Note that for vanishingly small changes in physical position $\Delta \ell \rightarrow 0$ , $\Delta z$ is locally constant ($\Delta z = \frac{\partial z}{\partial \ell}\Delta \ell$ by Taylor's theorem).}, this expression indicates that SFD will fail to be identified only if the demeaned derivative of a regressor with respect to the common cause ($\frac{\partial x}{\partial z}-\overline{\frac{\partial x(z_{i-1})}{\partial z}}$) is correlated with the derivative of the omitted variable with respect to the common cause ($\frac{\partial c}{\partial z}$) across the values of $z$ in the sample. In order for such a violation to occur, the second derivatives of $\mathbf{x}(z)$ and $\mathbf{c}(z)$ must correspond across the support of $z$, thereby generating correlation between the first derivatives in Eq. ((ref)). Such a situation is shown in Panel c of Figure (ref), where Eq. ((ref)) is likely violated. In this case, as space is traversed and $z_i$ varied, changes in $\Delta\mathbf{x}-\overline{\Delta\mathbf{x}}$ and changes in $\Delta\mathbf{c}$ will occur at exactly the same positions, leading to correlation in spatial first differences of these variables. Nonetheless, in most contexts, we believe scenarios similar to Panels a and b of Figure (ref) are generally more likely for most common causes one might postulate.
In practice, to satisfy Eq. ((ref)), it is generally sufficient for the curvature of $\mathbf{x}(z)$ and $\mathbf{c}(z)$ to differ over the range of $z$ plausibly contained within a sample, not everywhere in $z$. Of course, it is also possible for Eq. ((ref)) to be violated if the square of the first difference in $z$ is correlated with the product of the two derivatives, although we think such conditions are relatively exotic for most applied settings.
\paragraph{Rotation of coordinates}
Because $\mathbf{c}$ is unobserved, it is impossible to directly test whether Eq. ((ref)) holds. Nonetheless, we are able to offer a practical indirect check that is implementable using the data. We suggest “cross-checking” the identifying assumption of SFD across multiple implementations of the estimator. In Figure (ref), we pointed out that in two-dimensional environments, SFD could be implemented in the East-West direction or the North-South direction, providing two separate estimates of $\hat{\beta}_{SFD}$ that can be compared to one another for consistency. Results that match would suggest that either the East-West and North-South versions of Eq. ((ref)) both are true, or they both fail and somehow generate bias of similar structure despite exploiting different sources of variation in $\Delta\mathbf{x}$. We point out that in many environments, this intuition can be further generalized by noting that if Eq. ((ref)) holds in general, an econometrician should be able to estimate SFD by taking differences along an axis that has been rotated in space by an arbitrary angle $\theta$ and the resulting estimate $\hat{\beta}_{SFD}(\theta)$ should be relatively invariant across $\theta$. In our analysis of US maize yields below, we demonstrate this test for 180 estimates of key parameters as the coordinate system we use is rotated through $\theta=-89^\circ$ to $\theta=90^\circ$ by $1^\circ$ increments.
\paragraph{Spatial Double Differences}
Another indirect check that is implementable using the data and which requires different identifying assumptions involves taking higher order differences. Differencing Eq. ((ref))
we obtain the double difference of Eq. ((ref)) where $\Delta^2 \mathbf{x}_i = \Delta \mathbf{x}_i - \Delta \mathbf{x}_{i-1}$. Thus we propose the “Spatial Double Differences" (SDD) cross-sectional regression model
which can be estimated via OLS and will be unbiased so long as
which is a distinct restriction that is neither implied by nor a result of Eq. ((ref)). We suggest that in cases where one is uncertain that Eq. ((ref)) is true, an econometrician might estimate Eq. ((ref)) and compare whether $\hat{\beta}_{SDD}=\hat{\beta}_{SFD}$. If these estimates are very close then, similar to the comparison with a rotated coordinate system, in order for $\hat{\beta}_{SFD}$ to be biased by a failure of Eq. ((ref)), one must postulate a structure for $\mathbf{c}$ such that $\hat{\beta}_{SDD}$ is identically biased by the failure of Eq. ((ref)). If such a situation is deemed unreasonable, then it is likely that both Eq. ((ref)) and Eq. ((ref)) are true. Given sufficient data, simultaneous failure of these conditions but identical SFD and SDD estimates is difficult to achieve in practice, leading us to argue that this is a strong test. However, in many data environments, it maybe be challenging to estimate $\hat{\beta}_{SDD}$ because the variation in double-differenced variables is likely to be noisy, which can lead to imprecise estimates of $\hat{\beta}_{SDD}$ and attenuation bias. Nonetheless, we demonstrate implementation of this test below in our analysis of US maize yields.
The core benefit of the SFD research design is elimination of omitted variables bias induced by spatial correlations between regressors and omitted variables in cross section. The magnitude of this benefit is described by Eq. ((ref)), computed by comparing $\hat{\beta}_L$ and $\hat{\beta}_{SFD}$. In general, this difference will be dominated by the unobservable $\mathbf{W}_1$ term. In this section, we provide simple examples in one-dimensional space that demonstrate the magnitude of this benefit by comparing cross-sectional regression results using the standard levels model and SFD. We first consider an idealized simulation using synthetic data where the structure of the omitted variable bias is known exactly. Then we consider estimates for the returns to schooling in census blocks along 10th Avenue in New York and I-90 in Chicago, where omitted variables are not known but plausible values of $\beta$ have been well-documented and replicated in prior analyses, providing a benchmark for comparison.
In this simple simulation, we generate synthetic data to compare the performance of SFD to that of levels in the presence of a single, known omitted variable. This exercise demonstrates some conditions under which the SFD estimator performs well. The data generating process is as follows. Let $i = 1, \dots, 1000$ index evenly spaced observations along a line (note that here, $i = \ell_i$). We are interested in the outcome variable $y$ that is determined by $x$, which is observed, and $c$, which is not observed: $$ y_i = x_i\beta + c_i\gamma + \epsilon_i \label{Sim_Model} $$ where $\beta = \gamma = 1$ and $\epsilon_i \sim N(0,1)$. In order to allow us to smoothly vary the degree of spatial correlation between $x$ and $c$, we exploit sinusoidal functions (in degrees, not radians) and vary their wavelength. Specifically, we generate $$ x_i = \mathrm{sin}(i) + \delta_i \phi; \label{Sim_X} $$ $$ c_i = \mathrm{sin}\Big( \frac{360 i }{\lambda}\Big) + \eta_i \phi ; \label{Sim_C} $$ where $\delta_i$ and $\eta_i$ are disturbance terms that are both independently distributed $N(0,1)$. Throughout the simulation, the expected value of $x$ completes one cycle every $360$ observations, so its wavelength is $360$. The wavelength of $c$ is controlled by $\lambda$, such that $x$ and $c$ are most highly correlated when $\lambda = 360$. The noise terms, $\delta_i$ and $\eta_i$, are amplified by the parameter $\phi$, which we will also vary. We run 1,000 repetitions for each parameterization of the problem, defined by the values of $\lambda$ and $\phi$. The three subplots to the left in Figure (ref) show the explanatory variable $x$, the omitted variable $c$, and the outcome variable $y$ for the parameterization $\lambda = 360$ when $\phi = 0$ (panel a), $\phi = 0.04$ (panel b), and $\phi = 0.5$ (panel c).
Since $c$ is unobserved, we estimate $\beta$ (true value = 1) by regressing $y$ on $x$ and intentionally omit $c$ from the model. For each simulation we estimate the levels model $$ y_i = \hat{\alpha}_1 + x_i\hat{\beta}_L + \hat{\epsilon}_i \label{Sim_Levels} $$ and the SFD model $$ \Delta y_i = \hat{\alpha}_2 + \Delta x_i\hat{\beta}_{SFD} + \hat{ \Delta \epsilon}_i \label{Sim_SFD} $$ and record $\hat{\beta}_L$ and $\hat{\beta}_{SFD}$.
The results for three values of $\phi = \{0, 0.04, 0.5\}$ and integer values of $\lambda = \{1, \dots, 800\}$ are displayed in the right panels of Figure (ref). Each subplot displays the coefficient estimates as a function of $\lambda$, with the levels estimates shown in orange and the SFD estimates shown in blue. The lightly shaded areas show the inner 95% range of estimates over the 1,000 simulations, while the darker lines show the average across all 1,000 estimates. In the cases with no noise ($\phi = 0$; panel a), the mean $\hat{\beta}_{SFD}$ is very near to the mean $\hat{\beta}_{L}$ because $x$ and $c$ are perfectly correlated. Both estimators perform well for small values of $\lambda$, since at these values $x$ and $c$ have a low degree of spatial correlation. However, the two estimators perform poorly for $\lambda = 360$ because the observed variable and the omitted variable are perfectly correlated across space. In this situation, all the influence of omitted variable $c$ is picked up by $\hat{\beta}$ in both levels and SFD, creating large biases.
However, with even a minute amount of noise, $\Delta x$ and $\Delta c$ become uncorrelated in expectation and the bias becomes considerably smaller for the SFD estimate than for the levels estimate. This result is demonstrated in panel b of Figure (ref). When $\phi = 0.04$ (panel b), the bias of $\hat{\beta}_{SFD}$ is essentially zero, whereas the absolute bias for $\hat{\beta}_L$ is 0.08, on average across $\lambda$. An important trade-off in this case is that the efficiency of $\hat{\beta}_{SFD}$ relative to $\hat{\beta}_{L}$ is low ($\frac{Var(\hat{\beta}_{L})}{Var(\hat{\beta}_{SFD})}=0.08$). However, as the amplitude of noise increases, SFD remains less biased than levels, but its variance declines substantially relative to the levels estimator. This relative decline occurs because the presence of noise in $x$ increases the variation that can be exploited in first differences, reducing $(\Delta x' \Delta x)^{-1}$. When the quantity of noise is modest ($\phi = 0.5$; panel c), the bias of $\hat{\beta}_{SFD}$ is less than 0.2% that of $\hat{\beta}_{L}$, on average across $\lambda$, and the relative efficiency of $\hat{\beta}_{SFD}$ to $\hat{\beta}_{L}$ is much better at 0.5.
The difference between the SFD estimate and the levels estimate is most striking when the variable of interest and the omitted variable are highly spatially correlated, for the cases with non-zero $\phi$. While the levels estimator always performs badly when $\lambda$ is in the neighborhood of 360, such that $x$ and $c$ are perfectly in phase (as shown in the small subplots), the SFD estimator is remarkably unbiased, performing no worse than for other $\lambda$. For example, when $\lambda = 360$ and $\phi = 0.5$, the average value of $\hat{\beta}_{SFD}$ is $1.00$ (95% interval: $0.84 \leq \hat{\beta}_{SFD} \leq 1.16$), while the average value of $\hat{\beta}_{L}$ is $1.67$ (95% interval: $1.59 \leq \hat{\beta}_{L} \leq 1.75$).
This simple exercise illustrates the conditions under which the SFD estimator performs well. When the regressor of interest and the omitted variable are highly spatially correlated, the SFD estimator is remarkably less biased than the levels estimate. When these two variables are not as strongly correlated across space, the bias in SFD is no worse than levels. There is a tradeoff between bias and efficiency though, so the SFD estimate tends to have a larger variance. Nevertheless, so long as there is sufficient orthogonal variation in the independent variables, the efficiency of the SFD estimator may be comparable to that of the levels estimator. Next, we compare the results from the two estimators in a one-dimensional empirical example where plausible parameter values have been previously established.
In real world contexts, we do not observe omitted variables, making it impossible to know for certain how close $\hat{\beta}_{SFD}$ is to $\beta$. But what we can do, at least for an illustrative example, is estimate $\hat{\beta}_{SFD}$ for a relationship that is sufficiently well studied that we have some sense of the true value for $\beta$. In the following example, we demonstrate how SFD is implemented in one-dimensional space by conducting a simple analysis similar to the returns to schooling example discussed above, examining census blocks along the longest roads in Manhattan and Chicago. We then compare the $\hat{\beta}_{SFD}$ and $\hat{\beta}_{L}$ estimates with these samples to previous estimates of the returns to education.
We obtain data on average wages and average years of education from the 2010 American Community Survey 5-year estimates for New York City and Chicago ACSdataset. In New York, we produce a sample of 53 adjacent observations by following 10th Avenue from the lower to the upper tip of Manhattan and recording the census tracts along this path, as depicted in Figure (ref) (panel 1a). We follow the same procedure in Chicago, tracing Interstate-90 from the northwestern to the southeastern corner of the city and obtaining 54 sequential observations (panel 2a). We arrange the observations in order, such that we have a one-dimensional sequence of adjacent census blocks. Panels 1b and 2b of Figure (ref) show the levels of log weekly wages (grey) and years of schooling (blue) that would be observed driving north along 10th Avenue in Manhattan and southeast along I-90 in Chicago, respectively. For comparison, Panels 1c and 2c show the spatial first differences in log weekly wages (grey) and years of schooling (blue) along the same paths.
Indexing census tracts by $i$ according to their position along these two roads, we estimate the effect of years of education on wages via OLS, intentionally omitting all covariates. The levels model is
and the SFD model is
where $\hat{\alpha}_1$ is the intercept in the levels model and $\hat{\alpha}_2$ is the intercept in the SFD model.\footnote{We estimate returns to schooling on this dataset using the semi-parametric model proposed by Robinson (1998) in Appendix (ref).}
The results are shown in Table (ref). The semi-elasticities estimated using levels are 0.179 in Manhattan and 0.125 in Chicago, suggesting that each additional year of schooling increases wages by 18% and 12.5% in these cities, respectively. These values are larger than almost all previous estimates of the return to education. Card (2001) reports 17 previous estimates of the return to education in the United States, which range from 0.052 to 0.132 card2001estimating. The levels estimate in Manhattan is much larger than all of these estimates, and the estimate in Chicago is larger than all but one. In contrast, the coefficients we estimate using SFD, 0.089 in New York and 0.072 in Chicago, are in the center of the distribution of previous estimates. For comparison, the OLS and IV estimates from staiger1997instrumental, the largest study reviewed in Card (2001), are also reported in Table (ref).
The difference between the levels and SFD estimates reported in Table (ref) is the magnitude of the omitted variables bias that is eliminated from the cross section by differencing out spatial histories, as shown in Eq. ((ref)). As displayed clearly in Figure (ref), panels 1b and 2b, the correlation in spatial histories ($\mathbf{W}_1$) is large. The positive sign of this difference seems sensible when considering which omitted variables are likely to generate bias in the levels estimate. For example, if households that have high levels of education tend to live in areas with more Whites and Whites tend to earn more than other races, then race will be a spatially correlated omitted variable that will lead to upward bias in the levels estimate.
This analysis demonstrates that we recover well established estimates of the return to education using SFD even when {all} covariates are intentionally omitted from the analysis. In two different cities, with a modest sample size determined only by the location of the longest road, we produced estimates of the return to education that appear precise and match previous estimates.
We now develop new cross-sectional estimates for the effect of long-term climate and average soil conditions on maize yields in US counties, relationships notoriously confounded by unobserved heterogeneity, using SFD. Further, this rich data set allows us to explore three additional points: (i) implementation of SFD using irregular (non-gridded) two-dimensional data; (ii) systematic evaluation of the research design's vulnerability to omitted variables bias using an “extreme bounds" analysis leamer1985sensitivity; and, (iii) implementation of two novel robustness tests unique to SFD, the rotation of the coordinate system and SDD, which enable us to check the research design's underlying assumptions.
Understanding the impact of soil and climatic endowments on economic outcomes has become an important area of research in development economics, economic geography, and environmental economics, motivated by strong cross-sectional patterns of development and, more recently, climate change bloom1998geography, acemoglu2000colonial, easterly2003tropics, nordhaus2006geography, carleton2016social. However, analyses of economic responses to long-run geographic conditions, such as average temperature, are potentially affected by omitted variables bias since comparing outcomes in areas with different conditions directly (e.g. hot vs. cold locations) has long employed cross-sectional research designs. In this literature, unobserved heterogeneity has been identified as a key issue in particular, as climate variables are thought to be spatially correlated with other variables (e.g. ruggedness, political institutions) that may not be observed but likely affect economic outcomes deschenes2007economic, hsiang2016climate. As with other cross-sectional contexts, there is no systematic way to determine whether a key variable has been omitted. Nonetheless, prior analyses have worked to address the omitted variables issue by saturating the levels model with covariates nordhaus2006geography or assuming confounding factors are orthogonal to climate acemoglu2000colonial.
An alternative widely-used approach is to estimate a panel regression model that exploits random inter-temporal variation in environmental variables and includes location-specific fixed effects, which partial-out unobserved cross-sectional differences that influence the outcome of interest. In the context of climate effects, this approach has been used to estimate the effect of high-frequency weather variation on seasonal crop yields deschenes2007economic, schlenker2009nonlinear, lobell2011climate, auffhammer2014empirical, marginal effects that are then convolved with climate-specific weather distributions to reconstruct estimates for the effect of climate. However, this approach has remained contentious because it is unclear whether marginal effects of weather can be assumed similar to the marginal effects of changing the climate once populations have fully adapted hsiang2016climate. Furthermore, this approach cannot be applied to understand the effect of other environmental factors, such as soil quality, that exhibit little or no inter-temporal variation within a location on practical time scales. To our knowledge, the long-run effects of time-invariant natural factors, such as soil quality, remain essentially unstudied empirically in modern economics because the challenge of unobserved heterogeneity has been insurmountable.\footnote{hornbeck2012nature is an important and notable exception.}
Here, we demonstrate that SFD presents an appealing “third path”, allowing us to estimate the impact of long-run geographic conditions while eliminating the effects of unobserved heterogeneity. Our estimates can be thought of as the effect of long-run climate and soil quality net of adaptation and all long-run adjustments, in the sense originally articulated by mendelsohn1994impact. To evaluate the research design's performance in the presence of omitted variables, we again compare SFD estimates to standard levels estimates when covariates are intentionally and systematically withheld from the regression model.
The remainder of this section is organized as follows. First, we introduce our data and cross-sectional specification. Second, we estimate the effect of climate and soil on maize yields using the levels research design. Third, we demonstrate how SFD can be implemented with irregular units in two-dimensional space and estimate these effects via SFD, recovering new estimates for the effect of long-run climate and soil on yields net of unobserved heterogeneity. Fourth, we systematically withhold covariates from both regression models to compare how the two research designs perform in the presence of omitted variables. Lastly, we conduct two internal robustness checks for the SFD approach: a continuous rotation of the coordinate system and SDD.
We obtained the data on annual county-level maize yield, temperature, and rainfall for the years 1950-2005 used in schlenker2009nonlinear and the vector of five soil characteristics used in schlenker2006impact. The soil characteristic include the minimum permeability (inches per hectare), average water capacity (inches per inch), soil erodibility factor (0.02 for the least erodible soils to 0.64 for the most erodible), percent clay content (%), and percent of high-class top soil (%). Following these authors, we limit our analysis to the balanced panel of counties east of the 100th meridian. To demonstrate the benefits of SFD in cross-sectional data, we average the weather and yield data over the 56-year period, so there is only one observation per county. We then match these datasets with the soil data, which is already a cross section, to create a final cross section. Importantly, our long-term averages of weather (temperature and rainfall) can be viewed as measures of local climate defined in earlier work by mendelsohn1994impact. To the best of our knowledge, effects of cross-sectional variation in soil conditions have not been a focus of prior study, perhaps in part due to concerns of unobserved heterogeneity, although changes in soil quality over time have been analyzed hornbeck2012nature.
We employ a log-linear regression model with maize yield as the dependent variable. There are nine explanatory variables describing seven environmental conditions, since nonlinear effects of temperature and precipitation are each described by two variables.
Using data from schlenker2009nonlinear, we represent the non-linear effect of temperature on maize yields using a formulation that reflects a piecewise-linear spline in hourly temperatures during the growing season. This approach measures the amount of time a crop is exposed to various temperatures at high temporal resolution, but collapsed so that it may be matched to yield data that is collected after longer intervals of time over which crop growth occurs. Our specification allows the effect of hourly temperature to have a different marginal effect depending on whether the temperature is above or below 29\degree C. Degree-days below 29\degree C is a variable whose coefficient captures the effect on end-of-season yields from a marginal 24 hour period at temperatures between 0 and 29\degree C. The coefficient for the degree-days above 29\degree C variable describes the effect of a marginal 24 hour period at temperatures above 29\degree C. schlenker2009nonlinear established that large declines in maize yield occur for temperatures above 29\degree and the effects of hourly heat appear to be additively separable, motivating this specification.\footnote{Denoting temperature in Celsius as $T_{h}$ for each hour $h$ in growing season year $y$, these two variables are constructed:
These variables thus summarize hourly thermal exposure integrated over time across these two thermal ranges in units of degree-days. }
Our explanatory variables also include linear terms in the five soil characteristics described above and linear and quadratic terms in total growing-season precipitation (March to August, measured in mm). We note that the SFD approach is well-suited to capture non-linear effects, in this case those of rainfall and temperature, so long as the terms that describe the nonlinearity are computed prior to differencing.\footnote{To see this, note that we construct the SFD estimator for the model $y_i = \beta_1 p_i + \beta_2 p_i^2$ by writing $\Delta y_i= y_i - y_{i-1} = \beta_1 (p_i - p_{i-1})+ \beta_2 (p_i^2 - p_{i-1}^2)=\beta_1 \Delta p_i + \beta_2 \Delta p_i^2$. The coefficient $\beta_2$ maintains its interpretation even after differencing.}
To evaluate the performance of the levels and SFD models in the presence of omitted variables, we initially employ seven different specifications, altering which covariates are included in the regression. Specification (1) includes only temperature variables, (2) includes only rainfall variables, and (3) includes only soil variables. Specification (4) includes temperature and rainfall variables, (5) includes temperature and soil variables, and (6) includes rainfall and soil variables. Specification (7) includes all variables. For both temperature and precipitation, across all models, the two variables describing each (e.g. precipitation and precipitation-squared) are either included together or omitted together. By intentionally withholding known covariates from specifications (1)-(6), we mimic situations in which some (possibly unknown) variables are omitted from the model. We then evaluate how the levels and SFD estimators perform in these cases compared to specification (7), the most saturated model.
First, we estimate the effect of environmental conditions on maize yields using a standard cross-sectional specification in levels, analogous to the models estimated by mendelsohn1994impact and schlenker2006impact. The estimate the model
where $y_i$ is average maize yield in county $i$, $\alpha_1$ is a constant, $\mathbf{t}_i$ is the vector containing degree-days below 29\degree C and degree-days above 29\degree C, $\mathbf{p}_i$ is the vector containing precipitation and precipitation-squared, $\mathbf{s}_i$ is the vector of soil variables, and $\epsilon_i$ are unexplained variations. Terms are withheld in specifications (1)-(6) and the full model is estimated in specification (7).
The results are displayed in Table (ref). In the levels research design, the parameter estimates generally have the same sign across models but exhibit highly inconsistent point estimates in almost all cases. For example, the estimated effects of moderate temperature ({degree-days below 29\degree C}) and extreme heat ({degree-days above 29\degree C}) change substantially between the specification where only temperature is included (1) and the specification where both temperature and precipitation are included (4). Indeed, the coefficient estimate for days with moderate temperatures is positive in (1) but negative in (4), and the estimate for extreme heat is three times larger in (1) than it is in (4). The estimated effects for the soil variables also change considerably across specifications. For instance, across the five soil characteristics, the average difference between the estimates in the specification with only soil controls (3) and with all variables (7) is 50% of the point estimate in specification (3), with a high of 97% (for percent clay) and a low of 17% (for minimum permeability). Finally, we find it is worrisome that the estimated coefficient for the soil erodibility factor is positive and significant in all models, at odds with the agronomic literature that finds higher soil erodibility generally decreases yields renard1997predicting.
These results highlight the vulnerability of the standard levels model to omitted variables bias. Indeed, withholding covariates generally leads to large changes in the magnitude of estimated effects. In the case of degree-days below 29\degree C the sign of estimates are inconsistent across models, and in the case of the soil erodibility the estimate has the incorrect sign and is “statistically significant" across all models. Even when we include all three sets of controls in specification (7), we cannot be certain that important variables are not still missing, and given the inconsistency of parameter estimates across observed specifications, it is not unreasonable to expect that the estimates may once again change substantially if we included additional covariates in the regression model oster2014unobservable.
Next, we estimate the effect of environmental conditions on maize yields using the SFD research design. To do this, we first must overcome a key challenge to implementing SFD with county-level data: administrative boundaries do not follow the regular lattice structure depicted in Figure (ref)b. There exist many adjacencies between counties that may be exploited in an SFD research design, but ensuring that each county is differenced from exactly one neighboring county in a sequence is no longer trivial. Equation ((ref)) does not define a specific way in which to define neighbors and many arrangements would be valid. Importantly, sampling of differences must be organized such that no observation is double counted. Here, we develop a generalizable approach that imposes a “panel-like” structure on irregularly shaped counties, thereby identifying sequences of adjacent counties without double counting them. Notably, however, there exist other valid algorithms for setting up SFD in two dimensional space that we do not explore here.
The basic procedure is depicted in Figure (ref). First, we overlay the US map with 50 sampling “channels." These are long and narrow regions, approximately 30 miles wide, spanning West-East slices of the country and defined by a northern and southern boundary.\footnote{The width of our sampling channels was chosen to match the average width (from North to South) of the counties in our sample.} Adjacent channels share a boundary. Beginning with the northernmost channel, all the counties that intersect with the channel are recorded as sequentially adjacent and ordered by the longitude of their centroid. Then we move south to the next channel and repeat this process. To assure that each county is only included in one channel, if a county has already been sampled in a preceding more northern channel, it is omitted from the remaining southern channels. The sequence of counties within each channel are thus treated like a sequence of observations within a “panel-like unit” of the regularly shaped observations shown in Figure (ref)b. Finally, differences are computed between the ordered adjacent counties and a cross-sectional regression is estimated in these first-differences. Notably, SFD can be implemented in any direction by rearranging how the channels are oriented when they are first overlaid. We begin our analysis by computing SFD in both the West-East direction and in the North-South direction for comparison.
Figure (ref) illustrates a single sequence of adjacent counties within one channel (panel a) derived using this procedure and compares two variables of interest in levels and SFD. Comparing levels (panel b) to SFD (panel c), one can see how the use of spatial first differences eliminates low-frequency correlations contained in the “spatial history" of the two variables. What remains is the variation we use to estimate $\hat{\beta}_{SFD}$. Note that discontinuities in the levels of these variables, such as those that occur at the borders of Missouri, do not systematically affect the data after differencing.
Using the differences computed between ordered adjacent neighbors, we estimate the effect of climate and soil on maize yields via the SFD model
The model now contains differences of terms in non-linear functions for temperature and precipitation, but the coefficients for these differenced terms maintain the same interpretation as in the levels model.
The results for the same set of specifications (1)-(7) are displayed in Table (ref), although each specification is now estimated twice, once in the West-East direction and once in the North-South direction. Using SFD, the coefficient estimates are extremely consistent across models, especially in comparison to the levels estimates in Table (ref). There are no sign reversals and the median difference between the coefficient estimate in the model with just one set of variables and the model with all three sets of variables is 12% of the point estimate in the model with one set, as opposed to 38% for the levels estimates. Additionally, the SFD estimates are essentially unchanged when calculated in the West-East and North-South directions. The median difference between the West-East estimate and the North-South estimate is 10% of the point estimate, and this difference is less than 30% for all but one variable (percent clay).
The SFD estimates are also all consistent with the agronomic literature. Days with moderate temperatures are estimated to significantly increase yields across all specifications, in contrast to the levels model which recovered this result in one of four specifications. As expected, days with extreme heat reduce yields across all specifications. Precipitation changes have a smaller impact, but continue to have an inverted U-shaped effect on yields. A higher average water capacity increases yields, a higher percentage of clay reduces yields, a lower minimum permeability (which indicates drainage problems) reduces yields, a higher soil erodibility factor is harmful (the levels model consistently indicated the opposite), and better soils (as measured by best soil class) are beneficial. The SFD research design was able to recover these results in all specifications, in both the West-East and North-South directions. The relative invariance of all coefficient estimates across all models provides us with modest confidence that our results are robust to other important covariates that may not be observed even in the most saturated specification.
Our estimates for the constant term, $\hat{\alpha}_2$ are generally positive and significant in the West-East model and have a negative sign in the North-South model. In the SFD model, the constant term is interpretable as average trend in space as one moves in the direction of differencing, conditional on changes in the covariates. Thus a positive constant term in the West-East model indicates that yields are increasing as one moves from West to East (since we have subtracted from each observation values from its neighbor to the West). Similarly, a negative constant term in the North-South model indicates that yields are decreasing as one moves from North to South.
One practical question that arises with empirical estimation is how to calculate standard errors for SFD estimates. In all cases, we expect residuals estimated via OLS to be negatively auto-correlated (at least first-order) due to the first-differencing procedure since sequential residuals $\Delta \epsilon_i = \epsilon_i - \epsilon_{i-1}$ and $\Delta \epsilon_{i+1} = \epsilon_{i+1} - \epsilon_{i}$ share the component $\epsilon_i$ and it enters positively in one instance and negatively in the other. In this context, we also reject the null hypothesis of homoskedastic disturbances using the Bruesh-Pagan test. To address these two issues, we calculate five different sets of standard errors: (i) Conley standard errors, (ii) Newey-West standard errors, (iii) standard errors clustered by channel (iv) bootstrapped standard errors,\footnote{For the bootstrapped standard errors, we resample at the observation-level after differencing between adjacent counties.} and (v) block-bootstrapped standard errors block-resampled by sampling channel. These standard errors are presented in Appendix (ref), along with the OLS standard errors for comparison. The magnitude of these various standard error estimates are comparable across all five methods, but the Conley, Newey-West, and block bootstrapped standard errors are generally slightly larger, suggesting some additional spatial autocorrelation in $\Delta \epsilon$ beyond immediate neighbors. Thus, in Table (ref), we report standard errors that account for spatial autocorrelation using Conley's approach.
The recent literature in this field has been particularly interested in the effect of climate on yields, motivated by efforts to understand the economic consequences of climate change mendelsohn1994impact, deschenes2007economic, schlenker2009nonlinear, burke2016adaptation, hsiang2016climate. Because the temperature-yield and precipitation-yield relationships are nonlinear and the coefficient estimates are difficult to interpret, we plot these relationships in Figure (ref). Results from the levels model are shown in orange and those from the SFD model are shown in blue. We include the estimates from both the specification with no controls ((1) for temperature, (2) for precipitation) and the model with a full set of controls (7). Panel a shows our estimated effects for temperature. Two features of the results stand out. First, the two SFD estimates are very near one another, despite the differences being computed in orthogonal directions and thus exploiting different variation in the independent variables. Second, the SFD estimates are remarkably similar to previous estimates of the effect of long-run trends in temperature on US maize yields. Our coefficient estimate for the variable degree-days above 29\degree C is $-0.0050$ when SFD are computed in the West-East direction and $-0.0048$ when SFD are computed in the North-South direction. Using long differences over the period 1980-2000, burke2016adaptation estimated this same coefficient to be $-0.0053$ in their specification with time-specific fixed effects and $-0.0044$ in their specification with state fixed effects.\footnote{We use burke2016adaptation as a benchmark against which to compare our results because it is the only study to estimate the effect of long-run climate on agricultural productivity that is plausibly robust to unobserved heterogeneity. The authors employ a “long differences" approach and model county-level changes in yields over time as a function of changes in temperature and precipitation, accounting for time-invariant unobservables at the county level and time-trending unobservables at the state level.} Holding all else equal, these findings imply\footnote{$-0.00499$ per degree day $\times10$ degree-days $= -0.05$, i.e. 5 log points.} that substituting one full day (24 hours) at 29\degree C temperature with a full day at 40\degree C results in a predicted end-of-season yield decline of approximately 5%. Panel b of Figure (ref) shows the estimated effect of precipitation on maize yields. Once again, the SFD estimates are near one another across different regression models and differencing directions. Also noteworthy is the result that the four different SFD estimates largely agree with one another while the two levels estimates differ substantively, both from one another and from the range of SFD estimates.
With the goal of systematically evaluating both research designs' vulnerability to omitted variables bias, we experiment further by withholding all possible combinations of covariates when estimating effects for each explanatory variable. For each environmental variable of interest (e.g. temperature), this procedure generates a total of 192 different specifications where the remaining six controls are systematically withheld in all possible combinations (e.g. precipitation, precipitation $+$ water capacity, precipitation $+$ water capacity $+$ percent clay, $\dots$).\footnote{The temperature variables (degree-days below 29\degree C and degree-days above 29\degree C) always appear together. Similarly, the precipitation-squared variable is only included when precipitation is included.} This type of analysis is similar to the “extreme bounds" analysis proposed by leamer1985sensitivity, and taken to its logical extreme by sala1997just, who ran nearly two million growth regressions using different combinations of 62 explanatory variables. The goal of the procedure, as we are using it, is to gain more general insight into the magnitude and distribution of the omitted variables bias that is eliminated from the cross-sectional levels regression by differencing, as described in Eq. ((ref)). Specifically, for each variable of interest, we compute all possible estimates for $\hat{\beta}_L$ and $\hat{\beta}_{SFD}$ and compare their relative stability across these specifications. Variations in $\hat{\beta}$ that occur when covariates change are interpreted as evidence of omitted variables bias, although it is unknown which specification is “correct."
Using this extreme bounds analysis, SFD dramatically outperforms estimation in levels. The distribution of estimated marginal effects across the 192 regression specifications for each variable is shown in Figure (ref). The variance of this distribution (averaged across variables) is 80% smaller when employing SFD as opposed to a cross section on levels, with a high of 100% (for best soil class) and a low of 47% (for water capacity). Furthermore, under the SFD model, the coefficients for days with moderate heat and the soil erodibility factor have the expected positive and negative signs, respectively, across all 192 specifications. This pattern does not hold for the levels model, where the coefficient for days with moderate heat is often negative and the coefficient for the soil erodibility factor always has the wrong sign (positive). Reinforcing our findings from above, these results suggest that including/omitting variables often leads to substantial changes in the levels estimates but has limited effect on SFD estimates in this context.\\
This example demonstrates how the SFD research design can be robust to unobserved heterogeneity. In the context of maize yields, there appears to be a large degree of low-frequency spatial correlation between the covariates, leading to large biases in $\hat{\beta}_L$ when control variables are withheld from the model. In contrast, SFD recovers estimates for all seven environmental factors that are essentially unchanged regardless of whether or not key variables are included in the model. We suspect that if an important covariate were still missing from the saturated regression model (7) and it was discovered and included in a new specification, it would not dramatically change the SFD estimates.
As a final exercise, we demonstrate two internal robustness checks made uniquely possible in the SFD research design: a rotation of the coordinate system and spatial double differences (SDD). As discussed above, these tests should fail if the orthogonality condition in Eq. ((ref)) is not true.
First, we conduct a sensitivity analysis that exploits the rotation of the coordinate system. After imposing the sampling grid on the US map with the channels arranged in the West-East direction, we rotate the map by an angle $\theta$ at 1\degree increments on its axis under the grid for $\theta = -89\degree$ to $90\degree$ (see Figure (ref)a).\footnote{Equivalently, one could instead choose to keep the position of the map fixed, and rotate the sampling channels by the angle $-\theta$.} We then estimate the nonlinear effect of temperature on maize yields using SFD with no controls at each $\theta$, producing 180 different estimates of specification (1). These 180 estimated marginal effects of extreme heat and moderate temperatures are shown in Figures (ref)b and (ref)c, respectively. Note that the estimated effects are highly consistent across the different coordinate rotations. Indeed, the variance of the coefficient estimate for degree-days above 29\degree C is $1.025$ for a coefficient estimate of $-5$, producing a coefficient of variation equal to 0.2 (Figure (ref)b). Similarly, degree-days below 29\degree C have a coefficient of variation equal to 0.19 (Figure (ref)c). The fact that the estimates are consistent across all sampling directions implies either the identifying assumption holds (i.e. Eq. (ref) is true) or that it fails but somehow generates a similar bias for each of the 180 different estimates, despite these estimates exploiting different sources of variation in the independent variables.
Second, we repeat the analysis using SDD, as described in Eq. ((ref)). The SDD estimates are displayed in Table (ref). These estimates are near the SFD estimates (with the exception of those for percent clay), with a median difference of 18% of the SFD coefficient estimate. While the percent difference between the SFD and SDD estimates for percent clay are large, the SDD estimates for this variable still lie within the 95% confidence interval of the SFD estimates. In this context, it is not surprising that the SFD estimates and SDD estimates differ somewhat since it is difficult to estimate $\hat{\beta}_{SDD}$ precisely with this modestly sized sample. $\hat{\beta}_{SDD}$ will almost certainly be a substantially more variable estimate than $\hat{\beta}_{SFD}$ in almost all environments as variation in $\Delta^2 \mathbf{x}$ is much smaller than variation in $\Delta \mathbf{x}$. In practice, SDD is also vulnerable to attenuation bias since a relatively larger fraction of variation in $\Delta^2 \mathbf{x}$ may be due to measurement error. Nonetheless, similar to the rotation test above, the stability of estimates across SFD and SDD suggests that in order for $\hat{\beta}_{SFD}$ to be biased by the failure of the identifying orthogonality condition (Eq. (ref)), the structure of the omitted variables must be such that $\hat{\beta}_{SDD}$ is similarly biased by the failure of an entirely different orthogonality condition (Eq. (ref)) for each variable.
We interpret the results of these robustness checks as strong evidence that the assumptions underlying the SFD research design are very likely to be valid in the context of these new empirical estimates.
In standard cross-sectional approaches to inference, it is well understood that the omission of unobservable covariates may lead to large biases in estimated effects. Due to this fact, cross-sectional approaches are often not considered reliable research designs for obtaining causal estimates in many disciplines when instrumental variables are unavailable. We propose SFD as a simple, general, and robust alternative when observations are organized and densely packed in space.
We highlight that the Local Conditional Independence assumption underlying SFD is conceptually similar to the assumptions exploited in several well-established research designs. These include the assumption that immediately sequential observations within a time series are comparable in event study designs, the assumption that sequential observations within a panel unit are comparable in differences-in-differences panel analyses, and the assumption that observations just above and just below a treatment discontinuity are comparable in regression-discontinuity designs. Indeed, the assumptions necessary for the SFD approach to be valid are so nearly identical to these other assumptions that it seems difficult to logically reject one without also rejecting the other. Importantly, however, SFD is not suitable for all contexts and requires judgement from the analyst about whether observational units are “dense enough" in space, a data constraint that informs the potential validity of the Local Conditional Independence assumption in practice.
We imagine that the application of SFD could be applied in a number of different geometries. We demonstrate the application of SFD in one-dimensional space, in two-dimensional gridded data, and in US counties by imposing a “panel-like" structure on the data. However, with irregular (non-gridded) data, other approaches could be taken. For example, in the two-dimensional space, one could difference in a spiral structure to generate a single sequence of adjacent observations (rather than the “channels" approach we explore here). One could also combine differences taken in both the West-East and North-South directions---which exploit different variation in the variables---to increase the amount of variation in the sample. The SFD design could also be applied in other contexts. For instance, one might implement SFD along a coastline, throughout an infrastructure network, or even vertically up and down the floors of a skyscraper. Indeed, it remains an open question how to optimize the research design in different geometries and to what extent the performance that we document here generalizes.
It is our hope that the SFD research design reopens closed doors in the analysis of cross-sectional data. In many fields of economics---such as environment, development, geography, health, industrial organization, labor, public, growth, trade, and urban---there are core questions that are fundamentally cross-sectional in nature. Historically, econometricians that seek to address these questions have generally had two options: to trust that unobserved heterogeneity does not confound the analysis or to employ cross-sectional instruments which rely on exclusion restrictions that cannot be tested. The SFD research design may offer yet another path.
\singlespacing
\singlespacing