# Spatial Differencing for Sample Selection Models with Unobserved
Heterogeneity

Alexander Klein Guy Tchuente
University of Kent School of Economics, e-mail<EMAIL_ADDRESS>Corresponding author: School of Economics, University of Kent, e-mail:
<EMAIL_ADDRESS>Address: Kennedy Building, Park Wood Road, Canterbury,
Kent, CT2 7FS, Tel:+441227827249.

( June 2020)

###### Abstract

This paper derives identification, estimation and inference results using
spatial differencing in sample selection models with unobserved heterogeneity.
We show that under the assumption of smooth changes across space of the
unobserved sub-location specific heterogeneities and inverse Mills ratio, key
parameters of a sample selection model are identified. The smoothness of the
sub-location specific heterogeneities implies a correlation in the outcomes.
We assume that the correlation is restricted within a location or cluster and
derive asymptotic results showing that as the number of independent clusters
increases, the estimators are consistent and asymptotically normal. We also
propose a formula for standard error estimation. A Monte-Carlo experiment
illustrates the small sample properties of our estimator. The application of
our procedure to estimate the determinants of the municipality tax rate in
Finland shows the importance of accounting for unobserved heterogeneity.

Keywords: Sample selection, Spatial difference, Unobserved heterogeneity.

##  1 Introduction

In linear models, spatial differencing has been used to deal with unobserved
omitted variables. The availability of geographical locations allowed
empirical papers to take advantage of the spatial dimension of the data and
control for various unobserved heterogeneity (e.g. Duranton, Gobillon, and
Overman (2011), Black (1999) or Holmes (1998)). In general, spatial
differencing offers an identification strategy in the situations when
researchers face cross-sectional data with unobserved heterogeneity _and_ lack
suitable instrumental variables. This paper extends spatial differencing to a
model with sample selection.

For economists, the question of omitted variables is a serious concern in the
context of nonexperimental data. The solution is straightforward when omitted
variables are simply the result of not including all relevant variables for
which data exist - we add such variables to the model to avoid the bias
induced by their omission. When omitted variables are unobserved, researchers
have basically three options: they can use (i) proxies, (ii) instrumental
variables, or (iii) differencing the data across time or space.

The proxies reduce the bias if they manage to capture the effect of the
omitted variables such that what remains is uncorrelated with the error term.
However, it is often the case that the proxies are imperfect, hence they may
still be related to the unobserved heterogeneity, or the error term if it
turns out to be endogenous, or it could also be irrelevant after controlling
for observed covariates. In such cases, the inclusion of the proxy will not
solve the bias problem, it may even exacerbate it.111 See Todd and Wolpin
(2003) for details on the use of proxy and Oster (2019) for a rigorous
treatment of the evaluation of the robustness to omitted variables. The second
solution - using a valid set of the instruments - may help alleviate the bias.
However, as discussed in Todd and Wolpin (2003), the “quasi-experimental"
local average treatment effect (LATE) obtained in the instrumental variable
model may not correspond to the _ceteris paribus_ and thus, may not correspond
to the deep structural parameter of interest. Lastly, panel data sets allow
researchers to control for the unobserved heterogeneity. They indeed help them
to identify the causal effect when, for example, time constant unobserved
heterogeneity might cause endogeneity problem and strong instruments
satisfying exclusion restrictions cannot be found. However, there might be
situations when such data sets are not available.

Our paper is a contribution to the literature identifying and estimating model
parameters in the presence of unobserved omitted variables. We propose an
identification strategy based on spatial differencing. As was discussed above,
this approach has been used in the context of linear regressions. However,
little is known about its performance in non-linear models. We extend spatial
differencing into this direction, specifically to the case of cross-section
data with sample-selection. We show that under justifiable assumptions on the
smoothness of the unobserved heterogeneity (i.e. spatially close individuals
have similar unobserved heterogeneity and the derivative of their Inverse
Mill’s ratio are similar), spatial differencing eliminates the unobserved
effects even in the presence of a nonlinear element - in our case Mill’s
ratio. The parameters of interest of our sample selection model are estimated
using a standard two-step approach of Heckman (1974) Heckman (1979). We derive
asymptotic properties and propose a correction of standard errors accounting
for the two-step nature of our estimation and spatial differencing. The
asymptotic behavior of the estimator reveals important properties of spatial
differencing that researchers would need to be cautious about. The new
estimator and the standard errors correction are easy to implement.

The intuitions for the model of sample selection with unobserved spatial
heterogeneity that we consider in this paper can be described as follows.
Suppose we have a cross-sectional data on municipalities which are organized
into larger geographical units called regions, and which have the authority to
set the levels of local taxation. Municipality tax rates must be at least as
high as the threshold set by a central government. As a result, municipalities
self-select into those with the tax rates at the threshold and those above it.
We are interested in what determines the municipalities’ tax rates. The tax
rate will depend on various socio-economic characteristics, such as age
composition of population and income, but also on amenities. These can depend
on the region where the municipalities are located: for example, regions with
natural landscapes might have different level and composition of amenities
than regions without them. We can control for them with region-specific
dummies. However there can be a considerable unobserved heterogeneity at
municipality level. Controlling for that with municipality-specific dummies
might not be an option, since we may quickly run out of degrees of freedom.
Therefore, we face a problem of a self-selected cross-sectional sample with
unobserved heterogeneity which we cannot fully control with dummies, and which
has two spatial dimensions: high-level which we call locations (in our example
regions), and low-level which we call sub-locations (in our example
municipalities).

Spatial differencing will eliminate the sub-location specific unobserved
heterogeneity. It will, at the same time, also induce a correlation in the
error terms. We take that correlation into account, and derive the asymptotic
properties of our estimator using similar arguments to those present in the
derivation of the asymptotic behavior of the clustered standard errors: the
number of locations goes to infinity and the size of location is assumed
random and bounded almost surely. We find that this result also extends to a
linear model without sample selection. This has important implications that
researchers need to be cautious about. Indeed, the consistency of estimator
applied to the spatially differenced data requires _(i)_ a large number of
locations, _(ii)_ a limited number of individuals in each location. Monte
Carlo simulations also suggest that it would be better if the number of
individuals in sub-locations were small as well. Before we continue, let’s
notice that locations in our model are equivalent to clusters and we use
’location’ and ’cluster’ interchangeably.

Since our estimator is derived for a clustered sample with unobserved
heterogeneity, this paper contributes to the literature on the selection
correction in panel data. In this literature, the main challenge is the
presence of individual-specific unobserved heterogeneity in both the outcome
and the selection equations. The existing solutions are based on either a full
model specification or on a differencing procedure. Wooldridge (1995) uses a
Mundlak approach to specify the individual-specific unobserved heterogeneity
in both equations. He also imposes a special functional form to the selection
mechanism. Kyriazidou (1997), on the other hand, does not impose strong
restriction of the selection equation functional form and uses a nonparametric
approach to difference out the unobserved fixed-effect. Rochina-Barrachina
(1999) similarly relies on differencing to identify the parameters of the
model, but she also imposes additional distributional assumptions to the
selection equation. Even if our problem has similarities with the selection
correction in panel data model literature, the main difference is that we are
observing a clustered cross-section. In each cluster, there is a finer common
sub-location specific unobserved heterogeneity shared by some individuals in
that cluster. This heterogeneity, however, is different from the cluster and
individual-specific ones studied in panel data models and implies a different
cluster asymptotic. Since the outcome of the individuals are not independent
in our model, while it is in the panel data case, our asymptotic results are
thus derived using a large number of cluster asymptotic with heterogenous,
random and bounded cluster size.

The clustered dependence created by the finer sub-location specific unobserved
heterogeneity relates our asymptotic discussion to the papers dealing with
clustering at variance level (see Wooldridge (2010) for a textbook treatment).
The asymptotic in that literature is derived using either a large or a fixed
number of clusters. The fixed number of cluster leads to non-normal asymptotic
and discussion about recent contributions can be found in Hansen and Lee
(2019). A large number of cluster asymptotic was first derived by White (1984)
and has been investigated by several authors allowing either fixed cluster
size or heterogeneous cluster. Recent developments include Hansen and Lee
(2019) who propose conditions on the relation between the cluster sample sizes
and the full sample in a regular asymptotic, or Djogbenou, MacKinnon, and
Nielsen (2019) who derive asymptotic with varying cluster sizes and carry out
a cluster wild bootstrap. Our results complement this literature by extending
the cluster asymptotic to a sample selection model.

We present an empirical application of our new estimator. We examine the
determinants of tax rates across four hundred and eleven Finish municipalities
spread across nineteen Finish regions. In 1999, Finish central government has
decided to raise the lower bound of the tax rate municipalities could set from
0.2% to 0.5%. This has created a sample selection mechanism which resulted
with more than half of the municipalities opting for 0.5% while the rest
charging higher tax rate. We use our spatial differencing estimator to control
for unobserved municipalities effect which can be correlated with the error
term, creating thus an endogeneity problem and rendering the standard sample
selection estimator biased and inconsistent. Our results clearly show the
presence of unobserved heterogeneity across municipalities, limitations of
using only region dummies (nineteen in our case) to fully control for
municipalities’ unobserved heterogeneity, and the importance of spatial
differencing to control for it.

The structure of the paper is as follows. First, we expand the spatial
differing method in the case of linear regression model to the case of sample
selection. Then we discuss identification assumptions, propose an estimation
procedure, and derive the estimator of the corrected standard errors. Lastly,
we conduct Monte Carlo simulations and present an empirical application of our
estimator.

##  2 Sample Selection Models with Spatial Correlation

In many economic applications, we are interested in estimating the following
regression equation:

$y_{ij}=x_{ij}^{\prime}\delta+\gamma_{j}+\gamma_{j\alpha}+\varepsilon_{ij}$
(1)

where $\displaystyle x_{ij}^{\prime}$ is a vector of exogenous controls
variables, $\displaystyle\gamma_{j}$ is location fixed effect,
$\displaystyle\gamma_{j\alpha}$ is a sub-location specific effect for sub-
location $\displaystyle\alpha$ which is at a finer spatial scale than location
$\displaystyle j$, and $\displaystyle\varepsilon_{ij}$ is the error
term.222The sub-location specific component, $\displaystyle\gamma_{j\alpha}$,
is a simplification for $\displaystyle\gamma_{j\alpha_{i}}$. We are implicitly
assuming that the sub-location specific effects are the same for all its
individuals.

