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.
117,487 characters · 22 sections · 43 citation commands
Using generalized estimating equations to estimate nonlinear models with spatial data
In empirical economic and social studies, there are many examples of discrete data which exhibit spatial or cross-sectional correlations possibly due to the closeness of geographical locations of individuals or agents. One example is the technology spillover effect. The number of patents a firm received shows correlation with that received by other nearby firms (E.g. bloom2013identifying). Another example is the neighborhood effect. There is a causal effect between the individual decision whether to own stocks and the average stock market participation of the individual's community (E.g.brown2008neighbors). These two examples involves dealing with discrete data. The first example is concerned with count data and the second one handles binary response data. Nonlinear models are more appropriate than linear models for discrete response data. With spatial correlation, these discrete variables are no longer independent. Both the nonlinearity and the spatial correlation make the estimation difficult.
In order to estimate nonlinear models, one way is to use maximum likelihood estimation (MLE). A full MLE specifies the joint distribution of spatial random variables. This includes correctly specifying the marginal and the conditional distributions, which impose very strong assumptions on the data generating processes. However, given a spatial data set, the dependence structure is generally unknown. If the joint distribution of the variables is misspecified, MLE is in general not consistent. One of the alternative MLE method is partial-maximum likelihood estimation (PMLE), which only uses marginal distributions. wang2013partial use a bivariate Probit partial MLE to improve the estimation efficiency with a spatial Probit model. Their approach requires to correctly specify the marginal distribution of the binary response variable conditional on the covariates and distance measures\footnote{ A sample of spatial data is collected with a set of geographical locations. Spatial dependence is usually characterized by distances between observations. A distance measure is how one defines the distances between observations. Physical distance or economic distance could be two options. Information about agents locations is commonly imprecise, e.g. only zip code is known. conley2007spatial deals with the inference problem when there exist distance errors. In this paper we assume there are no measurement errors in pairwise distances.} There are two concerns with wang2013partial. First the computation is already hard for a bivariate distribution. The multivariate marginal distribution of {a higher dimensional variable}, e.g., trivariate, is more computationally demanding; second it also requires the correct specification of the marginal bivariate distribution to obtain consistency. The bivariate marginal distribution of a spatial multivariate normal distribution is bivariate normal, thus the bivariate Probit model can be derived. But there are other distributions whose marginal distribution is not the same anymore. For example, the marginal distribution of a multivariate Logit is not logistic. If the partial likelihood is misspecified, the estimation of the mean parameters could be not consistent. With less distributional assumptions, the quasi-maximum likelihood estimation (QMLE) can also be used to estimate nonlinear models. Using a density that belongs to a linear exponential family (LEF), QMLE is consistent if we correctly specify the conditional mean while other features of the density can be misspecified (gourieroux1984pseudo). lee2004asymptotic derives asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models by allowing not assuming normal distributions. In a panel data case, pooled or partial QMLE (PQMLE) which ignores serial correlations is consistent under some regularity conditions (wooldridge2010econometric).
We further relax distributional assumptions than those required in bivariate partial MLE as in wang2013partial. Suppose we only assume correct mean function and one working variance covariance matrix \footnote{The true variance covariance matrix is generally unknown. By specifying a working variance covariance matrix, one can capture some of the correlation structure between observations.} which may not be correct. Using QMLE in the LEF, we can consistently estimate the mean parameters as well as the average partial effects. The generalized estimating equations (GEE) approach is one of the QMLE methods. It is used in panel data models to account for serial correlation and thus get more efficient estimators. A generalized estimating equation is used to estimate the parameters of a generalized linear model with a possible unknown correlation between outcomes (liang1986longitudinal). Parameter estimates from the GEE are consistent even when the variance and covariance structure is misspecified under mild regularity conditions. This is quite related to a different terminology, composite likelihood. varin2011overview provide a survey of developments in the theory and application of composite likelihood. The motivation for the use of composite likelihood is usually computational, to avoid computing or modelling the joint distributions of {high dimensional} random processes. {One can find many related reference in the literature, such as bhat2010comparison.} As a special case of composite likelihood methods, one way is to use partial conditional distribution, and maximize the summand of log likelihoods for each observation. It assumes a working independence assumption, which means that the estimators are solved by ignoring dependence between individual likelihoods. The parameters can be consistently estimated if the partial log likelihood function satisfies certain regularity assumptions. However, a consistent variance estimator should be provided for valid inference\footnote{ Ignoring dependence in the estimation of parameters will result in wrong inferences if the variances are calculated in the way that independence is assumed. Dependence should be accounted for to the extent of how much one ignores it in the estimation.}. When there exists spatial correlation, the pooled maximum likelihoods (composite likelihoods) can be considered as misspecified likelihoods because of the independence assumption.
Generalized least squares (GLS) could be used to improve the estimation efficiency in a linear regression model even if the variance covariance structure is misspecified. { lu2017quasi propose a quasi-GLS method to estimate the linear regression model with an spatial error component. By first estimating the spatial parameter in the error variance and then using estimated variance matrix for within group observations, the quasi-GLS is computationally easier and would not loose much efficiency compared to GLS. } Similarly, the multivariate nonlinear weighted least squares estimator (MNWLS), {see Chapter 12.9.2 in wooldridge2010econometric}, is essentially a GLS approach applied in nonlinear models to improve the estimation efficiency.
{It is worth noting that the GEE approach discussed in this paper is a two-step method, which is essentially a special MNWLS estimator that uses a LEF variance assumption and a possibly misspecified working correlation matrix in the estimation. The GEE approach was first extended to correlated data by liang1986longitudinal, which propose a fully iterated GEE estimator in a panel data setting. In addition, zeger1986longitudinal fit the GEE method to discrete dependent variables. The iterated GEE method has solutions which are consistent and asymptotically Gaussian even when the temporal dependence is misspecified. The consistency of mean parameters only depends on the correct specification of the mean, not on the choice of working correlation matrix. GEE used in nonlinear panel data models and system of equations is supposed to obtain more efficient conditional mean parameters with covariance matrix accounting for the dependency structure of the data. In this paper, we apply a similar idea to grouped spatial data. We use the PQMLE as the initial estimator for the two-step GEE and study the efficiency properties of a two-step GEE estimator and expect that GEE can give more efficient estimators compared to PQMLE.}
Moreover, we demonstrate theoretically how to use our GEE approach within the QMLE framework in a spatial data setting to obtain consistent estimators. We give a series of assumptions, based on which QMLE estimators are consistent for the spatial processes. To derive the asymptotics for the GEE estimator we have to use a uniform law of large numbers (ULLN) and a central limit theorem (CLT) for spatial data. These limit theorems are the fundamental building blocks for the asymptotic theory of nonlinear spatial M-estimators, for example, maximum likelihood estimators (MLE) and generalized method of moments estimators (GMM) (Jenish2012). conley1999gmm makes an important contribution toward developing an asymptotic theory of GMM estimators for spatial processes. He utilizes bolthausen1982central CLT for stationary random fields. Jenish2009,Jenish2012 provide ULLNs and a CLTs for near-epoch dependent spatial processes. Using theorems in Jenish2009, Jenish2012, one can analyze more interesting economic phenomena. It should be noted that although GEE can be considered as a special case of M-estimation, we have carefully checked how the near-epoch dependence property of the underlying processes is translated to our responses and the partial sum processes involved in proving the asymptotics of the estimation. Our setup is different from the literature as it is with a grouped estimation structure. Finally, we have provided a consistency proof of the proposed semiparametric estimator of the variance covariance matrix.
{We contribute to the literature in three aspects. First, we propose a simple method which uses less distributional assumptions by only specifying the conditional mean for spatial dependent data. The method is computationally easier by dividing data into small groups compared to using all information. We model the spatial correlation as a moving average (MA) type in the underlying innovations instead of the spatial autoregressive (SAR) model in the dependent variable. Second, we proved the theoretical property of our estimator by applying ULLN and CLT in Jenish2009,Jenish2012 to the GEE estimator with careful checking the hyper assumptions. Third, we emphasize the possible efficiency gain from making use of spatial correlation from our simulation study, and we demonstrate how to use GEE with two types of data: count and binary response. }
In Section (ref), the GEE methodology in a QMLE framework under the spatial data context is proposed. In Section (ref), we look in detail at a Poisson model and Negative Binomial II model for count data with a multiplicative spatial error term. We further study a Probit model for binary response data with spatial correlation in the latent error term. In Section (ref), a series of assumptions are given based on Jenish2009, Jenish2012 under which GEE-estimators are consistent and have an asymptotic normal distribution. The asymptotic distributions for GEE for spatial data are derived. Consistent variance covariance estimators are provided for the nonlinear estimators. Section (ref) contains Monte Carlo simulation results which compare efficiency of different estimation methods for the nonlinear models explored in the previous section. Section (ref) contains an application to study the determinants of the inflow FDI to China using city level data. The technical details are delegated to Section (ref).
Unlike linear models, a very important feature of nonlinear models is that estimators cannot be obtained in a closed form, which requires new tools for asymptotic analysis: uniform law of large numbers (ULLN) and a central limit theorem (CLT). Jenish2009 develop ULLN and CLT for $\alpha $-mixing random fields on unevenly spaced lattices that allow for nonstationary processes with trending moments. But the mixing property can fail for quite a few reasons, thus we adopt the notion of near-epoch dependence (NED) as in Jenish2012 which refers to a generalized class of random fields that is "closed with respect to infinite transformations." We consider spatial processes located on a unevenly spaced lattice $D \subseteq \mathbb{R}^{d}, d \geq 1$. The space $\mathbb{R}^{d}$ is endowed with the metric $\rho(i, j) = max_{1 \leq l \leq d}|j_{l} - i_{l}|$ with the corresponding norm $|i|_{\infty} = max_{1 \leq l \leq d}|i_{l}|$, where $i_{l}$ is the $l$-th component of $i$. The distance between any subsets $U, V \in D$ is defined as $\rho(U, V) = \inf \{\rho(i, j): i\in U \text{ and } j\in V \}$. Further, let $|U|$ denote the cardinality of a finite subset $ U\subseteq D$. The setting is illustrated in Jenish2009,Jenish2012.
Let $Z = \{Z_{n, i}, i \in D_{n}, n \geq 1\}$ and $\varepsilon = \{\varepsilon_{n, i}, i \in T_{n}, n \geq 1\}$ be triangular arrays of random fields defined on a probability space $(\Omega, \mathscr{F}, P)$ with $D_{n}\subseteq T_{n} \subseteq D$ where $D$ satisfies A.1). The cardinality of $D_{n}$ and $T_{n}$ satisfy $\displaystyle \lim_{n \rightarrow \infty} |D_{n}|\rightarrow \infty, \displaystyle \lim_{n \rightarrow \infty} |T_{n}|\rightarrow \infty$. For any vector $v \in R^p$, $|v|_2$ denotes the $L_2$ norm of $v$. For any $n \times m$ matrix $A$ with element $a_{ij}$, denote $|A|_1 = \displaystyle \max_{1\leq j\leq m} \displaystyle \sum^{n}_{i=1} |a_{ij}|$ and $|A|_{\infty} = \displaystyle \max_{1\leq i\leq n} \displaystyle \sum^{m}_{j=1} |a_{ij}|$, $|A|_2$ denotes the 2-norm. For any random vector $X$, denote $\| X_{n,i} \|_{p} = (\mathop{\mbox{\sf E}}|X_{n,i}|^{p})^{1/p}$ as its $L_{p}$-norm, where the absolute $p$th moment exists. We brief $\| X_{n,i} \|_{2} $ as $\| X_{n,i} \|$. Let $\mathcal{F}_{n,i}(s) = \sigma(\varepsilon_{n,j}: j \in D_{n}, \rho(i, j) \leq s)$ as the $\sigma$- field generated by random vectors $\varepsilon_{n, j}$ located within distance $s$ from $i$. Given two sequences of positive numbers $x_n$ and $y_n$, write $x_n\lesssim y_n$ if there exists constant $C>0$ such that $x_n/y_n\leq C$, also we can write $x_n = {\mathcal{O}}(y_n)$. A sequence $x_n$ is said to be ${\scriptstyle{\mathcal{O}}}(y_n)$ if $x_n/y_n \to 0,$ as $n \to \infty$. In a similar manner, The notation, $X_n={\mathcal{O}}_p(a_n)$ means that the set of values $X_n/a_n$ is stochastically bounded. That is, for any $ \varepsilon > 0$, there exists a finite M > 0 and a finite N > 0 such that, $P(|X_n/a_n| > M) < \varepsilon, \forall n > N$. $|.|_a$ is the elementwise absolute value of a matrix $|A|_a$. $a\vee b $ is $ \max(a,b).$
We will present the $L_{2}$-NED properties of a random field $Z$ on some $\alpha$-mixing random field $\varepsilon$. The definition of the $\alpha$-mixing coefficient employed in the paper are stated as following.
Note that we suppress the dependence on $n$ from now on for the triangular array. Let $\left \{ \left(\mathbf{x}_{i},y_{i}\right) ,i=1,2,...,n\right \} $, where $\left( \mathbf{x}_{i},y_{i}\right) $ is the observation at location $s_{i}.$ $ \mathbf{x}_{i}$ is a row vector of independent variables which can be continuous, discrete or a combination. The dependent variable $y_{i}$ can be continuous or discrete. Let $\left( \mathbf{x}_{g},\mathbf{y} _{g}\right) $ be the observations in group $g$ and $B_{g}$ is the associated set of locations within the group $g$. We will focus on the case of a discrete dependent variable, a binary response and a count. Let $\theta \in \mathbf{R}^p$, $\gamma \in \mathbf{R}^q$ and $\mathbf{\theta \in \Theta, \gamma \in \Gamma}$, where $\mathbf{\Theta\times \Gamma} $ is a compact set, and $(\theta^0, \gamma^0)$ is the true parameter value.
{The GEE methodology proposed in equations (6) and (7) in liang1986longitudinal is an iterated approach to estimate the mean parameters. We simplify the procedure using a two-step method by first estimate the working correlation matrix and then apply MWNLS.} In the following, we write the GEE methodology in the group level notation. Groups are divided according to geographical properties or other researcher defined economic (social) relationships. Our asymptotic analysis is based on large number of groups $g = 1, \cdots, G$. The notation $D_{G}$ indicates the lattice containing group locations, each group location is denoted as vectorizing the elements in $B_g$. Let the total number of groups be $|D_G| = G,$ while the total number of observations is still $|D_n| = n.$ Let $L_{g}$ be the number of observations in group $g$. For simplicity assume $L_{g}=L,$ for all $g$. Let $\left \{ \left( \mathbf{x }_{g},\mathbf{y}_{g}\right) \right \} $ be the observations for group $g$, where $\mathbf{x}_{g}$ is an $L\times p$ matrix and $\mathbf{y}_{g}$ is an $ L\times 1$ vector$.$ There are two extreme cases of the group size. The first case is when the group size is $1$, the resulting estimator is the usual PQMLE estimator, which means we ignore all of the pairwise correlations. The second case is when the group size is $n$, which means we are using all the pairwise information. If the group size is not equal to $1$ or $n$, the estimation is actually a "partial" QMLE. By "partial", we mean that we do not use full information, but only the information within the same groups. Note that we work with the case with number of groups $G\to \infty$ in our theory, { while the groupsize $L$ is assumed to be fixed.}
Assume that we correctly specify conditional mean of $\mathbf{y}_{g},$ that is, the expectation of $\mathbf{y}_{g}$ conditional on $\mathbf{x}_{g}$ is
Assume the conditional variance-covariance matrix of $\mathbf{y}_{g}$ is $\mathbf{W}^*_{g}$ which is unknown in most cases, where $\mathbf{W}_g \stackrel{\mathrm{def}}{=} \mathop{\mbox{Cov}}(y_g,y_g|\mathbf x_g)= \mathop{\mbox{\sf E}}(\mathbf{y}_{g}\mathbf{y}_{g}^{\top}|\mathbf x_g) - \mathop{\mbox{\sf E}}(\mathbf{y}_{g}|\mathbf x_g)\mathop{\mbox{\sf E}}(\mathbf{y}_{g}|\mathbf x_g)^{\top}$. Usually we parameterize a corresponding weight matrix $\mathbf{W}_g$ by $\mathbf{W}_g( \theta,\gamma)$, where $\theta \in \Theta \subset \mathbf{R}^q$ and $\gamma \in \Gamma \subset \mathbf{R}^p$ as a nuisance parameter involved only in the estimation of the variance covariance matrix. {In practice, we usually preestimate $\gamma$ and thus it is replaced by a consistent estimate of $\hat{\gamma}$, then $\mathbf W_g$ is denoted as $\mathbf{W}( \theta, \hat{\gamma})$.}
The objective function for group $g$ and the whole sample are given as follows:
where {$M_{G}$ is a scaling constant defined in A.5) in section (ref).}
Theoretically, an estimator of $\theta^0, \gamma^0$ is given by
In practice a GEE estimator is obtained by a two-step procedure, where the first step is to estimate the nuisance parameter $\gamma$ and the second step is to have the parameter $\theta$ estimated with the plug-in estimator $\hat{\gamma}$ from step 1.
{Because this only uses the groupwise information, it actually is a "quasi" or "pseudo" MWNLS. }The quasi-score equation, which is the first order condition for GEE, is defined as follows:
where $\nabla _{\theta }\mathbf{m}_{g}\left( \mathbf{\theta }\right) $ is the gradient of $\mathbf{m}_{g}\left( \mathbf{\theta }\right) .$ $M_G$ is defined as the scaling constant in A.5) in section (ref). The GEE estimator $(\mathbf{\hat{\theta},\hat{\gamma}}) = \mbox{argzero}_{\theta \in \Theta, \gamma \in \Gamma} \mathbf{S}_{G}\left( \theta, \mathbf{{\gamma}}\right) .$
Denote the population version of loss as $\mathbf{S}_{\infty}\left( \theta,\gamma\right) = \mbox{lim}_{G\to \infty} \mathop{\mbox{\sf E}} \mathbf{S}_{G}\left(\theta,\gamma\right),$ and \\$Q_{\infty}(\theta,\gamma) = \lim_{G\to\infty} {(G M_G)}^{-1}\sum_g \mathop{\sf E} q_{g}(\theta,\gamma) .$ Thus the true parameter $( \theta^0,\gamma^0) \\= argzero_{\theta \in \Theta, \gamma \in \Gamma } \mathbf{S}_{\infty}\left( \mathbf{\theta}, \mathbf{\gamma} \right) = argmin_{ \theta \in \Theta, \gamma \in \Gamma }Q_{\infty}( \theta,\gamma).$
Frequently we restrict our attention to the exponential family, which embraces many frequency encountered distributions, such as Bernoulli, Poisson and Gaussian, etc.
Now we write this estimation in a QMLE framework. We suppress the parameter $\gamma$ for a moment. Assume the probability density function $f\left( \mathbf{y}_g|\mathbf{x}_g;\mathbf{\theta }\right) $ is in the LEF.({ See details in Appendix (ref) }.)
Without accounting for the spatial covariance, one characterization of QMLE in LEF is that the individual score function has the following form:
where $\nabla m_{i}\left( \mathbf{x}_{i};\mathbf{\theta } \right) $ is the $1\times p$ gradient of the mean function and $v_{i}\left( m_{i}\left( \mathbf{x}_{i},{D}_{n};\mathbf{\theta }\right) \right) $ is the {conditional} variance function associated with the chosen LEF density. For Bernoulli distribution, $v_{i}\left( m_{i}\left( \mathbf{x}_{i};\mathbf{ \theta }\right) \right) =m_{i}\left( \mathbf{x}_{i};\mathbf{\theta }\right) \left( 1-m_{i}\left( \mathbf{x}_{i};\mathbf{\theta }\right) \right) ,$ and for Poisson distribution, $v_{i}\left( m_{i}\left( \mathbf{x}_{i};\mathbf{ \theta }\right) \right) =m_{i}\left( \mathbf{x}_{i};\mathbf{\theta }\right) . $ Note that ((ref)) gives a consistent estimator but is not likely to be the most efficient estimator as it ignores the possible spatial correlations between observations. However, it accounts for possible heteroscedasticity.
We write the quasi-score function for a group. Let $\mathbf{v} _{g}\left( \mathbf{m}_{g}\left( \mathbf{x}_{g};\mathbf{\theta }\right) \right) $ be the conditional variance covariance matrix for group $g$. Then score involved in the estimation is denoted as
where
We specify a more general form of variance $\mathbf{v}_{g}\left( \theta \right)$ with the dependency of the nuisance parameter $\gamma$. The conditional mean vector is correctly specified for each individual E$\left( y_{i}|\mathbf{x}_{i}\right) =m_{i}\left( \mathbf{x}_{i}; \mathbf{\theta }^{0}\right) .$ Thus for each group, $ \mathbf{m}_{g}\left( \mathbf{x}_{g};\mathbf{\theta }^0\right) =\mathop{\mbox{\sf E}} \left( \mathbf{y}_{g}|\mathbf{x}_{g}\right) .$ Let ${s}_{g}\left( \theta,\mathbf{{\gamma}}\right) $ denote the $p\times 1$ vector of score for group $g$. Let $h_{g}\left( \mathbf{ \theta, {\gamma}}\right) $ be the $p\times p$ matrix of Hessian for group $g$. The score function for $Q_{G}\left( \mathbf{\theta,{\gamma} }\right) $ can be defined as $\mathbf{S}_{G}\left( \theta, {\gamma}\right) $ and the Hessian can be defined as $\mathbf{H}_{G}\left( \mathbf{ \theta, {\gamma} }\right) .$ The score function for GEE can be written as
and the Hessian is
where $\mbox{Vec}$ is denoted as the vectorization of a matrix $A$.
In this subsection, we demonstrate one way to find an estimator for $\gamma$ involved in $\mathbf{W}_{g}(\theta, \gamma).$ $\mathbf{W}_{g}(\theta,\gamma)$ can be written as
where $\mathbf{V}_{g}$ is the $L\times L$ diagonal matrix that only contains variances of $\mathbf{y}_g - \mathbf{m}_g(\mathbf x_g, \theta^0)$ and $\mathbf{R}_{g}$ is the $L\times L$ correlation matrix for group $g$.
Let
where the $l$th element on the diagonal is $v_{gl}=\mathrm{Var}(\mathbf{y} _{gl}\mathbf{|x}_{gl})$ in group $g,$ $\mathbf{y}_{gl}$ is the $l$th element in the vector $\mathbf{y}_{g}$ and $\mathbf{x}_{gl}$ is the $l$th row in $ \mathbf{x}_{g}$. And
Let $d_{glm}$ be the distance between the $l$th and the $m$th observations in group $g$. An example of a parametrization of the correlation i.e. the $l, m$th, $l\neq m,$ element of $\mathbf{R}_{g},$ as in cressie1992statistics is
where the spatial correlation parameters $\mathbf{\gamma =}\left( b,c,\rho \right) ,$ $b\geq 0,c\geq 0,\rho \geq 0,$ and $b+c\leq 2.\footnote{ See \cite{cressie1992statistics} p.61 for more examples.}$Set $b=c=1$ without loss of generality. Then
Although the above specification does not represent all the possibilities, it at least provides a way of how to parameterize the spatial correlation, and therefore the basis for testing spatial correlation.
The following provides a way to estimate $\mathbf{\gamma }$. Let $\mathbf{ \check{\theta}}$ be the first-step PQMLE estimator. $\check{u} _{i}=y_{i}-m_{i}\left( x_{i};\mathbf{\check{\theta}}\right) $ are the first-step residuals. $\check{v}_{i}=v\left( m_{i}\left( \mathbf{x}_{i}; \mathbf{\check{\theta}}\right) \right) $ is the fitted variance of individual $i$ corresponding to the chosen LEF density. Let $\check{r }_{i}=\check{u}_{i}/\sqrt{\check{v}_{i}}$ be the standardized residual. Let $ \mathbf{\check{r}}_{g}=$ $\left( \check{r}_{g1},\check{r}_{g2},...,\check{r} _{gL}\right) ^{\top }.$ Then $\mathbf{\mathbf{\check{r}}_{g}\mathbf{\check{r }}_{g}}^{\top }$ is the estimated sample correlation matrix for group $g$. Let $\mathbf{e}_{g}(\check{\theta})$ be a vector containing $L(L-1)/2$ different elements of the lower (or upper) triangle of $\mathbf{\mathbf{\check{r}}_{g}\mathbf{ \check{r}}_{g}}^{\top },$ excluding the diagonal elements. Let $\mathbf{z} _{g}(\gamma)$ be the vector containing the elements in $\mathbf{R}_{g}$ corresponding to the same entries of elements in $\mathbf{\mathbf{\check{r}} _{g}\mathbf{\check{r}}_{g}}^{\top }$. We can follow prentice1988applications, who provides one way to find a consistent estimator for $\mathbf{\gamma }$ by solving:
The setup of nonlinear models with spatial data varies with different models. For each model, we need to incorporate the spatial correlated term in an appropriate way. In this Section, we will demonstrate how we incorporate the spatial correlated error term in two types of discrete data and how to use a GEE procedure to estimate the nonlinear models. The first example is for count data and the second one is for binary response data.
A count variable is a variable that takes on nonnegative integer values, such as the number of patents applied for by a firm during a year. bloom2013identifying studies spillover effects of R&D between firms in terms of firm patents. Other examples include the number of times someone being arrested during a given year. Count data examples with upper bound include the number of children in a family who are high school graduates, in which the upper bound is number of children in the family (wooldridge2010econometric).
We first model the count data with a conditional Poisson density, $f\left( y| \mathbf{x}\right) =\exp \left[ -\mu \right] \mu ^{y}/y!,$ where $y!=1\cdot 2\cdot ...\cdot \left( y-1\right) \cdot y$ and $0!=1.$ $\mu $ is the conditional mean of $y.$ The Poisson QMLE requires us only to correctly specify the conditional mean. A default assumption for the Poisson distribution is that the mean is equal to the variance. Note that even if $y_{i}$ does not follow the Poisson distribution, the QMLE approach will give a consistent estimator if you use the Poisson density function and a correctly specified conditional mean (gourieroux1984pseudo). Moreover, $ y_{i}$ even need not to be a count variable. The most common mean function in applications is the exponential form:
When spatial correlation exists, we can characterize count data model with a multiplicative spatial error. silva2006log use the Poisson pseudo-maximum-likelihood (PPML), which is the Poisson QMLE in this paper, to estimate the gravity model for trade. They argue that constant elasticity models should be estimated in their multiplicative form, because using a log linear model can cause bias in coefficient estimates under heteroskedasticity. Now we further consider the Poisson regression model with spatial correlation in the multiplicative error,
where $v_{i}$ is the multiplicative spatial error term. Let $\mathbf{v}$ equal $\left( v_{1},v_{2},...,v_{n}\right) ^{\top }.$ (Note that for this example we treat location $i$ as an one dimensional object.) This model is characterized by the following assumptions:
(1) $\{(\mathbf{x}_{i},v_{i}),i=1,2,...,n\}$ is a mixing sequence on the sampling space $D_n$, with mixing coefficient $\alpha $.
(2) $\mathop{\mbox{\sf E}}\left( y_{i}|\mathbf{x}_{i},v_{i}\right) =v_{i}\exp \left( \mathbf{x}_{i}\mathbf{\beta }_{0}\right) .$
(3) $y_{i},y_{j}$ are independent conditional on $\mathbf{x}_{i},\mathbf{x} _{j},v_{i},v_{j},i\neq j.$
(4) $v_{i}$ has a conditional multivariate distribution, $\mathop{\mbox{\sf E}}\left( v_{i}|\mathbf{x}_{i}\right) =1$. $\mathrm{Var}\left( v_{i}|\mathbf{x} _{i}\right) =\tau ^{2},$ $\mathrm{Cov}\left( v_{i},v_{j}|\mathbf{x}_{i}, \mathbf{x}_{j}\right) =\tau ^{2}\cdot c\left( d_{ij},\rho \right) ,$ where $ c\left( d_{ij},\rho \right) $ is the correlation function of $v_{i}$ and $ v_{j}.$
Under the above assumptions, and again conditional on $D_n$ is suppressed, we can integrate out $v_{i}$ by using the law of iterated expectations.
If $x_{j}$ is continuous, the partial effects on $\mathop{\mbox{\sf E}}\left( y_{i}| \mathbf{x}_{i},D_n\right) $ is $\exp \left( \mathbf{x}_{i}\mathbf{ \beta }_{0}\right) \beta _{j}.$ If $x_{j}$ is discrete the partial effects is the change in $\mathop{\mbox{\sf E}}\left( y_{i}|\mathbf{x}_{i},\mathbf{D} _{n}\right) $ when, say, $x_{K}$ goes from $a_{K}$ to $a_{K}+1$ which is
The pooled QMLE gives a consistent estimator for the mean parameters, which solves:
Its score function is
Since this estimator does not account for any heteroskedasticity or spatial correlation, a robust estimator for the asymptotic variance of partial QMLE estimator is provided as follows,
where $k\left( d_{ij}\right) $ is a kernel function depending on the distance between observations $i$ and $j$.
Moreover, a very specific nature of the Poisson distribution is that we can write down the conditional variances and covariances of $y$:
The conditional variance of $y_{i}$ given $\mathbf{x}_{i}$ is a function of both the level and the quadratic of the conditional mean. The traditional Poisson variance assumption is that the conditional variance should equal the conditional mean. That is, $\mathrm{Var}\left( y_{i}|\mathbf{x}_{i}\right)=\exp \left( \mathbf{x}_{i}\mathbf{\beta }_{0}\right) .$ The Poisson GLM variance assumption is $\mathrm{Var}\left( y_{i}| \mathbf{x}_{i}\right) =\sigma ^{2}\exp \left( \mathbf{x}_{i}\mathbf{\beta } _{0}\right) $ with an overdispersion or underdispersion parameter $\sigma ^{2}$, which is a constant. Obviously, there is over-dispersion in ((ref)) since $\exp \left( 2\mathbf{x}_{i}\mathbf{\beta }_{0}\right) \cdot \tau ^{2}\geq 0,$ and the over-dispersion parameter is $1+\exp \left( \mathbf{x}_{i}\mathbf{\beta }_{0}\right) \cdot \tau ^{2}$, which is changing with $\mathbf{x}_{i}$. This does not coincide with Poisson variance assumption and the GLM variance assumption. What is more, the conditional covariances can be written in the following form,
In the group level notation,
Let $\mathbf{W}_{g}$ be the variance-covariance matrix for group $g$ evaluated at the true value $\beta_0,\rho_0$. The variance of the $l$th element in group $g$ is
and the covariance of the $l$th and $m$th elements in group $g$ is
Here $\mathbf{\gamma =}\left( \tau ^{2},\rho \right) ^{\top }$ and $ \mathbf{\hat{\gamma}=}\left( \hat{\tau}^{2},\hat{\rho}\right) ^{\top }$ is an estimator for $\mathbf{\gamma }$. Let $\mathbf{\check{\beta}}_{\mathrm{ PQMLE}}$ be the partial QMLE estimator in the first step. Then the elements in $\mathbf{W}_{g}$ can be estimated as
Based on the conditional distribution, the first order conditions for GEE is:
$\mathbf{\hat{\beta}}_{\mathrm{GEE}}$ is consistent and follows a normal distribution asymptotically by Theorem (ref) and (ref). We will brief $\mathbf{W}_{g}^{-1}\left( \mathbf{\hat{ \gamma},\hat{\theta}}\right)$ as $\hat{\mathbf{W}}_{g}^{-1}$ in the following text. The variance estimator for the asymptotic variance that is robust to misspecification of spatial correlation is:
where $k(d_{gh})$ is a kernel function depending on the distances between groups. The distances could be the smallest distance between two observations belonging to different groups.
The pivotal parameters, $\tau ^{2}$ and $\rho ,$ can be estimated using the Poisson QMLE residuals. Let $\check{u}_{i}^{2}=\left[ y_{i}-\exp \left( \mathbf{x}_{i}\mathbf{\check{\beta}}_{\mathrm{QMLE}}\right) \right] ^{2}$ be the squared residuals from the Poisson QMLE. Based on equation ((ref)), $\tau ^{2}$ can be estimated as the coefficient by regressing $ \check{u}_{i}^{2}-\exp \left( \mathbf{x}_{i}\mathbf{\check{\beta}}_{\mathrm{ QMLE}}\right) $ on $\exp \left( 2\mathbf{x}_{i}\mathbf{\check{\beta}}_{ \mathrm{QMLE}}\right) .$ The situation to estimate $\rho $ depends on the specific form of $c\left( d_{ij},\rho \right) $. We would like to assume a structure, though it might be wrong, to approximate the true covariance. For example, suppose the covariance structure of $e_{i}$ and $e_{j}$ is $\exp \left( \frac{\rho }{d_{ij}}\right) -1,$ and the correlation structure is $ c\left( d_{ij},\rho \right) =\frac{\exp \left( \frac{\rho }{d_{ij}}\right) -1 }{\mathrm{e}-1},$ then an estimator for $\rho $ is:
Then $\mathbf{\hat{W}}_{g}$ is obtained by plugging $\hat{\tau}^{2}$ and $ \hat{\rho}$ back in the variance-covariance matrix. We can also directly calculate $\hat{\rho}$ as
Since the conditional variances and covariances can be written in a specific form, we would consider NegBin II model of cameron1986econometric as a more appropriate model. The NegBin II model can be derived from a model of multiplicative error in a Poisson model. With an exponential mean, $y_{i}| \mathbf{x}_{i},v_{i},D_n\sim $Poisson$\left[ v_{i}\exp \left( \mathbf{x}_{i}\mathbf{\beta }_{0}\right) \right] $. Under the above assumptions for Poisson distribution, with the conditional mean ((ref) ) and conditional variance ((ref)), $y_{i}|\mathbf{x}_{i}$ is shown to follow a negative binomial II distribution. It implies overdispersion, but where the amount of overdispersion increases with the conditional mean,
Now the log-likelihood function for observation $i$ is
where $\Gamma \left( \cdot \right) $ is the gamma function defined for $r>0$ by $\Gamma \left( r\right) =\int_{0}^{\infty }z^{r-1}\exp \left( -z\right) dz$. For fixed $\tau ^{2}$, the log likelihood equation in ((ref)) is in the exponential family; see gourieroux1984pseudo. Thus the negative binomial QMLE using ((ref)) is consistent under conditional mean assumption only, which is the same as the Poisson QMLE. Since the negative binomial II likelihood captures the nature of the variance function, it should deliver more efficient estimation when the data generating process is correctly specified, although the spatial correlation is not accounted. Again, we can use a GEE working correlation matrix to account for the spatial correlation.
The Probit model is one of the popular binary response models. The dependent variable $y$ has conditional Bernoulli distribution and takes on the values zero and one, which indicates whether or not a certain event has occurred. For example, $y=1$ if a firm adopts a new technology, and $y=0$ otherwise. The value of the latent variable $y^{\ast }$ determines the outcome of $y$.
Assume the Probit model is
We do not observe $y_{i}^{\ast }$; we only observe $y_{i}.$ Let $\Phi \left( \cdot \right) $ be the standard normal cumulative density function ( CDF), and $\phi $ be the standard normal probability density function ( PDF). Assume that the mean function $m_{i}\left( \mathbf{x}_{i}; \mathbf{\beta }\right) \equiv $ $\mathop{\mbox{\sf E}}\left( y_{i}|\mathbf{x}_{i}, D_n\right) =\Phi \left( \mathbf{x}_{i}\mathbf{\beta }\right) $ is correctly specified. $e$ is the spatial correlated latent error. Let $ \mathbf{e}=\left( e_{1},e_{2},...,e_{n}\right) ^{\top }$. For example, pinkse1998contracting use the following assumption of $\mathbf{e}$:
where $\mathbf{\varepsilon =}\left( \varepsilon _{1},\varepsilon _{2},...,\varepsilon _{n}\right) $ which has a standard normal distribution. $W$ is a $n\times n$ weight matrix with zeroes on the diagonal and inverse of distances off diagonal. $\rho $ is a correlation parameter. We can see $e$ can be written as a function of $\varepsilon ,$
Thus the conditional expectation of $\mathbb{e}$ is zero. The variance covariance matrix of $\mathbb{e}$ is
If we assume that $e|x$ has a multivariate normal distribution with mean zero and variance matrix specified in ((ref)). Thus a much simpler specification is to directly model $e|x$ as a multivariate distribution. Different from the usual multivariate distribution\footnote{ A multivariate normal distribution usually specifies the mean vector and correlation matrix. The correlations do not depend on the pairwise distance between two variables.}, the covariances of $e$ should depend on the pairwise distances $d_{ij}$. We also let the covariances depend on a parameter $\rho $. The above equation can be written in a conditional mean form:
It is very natural to write the variance function for a Bernoulli distribution,
We are interested in the partial effects of $x$ to $y$. For a continuous $ x_{K}$ the partial effect is
For a discrete $x_{K}$, the partial effects when $x_{K}$ changes from $a_{K}$ to $a_{K}+1$ is
A simple one-step estimation is the pooled Bernoulli quasi-MLE (QMLE ), which is obtained by maximizing the pooled Probit log-likelihood. The log likelihood function for each observation is
Let $\check{u}_{i}=y_{i}-\Phi \left( \mathbf{x}_{i}\mathbf{\check{\beta}} \right) ,i=1,2,...,n$ be the residuals from the partial QMLE estimation. At this stage, a robust estimator for the asymptotic variance of $ \mathbf{\check{\beta}}_{\mathrm{PQMLE}}$ can be computed as follows:
where $k\left( d_{ij}\right) $ is the kernel weight function that depends on pairwise distances. This partial QMLE and its robust variance-covariance estimator provides a legitimate way of the estimation of the spatial Probit model.
We use partial QMLE as a first-step estimator. An estimator for the working variance matrix for each group is
And assume the working correlation function for $l$th and $m$th elements in group $g$ is
For example, suppose that
Let $\check{u}_{i}$ be the partial QMLE residual and $\hat{r}_{i}=\check{u} _{i}/\sqrt{\check{v}_{i}}$, for $i=1,2,...,n,$ be the standardized residuals. $\mathbf{\hat{C}}_{ij}$ equals the sample correlation of $\check{u }_{i}/\sqrt{\check{v}_{i}}$ and $\check{u}_{j}/\sqrt{\check{v}_{j}}$. Using the correlations within groups, one estimator of $\rho $ is
for $l<m.$
The second-step GEE estimator for $\mathbf{\beta }$ is
The first order condition is
$\mathbf{\hat{\beta}}_{\mathrm{GEE}}$ is consistent and follows a normal distribution asymptotically by Theorem (ref). $ \mathbf{\hat{\beta}}$ is consistent even for misspecified spatial correlation structure $\mathbf{\hat{W}}_{g}$. The asymptotic variance estimator that is robust to misspecification of spatial correlation is:
where $k(d_{gh})$ is a kernel function which depends on the distances between groups.
An alternative approach is to specify the specific distributions of the multivariate normal distribution of the latent error, and then find the estimator for the spatial correlation parameter for the latent error within a MLE framework. For example, see wang2013partial.
In this section, we provide the assumptions and results on the theoretical properties our GEE estimation.
From now on we work with group level asymptotics. Define the field $\tilde{\varepsilon} = \{\varepsilon_g: g \in 1, \cdots, G\}$ with grouped observations. First of all suppose that $D_n$ is divided by $G$ blocks with $\cup^G_1 B_g = D_n \subset T_n$, and the group level lattice is denoted as $D_G$. Define the distance between two groups $g,h$ as $\rho(g,h) = \mathbf{min}_{i \in B_g, j \in B_h} \rho(i,j).$ And the $\alpha-$ mixing coefficient between two union of groups for $U = \{g_1,\cdots, g_L\}$, $V = \{h_1, \cdots, h_M\}$, $\rho(U,V) = \mathbf{min}_{l \in 1 \cdots L,m \in 1,\cdots, M} \rho(g_l,h_m)$ is thus $\tilde{\alpha}(u, v, r) = \tilde{\alpha}(L\leq u, M \leq v, \rho(U,V) \geq r) = \sup_{L\leq u, M \leq v, \rho(U,V) \geq r} \alpha(\sigma(U), \sigma(V)) $. If the group size are the same, i.e. $L$, then the mixing coefficients of the grouped observations have the following relationship with respect to it in the original field $\tilde{\alpha}(u, v, r) = \alpha(uL, vL, r).$ We can assume $\tilde{\alpha}(u, v, r) = (uL+ vL)^{\tau} \hat{\alpha}(r)$.
Assume that $L^{\tau}\hat{\alpha}(r) \to 0$ as $r \to \infty,$ and $\tilde{\varepsilon}$ would maintain the $\alpha-$ mixing property. Define the ball around group $g$ with radius $s$ to be $\mathcal{F}_g(s) = \sigma\{\cup_{h: \rho(g,h)\leq s} B_h\}.$
Remark A.8) can be implied from positive definiteness of $\mathbf{W}_g(\theta, \gamma)$ and the same identification assumption $\liminf_{G\to \infty}\mbox{inf}_{\theta \in \Theta: \nu(\theta, \theta_0) \geq \varepsilon }Q'_{G} (\theta, \gamma)> c_0>0$ on $Q'_{G} (\theta, \gamma)\stackrel{\mathrm{def}}{=} \frac{1}{M_{G}|D_G|} \sum_{g\in |D_G|} \mathop{\mbox{\sf E}}\left[ \mathbf{y}_{g}-\mathbf{m} _{g}\left( \mathbf{x}_{g};\mathbf{\theta }\right) \right]^{\top }\left[ \mathbf{y}_{g}-\mathbf{m} _{g}\left( \mathbf{x}_{g};\mathbf{\theta }\right) \right]$. As it can be seen that with probability $1- {\scriptstyle{\mathcal{O}}}_p(1)$ \\$\liminf_{G\to \infty}inf_{\theta \in \Theta: \nu(\theta, \theta_0) \geq \varepsilon }Q_{G}\left( \mathbf{\theta }, \mathbf{\gamma}\right)> \liminf_{G\to \infty}inf_{\theta \in \Theta: \nu(\theta, \theta_0) \geq \varepsilon }\lambda_{min} \{\mathbf{W}_g(\theta, \gamma)\} Q'_{\infty} (\theta, \gamma),$ where $\lambda_{min} \{\mathbf{W}_g(\theta, \gamma)\}$ is the minimum eigenvalue of the matrix $\lambda_{min} \{\mathbf{W}_g(\theta, \gamma)\}$. If we assume that with probability $1-{{\mathcal{O}}}_p(1),$ $\lambda_{min} \{\mathbf{W}_g(\theta, \gamma)\}>c$ where $c$ is a positive constant. We now comment on assumptions, Condition A.2) is concerning the $L_2$ NED property of our data generating processes. A.3) and A.4) are the standard regularities assumptions. A.5) is a few moment assumptions on the statistical objects involved in the estimation. A.6) is the mixing coefficients restrictions after grouping observations. A.7) is again moment conditions on the elementwise Hessian matrices. A.8) is a condition on identification of our estimator. Given the assumptions, we can provide the consistency property of our estimation.
Theorem (ref) indicates the consistency of the estimation as long as the number of groups tends to infinity. The proof is in the Appendix. To prove further the asymptotic normality of the estimation we need in addition the following assumptions.
A.9) is concerning the the pre-estimation of the nuisance parameter $\gamma$, and A.10), A.11) are two standard assumptions on the regularities of the estimation. Note that $\mathbf{S}_{G}\left( \mathbf{ \hat{\theta},\hat{\gamma}} \right) = {\scriptstyle{\mathcal{O}}}_p(1) = 0$ if $\hat{\theta}, \hat{\gamma}$ lies in the interior point of the parameter space. In the following, we verify that with our proposal of estimating $\hat{\gamma}$ in ((ref)) in Section (ref) , we will achieve A.9).
$\mathbf{H}_{\infty} \stackrel{\mathrm{def}}{=} \lim_{G\to \infty} \mathop{\mbox{\sf E}}\mathbf{H}_G(\theta^0,\gamma^0),$ where $\mathbf{H}_G(\theta^0,\gamma^0)$ is defined in equation ((ref)). It is not surprising to see that our estimation will be asymptotically normally distributed, with a variance covariance matrix of a sandwich form $AV(\theta^0)$, which involves the Hessian. The rate of convergence is shown to be $\sqrt{G}$.
In this subsection, we propose a semiparametric estimator of the asymptotic variance in Theorem (ref), and prove its consistency. The estimation is tailored to account for the spatial dependency of the underlying process. This facilitates us to create a confidence interval for our estimation.
First let
where $\nabla \mathbf{\hat{m}}_{g}\equiv \nabla \mathbf{\hat{m}}_{g}\left(\mathbf{\hat{\theta}}\right) ,$ $\mathbf{\hat{W}}_{g}\equiv \mathbf{\hat{W}}_{g}(\mathbf{\hat{\gamma}}, \hat{\theta})$.
The estimator of $\mathrm{AV}\left( \mathbf{\theta}^0\right)$ which is robust to misspecification of the variance covariance matrix is
where $k\left( d_{gh}\right) $ is the kernel function depending on the distance between group $g$ and $h$, i.e. $\rho(g,h)$, and a bandwidth parameter $h_g$. As noted in kelejian2007hac, there are many choices for the kernel functions, such as rectangular kernel, Bartlett or triangular kernel, etc. In particular, without loss of generality, we can choose the Bartlett kernel function $k\left( d_{gh}\right) = 1- \rho(g,h)/h_g$, for $\rho(g,h) < h_g$ and $k\left(g,h\right) =0$ for $ \rho(g,h) \geq h_g$. Further, we can obtain the average partial effects (APE) of interest and carry on valid inference.
We now list the assumptions needed for the consistency of estimator of $AV(\theta^0)$.
B.1) is an assumption for decomposing the difference between the residuals and the true error, as in kelejian2007hac. B.2) is about the moment bound and B.3) is on property of the kernel function. B.4) constrains on the spatial dependence coefficients and the bandwidth length. We provide in the following theorem the consistency of the $\widehat{\mathrm{AV}}\left( \mathbf{\hat{\theta}}\right)$. It is worth noting that we prove an elementwise version of the consistency, and the results below can be verified equivalently in any matrix norm, as we consider fixed dimension parameter.
In this section, we use Monte Carlo simulations to investigate the finite sample performances of our proposed GEE approach with groupwise data compared to the partial QMLE. We simulated count data and binary response data separately. We show that our GEE method is very critical for improving the efficiency of our estimation.The simulation mechanism is described as follows.
We use sample sizes of 400 or 1600. We sample observations on a lattice. For example, for sample size of 400, the sample space is a $20\times 20$ square lattice. Each observation resides on the intersections of this lattice. The locations for the data are $\{(r,s):r,s=1,2,...,20\}$. The distance $d_{ij}$ between location $i$ and $j$ is chosen to be the Euclidean distance. Suppose $A(a_{i},a_{j})$ and $B(b_{i},b_{j})$ are the two points on the lattice; their distance $d_{ij}$ is $\sqrt{ (a_{i}-b_{i})^{2}+(a_{j}-b_{j})^{2}}$. The spatial correlation is based on a given parameter $\rho $ and $d_{ij}$. The data are divided into groups of 4 and the number of groups are set to be 100 for sample size 400. Similarly, for the sample size of 1600, we use a $40\times 40$ lattice. We still use sample size of 4 in each group and there are 400 groups in total. For simplicity, we keep the pairwise distances in different groups the same.
In the count data case, for a Poisson distribution the variances and covariances of the count dependent variable can be written in closed forms given the spatial correlation in the underlying spatial error term. That is, by knowing the correlations in the spatial error term, we can derive the correlations in the count dependent variable as shown in ((ref)) and ((ref)). Consider the following spatial count data generating process: 1. $ v_{i}$ is simulated as a multivariate lognormal variable with E$ \left( v_{i}\right) =1$, exponentiating an underlying multivariate normal distribution using with correlation matrix $W$. Let $a_{i}$ be the underlying multivariate normal distributed variable. Then $v_{i}=\exp \left( a_{i}\right) $ follows a multivariate lognormal distribution. We describe the underlying spatial process in Case 1, 2, and 3 as three special cases to demonstrate different spatial correlations. 2. The coefficient parameters and explanatory variables are set as follows: $\beta _{1}=0.5,\beta _{2}=1,\beta _{3}=1,\beta _{4}=1.$; $ x_{2}\sim $N$\left( 0,0.25\right), x_{3}\sim\mathrm{Uniform}\left( 0,1\right), x_{5}\sim $N$\left( 0,1\right), x_{4}=1[x_{5}>0].$ 3. The mean function for individual i is $m_{i}=v_{i}\exp \left( \beta _{1}+\beta _{2}x_{2}+\beta _{3}x_{3}+\beta _{4}x_{4}\right) ;$ 4. Finally we draw the dependent variable from the Poisson distribution with mean $m_{i}$: $y_{i}\sim\mathrm{Poisson}\left( m_{i}\right) .$ Specifically, the underlying spatial error $a_{i}$ has the following three cases.
Case 1. $a_{i}=\left( I-\rho W\right) ^{-1}e_{i},e_{i}\sim $N$\left( 0,1\right);$ $W$ is the matrix with $W_{g}$ on the diagonal, $g=1,2,...,G.$ Other elements in $W$ are equal to zero. For group size equal to four,
Case 2. $a_{i}=\left( I-\rho W\right) ^{-1}e_{i},e_{i}\sim $N$\left( 0,1\right);$ $W$ is the matrix with $W_{g}$ on the diagonal, $g=1,2,...,G.$ The $(l,m)$th element in $W_{g}$, $W_{g\_lm}=$\ \ $\frac{\rho }{6\ast d_{g_{\_}lm}},$ $\rho =0,0.5,1,1.5,l\neq m;$ $W_{g_{\_}lm}=0,l=m$ for group $g$. Correlations are zero if observations are in different groups. For group size equal to four,
Case 3. In this case, the DGP has the following differences from Case 1 and Case 2. $a_{i}$ is simulated as a multivariate lognormal variable by exponentiating an underlying multivariate normal distribution N$ \left( -\frac{1}{2},1\right) $ using with correlation matrix $W$. $W_{ij}=$\ \ $\frac{\rho }{d_{ij}},$ $\rho =0,0.2,0.4,0.6,i\neq j;$ $W_{ii}=1;$ $ i,j=1,2,...,N.$ The underlying normal distribution implies that $v_{i}$ follows a multivariate lognormal distribution with E$\left( v_{i}\right) =1$. We set $\beta _{1}=-1,\beta _{2}=1,\beta _{3}=1,\beta _{4}=1.$ $ x_{2}$ follows a multivariate normal distribution N$\left(0,W\right) ;$ In this case, the data has general spatial correlations for each pair of observations if $\rho \neq 0.$
Table (ref), Table (ref) and Table (ref) show three cases of simulation results with 1000 replications with two different samples and group sizes: (1) $N=400,$ $ G=100,$ $L=4$ (2) $N=1600,$ $G=400,$ $L=4$. There are four estimators, Poisson partial QMLE estimator, \ Poisson GEE, Negative Binomial II (NB II) partial QMLE, and NB II GEE. For simplicity, we use an exchangeable working correlation matrix for GEE estimators. We can see that, first as spatial correlation increases the GEE methods has smaller standard deviations than QMLE. Second, when there is little spatial correlation, GEE does not increase much finite sample bias due to accounting for possible spatial correlation.
In Case 1, when there is no spatial correlation, the Poisson QMLE should be as efficient as GEE asymptotically. We can see that when $\rho =0,$ the coefficient estimates and their standard deviations of Poisson QMLE and GEE are pretty close, which means that there is little finite sample bias due to accounting for possible spatial correlation when there is actually no spatial correlation. The standard deviations for the estimated coefficients of Poisson QMLE and GEE are almost the same. The standard deviation of $\hat{\beta}_{2}$ equals $ 0.259$ for Poisson QMLE and $0.260$ for Poisson GEE when $\rho =0$ for a sample size of 400. As $\rho $ grows larger. the GEE estimator shows more and more efficiency improvement over the partial QMLE. For example, for a sample size of 400, when $\rho =1,$ the standard deviation of $\hat{\beta} _{2}$ equals $0.267$ for Poisson QMLE and $0.259$ for Poisson GEE. When $ \rho =1.5,$ the standard deviation of $\hat{\beta}_{2}$ equals $0.320$ for Poisson GEE and $0.302$ for Poisson PQMLE. The NB II GEE also has some improvement over NB II PQMLE. When $\rho =1,$ the standard deviation of $ \hat{\beta}_{2}$ equals $0.234$ for NB II PQMLE and $0.226$ for NB II GEE. When $\rho =1.5,$ the standard deviation of $\hat{\beta}_{2}$ equals $0.276$ for NB II PQMLE and $0.261$ for NB II GEE. When sample size increases from 400 to 1600, we see the similar scenarios. Case 2 and Case 3 have shown similar efficiency results for the GEE estimators.
{
}
{\tiny
}
For the Probit model, the correlations of latent normal errors result in correlations of binary response variables, but we cannot easily find the specific form of the conditional variances and covariances for the binary dependent variables. The correlations in latent error do not reflect the exact correlations in the binary dependent variables. Consider the following cases of data generating process. 1. The latent variable $y^{\ast }=\beta _{1}+\beta _{2}x_{2}+\beta _{3}x_{3}+\beta _{4}x_{4}+e_{4},$ where $e_{4}$ is the latent spatial error term, and the parameters are set to be $\beta _{1}=\beta _{2}=\beta _{3}=\beta _{4}=1.$ Then the binary dependent variable is generated as $y_{i}=1$ if $y_{i}^{\ast }\geq 1.5$ and $y_{i}=0$ if $ y_{i}^{\ast }<1.5.$ The explanatory variables are set as follows: $\ x_{1}=1; $ $x_{2}\sim \mathrm{N}\left( 1,1\right) ;$ $\ x_{3}=0.2x_{2}-1.2e_{1},e_{1}\sim \mathrm{N}\left( 0,1\right) ;$ $\ x_{5}=0.2x_{2}+0.2x_{3}+e_{2},e_{2}\sim \mathrm{N}\left( 0,1\right) ;$ $ x_{4}=1\left[ x_{5}>0\right] .$ We consider two cases of latent spatial error terms and the corresponding binary response variables are generated as follows.
Case 1. The vector of spatial error $\mathbf{e}_{4}=\left( I-\rho W\right) ^{-1}\mathbf{e}_{3},e_{3}\sim \mathrm{N}\left( 0,1\right) ,$ where $\rho =0,0.5,1,1.5$ respectively. $W$ is the matrix with $W_{g}$ on the diagonal, $ g=1,2,...,G.$ Other elements in $W$ are equal to zero. In this case, only individuals within a group are correlated. For group size equal to four, $ W_{g}$ is the same as in ((ref)) in Case 1 for count data.
Case 2. The latent spatial error $e_{4}\sim \mathrm{MVN}(0,$W$ ),$ that is, $e_{4}$ follows a standard multivariate normal distribution with expectation zero and $W$ is $N\times N$ correlation matrix. $W_{ij}=\frac{\rho }{d_{ij}},$ $\rho =0,0.2,0.4,0.6,i\neq j;$ $W_{ii}=1;$ $ i,j=1,2,...,N.$ $W$ is the same as in ((ref)) in Case 3 for count data. Therefore, the data has general spatial correlations for each pair of observations if $\rho \neq 0.$
In the simulation, two estimators are compared, the Probit partial QMLE estimator, and the Probit GEE estimator with an exchangeable working correlation matrix. We show two cases of the simulation: (1) N=400, G=100, L=4; 2) N=1600, G=400, L=4. The replication times are 1000. The simulation results for Case 1 and Case 2 are in Table (ref) and Table (ref) separately. We find the following results.
First, in both cases, the GEE estimator is less biased than the partial QMLE estimator. For example, for N=400, in Case 1 when $\rho =1,\hat{\beta} _{2}$ equals 1.252 for QMLE and 1.203 for GEE. In Case 2 when $\rho =0.6, \hat{\beta}_{2}$ equals 1.148 for QMLE and 1.090 for GEE. Second, the GEE estimator has some obvious efficiency improvement over partial QMLE. For example, in case 1 when $\rho =1$, the standard deviation of $\hat{\beta} _{2} $ equals $0.280$ for QMLE and $0.173$ for GEE. In Case 2 when $\rho =0.6,$ the standard deviation of $\hat{\beta}_{2}$ equals 0.270 for QMLE and 0.164 for GEE for a sample size of 400. Third, when we increase the sample size to 1600 and number of groups to 400 correspondingly, the same scenario applies. What is more, the bias and especially standard deviations for both the Probit QMLE and GEE reduces. For example, for N=1600, in Case 1 when $\rho =1,$ the standard deviations of $\hat{\beta}_{2}$ reduce to 0.121 for QMLE and 0.081 for GEE.
In the empirical FDI literature, the gravity equation specification was initially adopted from the empirical literature on trade flows. The gravity equation has been widely used and extended in international trade since tinbergen1962analysis. anderson2003gravity specify the gravity equation as
where $T_{ij}$ is the trade flows between country $i$ and country $j$. $ T_{ij}$ is proportional to the product of the two countries' GDPs, denoted by $Y_{i}$ and $Y_{j}$, and inversely proportional to their distance. $ D_{ij} $ broadly represents trade resistance. Let $\eta _{ij}$ be a stochastic error that represents deviations from the theory. As a tradition in the existing literature, by taking the natural logarithms of both sides and adding other control variables represented by $Z_{ij}$, the log-linearized equation is:
For the above equation, a traditional estimation approach is to use ordinary least squares (OLS). However, there are two problems with the OLS estimation of the log linearized model. First, $T_{ij}$ must be positive in order to take the logarithm. A transformation of $\log(T_{ij}+1)$ can solve the problem of logarithm but it is not clear how to interpret the estimation results with respect to the original values. Second, the estimation heavily depends on the independence assumption of $\eta _{ij}$ and explanatory variables, which means the variance of $\eta _{ij}$ cannot depend on the explanatory variables. Because of taking the logarithm, only under very specific conditions on $\eta _{ij}$ is the log linear representation of the constant-elasticity model useful as a device to estimate the parameters of interest (silva2006log). Jensen's inequality implies that ($ \mathop{\mbox{\sf E}}\log Y $) is smaller than $\log \mathop{\mbox{\sf E}}(Y)$, thus log-linearized models estimated by OLS as elasticities can be highly misleading in the presence of heteroscedasticity. If the variance of $\eta _{ij}$ is dependent on the explanatory variables, ordinary least squares is not consistent any more.
We adopt this specification and augment it to the inflow FDI to cities of China. and use nonlinear estimation method, the GEE estimation. The estimating equation is specified as follows
where $FDI_{i}$ is the inflow FDI in actual use for city $i$, $X_{i}$ represents all explanatory variables. The control variables includes city level GDP, GDP per capita, the average wage, the government expenditure to science, and whether the city is on the border. We collect data of inflow FDI to 287 cities in 31 provincial administrative regions in 2007 in mainland China from the website of Development Research Center of the State Council of P. R. China \footnote{The website of Development Research Center of the State Council of P. R. China is www.drcnet.com.cn}. Three cities, Jiayuguan (Gansu Province), Dingxi (Gansu Province) and Karamay (Xinjiang Province), are dropped because of missing data on FDI. Thus we are using 284 cities in total. We collect the latitudes and longitudes of the center of each city using Google map and calculated the geographical distance matrix between cities. The city center is defined as the location of the city government. We use provinces as natural grouping so there are 31 groups. Each group has one to twenty cities. The descriptive statistics are in Table (ref). The grouping information is in Table (ref).
For comparison, we also provide the OLS estimates of the log-linearized model:
The log linearized model suffers from two main problems, first the dependent variable cannot take log if it is zero; second as mentioned in Silva and Tenreyro (2006) the log linearization can cause bias in parameter estimates if there exists heteroskedasticity in the error term $u_{i}$.
To estimate the equation for FDI, we use OLS, Poisson QMLE, Poisson GEE with the exchangeable working matrix, NB QMLE, NB GEE with the exchangeable working matrix. In Table (ref) the results show advantage of Poisson GEE estimation. All estimation results verifies the positive effect of GDP and GDP per capita in the gravity equation for FDI. These estimates are all significant at the 1% level. What is more, the standard error of GDP and GDP per capita for Poisson GEE is smaller than that for Poisson QMLE, which is smaller than that for OLS. The Poisson regression has significant results on the explanatory variables, log(wage), log(sciexp) and border, which are not significant in the OLS regression. The local average wage has a negative effect on inflow FDI to this city. Compared to other estimation methods, the Poisson GEE estimates on log(wage) is the most significant, at 1% level. It means that when the average wage increase by 1%, the inflow FDI would decrease by about 1%, which could due to the inhabiting effect of labor cost. Similarly, when local government increase science expenditure by 1%, the inflow FDI would increase by about 0.3%, which is shown by Poisson QMLE and Poisson GEE, and in which case the Poisson GEE estimate has smaller standard error than Poisson QMLE, which are 0.102 and 0.110 respectively.