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.
125,111 characters · 8 sections · 91 citation commands
Efficient closed-form estimation of large spatial autoregressions
Keywords: Spatial autoregression, efficiency, many parameters, networks
JEL Classification: C21, C31, C33, C36
\allowdisplaybreaks
Spatial autoregressive (SAR) models, introduced by cliff1973spatial, are popular tools for modelling cross-sectionally dependent economic data. The pre-eminent feature of such models is the presence of one or more `spatial weight' matrices, which parsimoniously capture the dependence between units in the sample. Such dependence need not be geographic in nature, indeed the spatial weight matrix is known by other terms such as `adjacency matrix', `network link matrix' and `sociomatrix'. For $n\times 1$ vectors $y_n$ and $u$ of responses and unobserved disturbances, respectively, and an $n\times k$ covariate matrix $X_n$, the SAR model is
where the elements of the $n\times n$ spatial weight matrices $W_{in}$ are inverse economic distances and $\lambda_{0n}=\left(\lambda_{01n},\ldots,\lambda_{0pn}\right)'$ and $\beta_{0n}$ are unknown parameter vectors. Subscripting with $n$ permits treatment of triangular arrays, an important issue for spatial models in general (see Robinson2011), and for SAR models even more so due to various normalizations of the $W_{in}$ that make it $n$-dependent. This paper justifies computationally straightforward estimation for the parameters of ((ref)) with the same asymptotic properties as pseudo maximum likelihood estimates.
SAR models allow dependence to occur across a very generalized notion of space: so long as a mapping exists between every pair of individuals to the real line a spatial weight matrix may be constructed. The flexible nature of the SAR model means that it may be used to model a very wide range of phenomena. Thus it has found application in many fields of economics such as development economics (case1991spatial, Helmers2014), industrial organization (pinkse2002), trade (Conley2003) and peer effects (Hsieh2018), to name only a few examples. Another frequently used approach to model cross-sectional dependence is the `common factor' technique, see e.g. Chudik2015 for a review.
Estimation of SAR models has long been considered in the regional science literature, see e.g. anselin1988spatial. Rigorous asymptotic theory for instrumental variables (IV) estimation was initially provided by kelejian1998generalized, leading to the present flourishing theoretical literature. lee2002consistency studied ordinary least squares (OLS) estimation of SAR models, stressing the need for lack of sparsity in the spatial weight matrix to establish desirable asymptotic properties such as consistency and efficiency. This was followed by lee2004asymptotic, a seminal contribution that provided a taxonomical asymptotic theory for Gaussian pseudo maximum likelihood estimates (PMLE) of SAR models. Recently Kuersteiner2013,Kuersteiner2020 have provided general theory for such models in a panel data setting.
The flexible nature of SAR modelling is further embellished by the seamless ability to integrate more than one spatial weight matrix in the model ((ref)), thus permitting simultaneous connections between units across a number of channels. This is an accurate representation of typical economic situations, e.g. countries are `connected' by both geographical proximity as well as trade ties. Furthermore, in many economic settings the sample partitions naturally into $p$ clusters or groups, leading to block diagonal structure for the spatial weight matrix $V_n=diag\left(V_{1n},\ldots,V_{pn}\right)$, where $V_{in}$ is $m_i\times m_i$ and $\sum_{i=1}^p m_i=n$. To permit the modelling of heterogenous spillover effects across clusters, one may take $W_{in}$ to be the $n\times n$ block diagonal matrix with the $m_i\times m_i$ dimensional $i$-th diagonal block given by $V_{in}$. This approach has been suggested by Gupta2013,Gupta2018. Here and more generally, the specification ((ref)) is termed a `higher-order' SAR model if $p>1$, see e.g blommestein1983specification, Lee2010, Li2017, Han2017, Kwok2019.
In the study of higher-order SAR models, Gupta2013,Gupta2018 have suggested that $p,k$ be allowed to diverge slowly to infinity as functions of sample size. The motivation for such generality is typically threefold: first, it is desirable to permit a richer model as the sample size permits. Second, clustered data as mentioned in the previous paragraph naturally imply asymptotic regimes with increasing $p$. For instance, when $m_i=m$ for each $i=1,\ldots,p$, we have $n=mp$ and the results of lee2004asymptotic imply that $p\rightarrow\infty$ is necessary for consistent estimation, analogous to the problems created in the spatial statistics literature by `infill asymptotics', see e.g. Lahiri1996. Finally, a theory that allows the model dimension to grow with sample size provides a more incisive analysis of large models in practice, much as typical asymptotic theory with a fixed parameter space itself can be thought as providing an approximation in finite samples.
The estimation of such increasing-order SAR models has been studied by Gupta2013,Gupta2018 using IV, OLS and PMLE approaches. The first two methods have the advantage of being in closed-form, while even for $p=1$ PMLE (in)famously requires grid search and the inversion of an $n\times n$ matrix in every iteration, leading to many ingenious solutions for faster computation, see e.g. Ord1975 and Pace1997. The computational cost of PMLE in SAR type models is particularly salient as data sets increase in size, as stressed by Zhu2020. Modern network data sets are amenable to modelling via SAR techniques and can feature, or accommodate, large parameter spaces but computation remains a serious challenge. Han2020 provide a discussion of the problems and propose a Bayesian solution.
These problems are naturally exacerbated if $p>1$, with grid search requiring more iterations to converge and each iteration requiring inversion of an $n\times n$ matrix, as well as risk of convergence to local optima. Furthermore, the requirement of a compact parameter space for $\lambda_{0n}$ can severely restrict the admissible parameter values (see Gupta2018). On the other hand, under Gaussianity the PMLE becomes the MLE and is efficient. This property is shared by OLS, but under rather delicate and specific conditions even for $p=1$ (see lee2002consistency). Thus the IV/OLS and PMLE approaches each have their advantages and it is desirable to combine the positive properties of both.
One method of obtaining closed-form estimates with the same asymptotic covariance matrix as a target estimate is to use Newton-type iterations commencing from an initial consistent estimator that is straightforward to compute. The approach dates back at least to Fisher1925 and LeCam1956. It enjoys the added attraction of avoiding a potentially complicated consistency proof for an implicitly defined estimate, as well as the compactness assumptions this typically entails. As a result, the technique has been used in a vast variety of settings, see e.g. Rothenberg1964 (simultaneous equations), Hartley1965 (nonlinear least squares), Janssen1985 ($M$-estimation), Rothenberg1984 (generalized least squares), Hualde2011, Kristensen2006, Robinson2005a (time series and adaptive estimation), Andrews1997 (generalized method of moments), Kasahara2008, Kristensen2017 (structural estimation), DeLuca2018 (generalized linear models) and Frazier2017 (efficient two-step estimation), to name just a few.
In this paper we use IV and OLS estimates as initial estimates to form a single Newton-step asymptotic approximation to the Gaussian PMLE with $p=p_n$ and $k=k_n$ allowed to diverge as functions of $n\rightarrow\infty$. The approach has been studied in the case of fixed-dimensional SAR models by robinson2010efficient and Lee2013, but the previous discussion hints at its particular usefulness when considering large models. One avoids grid search over a high-dimensional parameter space, compactness assumptions on this space and the inversion of large ($n\times n$) matrix for every search iteration, as well as various headaches related to convergence and local optima. When commencing from IV estimates, this leads to closed-form efficient estimates under Gaussianity. As suggested by the results of lee2002consistency and Gupta2013, commencing iteration from OLS preserves the efficiency property. However, we show that the Newton step approach cancels out certain terms of large stochastic order that allows for weaker rate conditions than those imposed in these papers.
In a simulation study, we demonstrate that the Newton step can lead to much improved estimates in finite samples, both in terms of bias and efficiency. While a single step is sufficient to establish desirable asymptotic properties, in our simulation study we also explore the finite sample implications of additional Newton steps, reporting results with up to six iterations. We find large finite sample gains in both bias and mean squared error that are robust to heavy tailed error distributions. We also observe fast convergence of iterations, which conforms to extant theoretical observations. The gains are particularly notable when the parameter space and sample size is large, a situation in which PMLE becomes computationally onerous. In a small illustration with real world data, we show that the estimates work well in practice and lead to more precise results.
We collect some frequently used notation here for the convenience of the reader. For a generic matrix $A$ denote $\left\Vert A\right\Vert=\left(\overline{\eta}\left(A'A\right)\right)^{\frac{1}{2}}$, with $\overline\eta(\cdot)$ and $\underline\eta(\cdot)$ denoting the largest and smallest eigenvalues, respectively, of a symmetric positive semidefinite matrix. Note that if $A$ is a vector then $\left\Vert A\right\Vert$ is simply its Euclidean norm. Let $\left\Vert A\right\Vert_R$ denote the maximum absolute row sum norm of $A$. For any parameter $\tau$, function $f(\tau)$ and generic estimate $\check\tau$, we will write $\check f\equiv f\left(\check\tau\right)$. We denote true parameter values with $0$ subscript and suppress the argument for a quantity evaluated at a true parameter value, i.e. $f\left(\tau_0\right)\equiv f$.
The ($-2/n$ times) log pseudo Gaussian likelihood function for model ((ref)) at any admissible point $\theta=(\lambda',\beta')'$ is given by
where $S_n(\lambda)=I_n-\sum_{i=1}^{p_n} \lambda_{in}W_{in}$, with $I_n$ denoting the $n\times n$ identity matrix. If $S_n$ is invertible, ((ref)) admits the reduced form $y_n=S_n^{-1}X_n\beta_{n}+S_n^{-1}u$, and we define $R_n=A_n+B_n,$ where $A_n = (G_{1n}X_n\beta_{n},\ldots,G_{{p_n}n}X_n\beta_{n}),B_n = (G_{1n}u,\ldots,G_{{p_n}n}u),$ $G_{in}(\lambda)=W_{in}S_n^{-1}(\lambda)$, $i=1,\ldots,p_n$, and so $R_n=(W_{1n}y_n,\ldots,W_{{p_n}n}y_n)$.
Defining $\mathcal{R}_n^{y}\left(\theta\right)=R_n\lambda_n+X_n\beta_n-y_n$, the derivative of ((ref)) at any admissible $\left(\theta,\sigma^2\right)$ is
where $ \varphi_n \left(\theta,\sigma^2\right)=2\sigma^{-2}n^{-1} \left(\sigma^{2}tr G_{1n}(\lambda)+y_n'W_{1n}'\mathcal{R}_n^{y}\left(\theta\right),\ldots,\sigma^{2}tr G_{p_nn}(\lambda)+y_n'W_{p_nn}'\mathcal{R}_n^{y}\left(\theta\right)\right)'. $ Because $\mathcal{R}_n^{y}=-u$, denoting $ \phi_n = \sigma_0^{-2}n^{-1} \left( \sigma_0^2 tr C_{1n}-u'C_{1n}u ,\ldots,\sigma_0^2 tr C_{p_n}-u'C_{p_nn}u \right)'$ with $C_{in}=G_{in}+G_{in}'$, we obtain
with $t_{n} = n^{-1} \left[A_n,X_n\right]'u$. The Hessian at any admissible point in the parameter space is
where $P_{ji,n}(\lambda_{n})$ is the $p_n\times p_n$ matrix with $(i,j)$-th element given by $tr \left(G_{jn}(\lambda_{n})G_{in}(\lambda_{n})\right)$.
Let $Z_n$ be an $n\times r_n$ matrix of instruments, with $r_n\geq p_n$, and define the IV and OLS estimates as
respectively, with $\hat{Q}_n=\hat{K}_n'J_n^{-1}\hat{K}_n$, $ \hat{K}_n=n^{-1} \left[Z_n, X_n\right]'[R_n,X_n],\hat{k}_n=n^{-1} \left[Z_n , X_n\right]'y_n,\;J_n=n^{-1} \left[Z_n ,X_n\right]'\left[Z_n,X_n\right], $ and $ \hat{L}_n=n^{-1} \left[R_n, X_n\right]'[R_n,X_n],\hat{l}_n=n^{-1} \left[R_n,X_n\right]'y_n. $ Define the respective `one-step' estimates $\hat{\hat{\theta}}_{n}$ and $\tilde{\tilde{\theta}}_{n}$ by the following equations
We observe that other initial estimates, such as the GMM estimates of Kelejian1999 and Lee2007, can also be used. However we choose initial estimates that are available in closed form for computational ease. While consistent initial estimates are needed to obtain a desirable asymptotic theory, even in the fixed-dimension parametric case these are permitted to be $n^\psi$-consistent, where $\psi<1/2$, see Robinson1988 and references therein.
While our theorems below establish desired asymptotic properties for the one step estimates, from a practical point of view more iterations may be desirable. In fact, these also improve the statistical rate of convergence to the target PMLE, yielding an even faster statistical counterpart to the famous quadratic numerical rate of convergence of Newton estimates, see for example Theorem 2 of Robinson1988 and p. 312-313 of Ortega1970. We examine this issue in more detail in the next section and also the Monte Carlo study.
The following assumptions are discussed in lee2002consistency,lee2004asymptotic, and Gupta2013,Gupta2018, amongst other spatial papers in which they are routinely employed. These conditions are by no means the weakest possible set, but we opt for tractability to convey the main message especially in view of the large number of spatial parameters involved. For example, stochastic regressors can be easily accommodated but complicate the notation.
Let $\Psi_n$ be an $s{\times}(p_n+k_n)$ matrix of constants with full row-rank. The claims of the following theorems also hold when $p_n$ and $k_n$ are fixed, but we state and prove the results for the more challenging case when these diverge.
In the `just identified' case $p_n=r_n$, condition ((ref)) is implied by ((ref)). Theorem (ref) $(i)$ shows that the one-step estimate asymptotically achieves the efficiency bound noted by lee2002consistency. On the other hand, Theorem (ref) $(ii)$ yields the same distributional result as for the OLS estimate (Theorem 4.3 of Gupta2013). This should come as no surprise since lee2002consistency has already established the efficiency of OLS under suitable conditions. Nevertheless, Theorem (ref) $(ii)$ imposes weaker conditions on the relative rates of $h_n$ and $n^{\frac{1}{2}}$ than those extant in the literature.
Indeed, for their result, Gupta2013 assumed ${n^{\frac{1}{2}}p_n^{\frac{1}{2}}}/{h_n}\rightarrow 0$ as compared to our ${n^{\frac{1}{2}}p_n^{\frac{5}{2}}}/{h_n^3}\rightarrow 0$. The latter is a quantity of smaller order as $\left({n^{\frac{1}{2}}p_n^{\frac{5}{2}}}/{h_n^3}\right)/\left({n^{\frac{1}{2}}p_n^{\frac{1}{2}}}/{h_n}\right)=(p_n/h_n)^2\rightarrow 0$. For fixed $p_n$ and $k_n$, our asymptotic normality result relies only on $n^{\frac{1}{2}}/h_n^{3}\rightarrow 0$, as $n\rightarrow\infty$. This is a weaker requirement as compared to lee2002consistency, who assumed $n^{\frac{1}{2}}/h_n\rightarrow 0$ as $n\rightarrow\infty$. The reason for these favourable outcomes is the cancellation of higher order terms when using the one-step approximation. The key difference is in the rates $\left\Arrowvert n^{-1}\left[B_n, 0\right]'u\right\Arrowvert= \mathscr{O}_{\mathpzc{p}}\left({p_n^{\frac{1}{2}}}/{h_n}\right)$ and $\left\Arrowvert\phi_n\right\Arrowvert=\mathscr{O}_{\mathpzc{p}}\left({p_n^{\frac{1}{2}}}/{n^{\frac{1}{2}}h_n^{\frac{1}{2}}}\right),$ the latter being sharper since $n/h_n\rightarrow\infty$ as $n\rightarrow\infty$.
To more transparently illustrate the implications of our weaker rate conditions, consider data collected in a `farmer-district' type of environment, such as in case1991spatial. Suppose that there are $D$ districts, each containing $m$ farmers, so that $n=Dm$, and $D,m\rightarrow\infty$ simultaneously. There is independence across districts, but equal dependence within districts, yielding $h_n=m-1$ (see lee2002consistency,lee2004asymptotic for a more detailed discussion). Then, with fixed $p_n$, lee2002consistency and Gupta2013 required $n^{\frac{1}{2}}/h_n=D^{\frac{1}{2}}/m^{\frac{1}{2}}=o(1)$, while our condition imposes $n^{\frac{1}{2}}/h_n^3=D^{\frac{1}{2}}/m^{\frac{5}{2}}=o(1)$. Thus, our condition permits $D$ to grow much faster as we only need $D^{\frac{1}{5}}=o(m)$ as compared to $D=o(m)$. The author thanks an anonymous referee for suggesting this illustration. We note that robinson2010efficient obtained asymptotic normality, indeed efficiency, in a semiparametric setup with $p_n=1$ requiring only $h_n\rightarrow\infty$ if the disturbances are symmetrically distributed or the weight matrix is symmetric. This condition would likely need to be suitably amended as $p_n\rightarrow\infty$.
If $h_n$ is bounded as $n\rightarrow\infty$, a more complicated analysis is required to establish that one-step estimates achieve the PMLE asymptotic covariance matrix, because the information equality does not hold asymptotically. Denote $\mu_l=\mathbb{E}\left(u_i^l\right)$ for natural numbers $l$, and introduce, with $i,j=1,\ldots,p_n$, the $p_n \times p_n$ matrix $\Omega_{\lambda\lambda,n}$ with $(i,j)$-th element $ \frac{4\mu_3}{n\sigma_0^4}\sum_{r=1}^n c_{rr,in}b_{r,jn}X_n\beta_{0n}+\frac{\left(\mu_4-3\sigma_0^4\right)}{n\sigma_0^4}\sum_{r=1}^n c_{rr,in}c_{rr,jn} $ and the $k_n \times p_n$ matrix $\Omega_{\lambda\beta,n}$ with $i$-th column $ \frac{2\mu_3}{n\sigma_0^4}\sum_{r=1}^n c_{rr,in}x_{r,n} $ where $c_{pq,in}$ is the $(p,q)$-th element of $C_{in}$, $b_{jn}=G_{jn}X_n\beta_{0n}$ with $t$-th element $b_{t,jn}$ ($j=1,\ldots,p_n$ and $t=1,\ldots,n$) and $x_{p,n}$ is the $p$-th column of $X_n'$. Define
Then $ \mathbb{E}\left(\xi_n\xi_n'\right)=n^{-1} \left(2\Xi_n+\Omega_n\right), $ where
When $h_n$ is bounded OLS cannot be consistent (see lee2002consistency), so the following theorem considers only initial IV estimates.
The rate condition ((ref)) can simplify depending on the value of $\delta$, i.e. the order of the finite moments assumed for $u_i$. As $\delta$ grows larger, the last term in the rate condition becomes redundant, indeed the numerator therein tends to $p_n^2k_n^2$ as $\delta\rightarrow\infty$, which is evidently dominated by the numerator of the other rate restriction. In the `farmer-district' setting discussed earlier, we have bounded $h_n=m-1$ in this case. To further illustrate the rate condition, suppose that we are in the just identified case $p_n=r_n$. Then ((ref)) requires $p_n^5k_n^3+p_n^4k_n^4+p_n^3k_n^5+\left(p_nk_n\right)^{2+\frac{8}{\delta}}=o(n)$. Then the term involving $\delta$ dominates the other three if $\delta\leq 8/3$.
As indicated earlier, further iterations on the Newton step can improve the rate of statistical convergence to the target as well as finite sample properties. To see this, let $\hat{\hat\theta}^{\ell}_n$ be the $\ell$-th Newton iteration towards the PMLE $\check\theta_n$. By Theorem 2 of Robinson1988, $\left\Vert\check\theta_n- \hat{\hat\theta}^{\ell+1}_n\right\Vert=\mathscr{O}_{\mathpzc{p}}\left(\left\Vert\check\theta_n- \hat{\hat\theta}_n\right\Vert^{2^\ell}\right)$, an identical bound holding also for $\tilde{\tilde\theta}^{\ell+1}_n$. A factor that depends on $\ell$ is suppressed in the stated stochastic bound, indicating that this is not uniform in $\ell$. Because the results of Gupta2018 and this paper show that one-step Newton estimates and $\check\theta_n$ are $n^{1/2}/\left(p_n+k_n\right)^{1/2}$-consistent, we have \[ \left\Vert\check\theta_n- \hat{\hat\theta}^{\ell+1}_n\right\Vert=\mathscr{O}_{\mathpzc{p}}\left(\left(n/\left(p_n+k_n\right)\right)^{-2^{\ell-1}}\right),\;\;\;\;\;\left\Vert\check\theta_n- \tilde{\tilde\theta}^{\ell+1}_n\right\Vert=\mathscr{O}_{\mathpzc{p}}\left(\left(n/\left(p_n+k_n\right)\right)^{-2^{\ell-1}}\right), \] thus yielding the rate at which the iterations approximate the target estimate in a statistical sense, pointwise in $\ell$.
We examine finite-sample performance of $\hat{\hat{\theta}}_n$ in this section, since the IV case entails a change in limiting distribution due to the Newton step and OLS requires divergent $h_n$ to be consistent. Following das2003 and the design in Gupta2013, define $W^*_{in}$ as the symmetric circulant matrix with first row
and take $W^c_{in}={\left\Vert W^*_{in}\right\Vert^{-1}} W^*_{in},$ where $\left\Vert W^*_{in}\right\Vert=\overline{\eta}\left(W^*_{in}\right)=2i,$ because $W^*_{in}$ is a symmetric, circulant matrix (see e.g. Davis1979 p. 73). Thus $W^c_{in}$ is also a symmetric circulant matrix with first row given by $w^*_{1j,in}/2i$. This is an example of spatial weight matrices with bounded $h_n$.
We now dispense with some $n$ subscripts for brevity. Our design generates $y=S^{-1}(X\beta+u)$ for sample sizes $n=200,400,800$ and $k=2$, with elements $x_{j1}$ and $x_{j2}$ of $X$ generated as iid replicates from a $U(0,1)$ distribution, $j=1,\ldots,n$. We generate the disturbance $u$ using two different distributions: $N(0,1)$ and $t_6$. PMLE becomes MLE under the first, while the second has heavier tails. Our experiments take $p=2,4,6$ for each of the described designs. We use a design with weights matrices given by $W^c_{i}$, $i=1,\ldots,p$. Finally, we set $\beta_1=1$, $\beta_2=0.5$ and $p=2:\lambda_1=0.4,\lambda_2=0.5;p=4: \lambda_1=0.3,\lambda_i=0.2,i=2,3,4;p=6: \lambda_i=0.15,i=1,\ldots,6.$ The choices of $\lambda_i$ satisfy the sufficient condition $\sum_{i=1}^p\left\vert\lambda_i\right\vert<1$ for invertibility of $S$.
With the aim of comparing initial IV estimates and MLE to Newton-step estimates, we first report three statistics: Monte Carlo mean, Monte Carlo mean squared error, and relative root Monte Carlo mean squared error, the latter being a straightforward ratio of the root MSE for IV and the iterated estimate. We also examine the use of more than one iteration in finite samples, and for this recall the notation $\hat{\hat\theta}^{\ell}_n$ for the $\ell$-th Newton iteration. Our results are reported for $\ell=1,3,6$. The set of instruments that we use for our initial estimates are the linearly independent columns of $Z=\left(W_{1}^cX,\ldots,W_p^cX,X\right)$.
In Tables (ref) and (ref), we report the Monte Carlo mean of our estimates for standard normal and $t_6$ errors, respectively. For standard normal errors, we notice that the initial IV estimate can be heavily biased but Newton iterations improve matters, sometimes spectacularly. Indeed, for $p=6$ and $n=200$ the performance of $\hat\theta_n$ can be appalling, with $\hat\lambda_5<0$. However after six Newton steps this has improved to 0.1216 and even three iterations lead to a significant improvement. The reduction of bias from Newton iterations is not a universal feature, however broadly speaking the Newton steps reduce bias in the estimates, even for smaller values of $p$. As the sample size increases the iterations converge substantially, with little to choose typically between $\hat{\hat\theta}^3_n$ and $\hat{\hat\theta}^6_n$ for $n=800$. However for $n<800$, we notice that three iterations usually do the job quite satisfactorily, especially when $p<6$.
For $t_6$ errors, Table (ref) paints a similar picture to Table (ref). Once again, the noticeable `rogue' estimate is for $\lambda_5$ when $p=6$ and $n=200$. Considering that all our simulations start from the same seed, this outlier may possibly be attributed to a bad draw. As in the normal errors case, results are quite stable for larger $n$ and smaller $p$, and typically show bias reduction due to Newton steps and near convergence after three iterations.
Tables (ref) and (ref) report mean squared error (MSE) for the IV estimates and iterated estimates with $N(0,1)$ and $t_6$ errors, respectively. As may be expected, MSE is very high for designs that combine the largest values of $p$ with the smallest values of $n$. The efficiency improvement due to the Newton step is apparent, with iterations leading to very clear improvements (i.e. reductions) in MSE. These gains can be spectacular in many cases, for example for the $\lambda_i$ estimates when $p=6$ and $n=800$. These patterns of improvement with iteration are similar for both error distributions but the magnitude of MSE is generally much larger for $t_6$ errors, which features heavier tails than the normal distribution.
In Tables (ref) and (ref), we report the ratio of the Monte Carlo root mean squared error of $\hat\theta_n$ to that of $\hat{\hat\theta}^\ell_n$, $\ell=1,3,6$, abbreviating this quantity to RRMSE. An RRMSE of two indicates that the RMSE of the IV estimate is twice that of the Newton iteration it is being compared to. Our results in Table (ref) show that Newton iterations can lead to tremendous finite sample gains in MSE. These gains are present in 100% of the cases considered, but are generally larger for the spatial parameters $\lambda_i$ than the regression parameters $\beta_i$.
We discuss the spatial parameter estimates first. Note that for greater sample sizes we have greater MSE gains, often the gains more than doubling from $n=200$ to $n=800$, and sometimes even tripling. As observed for the means in Table (ref), there is usually not much to choose from between the third and sixth iterations. With and $n=800$ we nearly always obtain Newton estimates with RMSE a quarter of that for IV, and occasionally even a fifth of the IV RMSE. In most cases three iterations are enough to achieve these superb gains.
These patterns for the $\lambda_i$ qualitatively repeat themselves when the errors are $t_6$, as seen in Table (ref). In this case when $n=800$ we achieve RMSE improvements over IV of a factor of 2.15 always when three iterations are carried out, with factors of three commonly seen and one case with nearly a fourfold improvement. The factors of efficiency improvement that we observe in our results can dominate similar precedents in other settings. Indeed, the greatest relative root MSE improvement that Robinson2005a finds in his fractional time series setting is $\sqrt{1/0.23}=2.085$ (see Table 4 of that paper).
Moving to the estimates of the regression parameters $\beta_1$ and $\beta_2$, in both Tables (ref) and (ref) we see almost universal improvement over IV. The exceptions are four cases out of a total of 54 in Table (ref), for the $t_6$ case. These RMSE gains are not as spectacular as for the $\lambda_i$, but are generally noticeably large as both $n$ and $p$ increase. Indeed, for $n=800$ we observe that the RMSE for the IV estimate can sometimes be almost one and a half times are large as the Newton iterations when $p=6$ and $n=800$. For $n\geq 400$, IV performs worse than the Newton iterations almost uniformly (there are only two exceptions for $t_6$ errors) over both $\beta_1$ and $\beta_2$, the values of $p$, the number of iterations and the error distribution. Thus there is evidence of the usefulness of Newton iterations even for the regression parameters, albeit the gains are greater for the spatial parameters.
Finally, we also present the RRMSE of MLE (denoted $\mathring\theta_n$) to our proposed iterated estimates in Table (ref) for $N(0,1)$ errors. Naturally, we anticipate MLE to outperform iterated IV estimates for smaller sample sizes and, because our iterations target the MLE limiting covariance matrix, a reasonable aim is to approach the RMSE of the MLE as $n$ grows larger. Indeed, we find that this is the case. Recall that our estimates are designed to approximate but not outperform MLE: the main focus of the paper is computational simplicity. Our estimates are available in closed form and can be computed much faster than those requiring grid search and inversion of an $n\times n$ matrix. Thus, approaching the MLE in RMSE as $n$ grows is an encouraging and desirable property of our estimates. Finally, we observe that the RMSE of $\mathring\beta_n$ is much closer to the iterated estimates than is the case for $\mathring\lambda_n$. For the latter, larger sample sizes are needed for the RRMSE to approach unity.
In this section we explore the performance of the Newton step estimates when the number of neighbours diverges with sample size, i.e. $h_n\rightarrow\infty$. This design, with diverging $h_n$, also allows us to study the performance of iterations on OLS starting values. For each $i=1,\ldots, p$, we generate a $n\times n$ matrix $W_{in}^*$ as $w^*_{rs,in}=\Phi \left(-d_{rs,i}\right) I\left(c_{rs,i}<n^{1/3}/100\right)$ if $r\neq s$, and $w^*_{rr,in}=0$, where $\Phi (\cdot )$ is the standard normal cdf, $d_{rs,i}\sim$iid $U[-3,3]$, and $c_{rs,i}\sim$iid $U[0,1]$. This construction generates $W_{in}^*$ with approximately $n^{1/3}\%$ (up to closest integer) nonzero elements. These $W_{in}^*$ are then symmetrized and normalized by spectral norm to ensure stability, yielding the final set of $W_{in}$ that we employ. The remaining design details are as in the previous subsection. To conserve space, we report results only for $N(0,1)$ errors.
Tables (ref) and (ref) display the Monte Carlo mean of $\hat\theta_n$, $\hat{\hat\theta}^{1}_n$, $\hat{\hat\theta}^{3}_n$, $\tilde\theta_n$, $\tilde{\tilde\theta}^{1}_n$ and $\tilde{\tilde\theta}^{3}_n$. Convergence of iterations is achieved after three Newton steps, so we do not report the sixth iteration as in the previous subsection. In fact, convergence is practically fully achieved by just a single iteration with the IV starting values $\hat\theta_n$, as Table (ref) indicates. Examining Table (ref) suggests that a third iteration has more influence for OLS starting values, but modestly so. Tables (ref) and (ref) report MSE for the same sets of estimates and we find a similar pattern: for IV starting values one iteration seems to do the job and reduces MSE. On the other hand, for OLS starting values the first iteration increases MSE but the third iteration reduces it, following which performance is stable and so we do not report further iterations.
In Table (ref) we report RRMSE of the estimates studies above. We notice that IV estimates improve in MSE with a single Newton step, and subsequent iterations do not help much, because convergence is achieved. On the other hand, when starting with OLS values $\tilde\theta_n$, further iterations are beneficial and yield more efficient estimates. Convergence is completely achieved after three iterations in this case. We also find that Newton steps, whether they commence from $\hat\theta_n$ or $\tilde\theta_n$, give greater efficiency gains for the spatial parameters $\lambda_i$ rather than the regression coefficients $\beta_i$. This matches the results in the previous subsection. Because the $\lambda_i$ correspond to the potentially endogenous spatial lags $W_{in}y$, we might expect initial estimates of these to have greater potential for improvement compared to the $\beta_i$.
In this design, we confirm the robustness of our findings to heteroskedasticity in the error distribution. We generate the errors using multiplicative heteroskedasticity via the regressors, and report only the bounded $h_n$ weight matrices of Section (ref) and designs with Gaussian errors to conserve space. Specifically, we employ a $N(0,h_{jn})$ distribution for the errors, where $h_{jn}=n\left(\sum_{r=1}^n\left(\left\vert x_{r1}\right\vert+\left\vert x_{r2}\right\vert\right)\right)^{-1}\left(\left\vert x_{j1}\right\vert+\left\vert x_{j2}\right\vert\right)$, see Liu2015 and also Lin2010. Monte Carlo mean, mean squared error and RRMSE of IV estimates to iterated Newton step estimates are presented in Tables (ref)-(ref). We find the same qualitative patterns as were observed for the homoskedastic designs presented earlier, with the `rogue' IV estimate for $\lambda_5$ appearing again because we start our simulations from the same seed. As far as quantitative results are concerned, the improvements due to the Newton step are generally smaller than the homoskedastic case but still substantial.
In this small empirical illustration we show that the Newton step estimates perform well in practice and can lead to more precise estimation. The example is based on kolympiris2011spatial (KKM), and is also studied in Gupta2013. KKM seek to model the venture capital funding (provided by venture capital firms (VCFs)) for dedicated biotechnology firms (DBFs) with a SAR model. The hypothesis is that the level of VC funding for a DBF increases with the number of VCFs located in close proximity. Denoting by $d_{lk}$ the distance in miles between the $l$-th and $k$-th DBFs, we estimate
where $W_{i}^b$ is the (row-normalised) weight matrix having off-diagonal $(l,k)$-th element equal to 1 if $i-1< d_{lk}\leq i$, $i=1\ldots,p,$ and if $d_{lk}=0$ for $i=1$. Thus the matrices are based on each one of $p$ sequential 1-mile rings from the origin DBF. $y$ is the vector of natural logs of the amount of VC funding (million \$) received by each of $n=816$ DBFs.
We first focus on estimates of the main parameters of interest $\lambda_i$ in ((ref)). We estimate ((ref)) with $p=2,4,6$ using initial IV and the Newton-step estimates that we have justified theoretically. We only report the Newton-step for a single iteration as convergence is achieved. Like Gupta2013, we find that only $\lambda_1$ and $\lambda_2$ are statistically significant at the 1% level, and the magnitude of our parameter estimates is also close to their findings, with our results reported in Table (ref). The table reports $t$ statistics in parentheses. In square brackets we report for each parameter estimate the ratio of IV standard error to Newton-step standard error, and find that this difference can be as great as 12.53%. Thus the iteration scheme we propose can lead to more accurate inference in practice as the estimates are more precise.
As far as the $\beta_i$ are concerned, our simulations generally show that the efficiency gains are smaller for these as compared to the $\lambda_i$. Table (ref) reports standard error ratios and absolute t-statistics for exclusion tests and confirms this. Indeed, all standard error ratios are very close to unity and the t-statistics are practically identical. We note that our proposed iteration does not make the estimation precision of the $\beta_i$ worse and improves the estimation precision of the $\lambda_i$, leading to an overall improvement in estimation quality.
We give a very brief description of the explanatory variables in $X_n$ and refer the reader to KKM for details. The covariates include the number of proximate VCFs and DBFs to capture the effects of being in areas of high VCF or DBF concentration. Firm-specific characteristics include the distance from each DBF to its funding VCFs, , the average age of each funding VCF, exposure of VCFs through syndication and an indicator for foreign VCF investment. Variables controlling for DBF-specific factors include firm age, dummies for receiving a grant and being in an R&D tax credit state, a cost of business index for the DBF's home state, distance to the closest university and the number of non-biotech establishments in the DBF's zip code. Two further variables recognize that additional factors can affect the cost of doing business in ways that influence the VC funding levels of a given DBF.
\setcounter{table}{0}