Examples of application of this model can be found in the estimation of the
fertilizer effect on wheat crops yield in farms growing multiple crops, or the
effect of local taxation on the growth of firms. The crop-yields depend on the
soil quality of location $\displaystyle\gamma_{j}$ (e.g. a village), but also
on the sub-location specific soil composition (e.g. a farm in the village)
Collins, Alva, Boydston, Cochran, Hamm, McGuire, and Riga (2006). Similarly,
the impact of local taxation on the growth of firms may vary by county but
also sub-locations such as neighborhoods as in Duranton, Gobillon, and Overman
(2011). We can control for $\displaystyle\gamma_{j}$ with location dummy
variables. However, they might not be enough to capture all unobserved
heterogeneity related to location $\displaystyle j$ as there can be
considerable heterogeneity at finer spatial scale of sub-locations: using the
example above, the firms are located in various neighborhoods
$\displaystyle\alpha$ which are sub-locations of location $\displaystyle j$.
Furthermore, standard location fixed effect $\displaystyle\gamma_{j}$ relies
upon an arbitrary specification of the comparison neighborhood group, as
pointed out by Gibbons and Machin (2003), making it an imperfect control for
sub-location specific effect $\displaystyle\gamma_{j\alpha}$. If
$\displaystyle\gamma_{j\alpha}$ is correlated with $\displaystyle x_{ij}$, OLS
estimate of $\displaystyle\delta$ will be biased. In the absence of suitable
instrumental variables for $\displaystyle x_{ij}$, the spatial differencing
offers a solution by differencing out the unobserved sub-location specific
effects $\displaystyle\gamma_{j\alpha}$.

Duranton, Gobillon, and Overman (2011), Black (1999) or Holmes (1998) use
spatial differencing in the case of linear models to solve endogeneity
problems arising from unobserved sub-location effect
$\displaystyle\gamma_{j\alpha}$. They take advantage of the fact that for
sufficiently small distances between sub-locations, their specific effect
$\displaystyle\gamma_{j\alpha}$ changes smoothly across space, allowing thus
to difference them out. This corresponds to the following assumption.

Assumption I1: The sub-location specific unobservable effect is homogenous in
the a neighborhood of the individual ie
$\displaystyle\Delta_{d}\gamma_{j\alpha}=0$ for $\displaystyle d$ small
enough.

In several economic models, in addition to the sub-location specific fixed
effects $\displaystyle\gamma_{j\alpha}$, the outcome of interest is not
observed for the selected sub-sample. The selection can be the result of the
decision of individuals or the researcher. The presence of sample selection
introduces nonlinearity to the model (1).

We specify the model with sample selection as follows. Consider two latent
dependent variables $\displaystyle y^{*}_{1ij}$, and $\displaystyle
y^{*}_{2ij}$ in a cross-section which follow a regular linear model for
individual $\displaystyle i$ in a location $\displaystyle j$:

$\displaystyle
y_{1ij}^{\ast}=z_{ij}^{\prime}\beta+\theta_{j\alpha}+\theta_{j}+\varepsilon_{1ij}$
\- selection equation,

$\displaystyle
y_{2ij}^{\ast}=x_{ij}^{\prime}\delta+\gamma_{j\alpha}+\gamma_{j}+\varepsilon_{2ij}$
\- outcome equation.

Individual error terms are $\displaystyle\varepsilon_{1ij}$ and
$\displaystyle\varepsilon_{2ij}$ ; $\displaystyle\theta_{j\alpha}$ and
$\displaystyle\gamma_{j\alpha}$ are sub-location specific effects for a sub-
location $\displaystyle\alpha$ in location $\displaystyle j$, affecting the
selection and the outcome equation respectively. The exogenous characteristics
$\displaystyle x_{ij}$ affect the outcome. They could be correlated with
$\displaystyle\gamma_{j\alpha}+\gamma_{j}$ but not with
$\displaystyle\varepsilon_{1ij}$ and $\displaystyle\varepsilon_{2ij}$. The
variables $\displaystyle z_{ij}$ are exogenous variables determining
selection, they can be a subset of $\displaystyle x_{ij}$. However, for
identification purposes, some elements of $\displaystyle z_{ij}$ are assumed
to be absent from $\displaystyle x_{ij}$.

Assumption I2: $\displaystyle\varepsilon_{1ij}$ and
$\displaystyle\varepsilon_{2ij}$ are independent identically distributed
normal random variables for all $\displaystyle i,j$.
The outcome is modelled in the form of a truncated sample selection model and
is represented by equation (2).

