title: Efficient closed-form estimation of large spatial autoregressions
authors: Abhimanyu Gupta
text: 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  equationsection
  tablesection
  figuresection
  
  
  
   
  
  
  
  
  
  
  
   
    
    .eps,.bmp,.png
    .emfbmp
  
    
    
    
  
  
  
  
  
  theoremTheorem[section]
  lemmaLemma[section]
  proposition[theorem]Proposition
  corollary[lemma]Corollary
  conjecture[theorem]Conjecture
  remark
  assumptionAssumption
  
  definition[1][Definition]
  
    #1
  example[1][Example]
  
    #1
  remark[1][Remark]
  
    #1
  
  
  
  
  
  
  
  
  
  OT1pzcmit
  
  
  
  
  
  
  
  
  
  
  
  
  
  acknowledgement[theorem]Acknowledgement
  algorithm[theorem]Algorithm
  axiom[theorem]Axiom
  case[theorem]Case
  claim[theorem]Claim
  conclusion[theorem]Conclusion
  condition[theorem]Condition
  
  
  
  
  
  
  
  
  
  
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  =9999        
  =9999
  
  
  
  
  
  
  
  
  
  Efficient closed-form estimation of large spatial autoregressionsI am grateful to the associate editor, three referees and Peter Robinson for very helpful comments that improved the paper significantly, and to Christos Kolympiris for permission to use the data in the empirical illustration.
      Abhimanyu GuptaDepartment of Economics, University of Essex, Wivenhoe Park, Colchester CO4 3SQ, UK. E-mail: a.gupta@essex.ac.uk. Research supported by ESRC grant ES/R006032/1. 
      December 30, 2023
  ===================================================================================================================================================================================================================================================================================================
  
  
  
  
  Newton-step approximations to pseudo maximum likelihood estimates of spatial autoregressive models with a large number of parameters are examined, in the sense that the parameter space grows slowly as a function of sample size. These have the same asymptotic efficiency properties as maximum likelihood under Gaussianity but are of closed form. Hence they are computationally simple and free from compactness assumptions, thereby avoiding two notorious pitfalls of implicitly defined estimates of large spatial autoregressions. When commencing from an initial least squares estimate, the Newton step can also lead to weaker regularity conditions for a central limit theorem than some extant in the literature. A simulation study demonstrates excellent finite sample gains from Newton iterations, especially in large multiparameter models for which grid search is costly. A small empirical illustration shows improvements in estimation precision with real data.
  
  
  Keywords: Spatial autoregression, efficiency, many parameters, networks
  
  JEL Classification: C21, C31, C33, C36 
  
  
  
  
  
  
  
  § INTRODUCTION
  
  Spatial autoregressive (SAR) models, introduced by <cit.>, 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× 1 vectors y_n and u of responses and unobserved disturbances, respectively, and an n× k covariate matrix X_n, the SAR model is 
  
      y_n=∑_i=1^pλ_0inW_iny_n+X_nβ_0n+u,
  
  where the elements of the n× n spatial weight matrices W_in are inverse economic distances and λ_0n=(λ_01n,…,λ_0pn)' and β_0n are unknown parameter vectors. Subscripting with n permits treatment of triangular arrays, an important issue for spatial models in general (see <cit.>), 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 (<cit.>, <cit.>), industrial organization (<cit.>), trade (<cit.>) and peer effects (<cit.>), to name only a few examples. Another frequently used approach to model cross-sectional dependence is the `common factor' technique, see e.g. <cit.> for a review.
  
  Estimation of SAR models has long been considered in the regional science literature, see e.g. <cit.>. Rigorous asymptotic theory for instrumental variables (IV) estimation was initially provided by <cit.>, leading to the present flourishing theoretical literature. <cit.> 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 <cit.>, a seminal contribution that provided a taxonomical asymptotic theory for Gaussian pseudo maximum likelihood estimates (PMLE) of SAR models. Recently <cit.> 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(V_1n,…,V_pn), where V_in is m_i× m_i and ∑_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× n block diagonal matrix with the m_i× m_i dimensional i-th diagonal block given by V_in. This approach has been suggested by <cit.>. Here and more generally, the specification (<ref>) is termed a `higher-order' SAR model if p>1, see e.g <cit.>, <cit.>, <cit.>, <cit.>, <cit.>. 
  
  In the study of higher-order SAR models, <cit.> 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,…,p, we have n=mp and the results of <cit.> imply that p→∞ is necessary for consistent estimation, analogous to the problems created in the spatial statistics literature by `infill asymptotics', see e.g. <cit.>. 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 <cit.> 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× n matrix in every iteration, leading to many ingenious solutions for faster computation, see e.g.  <cit.> and <cit.>. The computational cost of PMLE in SAR type models is particularly salient as data sets increase in size, as stressed by <cit.>. 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. <cit.> 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× n matrix, as well as risk of convergence to local optima. Furthermore, the requirement of a compact parameter space for λ_0n can severely restrict the admissible parameter values (see <cit.>). 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 <cit.>). 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 <cit.> and <cit.>. 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. <cit.> (simultaneous equations), <cit.> (nonlinear least squares), <cit.> (M-estimation),  <cit.> (generalized least squares), <cit.>, <cit.>, <cit.> (time series and adaptive estimation),  <cit.> (generalized method of moments), <cit.>, <cit.> (structural estimation), <cit.> (generalized linear models) and <cit.> (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 . The approach has been studied in the case of fixed-dimensional SAR models by <cit.> and <cit.>, 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× 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 <cit.> and <cit.>, 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 ‖ A‖=(η(A'A))^1/2, with η(·) and η(·) denoting the largest and smallest eigenvalues, respectively, of  a symmetric positive semidefinite matrix. Note that if A is a vector then ‖ A‖ is simply its Euclidean norm. Let ‖ A‖_R denote the maximum absolute row sum norm of A. For any parameter τ, function f(τ) and generic estimate τ̌, we will write f̌≡ f(τ̌). We denote true parameter values with 0 subscript and suppress the argument for a quantity evaluated at a true parameter value, i.e. f(τ_0)≡ f.
  
  
  § APPROXIMATIONS TO GAUSSIAN PMLE
  
  The (-2/n times) log pseudo Gaussian likelihood function for model (<ref>) at any admissible point θ=(λ',β')' is given by
  
      𝒬_n(θ,σ^2)=log(2πσ^2)-2/nlog|S_n(λ)|
       +  1/nσ^2(S_n(λ)y_n-X_nβ)'(S_n(λ)y_n-X_nβ),
   
  where S_n(λ)=I_n-∑_i=1^p_nλ_inW_in, with I_n denoting the n× n identity matrix.  If S_n is invertible, (<ref>) admits the reduced form y_n=S_n^-1X_nβ_n+S_n^-1u, and we  define R_n=A_n+B_n, where A_n = (G_1nX_nβ_n,…,G_p_nnX_nβ_n),B_n = (G_1nu,…,G_p_nnu),
  G_in(λ)=W_inS_n^-1(λ), i=1,…,p_n, and so R_n=(W_1ny_n,…,W_p_nny_n). 
  
  Defining ℛ_n^y(θ)=R_nλ_n+X_nβ_n-y_n, the derivative of (<ref>) at any admissible (θ,σ^2) is
  
      ξ_n(θ,σ^2)=(φ_n'(θ,σ^2), 2σ^-2ℛ_n^y'(θ)X_n)',
   
  where
  φ_n (θ,σ^2)=2σ^-2(σ^2tr G_1n(λ)+y_n'W_1n'ℛ_n^y(θ),…,σ^2tr G_p_nn(λ)+y_n'W_p_nn'ℛ_n^y(θ))'.
  Because ℛ_n^y=-u, denoting 
  ϕ_n   =  σ_0^-2( σ_0^2 tr C_1n-u'C_1nu ,…,σ_0^2 tr C_p_n-u'C_p_nnu  )' with C_in=G_in+G_in', we obtain
  
      ξ_n≡∂𝒬_n /∂θ=(ϕ_n ',0)'-2σ_0^-2t_n,
  
  with t_n  = [A_n,X_n]'u. The Hessian at any admissible point in the parameter space is
  
  
      H_n(θ_n,σ^2)≡∂^2𝒬_n(θ_n,σ^2)/∂θ∂θ'=([ 2/n P_ji,n(λ_n)+2/nσ^2R_n'R_n                 2/nσ^2R_n'X_n;                              ;                 2/nσ^2X_n'R_n                 2/nσ^2X_n'X_n ])
  
  where P_ji,n(λ_n) is the p_n× p_n matrix with (i,j)-th element given by tr (G_jn(λ_n)G_in(λ_n)). 
  
  
  
  
  
  
  
  
  
  
  
  
  
   Let Z_n be an n× r_n matrix of instruments, with r_n≥ p_n, and  define the IV and OLS estimates as 
  
      θ̂_n    =    Q̂_n^-1K̂_n'J_n^-1k̂_n,        σ̂^2_n=‖ y_n-(R_n,X_n)θ̂_n‖^2, 
      θ̃_n    =    L̂^-1_nl̂_n ,                    σ̃^2_n=‖ y_n-(R_n,X_n)θ̃_n‖^2,
   
  respectively, with
  Q̂_n=K̂_n'J_n^-1K̂_n, 
  K̂_n=[Z_n, X_n]'[R_n,X_n],k̂_n=[Z_n , X_n]'y_n, J_n=[Z_n ,X_n]'[Z_n,X_n], and 
  L̂_n=[R_n, X_n]'[R_n,X_n],l̂_n=[R_n,X_n]'y_n.
  Define the respective `one-step' estimates θ̂̂̂_n and θ̃̃̃_n by the following equations
  
      θ̂̂̂_n    =    θ̂_n-Ĥ_n^-1ξ̂_n, 
      θ̃̃̃_n    =    θ̃_n-H̃_n^-1ξ̃_n.
   
  We observe that other initial estimates, such as the GMM estimates of <cit.> and <cit.>, 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^ψ-consistent, where ψ<1/2, see <cit.> 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 <cit.> and p. 312-313 of <cit.>. We examine this issue in more detail in the next section and also the Monte Carlo study. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § ASYMPTOTIC PROPERTIES
  
  The following assumptions are discussed in <cit.>, and <cit.>, 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.
  
  u=(u_1,…,u_n)' has iid elements with zero mean and finite variance σ_0^2.
  
  
  For i=1,…,p_n, the elements of W_in are uniformly (1/h_n), where h_n is some positive sequence which may be bounded or divergent, but always bounded away from zero and such that n/h_n→∞ as . The diagonal elements of each W_in are zero. 
  
  
  S_n is non-singular for all sufficiently large n.
  
  
  
  ‖ S_n^-1‖_R, ‖ S_n'^-1‖_R, ‖ W_in‖_R and ‖ W'_in‖_R are uniformly bounded in n and i for all i=1,…,p_n and sufficiently large n.
  
   
  The elements of X_n are constants and are uniformly bounded in n, in absolute value, for all sufficiently large n.
  
  
  The elements of Z_n are constants and are uniformly bounded in absolute value, for all sufficiently large n.
  
  
  _n→∞η(J_n)<∞ and _n→∞η(K_n'K_n)>0.
  
  
  _n→∞η(J_n)>0 and _n→∞η(K_n'K_n)<∞.
  
  
  _n→∞η(L_n)<∞.
  
  
  _n→∞η(L_n)>0.
  
  
  𝔼( u_i^4)≤ C for i=1,…,n.
  
  Let Ψ_n be an s×(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. 
  
  
  
  
    (i) Let Assumptions <ref>-<ref> hold along with
  
      1/p_n+1/r_n+1/k_n+p_n^3k_n^4/n+p_n^3/2k_n/h_n→ 0   as ,
  
  and
  
      r_n^2/n bounded as .
  
  Then
  
      n^1/2/(p_n+k_n)^1/2Ψ_n(θ̂̂̂_n-θ_0n)d⟶ N(0,lim_σ_0^2/p_n+k_nΨ_n L_n^-1Ψ_n'),
  
  where the asymptotic covariance matrix exists, and is positive definite, by Assumptions <ref> and <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
    (ii) Let Assumptions <ref>-<ref> and <ref>-<ref> hold. Suppose also that
  
      1/p_n+1/k_n+p_n^3k_n^4/n+p_n^3/2k_n/h_n+n^1/2p_n^5/2/h_n^3→ 0.
  
  
  Then
  
      n^1/2/(p_n+k_n)^1/2Ψ_n(θ̃̃̃_n-θ_0n)d⟶ N(0,lim_σ_0^2/p_n+k_nΨ_n L_n^-1Ψ_n'),
  
  where the asymptotic covariance matrix exists, and is positive definite, by Assumptions <ref> and <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  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 <cit.>. On the other hand, Theorem <ref> (ii) yields the same distributional result as for the OLS estimate (Theorem 4.3 of <cit.>). This should come as no surprise since <cit.> 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^1/2 than those extant in the literature. 
  
  Indeed, for their result, <cit.> assumed n^1/2p_n^1/2/h_n→ 0 as compared to our n^1/2p_n^5/2/h_n^3→ 0. The latter is a quantity of smaller order as (n^1/2p_n^5/2/h_n^3)/(n^1/2p_n^1/2/h_n)=(p_n/h_n)^2→ 0. For fixed p_n and k_n, our asymptotic normality result relies only on n^1/2/h_n^3→ 0, as . This is a weaker requirement as compared to <cit.>, who assumed n^1/2/h_n→ 0 as . 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 n^-1[B_n, 0]'u= (p_n^1/2/h_n)
  and ϕ_n=(p_n^1/2/n^1/2h_n^1/2),
  the latter being sharper since n/h_n→∞ as . 
  
  To more transparently illustrate the implications of our weaker rate conditions, consider data collected in a `farmer-district' type of environment, such as in <cit.>. Suppose that there are D districts, each containing m farmers, so that n=Dm, and D,m→∞ simultaneously. There is independence across districts, but equal dependence within districts, yielding h_n=m-1 (see <cit.> for a more detailed discussion). Then, with fixed p_n, <cit.> and <cit.> required n^1/2/h_n=D^1/2/m^1/2=o(1), while our condition imposes n^1/2/h_n^3=D^1/2/m^5/2=o(1). Thus, our condition permits D to grow much faster as we only need D^1/5=o(m) as compared to D=o(m). The author thanks an anonymous referee for suggesting this illustration. We note that <cit.> obtained asymptotic normality, indeed efficiency, in a semiparametric setup with p_n=1 requiring only h_n→∞ if the disturbances are symmetrically distributed or the weight matrix is symmetric. This condition would likely need to be suitably amended as p_n→∞.
  
  If h_n is bounded as , 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 μ_l=𝔼(u_i^l) for natural numbers l, and introduce, with i,j=1,…,p_n, the p_n × p_n matrix Ω_λλ,n with (i,j)-th element
  4μ_3/nσ_0^4∑_r=1^n c_rr,inb_r,jnX_nβ_0n+(μ_4-3σ_0^4)/nσ_0^4∑_r=1^n c_rr,inc_rr,jn
  and the k_n × p_n matrix Ω_λβ,n with i-th column
  2μ_3/nσ_0^4∑_r=1^n c_rr,inx_r,n
  where c_pq,in is the (p,q)-th element of C_in, b_jn=G_jnX_nβ_0n with t-th element b_t,jn (j=1,…,p_n and t=1,…,n) and x_p,n is the p-th column of X_n'. Define
  
      Ω_n=([  Ω_λλ,n Ω_λβ,n';        ;  Ω_λβ,n       0 ]).
  
  Then 
  𝔼(ξ_nξ_n')=(2Ξ_n+Ω_n),
  where 
  
      Ξ_n=𝔼(H_n)=([ 2/n(P_ji,n+P_j'i,n+1/σ_0^2A_n'A_n)                    2/nσ_0^2A_n'X_n;                                   ;                    2/nσ_0^2X_n'A_n                    2/nσ_0^2X_n'X_n ]).
  
  When h_n is bounded OLS cannot be consistent (see <cit.>), so the following theorem considers only initial IV estimates.
  
  Let Assumptions <ref>-<ref> hold. Suppose that h_n is bounded away from zero and that there is a real number δ>0 such that
  𝔼 u_i^4+δ≤ C
  for i=1,…,n. In addition, assume that
  
      _η(2Ξ_n^-1+Ξ_n^-1Ω_nΞ_n^-1)<∞, _η(2Ξ_n^-1+Ξ_n^-1Ω_nΞ_n^-1)>0 and _η(Ξ_n)>0.
  
  Suppose also that 
  
      1/p_n+1/r_n+1/k_n+p_n^3k_n^2(p_n(r_n+k_n)+k_n^2)/n+(p_nk_n)^2+8/δ/n→ 0  as .
  
  Then
  
      n^1/2/(p_n+k_n)^1/2Ψ_n(θ̂̂̂_n-θ_0n)d⟶ N(0,lim_σ_0^2/p_n+k_nΨ_n (2Ξ_n^-1+Ξ_n^-1Ω_nΞ_n^-1)Ψ_n'),
  
  where the asymptotic covariance matrix exists, and is positive definite, by (<ref>).
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  The rate condition (<ref>) can simplify depending on the value of δ, i.e. the order of the finite moments assumed for u_i. As δ grows larger, the last term in the rate condition becomes redundant, indeed the numerator therein tends to p_n^2k_n^2 as δ→∞, 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+(p_nk_n)^2+8/δ=o(n). Then the term involving δ dominates the other three if δ≤ 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 θ̂̂̂^ℓ_n be the ℓ-th Newton iteration towards the PMLE θ̌_n.
  By  Theorem 2 of <cit.>, ‖θ̌_n- θ̂̂̂^ℓ+1_n‖=(‖θ̌_n- θ̂̂̂_n‖^2^ℓ), an identical bound holding also for θ̃̃̃^ℓ+1_n. A factor that depends on ℓ is suppressed in the stated stochastic bound, indicating that this is not uniform in ℓ. Because the results of <cit.> and this paper show that one-step Newton estimates and θ̌_n are n^1/2/(p_n+k_n)^1/2-consistent, we have
  
      ‖θ̌_n- θ̂̂̂^ℓ+1_n‖=((n/(p_n+k_n))^-2^ℓ-1),     ‖θ̌_n- θ̃̃̃^ℓ+1_n‖=((n/(p_n+k_n))^-2^ℓ-1),
   thus yielding the rate at which the iterations approximate the target estimate in a statistical sense, pointwise in ℓ.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § FINITE-SAMPLE PERFORMANCE OF NEWTON-STEP ESTIMATES
  
  
  
   §.§ Fixed number of neighbours (bounded h_n)
  We examine finite-sample performance of θ̂̂̂_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 <cit.> and the design in <cit.>, define W^*_in as the symmetric circulant matrix with first row
  
      w^*_1j,in={[                               0      if  j=1  or  j=i+2,…,n-i ;;                               1 if  j=2,…,i+1  or  j=n-i+1,…,n,;                                 ].
   
  and take W^c_in=‖ W^*_in‖^-1 W^*_in,
  where ‖ W^*_in‖=η(W^*_in)=2i, because W^*_in is a symmetric, circulant matrix (see e.g. <cit.> 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β+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,…,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,…,p. Finally, we set β_1=1, β_2=0.5 and p=2:λ_1=0.4,λ_2=0.5;p=4: λ_1=0.3,λ_i=0.2,i=2,3,4;p=6: λ_i=0.15,i=1,…,6. The choices of λ_i satisfy the sufficient condition ∑_i=1^p|λ_i|<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 θ̂̂̂^ℓ_n for the ℓ-th Newton iteration. Our results are reported for ℓ=1,3,6. The set of instruments that we use for our initial estimates are the linearly independent columns of Z=(W_1^cX,…,W_p^cX,X).
  
  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 θ̂_n can be appalling, with λ̂_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 θ̂̂̂^3_n and θ̂̂̂^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 λ_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 λ_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 θ̂_n to that of θ̂̂̂^ℓ_n, ℓ=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 λ_i than the regression parameters β_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 λ_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 <cit.>  finds in his fractional time series setting is √(1/0.23)=2.085 (see Table 4 of that paper).  
  
  Moving to the estimates of the regression parameters β_1 and β_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 λ_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≥ 400, IV performs worse than the Newton iterations almost uniformly (there are only two exceptions for t_6 errors) over both β_1 and β_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 θ_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× 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 β_n is much closer to the iterated estimates than is the case for λ_n. For the latter, larger sample sizes are needed for the RRMSE to approach unity.
  
  
   §.§ Growing number of neighbours (divergent h_n)
  
  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→∞. This design, with diverging h_n, also allows us to study the performance of iterations on OLS starting values. For each i=1,…, p, we generate a n× n matrix W_in^* as w^*_rs,in=Φ(-d_rs,i) I(c_rs,i<n^1/3/100) if r≠ s, and w^*_rr,in=0, where Φ (· ) is the standard
  normal cdf, d_rs,i∼iid U[-3,3], and c_rs,i∼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 θ̂_n, θ̂̂̂^1_n,  θ̂̂̂^3_n, θ̃_n, θ̃̃̃^1_n and θ̃̃̃^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 θ̂_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 θ̃_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 θ̂_n or θ̃_n, give greater efficiency gains for the spatial parameters λ_i rather than the regression coefficients β_i. This matches the results in the previous subsection. Because the λ_i correspond to the potentially endogenous spatial lags W_iny, we might expect initial estimates of these to have greater potential for improvement compared to the β_i.
   
  
   §.§ Heteroskedastic errors
  
  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(∑_r=1^n(| x_r1|+| x_r2|))^-1(| x_j1|+| x_j2|), see <cit.> and also <cit.>. 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 λ_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.
  
  
  
  
  § EMPIRICAL ILLUSTRATION
  
  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 <cit.> (KKM), and is also studied in <cit.>. 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
  
      y=∑_i=1^pλ_i W_i^by+Xβ+U,
  
   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≤ i, i=1…,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 λ_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 <cit.>, we find that only λ_1 and λ_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 β_i are concerned, our simulations generally show that the efficiency gains are smaller for these as compared to the λ_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 β_i worse and improves the estimation precision of the λ_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.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
   
  
         
   
  
   
   
  
  
   
  
   
  
   
  
   
  
   
  
  
  
   
  
  
  
       
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
  
  
   
  
  
  
  
  
  
   §.§ Proofs of theorems
  
  Write a_n=p_n+k_n, b_n=r_n+k_n, c_n=p_nk_n^2+k_n and τ_n=n^1/2/a_n^1/2
  . 
  
  
  
  
  
    (i) By the mean value theorem (<ref>) implies that
  
      θ̂̂̂_n-θ_0n    =    ( I_a_n-Ĥ_n^-1H_n )(θ̂_n-θ_0n)-Ĥ_n^-1ξ_n
          =    θ̂_n-θ_0n-Ĥ_n^-1H_n (θ̂_n-θ_0n)-Ĥ_n^-1ξ_n
  
  where H_n=∂^2𝒬_n(θ_n,σ̂_n^2)/∂θ∂θ'
  and θ_n^-θ_0n≤θ̂_n-θ_0n, with each row of the Hessian matrix evaluated at possibly different θ_n. The latter point is a technical comment that we take as given in the remainder of the paper whenever a mean-value theorem is applied to vector of values. For any s×1 vector α, we can use (<ref>) to write 
  
      τ_nα'Ψ_n(θ̂̂̂_n-θ_0n)     =    τ_nα'Ψ_nĤ_n^-1(Ĥ_n-H_n )(θ̂_n-θ_0n) 
          -    τ_nα'Ψ_nĤ_n^-1ξ_n,
  
  recalling that τ_n=n^1/2/a_n^1/2. The first term on RHS above has modulus bounded by
  τ_nαΨ_nĤ_n^-1Ĥ_n-H_nθ̂_n-θ_0n,
  where the second factor in norms is (a_n^1/2), the third is bounded for sufficiently large n by Lemma <ref>, by Lemma <ref> the fourth is (max{p_n^3/2b_n^1/2/n^1/2h_n,b_n^1/2c_n^1/2/n,p_n^1/2b_n^1/2/n^1/2h_n^1/2,b_n/n}) and the fifth is (b_n^1/2/n^1/2) by Theorem 3.1 of <cit.>.  We conclude that the first term on the RHS of (<ref>) is (max{p_n^3/2b/n^1/2h_n,b_nc_n^1/2/n,p_n^1/2b/n^1/2h_n^1/2,b_n^3/2/n}),
  which is negligible by (<ref>) and (<ref>) because
  p_n^3b_n^2/nh_n^2≤  C (p_n^3r_n^2+p_n^3k_n^2/nh_n^2)=(p_n^3/h_n^2r_n^2/n+p_n^3k_n^2/nh_n^2),
  b_nc_n^1/2/n≤  C (r_np_n^1/2k_n+p_n^1/2k_n^2/n^2)
  =(p_n^1/2k_n/n^1/2r_n/n^1/2+p_n^1/2k_n^2/n^2), p_nb_n^2/nh_n≤ C(p_nr_n^2+p_nk_n^2/nh_n)=(r_n^2/np_n/h_n+p_nk_n^2/n1/h_n), b_n^3/n^2≤ C(r_n^3+k_n^3/n^2). For the negligibility of the last term note that r_n^3/2/n=r_n^2/nr_n^-1/2.
  
  Thus we only need to find the asymptotic distribution of -τ_nα'Ψ_nĤ_n^-1ξ_n. We can write
  
      -τ_nα'Ψ_nĤ_n^-1ξ_n=2/σ_0^2τ_nα'Ψ_nĤ_n^-1t_n-τ_nα'Ψ_nĤ_n^-1ϕ_n.
  
  
  
  
  
  
  We have 
  𝔼ϕ_n^2  ≤∑_i=1^p_n𝔼(  tr C_in-σ_0^-2u'C_inu)^2 
   =  ∑_i=1^p_nvar(  u'C_inu)  = (p_n/nh_n),
   (see (A.20) in the proof of Theorem 3.3 and Lemma B.2 in <cit.>), so that 
  
      ϕ_n=(p_n^1/2/n^1/2h_n^1/2).
   
  Therefore the second term on the right of (<ref>) has modulus bounded by τ_n times
  
      αΨ_n Ĥ_n^-1 ϕ_n,
   
  where the second factor is (a_n^1/2), the third is bounded for sufficiently large n by Lemma <ref> and the last is (p_n^1/2/n^1/2h_n^1/2). Thus (<ref>) is (p_n^1/2a_n^1/2/n^1/2h_n^1/2) and the second term on the right of (<ref>) is (p_n^1/2/h_n^1/2) which is negligible by (<ref>).
  Then the asymptotic distribution required is that of
  
      2/σ_0^2τ_nα'Ψ_nĤ_n^-1t_n=∑_i=1^3Υ_in+τ_nα'Ψ_n L_n^-1t_n,
  
  
  Υ_1n  =  2/σ_0^2τ_nα'Ψ_nĤ_n^-1(Ĥ_n-H_n)H_n^-1t_n,
  Υ_2n  =  2/σ_0^2τ_nα'Ψ_nΞ_n^-1(H_n-Ξ_n)H_n^-1t_n,
  Υ_3n  =  τ_nα'Ψ_nL_n^-1(σ_0^2/2Ξ_n-L_n)(σ_0^2/2Ξ_n)^-1t_n.
  We will demonstrate that |Υ_in|=(1), i=1,2,3.
  
  
  
  
  First we observe that
  |Υ_1n|≤2/σ_0^2τ_n‖α‖‖Ψ_n‖‖Ĥ_n^-1‖‖Ĥ_n-H_n‖‖ H_n^-1‖‖ t_n ‖,
  
  
  
  
  where the second factor in norms is (a_n^1/2), the third and fifth are bounded for sufficiently large n by Lemma <ref>, the fourth is (max{p_n^3/2b_n^1/2/n^1/2h_n,b_n^1/2c_n^1/2/n,p_n^1/2b_n^1/2/n^1/2h_n^1/2,b_n/n}) by Lemma <ref>, and by (A.13) of <cit.> the last is (c_n^1/2/n^1/2) . 
  
  Then 
  |Υ_1n|=(max{p_n^3/2b_n^1/2c_n^1/2/n^1/2h_n,b_n^1/2c_n/n,p_n^1/2b_n^1/2c_n^1/2/n^1/2h_n^1/2,b_nc_n^1/2/n}),
  which is negligible by (<ref>) and (<ref>) because
  p_n^3b_nc_n/nh_n^2≤ C(p_n^4r_nk_n^2+p_n^4k_n^3/nh_n^2)=(r_n/n^1/2p_n^3/2k_n^2/n^1/2p_n^5/2/h_n^2+p_n^3k_n^3/np_n/h_n^2), b_nc_n^2/n^2≤ C(r_np_n^2k_n^4+p_n^2k_n^5/n^2)=(r_n/n^1/2p_n^2k_n^4/n^3/2+p_n^2k_n^4/nk_n/n), p_nb_nc_n/nh_n≤ C(r_np_n^2k_n^2+p_n^2k_n^3/nh_n)= 
  (r_n/n^1/2p_n^3/2k_n^2/n^1/2p_n^1/2/h_n+p_n^2k_n^3/nh_n) and b_nc_n^1/2/n=o(1) has been shown earlier.
  
  
  
  
  Next |Υ_2n|≤2σ_0^-2τ_n‖α‖‖Ψ_n‖‖ H_n^-1‖‖ H_n-Ξ_n‖‖Ξ_n^-1‖‖ t_n‖,
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  where the second factor in norms is (a_n^1/2), the third and fifth are bounded for sufficiently large n by Lemma <ref>, the fourth is (p_nk_n/n^1/2) by Lemma <ref> and the last is (c_n^1/2/n^1/2) as above. Then
  |Υ_2n|=(p_nk_nc_n^1/2/n^1/2)
  which is negligible by (<ref>) because p_n^2k_n^2c_n/n≤ Cp_n^3k_n^4/n.
  
  
  
  
  
  Similarly |Υ_3n|=(p_nc_n^1/2/h_n) by Lemma <ref>, which is negligible by (<ref>) because p_n^2c_n/h_n^2≤ Cp_n^3k_n^2/h_n^2. 
  Then we only need to find the asymptotic distribution of the last term term in (<ref>), but this is precisely the proof of Theorem 3.3 of <cit.>. Replicating those arguments leads to the theorem.
  
    (ii) In view of Lemmas <ref>, <ref> and <ref>, the theorem is proved exactly like Theorem <ref> (i), except for different orders of magnitudes of various expressions. In this case two of the orders will be different from the analogous ones considered in the the proof of Theorem <ref> (i). Indeed, the analogue of the bound for the first term in (<ref>) is 
  
      (n^1/2max{p_n^3/2c_n^1/2/n^1/2h_n,p_n^2/h_n^2,p_n^7/4c_n^1/4/n^1/4h_n^3/2,c_n/n}max{c_n^1/2/n^1/2,p_n^1/2/h_n})
          =    (max{π_1n,π_2n,π_3n,π_4n,π_5n,π_6n}),
  
  where
  π_1n  = p_n^3/2c_n/n^1/2h_n, π_2n  = p_n^2c_n^1/2/h_n^2, π_3n  = p_n^7/4c_n^3/4/n^1/4h_n^3/2,
  π_4n  = c_n^3/2/n, π_5n = n^1/2p_n^5/2/h_n^3, π_6n  = n^1/4p_n^9/4c_n^1/4/h_n^5/2. Now π_5n is assumed to tend to zero by (<ref>), while  the remaining π_in terms are also negligible by (<ref>) because  π_1n^2 =(p_n^5k_n^4/nh_n^2), π_2n^2 =(p_n^5k_n^2/h_n^4), π_3n^4 =((p_n^4k_n^3/n)(p_n^6k_n^3/h_n^6)), π_4n^2 =(p_n^3k_n^6/n^2), π_6n^2 =(n^1/2p_n^5k_n/h_n^5)=
  (π_5n(p_n^5/2k_n/h_n^2)). 
  
   The analogue for the bound on Υ_1n is of order 
  
      (n^1/2)(max{p_n^3/2c_n^1/2/n^1/2h_n,p_n^2/h_n^2,p_n^7/4c_n^1/4/n^1/4h_n^3/2,c_n/n})(c_n^1/2/n^1/2)
         =   (max{π_1n,π_2n,π_3n,π_4n}),
  
  which was shown to be negligible under the assumed conditions. All other bounds remain unchanged and will be also be negligible under (<ref>), as in the proof of Theorem <ref> (i).
  
  
  
  
  
  Proceeding as in the proof of Theorem <ref> (i), we can write
  
      τ_nα'Ψ_n(θ̂̂̂_n-θ_0n)     =    τ_nα'Ψ_nĤ_n^-1(Ĥ_n-H_n )(θ̂_n-θ_0n) 
          -    τ_nα'Ψ_n(Ĥ_n^-1-Ξ_n^-1)ξ_n-τ_nα'Ψ_nΞ_n^-1ξ_n.
  
   As in the proof of Theorem <ref> (i), the first term on the RHS above is negligible by (<ref>). Lemma <ref> (for bounded h_n) indicates that the second term on the RHS of (<ref>) is bounded in modulus by a constant times
  τ_n‖Ψ_n‖(‖ t_n‖+‖ϕ_n‖)(‖Ĥ_n-H_n‖+‖H_n-Ξ_n‖)
  which is 
  (n^1/2max{c_n^1/2/n^1/2,p_n^1/2/n^1/2h_n^1/2}max{p_n^3/2b_n^1/2/n^1/2h_n,b_n^1/2c_n^1/2/n,p_n^1/2b_n^1/2/n^1/2h_n^1/2,b_n/n,p_nk_n/n^1/2}),
  using (A.13) of <cit.>, (<ref>) and Lemmas <ref> and <ref> (i). This is negligible by (<ref>), in a similar way to the preceding proofs. Thus we need to establish the asymptotic distribution of -τ_nα'Ψ_nΞ_n^-1ξ_n,
  which is established under the assumed conditions in Theorem 3.4 of <cit.>.
  
  
  
   §.§ Lemmas
  
  
  In the subsequent lemmas the assumptions of the theorems that these are used to prove are taken to hold.
  
   (Lemma LS.4 of <cit.>, Supplementary Material) 
  ‖ B_n'A_n‖=‖ A_n'B_n‖=(n^1/2p_nk_n).
  
  
  
  (Lemma LS.4 of <cit.>, Supplementary Material)
  ‖ X_n'B_n‖=‖ B_n'X_n‖=(n^1/2p_n^1/2k_n^1/2).
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
  
  
  Ĥ_n-H_n and Ĥ_n-H_n are 
  (max{p_n^3/2b_n^1/2/n^1/2h_n,b_n^1/2c_n^1/2/n,p_n^1/2b_n^1/2/n^1/2h_n^1/2,b_n/n}).
  
  
  By the triangle inequality
  Ĥ_n-H_n≤Ĥ_n-H_n+H_n^^-H_n^, and again by the triangle inequality Ĥ_n-H_n is bounded by
  
      2/n‖ P_ji,n(λ̂_n)-P_ji,n‖  +  2/n|1/σ̂_n^2-1/σ_0^2|(‖ R_n'R_n‖ + 2‖ X_n'R_n‖+ ‖ X_n'X_n‖) .
  
  The first term in (<ref>) is bounded by
  
      {∑_i,j=1^p_n(2/n tr (G_jn(λ̂_n)G_in(λ̂_n)) - 2/n tr (G_jnG_in))^2}^1/2
  
   By the MVT,  we have tr (G_jn(λ̂_n)G_in(λ̂_n))= tr (G_jnG_in)+ζ_ij, n'(λ̂_n-λ_0n), where ζ_ij,n has elements tr(G_in(λ_n)G_sn(λ_n)G_jn(λ_n)
  +   G_sn(λ_n)G_in(λ_n)G_jn(λ_n)), s=1,…,p_n, and λ_n^-λ_0n≤λ̂_n^-λ_0n. Thus, the summands in (<ref>) are 4n^-2 (ζ_ij,n'(λ̂_n^-λ_0n))^2 ≤4n^-2 ‖ζ_ij,n‖^2‖λ̂_n^-λ_0n‖^2
  by Cauchy-Schwarz inequality, where the first factor in norms on the RHS is (p_nn^2/h_n^2) by Lemma LS.3 of <cit.>, supplementary material. The second factor is bounded by θ̂_n-θ_0n^2 =  (b_n/n) (see (A.6) of <cit.>),
  so we conclude that the summands in (<ref>) are ( b_np_n/nh_n^2) and therefore (<ref>) is (p_n^3/2b_n^1/2/n^1/2h_n) and it follows that so is the first term in (<ref>).
  By (A.7) of <cit.>, 
  
      |1/σ̂_n^2-1/σ_0^2|=(max{b_n^1/2c_n^1/2/n,p_n^1/2b_n^1/2/n^1/2h_n^1/2,b_n/n}),
  
  which handles the second factor in the second term in (<ref>). We shall now bound the terms inside the parentheses in the second term in (<ref>). These are (n) because n^-1/2‖ R_n‖=(1), n^-1/2‖ X_n‖=(1) and n^-1‖ X_n'R_n‖=(1), by Assumption <ref>. 
  From (<ref>), (<ref>), we conclude that
  
      Ĥ_n-H_n    =   (p_n^3/2b_n^1/2/n^1/2h_n)  + (max{b_n^1/2c_n^1/2/n,p_n^1/2b_n^1/2/n^1/2h_n^1/2,b_n/n})
  
  Similarly, it may be shown that H_n-H_n has the same order, whence the lemma follows.
  
  
   (Lemma B.2 of <cit.>)
  (i) H_n-Ξ_n=(p_nk_n/n^1/2) if  Assumption <ref> holds, 
  (ii) L_n-(σ_0^2/2)Ξ_n=(p_n/h_n).
  
  
  
   (Lemma B.3 of <cit.>) The following inequalities are satisfied: 
  Ĥ_n^-1≤ C H_n^-1≤ Clim_n→∞Ξ_n^-1≤ C(_n→∞η(L_n))^-1≤ C.
  If h_n does not diverge, the above result becomes 
  Ĥ_n^-1≤ C H_n^-1≤ C(_n→∞η(Ξ_n))^-1≤ C,
  if also _n→∞η(Ξ_n)>0.
  
  
   
  
  
  H̃_n-H_n and H̃_n-H_n are (max{c_n/n,p_n^2/h_n^2,p_n^3/2c_n^1/2/n^1/2h_n,p_n^7/4c_n^1/4/n^1/4h_n^3/2}).
  
  
  The proof is similar to that of Lemma <ref> and we only elaborate on the differences from that proof. In this case we need to bound
  
      2/n‖ P_ji,n(λ̃_n)-P_ji,n‖  +  2/n|1/σ̃_n^2-1/σ_0^2|(‖ R_n'R_n‖ + 2‖ X_n'R_n‖+ ‖ X_n'X_n‖) .
  
   In the OLS case we have σ̃_n^2-σ_0^2=(max{c_n/n,p_n/h_n^2,p_n^1/2c_n^1/2/n^1/2h_n})
  and θ̃_n-θ_0n=(max{c_n^1/2/n^1/2,p_n^1/2/h_n}),  from (A.23) and (A.21) of <cit.>, respectively. The first term in (<ref>) is then (max{p_n^3/2c_n^1/2/n^1/2h_n,p_n^2/h_n^2,p_n^7/4c_n^1/4/n^1/4h_n^3/2})
  while the second one is (max{c_n/n,p_n/h_n^2,p_n^1/2c_n^1/2/n^1/2h_n}).
  We may then argue in a similar way that the Hessian evaluated at the OLS estimate differs from its value at an intermediate point in norm by the same to conclude the proof.
  
  
   (Lemma B.3 of <cit.>) 
  H̃_n^-1≤ H_n^-1≤lim_n→∞Ξ_n^-1≤σ_0^2/2(_n→∞η(L_n))^-1≤ C.
  
  
  
  
  
  
  chicago
  
  
  
  