$y_{2ij}=\left\\{\begin{array}[]{c}y_{2ij}^{\ast}\hskip 5.69046ptif\hskip
5.69046pty_{1ij}^{\ast}>0\\\ -\hskip 5.69046ptif\hskip
5.69046pty_{1ij}^{\ast}\leq 0\end{array}\right.$ (2)

Let us consider the following conditions.

  1. Condition 1: $\displaystyle Cov[z_{ij},\theta_{j\alpha}+\theta_{j}+\varepsilon_{1ij}]=0$; $\displaystyle z_{ij}$ is exogenous

  2. Condition 2: $\displaystyle Cov[x_{ij},\gamma_{j\alpha}+\gamma_{j}+\varepsilon_{2ij}]=0$; $\displaystyle x_{ij}$ is exogenous

  3. Condition 3: and errors $\displaystyle(\varepsilon_{1ij},\varepsilon_{2ij})$ satisfy $\displaystyle\varepsilon_{2ij}=\rho\times\varepsilon_{1ij}+v_{ij}$ with $\displaystyle\varepsilon_{1ij}\sim\mathcal{N}(0,1)$ and independent of $\displaystyle v_{ij}$.

It is possible to consistently estimate $\displaystyle\delta$ by Tobit
regression under these three conditions.333Identification required an
exclusion restriction ie a variable that affects $\displaystyle y^{*}_{1ij}$
but not $\displaystyle y^{*}_{2ij}$. Otherwise, identification relies on the
nonlinearity of the inverse Mills ratio. In most applications, the Condition 1
and 2 are unlikely to hold because there is a possibility that, within a
location, there could be a sub-location specific omitted variable affecting
both the outcome and some observed characteristics of interest. Thus, it is
possible that $\displaystyle Cov[z_{ij},\theta_{j\alpha}+\theta_{j}]\neq 0$
and $\displaystyle Cov[x_{ij},\gamma_{j\alpha}+\gamma_{j}]\neq 0$. The
standard way to deal with the correlation between $\displaystyle x_{ij}$ and
$\displaystyle\gamma_{j\alpha}$ would be to find a suitable instrument for the
$\displaystyle x_{ij}$ and run a IV Tobit or IV two-stage Heckit.

The very local nature of the sub-location specific effect means that it is not
always evident to find a variable correlated with $\displaystyle x_{ij}$ and
uncorrelated with $\displaystyle\gamma_{j\alpha}$. The exclusion restriction
is likely to be violated and IV two-stage Heckit will yield inconsistent
estimates for $\displaystyle\delta$. Another option is to use the finer
location fixed effect and estimate the model using classic Heckman two-stage
procedure, but this will in practice lead to a proliferation of variable and
lose of degrees of freedom.

###  2.1 Identification via Spatial Differencing

This section investigates the application of this spatial differencing
technique to the case of cross-section sample selection models. We denote
$\displaystyle\Delta_{d}$ to be a spatial difference operator. One example is
a pair-wise difference operator which takes the difference between each
observation and another observation located at distance less than
$\displaystyle d$ from that observation. In a location $\displaystyle j$, with
individual $\displaystyle i$ and $\displaystyle k$ who are neighbours. The
pair-wise differencing of the variable $\displaystyle A$ is:

$\displaystyle\Delta_{d}A=A_{ij}-A_{kj}.$

Another example is the difference between the individual outcome and the
average outcome of his/her neighbourhood $\displaystyle\mathcal{N}_{id}$. This
operator is similar to the neighbourhood fixed effect operator, the difference
being that the neighbourhoods can overlap. We call this operator the fixed-
effect difference operator. Let $\displaystyle\mathcal{N}_{id}=\\{k,\hskip
5.69046ptin\hskip 5.69046ptneighbourhood\hskip 5.69046ptd\\}$, the sample size
of $\displaystyle\mathcal{N}_{id}$ is $\displaystyle N_{d}$, the differencing
is given by:

$\displaystyle\Delta_{df}A=A_{ij}-\frac{1}{N_{d}}\sum_{k\in\mathcal{N}_{id}}A_{kj}.$

A further possibility is to use a kernel as in Kyriazidou (1997) to weight
neighbour in $\displaystyle\mathcal{N}_{id}$ according to how far they are, in
term of observable characteristics. This operator is the kernel difference
operator.

$\displaystyle\Delta_{dK}A=A_{ij}-\sum_{k\in\mathcal{N}_{id}}\psi(i,k)A_{kj}.$

Where
$\displaystyle\psi(i,k)=\frac{1}{h_{N_{d}}}K\left(\frac{(z_{ij}^{\prime}-z_{kj}^{\prime})\beta+(x_{ij}^{\prime}-x_{kj}^{\prime})\delta}{h_{N_{d}}}\right)$,
$\displaystyle K$ is a kernel density function while $\displaystyle
h_{N_{d}}$, is a sequence of bandwidths. To illustrate our identification
strategy and for the asymptotic derivation, we use the pairwise spatial
difference operator, while for the empirical application and for the Monte
Carlo simulations, the fixed effect difference is used.

For the spatial difference operator $\displaystyle\Delta_{d}$,
$\displaystyle\Delta_{d}y_{2ij}=y_{2ij}-y_{2kj}$ with $\displaystyle k$ an
observation in the neighborhood $\displaystyle d$ of $\displaystyle i$. Let
$\displaystyle\xi_{ij}\equiv\\{x_{ij},z_{ij},y_{1ij}^{\ast}>0,\gamma_{id},\theta_{id}\\}$
with $\displaystyle\gamma_{id}=\\{\gamma_{kj}\hskip 5.69046ptwith\hskip
5.69046ptk\in\mathcal{N}_{id}\cup\\{i\\}\\}$ and
$\displaystyle\theta_{id}=\\{\theta_{kj}\hskip 5.69046ptwith\hskip
5.69046ptk\in\mathcal{N}_{id}\cup\\{i\\}\\}.$

$\displaystyle\displaystyle E[\Delta_{d}y_{2ij}|\xi_{ij},\xi_{kj}]$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle
E[y_{2ij}-y_{2kj}|\xi_{ij},\xi_{kj}]$ (3) $\displaystyle\displaystyle=$
$\displaystyle\displaystyle E[y_{2ij}|\xi_{ij}]-E[y_{2kj}|\xi_{kj}]$ (4)
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle
x_{ij}^{\prime}\delta+\gamma_{aj}+\gamma_{j}+\rho\lambda(z_{ij}^{\prime}\beta+\theta_{ja}+\theta_{j})$
$\displaystyle\displaystyle-$
$\displaystyle\displaystyle[x_{kj}^{\prime}\delta+\gamma_{j\alpha}+\gamma_{j}+\rho\lambda(z_{kj}^{\prime}\beta+\theta_{ja}+\theta_{j})]$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\Delta_{d}x_{ij}^{\prime}\delta+\Delta_{d}\gamma_{j\alpha}+\rho\Delta_{d}\lambda(z_{ij}^{\prime}\beta+\theta_{ja}+\theta_{j})$
(6)

where $\displaystyle\lambda(c)=\phi(c)/\Phi(c)$ is the inverse Mill’s ratio
while $\displaystyle\phi(c)$ and $\displaystyle\Phi(c)$ are respectively the
density and distribution function of a normal random variable with mean zero
and variance 1.

To go from Equation (3) to Equation (4) we use the linearity of expectation
and the mean independence of $\displaystyle y_{2ij}$ and $\displaystyle
y_{1kj}^{\ast}$ conditional on $\displaystyle\xi_{ij}$, as well as the mean
independence of $\displaystyle y_{2kj}$ and $\displaystyle y_{1ij}^{\ast}$
conditional on $\displaystyle\xi_{kj}$, since we have assumed in Assumption I2
that $\displaystyle\varepsilon_{1ij}$ and $\displaystyle\varepsilon_{2ij}$ are
$\displaystyle iid$. The separation of the conditional set,
$\displaystyle\xi_{ij}$ and $\displaystyle\xi_{kj}$, is possible because we
are working with cross-sectional data. Such separation of the conditional set
is not possible for panel data. Indeed, in the context of panel data with
individual effects and sample selection, when the differencing is used to
remove the fixed-effects, the conditional set cannot be separated as we have
done to move from Equation (3) to (4). For example, Kyriazidou (1997) has to
impose a “conditional exchangeability" assumption that is conditioned on the
variable related to the two periods used in differencing. In case of models
with censoring, Lee (2001) discusses conditions under which first-difference
can be applied, and applies the linear implication of the "conditional
exchangeability" assumption. In a similar context using first difference,
Rochina-Barrachina (1999) imposes a joint normality between the difference in
the error of the outcome equation and the error in the selections equation in
the two time periods.444See Dustmann and Rochina-Barrachina (2007) for a
review on selection correction in panel data models.

Estimating equation (6) presents two challenges for the identification of the
parameter of interest $\displaystyle\delta$ and the sample selection parameter
$\displaystyle\rho$: the sub-location specific difference
$\displaystyle\Delta_{d}\gamma_{j\alpha}$, and the sample selection term
$\displaystyle\rho\Delta_{d}\lambda(z_{ij}^{\prime}\beta+\theta_{j\alpha}+\theta_{j})$.
As for the sub-location specific difference
$\displaystyle\Delta_{d}\gamma_{j\alpha}$, under Assumption I1 and I2,
equation (6) becomes

$E[\Delta_{d}y_{2ij}|\xi_{ij},\xi_{kj}]=\Delta_{d}x_{ij}^{\prime}\delta+\rho\Delta_{d}\lambda(z_{ij}^{\prime}\beta+\theta_{j\alpha}+\theta_{j})$
(7)

These assumptions allow us to difference-out the sub-location specific
unobserved effect $\displaystyle\gamma_{j\alpha}$, a strategy that was applied
by Duranton, Gobillon, and Overman (2011).

As for the sample selection term
$\displaystyle\rho\Delta_{d}\lambda(z_{ij}^{\prime}\beta+\theta_{j\alpha}+\theta_{j})$,
we see that it depends on the unobservable sub-location specific and location
effects $\displaystyle\theta_{j\alpha}+\theta_{j}$. Because that sample
selection term is a nonlinear function, a simple spatial differencing will not
always work unlike the case of $\displaystyle\gamma_{ja}$. Therefore, the
following assumption helps us to deal with this challenge:

Assumption I3:
(i)

The sub-location specific unobservable selection effect is homogeneous in a
neighborhood of the individual i.e. $\displaystyle\Delta_{d}\theta_{ja}=0$ for
$\displaystyle d$ small enough.

(ii) The changes in the inverse Mill’s-Ratio in a neighborhood of the
individual i.e.

$\frac{\lambda(z_{ij}^{\prime}\beta+\theta_{j\alpha_{i}}+\theta_{j})-\lambda(z_{ij}^{\prime}\beta)}{\theta_{j\alpha_{i}}+\theta_{j}}=\lambda^{\prime}(c_{i})=\lambda^{\prime}(c_{k})=\frac{\lambda(z_{kj}^{\prime}\beta+\theta_{j\alpha_{k}}+\theta_{j})-\lambda(z_{kj}^{\prime}\beta)}{\theta_{j\alpha_{k}}+\theta_{j}}$
(8)

for $\displaystyle i$ and $\displaystyle k$ in a neighborhood $\displaystyle
d$ small enough, $\displaystyle\theta_{j\alpha_{i}}+\theta_{j}$ and
$\displaystyle\theta_{j\alpha_{k}}+\theta_{j}$ both different from 0,
$\displaystyle\lambda^{\prime}(.)$ is the first derivative the inverse Mill’s
ratio, $\displaystyle c_{i},$ and $\displaystyle c_{k}$ are, respectively, in
the intervals formed by $\displaystyle[z_{ij}^{\prime}\beta,\hskip
5.69046ptz_{ij}^{\prime}\beta+\theta_{j\alpha_{i}}+\theta_{j}]$ and
$\displaystyle[z_{kj}^{\prime}\beta,\hskip
5.69046ptz_{kj}^{\prime}\beta+\theta_{j\alpha_{k}}+\theta_{j}]$ such that
Equation (8) holds.

Assumption I3 (i) is similar to assumption I1. It seems plausible that if that
assumption holds for the outcome equation, it will hold true for the selection
equation as well.

Assumption I3 (ii) is novel and one of the contributions of this paper. It
assumes that if the exact Taylor approximation is applied on the individual
inverse Mill’s ratio for individuals $\displaystyle i$ and $\displaystyle k$
in the location $\displaystyle j$, the intermediate points $\displaystyle
c_{i}$ and $\displaystyle c_{k}$ should be similar. If the level of
nonlinearity of $\displaystyle\lambda(.)$ is low, then the assumption will
also hold. In the extreme case of local linearity of the inverse Mill’s ratio,
the Assumption 3 (ii) perfectly holds.

The combination of assumptions I3 (i) and I3 (ii) implies that

$\displaystyle\lambda(z_{ij}^{\prime}\beta+\theta_{j\alpha_{i}}+\theta_{j})-\lambda(z_{ij}^{\prime}\beta)=\lambda(z_{kj}^{\prime}\beta+\theta_{j\alpha_{k}}+\theta_{j})-\lambda(z_{kj}^{\prime}\beta).$

Thus,
$\displaystyle\Delta_{d}\lambda(z_{ij}^{\prime}\beta)=\Delta_{d}\lambda(z_{ij}^{\prime}\beta+\theta_{j\alpha}+\theta_{j})$

######  Theorem 1.

Let us consider the sample selection model presented in Equation 2. Under
assumptions I1 to I3 the parameters $\displaystyle\delta$ and
$\displaystyle\rho$ are identified.

Proof of Theorem 1
We have already shown that under the assumptions I1 and I2, we can obtain
Equation (7). Applying the assumption I3, to Equation (7) leads to the
following equation

$E[\Delta_{d}y_{2ij}|\xi_{ij},\xi_{kj}]=\Delta_{d}x^{\prime}_{ij}\delta+\rho\Delta_{d}\lambda(z_{ij}^{\prime}\beta).$
(9)

Thus, assumptions I1 to I3 are sufficient for the identification of
$\displaystyle\delta$ and $\displaystyle\rho$.

We have derived the results using the pairwise spatial difference operator.
However, the identification result holds for other spatial difference
operators as well. In the case of the average or kernel difference operator,
the conditioning in equation 9 is on $\displaystyle\xi_{kj}$ with
$\displaystyle k\in\mathcal{N}_{id}$ for the average difference operator and
$\displaystyle k$ is in the full sample for the kernel operator. Note that
under assumptions I1 and I3, any difference of the weighted average in a
neighborhood of the individual will enable us to remove the sub-location
specific effect. The conditional expectation presented in Equation (9) depends
on exogenous observable variables and parameters of interest.

###  2.2 Estimation and Asymptotic Properties

In this section, we present an estimation procedure and derive asymptotic
properties of the proposed estimator. The estimation procedure involves two-
steps. In the first step, probit model is estimated and the inverse Mill’s
ratio predicted. In the second step, a spatial difference operator differences
out both location and the sub-location specific unobserved heterogeneity. The
model is then estimated using an ordinary least square estimator. When we have
a sample of $\displaystyle N$ individuals, the estimation procedure is thus as
follows:

  * Step 1: Estimate $\displaystyle\beta$ by probit with location effect $\displaystyle\gamma_{j}$; and calculate $\displaystyle\hat{\lambda}_{i}=\lambda(z_{ij}^{\prime}\hat{\beta})$.

  * Step 2: Estimate $\displaystyle\delta$ and $\displaystyle\rho$ in the OLS regression

$\Delta_{d}y_{2ij}=\Delta_{d}x^{\prime}_{ij}\delta+\rho\Delta_{d}\lambda(z_{ij}^{\prime}\hat{\beta})+w_{ikj}.$
(10)

Since we used spatial differencing and
$\displaystyle\lambda(z_{ij}^{\prime}\hat{\beta})$ is estimated in the first
step, a particular structure of the variance-covariance matrix emerges.
Therefore, we also need to derive the correct estimator of standard errors
$\displaystyle w_{ikj}$ which we will do in section 2.3.

We will now show that the estimator obtained by the above procedure is
consistent and asymptotically normal. To derive the asymptotic properties we
use similar arguments as those used to derive the asymptotic properties of the
clustered standard errors. Specifically, the population size of each location
is assumed random and bounded almost surely, and the law of large numbers is
applied by letting the number of locations (clusters in case of clustered
standard errors) go to infinity.

We consider a generic matrix of spatial difference $\displaystyle\Delta$. The
matrix form notation of equation (10) can be expressed in as 555The variables
without subscript represent vector or matrices of all observation in the
sample.

$\Delta y_{2}=\Delta
x^{\prime}\delta+\rho\Delta\lambda(z^{\prime}\hat{\beta})+\Delta\eta$ (11)

where $\displaystyle\eta_{ij}$ are the same error as in standard sample
selection models.666We assume the notation that
$\displaystyle\lambda(z^{\prime}\hat{\beta})$ is a vector with typical element
$\displaystyle\lambda(z_{ij}^{\prime}\hat{\beta})$.  Let us denote
$\displaystyle\theta=(\delta,\rho)^{\prime}$ and $\displaystyle
W=[x^{\prime},\lambda(z^{\prime}\hat{\beta})]$. The simplified estimation
Equation (11) is

$\displaystyle\Delta y_{2}=\Delta W\theta+\Delta\eta$

and OLS estimator of $\displaystyle\theta$ is

$\hat{\theta}=[(\Delta W)^{\prime}\Delta W]^{-1}[(\Delta W)^{\prime}\Delta
y_{2}]$ (12)

The spatial nature of data implies that an observation $\displaystyle k$ with
$\displaystyle n$ neighbours may appear in several pairs. This induces
correlation in the error term $\displaystyle\Delta\eta$ for all $\displaystyle
n$ of these pairs because of the spatial differencing in the second step of
the estimation procedure. As a result, a particular structure of the
covariance matrix emerges, and we need to take that into account when
calculating the standard errors.

To proceed further, we need to introduce assumptions under which the
asymptotic properties of our estimator are derived.
Assumption E1: The sample is formed if $\displaystyle N$ individuals from the
population.
$\displaystyle(i)$ We observed $\displaystyle\\{x_{ij},z_{ij}\\}$ independent
and identically distributed random variable with $\displaystyle i=1,...,N$ and
$\displaystyle j=1,...,J$.
$\displaystyle(ii)$ The number of individuals in a location $\displaystyle j$,
$\displaystyle N_{j}$, is exogenous, random, identically distributed with
$\displaystyle N_{j}<n_{0}$ almost surely and $\displaystyle E(N_{j})<\infty$;
where $\displaystyle n_{0}$ is a scalar.
$\displaystyle(iii)$ The outcomes and the latent variables are independent
across location i.e. $\displaystyle j_{1}\neq j_{2}$ the variables
$\displaystyle y_{2ij_{1}}\bot y_{2ij_{2}}$ and $\displaystyle
y^{*}_{1ij_{1}}\bot y^{*}_{1ij_{2}}$.

An implication of assumption E1 $\displaystyle(i)$ in conjunction with
assumption I2 is that $\displaystyle\theta_{j}$ and $\displaystyle\gamma_{j}$
are $\displaystyle iid$. However, $\displaystyle within$ a location
$\displaystyle j$, there is a certain level of correlation among individuals
which operates through $\displaystyle\theta_{j\alpha_{i}}$ or
$\displaystyle\gamma_{j\alpha_{i}}$. This means that our assumptions restrict
how that within-location individual correlations occurs.

Assumption E1 $\displaystyle(ii)$ restricts the location size to be bounded
and implies that the number of locations has to grow to achieve a large sample
size in our asymptotic calculation. This assumption is similar to those held
in the literature of cluster samples asymptotic and it leads to a “large
number of cluster" asymptotic theory similar to the one discussed in
Wooldridge (2010), who assumes fixed cluster size. This assumption corresponds
to a specific case of the Assumption 1 in Hansen and Lee (2019), who allow for
different cluster size ranging from fixed to infinite. We have, however,
derived the asymptotic of our estimator under the more restrictive condition
of Assumption E1 ($\displaystyle ii$). The reason is that it can be proven
that under a joint asymptotic ($\displaystyle N,J\rightarrow\infty$),
Assumption 1 is equivalent to assuming that the size of the sample in each
location is bounded. If we instead allow for a sequential asymptotic where the
number of locations is fixed and the sample size goes to infinity, then there
exists at least one location with an infinite number of individuals and the
inequality used in the proof of Hansen and Lee (2009)’s Theorem 1 becomes
invalid.

To better illustrate our argument, let us consider the location sample size
proposed by Hansen and Lee (2019): $\displaystyle N_{j}=N^{\alpha}$ with
$\displaystyle 0\leq\alpha<1$; we can prove that $\displaystyle
1-\alpha=\frac{ln(J)}{ln(N)}$. If we allow for a joint asymptotic,
$\displaystyle\alpha$ is not define. If on the contrary we assume that the
number of locations $\displaystyle J$ is fixed, then, $\displaystyle\alpha$
goes to $\displaystyle 1$. In both cases, relaying on Hansen and Lee (2009)’s
Assumption 1 seems not enough to warrant the desire asymptotic regularities.

Assumption E2: $\displaystyle z^{\prime}$ and $\displaystyle W$ are full rank
column, with each element having up to its $\displaystyle 4^{th}$ moment.

######  Theorem 2.

We consider the sample selection model presented in Equation 2. Under
assumptions I1 to I3, E1 and, E2.
$\displaystyle(i)$ $\displaystyle\hat{\theta}\rightarrow^{p}\theta$ as
$\displaystyle N\rightarrow\infty$
$\displaystyle(ii)$
$\displaystyle\sqrt{N}(\hat{\theta}-\theta)\rightarrow^{d}\mathcal{N}(0,\Theta)$
with $\displaystyle\Theta=C\Gamma C^{\prime}$
where $\displaystyle C^{-1}=E((\Delta W_{ij})^{\prime}\Delta W_{ij}),$
$\displaystyle\Gamma=\rho^{2}E[(\Delta W_{ij})^{\prime}\Omega_{ij}\Delta
W_{ij}]+E[(\Delta W_{ij})^{\prime}\Delta e_{ij}\Delta e_{ij}(\Delta W_{ij})],$
and
$\displaystyle\Omega_{ij}=[\lambda^{\prime}(z_{ij}^{\prime}\beta)]^{2}z_{ij}^{\prime}V_{\beta}z_{ij}$
taking $\displaystyle V_{\beta}$ as the first step probit variance-covariance
matrix.

Proof of Theorem 2: In appendix.

It is important to notice that the same type of asymptotics should be used in
a linear model. In this respect, we complement Duranton, Gobillon, and Overman
(2011) who propose a correction for the standard errors, but do not discuss
the asymptotic properties of their estimators. Similarly, Black (1999) and
Holmes (1998) use spatial differencing, but do not account for the fact that
differencing will lead to a correlation between pairs where an individual is
present. Our asymptotic derivations do account for the presence of correlation
between pairs, and are valid not only for a model with but also without sample
selection (in our model, the absence of selection implies
$\displaystyle\rho=0$). They also have important practical implications: the
consistency of the estimator requires a large number of locations
$\displaystyle\gamma_{j}$, and a small number of individuals in each sub-
location $\displaystyle\gamma_{j\alpha}$

###  2.3 Estimator of Variance

This section derives a procedure estimate the variance-covariance of the
estimator in Equation (12) which has a particular structure arising from
$\displaystyle(i)$ spatial differencing and $\displaystyle(ii)$ Heckman’s two-
step estimation procedure.

We consider $\displaystyle B=\left[(\Delta W)^{\prime}\Delta W\right]^{-1}$
and $\displaystyle\Sigma=Var[(\Delta W)^{\prime}\Delta\eta]$ such that the
conditional variance-covariance matrix of $\displaystyle\hat{\theta}$ is

$Var(\hat{\theta})=B\Sigma B^{\prime}$

Note that

$\Sigma=(\Delta W)^{\prime}Var(\Delta\eta)(\Delta W)$

This means that we need a consistent estimator of $\displaystyle
Var(\Delta\eta)$ to compute correct standard error for
$\displaystyle\hat{\theta}$.

Let us consider that $\displaystyle Var(\Delta\eta)=V_{1}+V_{2}$ with

$\displaystyle\displaystyle V_{1}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\Delta Var(e)\Delta^{\prime}$
$\displaystyle\displaystyle=\rho^{2}$ $\displaystyle\displaystyle\Delta
R\Delta^{\prime}$

where $\displaystyle R$ a diagonal matrix of dimension $\displaystyle N$
(total number of observations), with $\displaystyle
d_{ij}=1-\lambda(z_{ij}^{\prime}\beta)[z_{ij}^{\prime}\beta+\lambda(z_{ij}^{\prime}\beta)]$
as the diagonal elements.

$\displaystyle\displaystyle V_{2}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\rho^{2}\Delta
Var\left[\lambda(z^{\prime}\hat{\beta})-\lambda(z^{\prime}\beta)\right]\Delta^{\prime}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\rho^{2}\Delta
DzV_{\beta}z^{\prime}D\Delta^{\prime}$

where $\displaystyle D$ is the square, diagonal matrix of dimension
$\displaystyle N$ with $\displaystyle 1-d_{ij}$ as the diagonal elements;
$\displaystyle z$ is the data matrix of selection equation; and $\displaystyle
V_{p}$ is the variance-covariance estimate from the probit estimation of the
selection equation.

######  Theorem 3.

We consider the sample selection model presented in Equation 2. Under
assumptions I1 to I3, E1 and, E2. The variance-covariance estimator of the
$\displaystyle\hat{\theta}$ is given by

$V_{twostep}=B(\Delta W)^{\prime}[\hat{V}_{1}+\hat{V}_{2}](\Delta
W)B^{\prime}$ (13)

where $\displaystyle\hat{V}_{1}=\hat{\rho^{2}}\Delta\hat{R}\Delta^{\prime}$
and
$\displaystyle\hat{V}_{2}=\hat{\rho}^{2}\Delta\hat{D}z\hat{V}_{\beta}z^{\prime}\hat{D}\Delta^{\prime}$
with all unknown parameters replaced by their estimates. Moreover, this is a
consistent estimator $\displaystyle Var(\hat{\theta})$.

Proof of Theorem 3:
The result holds by construction.

##  3 Monte Carlo Simulation

In this section we present the results of Monte Carlo simulations to (i) to
describe the behavior of the estimator proposed in this paper and (ii) offer
empirical guidance for applied research. Regarding the latter, we will pay a
close attention to the implication of assumption E1 $\displaystyle(ii)$
according to which it is important to have a large number of locations
relative to the number of individuals in the sub-locations. Monte Carlo
experiments will offer empirical guidance as to when the number of locations
is large enough.

The estimator developed in this paper is referred to as the “Sub-location
Differencing" and it accounts for sub-location specific effect
$\displaystyle\gamma_{j\alpha}$. To highlight its features, we compare it to
two other estimators. One ignores the presence of both
$\displaystyle\gamma_{j}$ and $\displaystyle\gamma_{j\alpha}$ and applies a
simple two step estimator with no spatial differencing - we call it “No-
Differencing" estimator. The other accounts only for the location fixed
effects $\displaystyle\gamma_{j}$ and we call it “Location Differencing"
estimator. For each estimator, the mean bias and the coverage rate for the 95%
confidence level test are reported in tables I to III.

The data is obtained using the following data generating process. We assume
that there are $\displaystyle J=20,30,100$ non-overlapping locations, each
location is divided into $\displaystyle s=2,4,8$ sub-locations. There are
$\displaystyle n_{j}=3,5,8,10$ individuals sharing the same sub-location. The
latent variables are $\displaystyle
y^{*}_{1ij}=z_{ij}\beta+\theta_{ijs}+\theta_{j}+\varepsilon_{1ij}$ and
$\displaystyle
y^{*}_{2ij}=x_{ij}\delta+\gamma_{ijs}+\gamma_{j}+\varepsilon_{2ij}$, where
$\displaystyle\theta_{ija}=10^{-5}j\times s$ and
$\displaystyle\gamma_{ijs}=5j\times s$ is the sub-location specific effect,
while $\displaystyle\theta_{j}=10^{-5}j$ and $\displaystyle\gamma_{j}=10j$ are
the location effects; for all $\displaystyle i$ and $\displaystyle j$,
$\displaystyle x_{ij}\sim\mathcal{N}(0,1)$ , $\displaystyle z_{ij}\sim U(0,1)$
each drawn independently; $\displaystyle\delta=1$, $\displaystyle\beta=0.2$.
The error terms in both equations for all $\displaystyle i$ and $\displaystyle
j$ are generated as follows:
$\displaystyle\varepsilon_{1ij}\sim\mathcal{N}(0,1)$,
$\displaystyle\varepsilon_{2ij}=\rho\varepsilon_{1ij}+v_{ij}$ where
$\displaystyle v_{ij}\sim\mathcal{N}(0,1)$ is independent of
$\displaystyle\varepsilon_{1ij}$ and $\displaystyle\rho=0.7$.

Table I: Simulation Results with 20 Locations Numb. Of sub-location | sub-location-size | Estimators | Mean bias | Coverage rate
---|---|---|---|---
2 | 3 | No-differencing | -0.079 | 95.8
Location Differecing | -0.393 | 74.2
Sub-location Differencing | -0.039 | 82.2
5 | No-differencing | -0.344 | 94.5
Location Differecing | -0.953 | 82.1
Sub-location Differencing | 0.011 | 87.0
8 | No-differencing | -0.173 | 94.4
Location Differecing | -1.307 | 90.7
Sub-location Differencing | 0.016 | 79.5
10 | No-differencing | 0.233 | 95.8
Location Differecing | -2.264 | 93.6
Sub-location Differencing | 0.067 | 93.6
4 | 3 | No-differencing | 0.369 | 95.3
Location Differecing | 0.117 | 88.1
Sub-location Differencing | 0.001 | 75.5
5 | No-differencing | 0.996 | 96.0
Location Differecing | 2.798 | 92.7
Sub-location Differencing | -0.019 | 82.5
8 | No-differencing | 0.085 | 94.5
Location Differecing | 5.824 | 95.2
Sub-location Differencing | -0.035 | 80.5
10 | No-differencing | -0.241 | 95.6
Location Differecing | -2.494 | 96.1
Sub-location Differencing | -0.066 | 88.2
8 | 3 | No-differencing | 0.833 | 94.7
Location Differecing | -0.310 | 94.9
Sub-location Differencing | -0.005 | 66.3
5 | No-differencing | -0.176 | 93.6
Location Differecing | 1.678 | 96.8
Sub-location Differencing | 0.011 | 76.8
8 | No-differencing | 0.013 | 95.4
Location Differecing | -0.269 | 99.0
Sub-location Differencing | 0.006 | 83.3
10 | No-differencing | 0.271 | 95.2
Location Differecing | -2.573 | 99.3
Sub-location Differencing | -0.016 | 86.1
Table II: Simulation Results with 30 Locations Numb. Of sub-location | sub-location-size | Estimators | Mean bias | Coverage rate
---|---|---|---|---
2 | 3 | No-differencing | -0.392 | 95.3
Location Differecing | 0.497 | 75.1
Sub-location Differencing | -0.006 | 78.4
5 | No-differencing | 0.095 | 95.2
Location Differecing | -0.965 | 84.5
Sub-location Differencing | 0.007 | 84.8
8 | No-differencing | -0.108 | 94.8
Location Differecing | 1.330 | 92.3
Sub-location Differencing | -0.053 | 85.8
10 | No-differencing | 0.043 | 94.5
Location Differecing | -2.052 | 95.8
Sub-location Differencing | -0.028 | 79.2
4 | 3 | No-differencing | -0.177 | 95.2
Location Differecing | -1.027 | 89.5
Sub-location Differencing | 0.004 | 71.3
5 | No-differencing | -0.227 | 95.3
Location Differecing | -1.387 | 92.8
Sub-location Differencing | -0.017 | 78.7
8 | No-differencing | 0.424 | 93.6
Location Differecing | -1.437 | 95.0
Sub-location Differencing | 0.012 | 84.2
10 | No-differencing | -0.156 | 95.1
Location Differecing | 0.348 | 96.8
Sub-location Differencing | 0.025 | 78.9
8 | 3 | No-differencing | 0.031 | 95.0
Location Differecing | -0.200 | 95.0
Sub-location Differencing | 0.010 | 67.4
5 | No-differencing | 0.279 | 94.4
Location Differecing | 2.118 | 97.0
Sub-location Differencing | -0.008 | 73.6
8 | No-differencing | -0.101 | 93.2
Location Differecing | 4.108 | 98.0
Sub-location Differencing | 0.005 | 80.8
10 | No-differencing | -1.369 | 95.6
Location Differecing | -1.866 | 99.5
Sub-location Differencing | -0.041 | 83.6
Table III: Simulation Results with 100 Locations Numb. Of sub-location | sub-location-size | Estimators | Mean bias | Coverage rate
---|---|---|---|---
2 | 3 | No-differencing | -0.702 | 94.9
Location Differecing | 1.017 | 74.0
Sub-location Differencing | 0.007 | 63.6
5 | No-differencing | -0.716 | 95.3
Location Differecing | -2.314 | 82.0
Sub-location Differencing | -0.001 | 72.2
8 | No-differencing | -0.089 | 94.6
Location Differecing | -5.000 | 90.5
Sub-location Differencing | -0.005 | 84.2
10 | No-differencing | 0.198 | 94.9
Location Differecing | -0.197 | 96.2
Sub-location Differencing | -0.050 | 85.3
4 | 3 | No-differencing | -0.096 | 94.5
Location Differecing | 0.858 | 89.1
Sub-location Differencing | -0.001 | 59.4
5 | No-differencing | -0.691 | 94.7
Location Differecing | 0.712 | 89.6
Sub-location Differencing | 0.001 | 64.0
8 | No-differencing | -0.177 | 96.7
Location Differecing | 1.216 | 96.9
Sub-location Differencing | 0.015 | 79.2
10 | No-differencing | -0.384 | 94.0
Location Differecing | -10.139 | 97.9
Sub-location Differencing | -0.075 | 79.6
8 | 3 | No-differencing | -1.018 | 95.2
Location Differecing | 0.598 | 94.1
Sub-location Differencing | 0.002 | 53.7
5 | No-differencing | 0.050 | 95.1
Location Differecing | 1.400 | 95.7
Sub-location Differencing | -0.001 | 61.0
8 | No-differencing | -1.093 | 95.7
Location Differecing | -10.195 | 98.1
Sub-location Differencing | 0.010 | 75.7
10 | No-differencing |  -0.642 | 94.7
Location Differecing | 6.82 | 98.9
Sub-location Differencing | 0.009 | 79.1

We summarize the main results of the simulations in the four points below but
in general, the “Sub-location Differencing" estimator has the smallest mean
bias and delivers a conservative coverage rate.777There is room for
improvement concerning our inference strategy. Cluster robust inference is
part of a large and growing literature and our work gives some insight as to
how diffrencing can be used in cross-sectional data. Future work will
investigate the importance of heteroscedasticity, and small sample procedures
such as bootstrap will be used to improve inference.

  1. 1.

As expected, the “No-Differencing" estimator has a larger mean bias in the
presence of spatial heterogeneity. This result holds for both small and large
numbers of locations as well as for few or many individuals having the same
sub-location specific unobserved heterogeneity
$\displaystyle\mathcal{N}_{id}$.

  2. 2.

The mean bias of the “Sub-location Differencing" estimator is smaller than
other estimators. It increases with the number of individuals in the sub-
locations, and decreases with the number of locations. For example, in a
sample of 600 individuals which are spread across 100 locations with 2 sub-
locations and 3 individuals in each sub-location, the mean bias is of 0.007.
However, for the same sample size but spread across 30 locations with 2 sub-
locations and 10 individuals in each sub-locations, the mean bias is
$\displaystyle-0.028$. This result is in line with our asymptotic derivations.

  3. 3.

For a fixed number of locations, the bias increases with the number of
individuals in sub-locations. The empirical consequence of this results is
that our estimator should be applied when the size of sub-locations is small.

  4. 4.

The coverage rate of the “Sub-location Differencing" estimator using the
variance-covariance estimator in Equation (13) and the normal asymptotic
distribution suggests a conservative coverage.

We now turn to an empirical application of our estimation strategy.

##  4 Empirical Application

This section shows the empirical importance of spatial differencing
methodology proposed in the previous sections. To illustrate the importance of
our estimator, we ask what determines tax rates set by regional governing
bodies.888There is a large literature which examines a range of factors
influencing local tax rates e.g. Charney (1983), Ashworth and Heyndels (1997),
Ross and Yinger (1999), Charlot and Paty (2007); Charlot, Paty, and Piguet
(2015), Crowley and Sobel (2011), Baskaran (2014), Buettner and von Schwerin
(2016). This question opens an important issue of identification since
circular causation or omitted variable bias leads to biased and inconsistent
estimators. We will use our spatial differencing method to examine the case of
changes in the Finish local property tax rate at the turn of the millennium.

Finland consists of 411 municipalities (in 1999) spread across 19 regions
which choose property tax rate within the limits set by the central
government. In 1999, the central government decided to raise the lower limit
for the year 2000 from 0.2% to 0.5%. This change created a probability mass of
municipalities at the lower bound: more than half of municipalities have a
taxation rate of 0.5%, making the data sample censored. We investigate what
affected municipalities tax rate in the year 2000.

We estimate the parameter of the outcome equation in the model represented as
in Equation (2). Specifically, the outcome variable is the level of general
property tax in a municipality $\displaystyle i$ in a region $\displaystyle
j$, and explanatory variables include municipality’s $\displaystyle i$ age
structure of the population, level of municipality’s income, received
subsidies, local income tax rate and a dummy for region $\displaystyle j$ in
which the municipality $\displaystyle i$ is located. The selection equation
determines whether the municipality sets its general tax rate at the mandatory
minimum of 0.5% or above and contains all the variables which are in the
outcome equation except for local income tax rate.

As illustrated in Equation (1), there can be an unobserved sub-location
specific effect operating at a finer spatial scale than region $\displaystyle
j$, in our case at the level of municipalities which region $\displaystyle j$
consists of. Indeed, municipalities’ tax level can depend not only on its
population size, income and subsidies received from the central government,
but it can also depend on the level of amenities in the municipality. It is
usually difficult to measure them. More importantly, even if we have a few
measures of amenities or their proxies, they might not be able to capture all
of them, leaving some amenities unobserved. In our case, unobserved amenities
can be correlated with municipalities’ population, income level, or the level
of subsidies which implies that not controlling for them will render the
estimates biased and inconsistent. Therefore, using the fact that two
municipalities from the same region sharing borders are neighbor, we use our
spatial differencing method to tackle this problem.

We estimate equation (2) with spatial differencing conducted as the difference
between municipality $\displaystyle i$ and the average of its neighbours.
Columns 1 and 2 present the results without spatial differencing, columns 3-6
with spatial differencing. Estimation is conducted with and without regional
dummies respectively, and with two different estimators of standard errors:
wild cluster bootstrap, and spatially-adjusted standard errors derived in
Section 2.3. Clustering of the standard errors is done at the level of region
$\displaystyle j$. Since there are only 19 regions, we use wild clustered
bootstrap procedure developed by Cameron, Gelbach, and Miller (2008), which
properties were studied by e.g. Davidson and MacKinnon (2010), MacKinnon
(2013), and MacKinnon and Webb (2017). Specifically, we use a recently
developed wild bootstrap package $\displaystyle boottest$ by Roodman, Nielsen,
MacKinnon, and Webb (2019) implemented in Stata.

We begin with discussing the results with spatial differencing: columns 3-6,
Columns 3 and 4 present the results when we use spatial differencing with the
standard errors calculated using the formula derived in Section 2.3. There is
only one significant variable: the share of population older than 75 in column
4. This is not surprising, since our estimator of variance-covariance matrix
is an asymptotic estimator, while the estimation is done on a sample with
small number of clusters (nineteen regions). Therefore we use wild-bootstrap
procedure which is known to be suitable for a small number of clusters. The
results using wild-bootstrapping are shown in columns 5 and 6 and we see a
considerable increase in the number of statistically significant results.

The comparison of columns 5 and 6 with columns 1 and 2 reveals the importance
of spatial differencing. Controlling for the sub-location specific unobserved
effect $\displaystyle\gamma_{j\alpha}$ by spatial differencing renders four
variables statistically significant: share of population younger than 15,
share of population older than 75, government grants, and income tax rate.
This is in contrast to columns 1 and 2 in which income tax rate is the only
significant variable. Not controlling for $\displaystyle\gamma_{j\alpha}$
leads to the omitted variable bias which, apart from rendering the estimates
inconsistent, inflates the standard errors and makes the estimates mostly
insignificant. Spatial differencing controls for this omitted variable bias,
which means that they do not ’end up’ in the error term and do not inflate the
standard errors. In addition to comparing estimates with and without spatial
differencing, it is also instructive to compare columns 5 and 6: spatial
differencing with and without regional dummies. We see that controlling for
the regional unobserved effect $\displaystyle\gamma_{j}$ does not help to
control for sub-location specific effects $\displaystyle\gamma_{j\alpha}$.
Indeed, the magnitude and the statistical significance of the estimates change
very little, and even the one variable that looses its statistical
significance after including regional dummies (municipality’s income) is only
marginally significant without these dummies.

Overall, our empirical analysis shows that controlling for the unobserved
municipality effects matters. Estimations which control for spatial unobserved
effects only at the regional suggest that the income tax rate is the only
determinant of general tax rate set by the municipalities. However, after
controlling for the sub-location specific unobserved effects, the tax rate
depends not only on the income tax rate but also on the age composition of
their population - the share of young as well as the share of elderly
population. These results thus indicate that spatial differencing is an
important tool to deal with omitted variable bias which often plagues
empirical studies on local taxation.

Table IV: Determinants of Municipality Taxation Rate | No Spatial Differencing | Spatial Differencing
---|---|---
Wild bootstrap | Spatially adj se | Wild bootstrap
No Reg. Dummies | Regional Dummies | No Reg. Dummies | Regional Dummies | No Reg. Dummies | Regional Dummies
| (1) | (2) | (3) | (4) | (5) | (6)
Population | -1.897 | -0.7047 | 0.1346 | 0.2448 | 0.1346 | 0.2448
| [-1.420] | [-1.0432] | [0.0489] | [0.2656] | [0.507] | [0.8793]
Share pop.<15 | 0.009 | -0.0015 | -0.0116 | -0.0113 | -0.0116** | -0.0113**
| [0.492] | [-0.0996] | [-0.0933] | [-0.1373] | [-2.536] | [-2.3921]
Share 61<pop.<74 | -0.0055 | -0.0054 | -0.0089 | -0.0092 | -0.0089 | -0.0092
| [-0.923] | [-0.343] | [-0.1106] | [-0.0721] | [-1.495] | [-1.4902]
Share pop.>75 | 0.0211 | 0.0049 | -0.0145 | -0.0140** | -0.0145* | -0.014*
| [1.027] | [0.4688] | [-0.1575] | [-1.8771] | [-1.794] | [-1.6746]
Income | 2.40E-07 | 6.48E-06 | 1.23E-05 | 0.00001 | 1.23E-05* | 1.1E-05
| [0.047] | [0.8092] | [0.0234] | [0.0213] | [1.924] | [1.6688]
Gov. grant | -1.7E-05 | 1.8E-05 | -1.78E-05 | -2.16E-05 | -1.78E-05 | -2.16E-05
| [-0.868] | [0.2555] | [-0.0458] | [-0.0925] | [-0.483] | [-0.5666]
Income tax rate | 0.0350*** | 0.0396** | 0.0482 | 0.0453 | 0.0482*** | 0.0453***
| [3.383] | [3.547] | [0.2775] | [0.2965] | [3.362] | [3.088]
Inverse Mill’s ratio | 0.3028 | 0.1731 | 0.0159 | 0.0034 | 0.0159 | 0.0034
| [1.505] | [1.1468] | [0.0346] | [0.0298] | [0.854] | [0.1525]
Constant | -0.5327 | -0.3601 | 0.0111 | -0.006 | 0.0111 | -0.006
| [-0.713] | [-0.598] | [1.540] | [-0.230] | [1.540] | [-0.199]
Observations | 403 | 403 | 273 | 273 | 273 | 273
Regional Dummies | NO | YES | NO | YES | NO | YES
Number of Dummies | 19 | 19 | 19 | 19 | 19 | 19
R-squared | 0.197 | 0.271 | 0.248 | 0.279 | 0.248 | 0.279
Source: see text; Note: t-statistics in brackets, *** p<0.01, ** p<0.05, *
p<0.1. The table reports regression results with municipality taxation rate as
the dependent variable. Selection equation excludes local income tax rate.
Wild bootstrapping with 999 interactions was used except for column 4 which
reports t-statistics calculated with the standard errors adjusted for spatial
differencing.

##  5 Conclusion

This paper has investigated a sample selection model with unobserved
heterogeneity at a very fine location level. It proposes spatial differencing
as an alternative identification strategy when instrumental variable and/or a
panel data are not available. We discuss the assumptions under which the
parameters of the model are identified. The estimation of the parameters is
done using the classic Heckman’s two-step estimation procedure. The
differecing and the two-step procedure lead to a novel estimator with
properties that are also relevant for spatial differencing in linear models.
To understand the behavior of the new estimator, we derive a cluster
asymptotic of the estimator. The derivation reveals two important implications
for its empirical implementation: $\displaystyle(i)$ the number of clusters
needs to be large for inference to be based on normal distribution.
($\displaystyle ii$) each cluster should have a bounded number of individuals.

Monte Carlo experiments show that accounting for sub-location specific
heterogeneity is crucial for identification. It also confirms the estimator’s
properties derived in our asymptotic. In particular, the estimator performs
better with the increasing number of locations, and fewer individuals in sub-
locations. In addition, ignoring sub-locations and applying spatial
differencing only to more aggregate geographical units subsuming sub-
locations, the mean bias is larger. The coverage rate of the test based on the
corrected standard error has an empirical coverage lower that the theoretical
one.

In the empirical application which looked at the determinants of municipal tax
rate, we show that using spatial differencing in combination with cluster
wild-bootstrap inference tools can be extremely useful. Indeed, the new
estimator reveals several determinants of the municipal tax rate that would
have been missed otherwise. The development of a bootstrap appropriate sample
selection models is left for future research.

## Appendix

Proof of Theorem 2:
The proof is written conditional on the set of number of individuals in the
locations. Thus, when $\displaystyle E(N_{j})$ is used, it can be considered
as a constant.

The substitution of the true value of $\displaystyle\Delta y_{2}$ in equation
(12) yields the following equality

$\hat{\theta}=\theta+[(\Delta W)^{\prime}\Delta W]^{-1}[(\Delta
W)^{\prime}\Delta\eta]$

Let us assume that $\displaystyle
y_{2ij}=x_{ij}^{\prime}\delta+\gamma_{j\alpha}+\gamma_{j}+\rho\lambda(z_{ij}^{\prime}\beta+\theta_{j\alpha_{k}}+\theta_{j})+e_{ij}$
with $\displaystyle E(e_{ij}|\xi_{ij})=0$. Thus, $\displaystyle\Delta
y_{2ij}=\Delta
x_{ij}^{\prime}\delta+\rho\Delta\lambda(z_{ij}^{\prime}\beta+\theta_{j\alpha}+\theta_{j})+\Delta
e_{ij}$. Under the identification assumptions I1 to I3 have

$\displaystyle\Delta y_{2ij}=\Delta
x_{ij}^{\prime}\delta+\rho\Delta\lambda(z_{ij}^{\prime}\beta)+\Delta e_{ij}.$

The second step regression equation is equivalent to

$\displaystyle\Delta y_{2ij}=\Delta
x_{ij}^{\prime}\delta+\rho\tilde{\Delta}\lambda(z_{ij}^{\prime}\beta)+\Delta[\rho(\lambda(z_{ij}^{\prime}\beta)-\lambda(z_{ij}^{\prime}\hat{\beta}))+e_{ij}]$

$\displaystyle\hat{\beta}$ is estimated by maximum likelihood probit in the
first step with variance-covariance matrix $\displaystyle V_{\beta}$. Given
that $\displaystyle\lambda(.)$ is twice differentiable, the continuous mapping
theorem implies that
$\displaystyle\lambda(z_{ij}^{\prime}\beta)-\lambda(z_{ij}^{\prime}\hat{\beta})$
goes to zero in probability and is asymptotically normal. If we assume that
$\displaystyle N_{1}$ is the full sample while $\displaystyle N$ is the
selected sample.999 We assume that $\displaystyle N/N_{1}\rightarrow 1$. We,
therefore, have

$\sqrt{N_{1}}(\lambda(z_{ij}^{\prime}\beta)-\lambda(z_{ij}^{\prime}\hat{\beta}))\rightarrow^{d}\mathcal{N}(0,\Omega_{ij})$
(14)

where
$\displaystyle\Omega_{ij}=[\lambda^{\prime}(z_{ij}^{\prime}\beta)]^{2}z_{ij}^{\prime}V_{\beta}z_{ij}$.

We are interested in the limiting distribution of
$\displaystyle\sqrt{N}(\hat{\theta}-\theta)$.

$\displaystyle\displaystyle\sqrt{N}(\hat{\theta}-\theta)$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle N[(\Delta
W)^{\prime}\Delta W]^{-1}\frac{1}{\sqrt{N}}[(\Delta W)^{\prime}\Delta\eta]$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\left[\frac{\sum_{i=1}^{N}(\Delta
W_{i})^{\prime}\Delta
W_{i}}{N}\right]^{-1}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}(\Delta
W_{i})^{\prime}\Delta\eta_{i}$

While $\displaystyle W_{i}$ are $\displaystyle iid$, $\displaystyle\Delta
W_{i}$ are not independent because an individual is allowed to appear in many
pairs. We will therefore have to use LLN and CLT for non-independent random
variables. The dependence structure is driven by the operator
$\displaystyle\Delta$. If $\displaystyle\Delta$ is such that each individual
appear only in one pair then the classical CLT and LLN could be applied.
However, if individuals are allowed to appear in several pairs, then we need
to apply CLT and LLN accounting for correlation.

$\displaystyle\displaystyle\frac{\sum_{i=1}^{N}(\Delta W_{i})^{\prime}\Delta
W_{i}}{N}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\frac{\sum_{j=1}^{J}\sum_{k=1}^{N_{j}}(\Delta
W_{kj})^{\prime}\Delta W_{kj}}{N}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\frac{1}{J}\sum_{j=1}^{J}\frac{1}{E(N_{j})}\sum_{k=1}^{N_{j}}(\Delta
W_{kj})^{\prime}\Delta W_{kj}$

Let us consider $\displaystyle
Y_{j}=\frac{1}{E(N_{j})}\sum_{k=1}^{N_{j}}(\Delta W_{kj})^{\prime}\Delta
W_{kj}$, these variables are $\displaystyle iid$, moreover, note that
$\displaystyle N=N_{1}+N_{2}+....+N_{J}=JE(N_{j})$. Under assumption E1, all
locations have a bounded maximum capacity of $\displaystyle N_{j}<n_{0}$ with
$\displaystyle n_{0}$ a scalar. Under the assumption that all second moments
of the variables in $\displaystyle W$ exist (Assumption E2).
$\displaystyle
H_{J}=\frac{1}{J}\sum_{j=1}^{J}\frac{1}{E(N_{j})}\sum_{k=1}^{N_{j}}Y_{j}$ is a
matrix. Thus, the law of large number apply to it if and only if it applies to
all is elements. Let $\displaystyle a_{j}$ be a typical element of the matrix
$\displaystyle\frac{1}{E(N_{j})}\sum_{k=1}^{N_{j}}Y_{j}$. Let $\displaystyle
t$ and $\displaystyle m$ be two variables from the set of variables forming
$\displaystyle W$. For example, we can consider $\displaystyle t=x_{1}$ the
first column of the random variable $\displaystyle x$. If $\displaystyle t\neq
m$ then,

$\displaystyle\displaystyle E|a_{j}|$ $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\frac{1}{E(N_{j})}\sum_{k=1}^{N_{j}}E|\Delta
t_{kj}\Delta m_{kj}|$ (15) $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\frac{1}{E(N_{j})}\sum_{k=1}^{N_{j}}E|(t_{kj}-t_{ij})(m_{kj}-m_{ij})|$
(16) $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\frac{4}{E(N_{j})}\sum_{k=1}^{N_{j}}E|t_{kj}m_{kj}|$
(17) $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\frac{4}{E(N_{j})}\sum_{k=1}^{N_{j}}\sqrt{E(|t_{kj}|^{2})E(|m_{kj}|^{2})}$
(18) $\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle M_{0}$ (19)

with $\displaystyle M_{0}$ a constant.

The result is obtained by using successively the triangular inequality, the
identical distribution of variable in $\displaystyle W$, the Cauchy-Schwars’s
inequality and the existence of moment up-to its fourth (which means that the
second moment exists). If $\displaystyle t=m$ we have,

$\displaystyle\displaystyle E|a_{j}|$ $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\frac{1}{E(N_{j})}\sum_{k=1}^{N_{j}}E|\Delta
t_{kj}\Delta t_{kj}|$ (20) $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\frac{1}{E(N_{j})}\sum_{k=1}^{N_{j}}E|(t_{kj}-t_{ij})^{2}|$
(21) $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\frac{2}{E(N_{j})}\sum_{k=1}^{N_{j}}E|t_{kj}t_{ij}|+E(t_{kj}^{2})$
(22) $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\frac{1}{E(N_{j})}\sum_{k=1}^{N_{j}}(E|t_{kj}|)^{2}+E(t_{kj}^{2})$
(23) $\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle M_{0}$ (24)

Thus, the LLN implies that

$\displaystyle\frac{\sum_{i=1}^{N}(\Delta W_{i})^{\prime}\Delta
W_{i}}{N}\rightarrow^{p}E((\Delta W_{ij})^{\prime}\Delta W_{ij})=C^{-1}.$

We can also show that

$\displaystyle\displaystyle\frac{1}{\sqrt{N}}\sum_{i=1}^{N}(\Delta
W_{i})^{\prime}\Delta\eta_{i}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\frac{\rho}{\sqrt{N}}\sum_{i=1}^{N}(\Delta
W_{i})^{\prime}\Delta(\lambda(z_{ij}^{\prime}\beta)-\lambda(z_{ij}^{\prime}\hat{\beta}))$
$\displaystyle\displaystyle+$
$\displaystyle\displaystyle\frac{1}{\sqrt{N}}\sum_{i=1}^{N}(\Delta
W_{i})^{\prime}\Delta e_{ij}.$

We consider $\displaystyle\Lambda_{j}=\sum_{k=1}^{N_{j}}(\Delta
W_{kj})^{\prime}\Delta(\lambda(z_{kj}^{\prime}\beta)-\lambda(z_{kj}^{\prime}\hat{\beta}))$
and $\displaystyle E_{j}=\sum_{k=1}^{N_{j}}(\Delta W_{kj})^{\prime}\Delta
e_{kj}$. Conditional on $\displaystyle\hat{\beta}$, $\displaystyle\Lambda_{j}$
are $\displaystyle iid$ random variables; $\displaystyle E_{j}$ are too. We
assume that the number of individuals in a group is $\displaystyle iid$ with
finite mean $\displaystyle E(N_{j})$.101010The application of the LLN implies
for consistency reason that $\displaystyle N/J\rightarrow^{p}E(N_{j}).$ Thus
$\displaystyle JE(N_{j})\approx N.$ Given all locations are assumed to be
disjoint,

$\displaystyle\displaystyle\frac{1}{\sqrt{N}}\sum_{i=1}^{N}(\Delta
W_{i})^{\prime}\Delta\eta_{i}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\frac{\rho}{\sqrt{N}}\sum_{j=1}^{J}\Lambda_{j}$
$\displaystyle\displaystyle+$
$\displaystyle\displaystyle\frac{1}{\sqrt{N}}\sum_{j=1}^{J}E_{j}.$

We have $\displaystyle E(E_{j})=0$ for each $\displaystyle j.$ Moreover,

$\displaystyle\displaystyle Var(E_{j})$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle E[\sum_{k=1}^{N_{j}}(\Delta W_{kj})^{\prime}\Delta
e_{kj}(\sum_{k=1}^{N_{j}}(\Delta W_{kj})^{\prime}\Delta e_{kj})^{\prime}]$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle
E[\sum_{k=1}^{N_{j}}(\Delta W_{kj})^{\prime}\Delta e_{kj}\Delta e_{kj}(\Delta
W_{kj})]$ $\displaystyle\displaystyle=$ $\displaystyle\displaystyle
E(N_{j})E[(\Delta W_{kj})^{\prime}\Delta e_{kj}\Delta e_{kj}(\Delta W_{kj})].$

Under Assumption E2, $\displaystyle Var(E_{j})$ is finite, because all
variables have up-to the fourth moments. Indeed, if we consider a typical
element of $\displaystyle E[(\Delta W_{kj})^{\prime}\Delta e_{kj}\Delta
e_{kj}(\Delta W_{kj})]$, form by the variables $\displaystyle t$ and
$\displaystyle m$,

$\displaystyle\displaystyle E[(t_{kj}-t_{ij})(m_{kj}-m_{ij})\Delta
e_{kj}(\Delta W_{kj})]$ $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle 4E[t_{kj}m_{kj}(\Delta e_{kj})^{2}]$
$\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle
4E|t_{kj}m_{kj}\Delta e_{kj}^{2}|$ $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle 4\sqrt[4]{E(|t_{kj}|^{4})E[(\Delta
e_{kj})^{2}]E(|m_{kj}|^{4})E[(\Delta e_{kj})^{2}]}$
$\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle M_{0}$

It should be noted that $\displaystyle E[(\Delta
e_{kj})^{2}]=2E(e_{kj}^{2})<\infty$.

Similarly, we can show that $\displaystyle E(\Lambda_{j})=0$, and

$\displaystyle\displaystyle Var(\Lambda_{j})$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle E[(\sum_{k=1}^{N_{j}}(\Delta
W_{kj})^{\prime}\Delta(\lambda(z_{kj}^{\prime}\beta)-\lambda(z_{kj}^{\prime}\hat{\beta})))(\sum_{k=1}^{N_{j}}(\Delta
W_{kj})^{\prime}\Delta(\lambda(z_{kj}^{\prime}\beta)-\lambda(z_{kj}^{\prime}\hat{\beta})))^{\prime}]$
(25) $\displaystyle\displaystyle=$ $\displaystyle\displaystyle
E[\sum_{k=1}^{N_{j}}(\Delta
W_{kj})^{\prime}\Delta(\lambda(z_{kj}^{\prime}\beta)-\lambda(z_{kj}^{\prime}\hat{\beta}))(\Delta
W_{kj})^{\prime}\Delta(\lambda(z_{kj}^{\prime}\beta)-\lambda(z_{kj}^{\prime}\hat{\beta}))^{\prime}]$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle E(N_{j})E((\Delta
W_{kj})^{\prime}\Delta(\lambda(z_{kj}^{\prime}\beta)-\lambda(z_{kj}^{\prime}\hat{\beta}))\Delta(\lambda(z_{kj}^{\prime}\beta)-\lambda(z_{kj}^{\prime}\hat{\beta}))(\Delta
W_{kj}))$ (26) $\displaystyle\displaystyle=$ $\displaystyle\displaystyle
E(N_{j})E[(\Delta W_{kj})^{\prime}\Omega_{kj}(\Delta W_{kj})]$ (27)

We need to show that $\displaystyle E[(\Delta
W_{kj})^{\prime}\Omega_{kj}(\Delta W_{kj})]$, with
$\displaystyle\Omega_{kj}=[\lambda^{\prime}(z_{kj}^{\prime}\beta)]^{2}z_{kj}^{\prime}V_{\beta}z_{kj}$
is finite. A typical element of this matrix is given by, $\displaystyle
E[(t_{kj}-t_{ij})\Omega_{kj}(m_{kj}-m_{ij})]$. We can show the following using
a Cauchy-Schwarz’s inequality.

$\displaystyle\displaystyle E[(t_{kj}-t_{ij})\Omega_{kj}(m_{kj}-m_{ij})]$
$\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle
4E[t_{kj}m_{kj}\Omega_{kj}]$ (28) $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle 4E[|t_{kj}m_{kj}\Omega_{kj}|]$
$\displaystyle\displaystyle\leq$ $\displaystyle\displaystyle
4\sqrt[4]{E(|t_{kj}|^{4})(E([\lambda^{\prime}(z_{kj}^{\prime}\beta)]^{2}z_{kj}^{\prime}V_{\beta}z_{kj}))^{2}E(|m_{kj}|^{4})}$

It remains to be proofed that $\displaystyle
E([\lambda^{\prime}(z_{kj}^{\prime}\beta)]^{2}z_{kj}^{\prime}V_{\beta}z_{kj})<\infty$.
The application of the Cauchy-Schwarz’s inequality implies,

$\displaystyle\displaystyle
E([\lambda^{\prime}(z_{kj}^{\prime}\beta)]^{2}z_{kj}^{\prime}V_{\beta}z_{kj})$
$\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\sqrt{E([\lambda^{\prime}(z_{kj}^{\prime}\beta)]^{4}E[(z_{kj}^{\prime}V_{\beta}z_{kj})^{2}]}$
(30) $\displaystyle\displaystyle\leq$
$\displaystyle\displaystyle\sqrt{E[(z_{kj}^{\prime}V_{\beta}z_{kj})^{2}]}<\infty$
(31)

This follows from noting that $\displaystyle|\lambda^{\prime}(.)|\leq 1$ and
the elements of $\displaystyle z$ have up-to their fourth moments.

The moment of a typical element $\displaystyle
E[(t_{kj}-t_{ij})\Omega_{kj}(m_{kj}-m_{ij})]<\infty.$

This proof that the variance is finite.

It is important to notice that conditional $\displaystyle W$,
$\displaystyle\sum_{i=1}^{N}(\Delta
W_{i})^{\prime}\Delta(\lambda(z_{ij}^{\prime}\beta)-\lambda(z_{ij}^{\prime}\hat{\beta}))$
and $\displaystyle\sum_{i=1}^{N}(\Delta W_{i})^{\prime}\Delta e_{ij}$ are
independent random variables. Therefore,

$\frac{1}{\sqrt{N}}\sum_{i=1}^{N}(\Delta
W_{i})^{\prime}\Delta\eta_{i}\rightarrow^{d}\mathcal{N}(0,\Gamma),$ (32)

where $\displaystyle\Gamma=\rho^{2}E[(\Delta W_{ij})^{\prime}\Omega_{ij}\Delta
W_{ij}]+E[(\Delta W_{ij})^{\prime}\Delta e_{ij}\Delta e_{ij}(\Delta W_{ij})].$

$\sqrt{N}(\hat{\theta}-\theta)\rightarrow^{d}\mathcal{N}(0,\Theta)$ (33)

with $\displaystyle\Theta=C\Gamma C^{\prime}$. This proves the asymptotic
normality of our two step estimator.

We have proven that under assumptions I1, I2, I3, E1 and E2,

$\displaystyle\frac{\sum_{i=1}^{N}(\Delta W_{i})^{\prime}\Delta
W_{i}}{N}\rightarrow^{p}E((\Delta W_{1})^{\prime}\Delta W_{1})=C^{-1}.$

Using similar arguments we can show that

$\displaystyle\frac{\sum_{i=1}^{N}(\Delta
W_{i})^{\prime}\Delta\eta_{i}}{N}\rightarrow^{p}E((\Delta
W_{1})^{\prime}\Delta\eta)=0.$

Which means that $\displaystyle\hat{\theta}$ is a consistent estimator of
$\displaystyle\theta.$ We have proven the estimator is both consistent and
asymptotically normal. This ends the proof of Theorem 2.

## References

  * (1)
  * Ashworth and Heyndels (1997) Ashworth, J., and B. Heyndels (1997): “Politicians’ preferences on local tax rates: an empirical analysis,” _European Journal of Political Economy_ , 13(3), 479–502.
  * Baskaran (2014) Baskaran, T. (2014): “Identifying local tax mimicking with administrative borders and a policy reform,” _Journal of Public Economics_ , 118, 41–51.
  * Black (1999) Black, S. E. (1999): “Do Better Schools Matter? Parental Valuation of Elementary Education,” _The Quarterly Journal of Economics_ , 114(2), 577–599.
  * Buettner and von Schwerin (2016) Buettner, T., and A. von Schwerin (2016): “Yardstick competition and partial coordination: Exploring the empirical distribution of local business tax rates,” _Journal of Economic Behavior & Organization_, 124, 178–201.
  * Cameron, Gelbach, and Miller (2008) Cameron, A. C., J. B. Gelbach, and D. L. Miller (2008): “Bootstrap-based improvements for inference with clustered errors,” _The Review of Economics and Statistics_ , 90(3), 414–427.
  * Charlot and Paty (2007) Charlot, S., and S. Paty (2007): “Market access effect and local tax setting: evidence from French panel data,” _Journal of Economic Geography_ , 7(3), 247–263.
  * Charlot, Paty, and Piguet (2015) Charlot, S., S. Paty, and V. Piguet (2015): “Does fiscal coopération increase local tax rates in urban areas?,” _Regional Studies_ , 49(10), 1706–1721.
  * Charney (1983) Charney, A. H. (1983): “Intraurban manufacturing location decisions and local tax differentials,” _Journal of Urban Economics_ , 14(2), 184–205.
  * Collins, Alva, Boydston, Cochran, Hamm, McGuire, and Riga (2006) Collins, H., A. Alva, R. Boydston, R. Cochran, P. Hamm, A. McGuire, and E. Riga (2006): “Soil microbial, fungal, and nematode responses to soil fumigation and cover crops under potato production,” _Biology and Fertility of Soils_ , 42(3), 247–257.
  * Crowley and Sobel (2011) Crowley, G. R., and R. S. Sobel (2011): “Does fiscal decentralization constrain Leviathan? New evidence from local property tax competition,” _Public Choice_ , 149(1-2), 5.
  * Davidson and MacKinnon (2010) Davidson, R., and J. G. MacKinnon (2010): “Wild bootstrap tests for IV regression,” _Journal of Business & Economic Statistics_, 28(1), 128–144.
  * Djogbenou, MacKinnon, and Nielsen (2019) Djogbenou, A. A., J. G. MacKinnon, and M. Ø. Nielsen (2019): “Asymptotic theory and wild bootstrap inference with clustered errors,” _Journal of Econometrics_ , 212(2), 393–412.
  * Duranton, Gobillon, and Overman (2011) Duranton, G., L. Gobillon, and H. G. Overman (2011): “Assessing the Effects of Local Taxation using Microgeographic Data,” _The Economic Journal_ , 121(555), 1017–1046.
  * Dustmann and Rochina-Barrachina (2007) Dustmann, C., and M. E. Rochina-Barrachina (2007): “Selection correction in panel data models: An application to the estimation of females’ wage equations,” _The Econometrics Journal_ , 10(2), 263–293.
  * Gibbons and Machin (2003) Gibbons, S., and S. Machin (2003): “Valuing English primary schools,” _Journal of Urban Economics_ , 53(2), 197–219.
  * Hansen and Lee (2019) Hansen, B. E., and S. Lee (2019): “Asymptotic theory for clustered samples,” _Journal of econometrics_ , 210(2), 268–290.
  * Heckman (1974) Heckman, J. (1974): “Shadow Prices, Market Wages, and Labor Supply,” _Econometrica_ , 42(4), 679–694.
  * Heckman (1979) Heckman, J. J. (1979): “Sample Selection Bias as a Specification Error,” _Econometrica_ , 47(1), 153–161.
  * Holmes (1998) Holmes, T. J. (1998): “The Effect of State Policies on the Location of Manufacturing: Evidence from State Borders,” _Journal of Political Economy_ , 106(4), 667–705.
  * Kyriazidou (1997) Kyriazidou, E. (1997): “Estimation of a Panel Data Sample Selection Model,” _Econometrica_ , 65(6), 1335–1364.
  * Lee (2001) Lee, M.-j. (2001): “First-difference estimator for panel censored-selection models,” _Economics Letters_ , 70(1), 43–49.
  * MacKinnon (2013) MacKinnon, J. G. (2013): “Thirty years of heteroskedasticity-robust inference,” in _Recent advances and future directions in causality, prediction, and specification analysis_ , pp. 437–461. Springer.
  * MacKinnon and Webb (2017) MacKinnon, J. G., and M. D. Webb (2017): “Wild bootstrap inference for wildly different cluster sizes,” _Journal of Applied Econometrics_ , 32(2), 233–254.
  * Oster (2019) Oster, E. (2019): “Unobservable selection and coefficient stability: Theory and evidence,” _Journal of Business & Economic Statistics_, 37(2), 187–204.
  * Rochina-Barrachina (1999) Rochina-Barrachina, M. E. (1999): “A new estimator for panel data sample selection models,” _Annales d’Economie et de Statistique_ , pp. 153–181.
  * Roodman, Nielsen, MacKinnon, and Webb (2019) Roodman, D., M. Ø. Nielsen, J. G. MacKinnon, and M. D. Webb (2019): “Fast and wild: Bootstrap inference in Stata using boottest,” _The Stata Journal_ , 19(1), 4–60.
  * Ross and Yinger (1999) Ross, S., and J. Yinger (1999): “Sorting and voting: A review of the literature on urban public finance,” _Handbook of regional and urban economics_ , 3, 2001–2060.
  * Todd and Wolpin (2003) Todd, P. E., and K. I. Wolpin (2003): “On the specification and estimation of the production function for cognitive achievement,” _The Economic Journal_ , 113(485), F3–F33.
  * White (1984) White, H. (1984): _Asymptotic theory for econometricians_. Academic press.
  * Wooldridge (1995) Wooldridge, J. M. (1995): “Selection corrections for panel data models under conditional mean independence assumptions,” _Journal of econometrics_ , 68(1), 115–132.
  * Wooldridge (2010) (2010): _Econometric analysis of cross section and panel data_. MIT press.