4 Maximum Likelihood Estimation: Theory

This chapter begins the study of estimation methods for spatial regression models. We focus first on maximum likelihood (ML), one of the earliest and most widely used approaches in spatial econometrics. The method accounts explicitly for the simultaneity generated by a spatially lagged dependent variable and, under a correctly specified probability model, provides a coherent framework for estimation and inference.

The chapter proceeds in three stages. We first explain why Ordinary Least Squares (OLS) is generally unsuitable for the Spatial Lag Model. We then derive the likelihood, score, concentrated likelihood, and Hessian for the Spatial Lag Model and the Spatial Error Model. Finally, we study the consistency and asymptotic distribution of the Gaussian quasi-maximum likelihood estimator (QMLE), drawing on the regular asymptotic framework of Lee (2004).

The principal asymptotic results developed in this chapter concern the regular case in which the information matrix is nonsingular and the estimator has the usual root-\(n\) rate. Lee’s broader theory also allows the spatial weights to become increasingly diffuse and shows that weaker rates may arise when the information about the spatial autoregressive parameter degenerates. We will explain that distinction after establishing the regular case.

4.1 Why OLS Fails in the Spatial Lag Model

Recall from Section 2.1.1 that the Spatial Lag Model contains the spatially lagged dependent variable \(\mathbf W_n\mathbf y_n\) as a regressor. Unlike an ordinary explanatory variable, this spatial lag is jointly determined with \(\mathbf y_n\). It is therefore generally correlated with the innovation vector, even when \(\mathbf W_n\) is fixed and the innovations are mutually independent.

4.1.1 A Zero-Mean Pure Spatial Lag Model

To isolate this simultaneity problem, consider the pure Spatial Lag Model without an intercept or additional explanatory variables: \[ \mathbf y_n = \rho_0\mathbf W_n\mathbf y_n + \boldsymbol\varepsilon_n, \tag{4.1}\] where \(\rho_0\) is the true spatial autoregressive parameter, \(\mathbf W_n\) is a nonrandom \(n\times n\) spatial weights matrix, and \(\mathbb E(\boldsymbol\varepsilon_n\mid\mathbf W_n)=\mathbf{0}\).

Define the spatial filter at the true parameter value by \[ \mathbf B_{n0} = \mathbf I_n-\rho_0\mathbf W_n. \tag{4.2}\] Assume that \(\mathbf B_{n0}\) is nonsingular. The reduced form is then \[ \mathbf y_n = \mathbf B_{n0}^{-1}\boldsymbol\varepsilon_n. \tag{4.3}\] This is the same reduced-form operation studied in Section 2.2, now written with an \(n\) subscript to make the sequence of spatial systems explicit.

Premultiplying Equation 4.3 by \(\mathbf W_n\) gives \[ \mathbf W_n\mathbf y_n = \mathbf G_n\boldsymbol\varepsilon_n, \tag{4.4}\] where \[ \mathbf G_n = \mathbf W_n\mathbf B_{n0}^{-1}. \tag{4.5}\] Thus, the regressor \(\mathbf W_n\mathbf y_n\) is itself a linear transformation of the innovation vector.

If OLS is applied to Equation 4.1 without an intercept, the estimator is \[ \widehat\rho_{\mathrm{OLS},n} = \left[ (\mathbf W_n\mathbf y_n)^{\top} (\mathbf W_n\mathbf y_n) \right]^{-1} (\mathbf W_n\mathbf y_n)^{\top}\mathbf y_n. \tag{4.6}\]

Substituting the structural equation \(\mathbf y_n=\rho_0\mathbf W_n\mathbf y_n+\boldsymbol\varepsilon_n\) into Equation 4.6 yields the exact sampling-error decomposition \[ \begin{aligned} \widehat\rho_{\mathrm{OLS},n}-\rho_0 &= \left[ (\mathbf W_n\mathbf y_n)^{\top} (\mathbf W_n\mathbf y_n) \right]^{-1} (\mathbf W_n\mathbf y_n)^{\top} \boldsymbol\varepsilon_n \\ &= \left( \boldsymbol\varepsilon_n^{\top} \mathbf G_n^{\top}\mathbf G_n \boldsymbol\varepsilon_n \right)^{-1} \boldsymbol\varepsilon_n^{\top} \mathbf G_n^{\top} \boldsymbol\varepsilon_n. \end{aligned} \tag{4.7}\]

Equation Equation 4.7 displays the central difficulty. Both the numerator and denominator are quadratic forms in the same innovation vector.

4.1.2 Endogeneity of the Spatial Lag

Suppose that, conditional on \(\mathbf W_n\), the innovations are independent with \[ \mathbb E(\varepsilon_{ni}\mid\mathbf W_n)=0, \qquad \mathbb E(\varepsilon_{ni}^2\mid\mathbf W_n)=\sigma_0^2. \]

Using the first-moment formula for quadratic forms from Lemma 3.10, we obtain \[ \begin{aligned} \mathbb E\left[ (\mathbf W_n\mathbf y_n)^{\top} \boldsymbol\varepsilon_n \mid\mathbf W_n \right] &= \mathbb E\left[ \boldsymbol\varepsilon_n^{\top} \mathbf G_n^{\top} \boldsymbol\varepsilon_n \mid\mathbf W_n \right] \\ &= \sigma_0^2\operatorname{tr}(\mathbf G_n^{\top}) \\ &= \sigma_0^2\operatorname{tr}(\mathbf G_n). \end{aligned} \tag{4.8}\] At the level of an individual spatial unit, \[ \mathbb E\left[ (\mathbf W_n\mathbf y_n)_i\varepsilon_{ni} \mid\mathbf W_n \right] = \sigma_0^2g_{n,ii}, \tag{4.9}\] where \(g_{n,ii}\) is the \(i\)th diagonal element of \(\mathbf G_n\).

The zero diagonal of \(\mathbf W_n\) does not generally imply a zero diagonal for \(\mathbf G_n\). When the Neumann-series condition in Lemma 2.4 holds, \[ \begin{aligned} \mathbf G_n &= \mathbf W_n (\mathbf I_n-\rho_0\mathbf W_n)^{-1} \\ &= \mathbf W_n + \rho_0\mathbf W_n^2 + \rho_0^2\mathbf W_n^3 + \cdots. \end{aligned} \tag{4.10}\]

Although \(\operatorname{diag}(\mathbf W_n)=\mathbf{0}\), the matrix \(\mathbf W_n^2\) commonly has nonzero diagonal entries. These entries represent feedback paths that leave a spatial unit, pass through another unit, and return to the original unit (see Section 1.5). Consequently, the innovation \(\varepsilon_{ni}\) affects neighboring outcomes and can feed back into the spatial lag used in the \(i\)th equation.

Endogeneity is not yet an exact bias formula

Equation Equation 4.8 establishes that the spatial lag is generally endogenous. It does not justify taking the random denominator in Equation 4.7 outside an expectation. Therefore, the nonzero expectation of the numerator alone does not provide an exact formula for \(\mathbb E(\widehat\rho_{\mathrm{OLS},n})\). The finite-sample bias depends on the joint distribution of the numerator and denominator. Inconsistency can instead be established from the probability limits of the two normalized quadratic forms.

4.1.3 Probability Limit of the OLS Estimator

The following proposition applies the quadratic-form Law of Large Numbers from the preceding chapter.

Proposition 4.1 (Probability Limit of OLS in the Pure Spatial Lag Model) For each \(n\), suppose that:

\(\boldsymbol\varepsilon_n\) has i.i.d. elements with mean zero, variance \(\sigma_0^2>0\), and finite fourth moment;
\(\mathbf W_n\) is nonrandom and \(\mathbf B_{n0}=\mathbf I_n-\rho_0\mathbf W_n\) is nonsingular;
the absolute row and column sums of \(\mathbf G_n=\mathbf W_n\mathbf B_{n0}^{-1}\) are uniformly bounded;
the limits \[ \tau_G = \lim_{n\to\infty} \frac{1}{n}\operatorname{tr}(\mathbf G_n) \] and \[ q_G = \lim_{n\to\infty} \frac{1}{n} \operatorname{tr} (\mathbf G_n^{\top}\mathbf G_n) >0 \] exist.

Then \[ \widehat\rho_{\mathrm{OLS},n} \overset{p}{\longrightarrow} \rho_0+\frac{\tau_G}{q_G}. \tag{4.11}\] Consequently, OLS is inconsistent whenever \(\tau_G\neq0\).

Proof. Start from the exact decomposition in Equation 4.7 and divide the numerator and denominator by \(n\): \[ \widehat\rho_{\mathrm{OLS},n}-\rho_0 = \frac{ n^{-1}\boldsymbol\varepsilon_n^{\top} \mathbf G_n^{\top} \boldsymbol\varepsilon_n }{ n^{-1}\boldsymbol\varepsilon_n^{\top} \mathbf G_n^{\top}\mathbf G_n \boldsymbol\varepsilon_n }. \tag{4.12}\]

Step 1: Probability limit of the numerator.

Uniform boundedness of the row and column sums of \(\mathbf G_n\) implies the same property for \(\mathbf G_n^{\top}\). Applying Corollary 3.3 with \(\mathbf A_n=\mathbf G_n^{\top}\) gives \[ \begin{aligned} \frac{1}{n} \boldsymbol\varepsilon_n^{\top} \mathbf G_n^{\top} \boldsymbol\varepsilon_n \overset{p}{\longrightarrow} \sigma_0^2 \lim_{n\to\infty} \frac{1}{n} \operatorname{tr}(\mathbf G_n^{\top}) = \sigma_0^2\tau_G. \end{aligned} \tag{4.13}\]

Step 2: Probability limit of the denominator.

By Lemma 3.6, the row and column sums of \(\mathbf G_n^{\top}\mathbf G_n\) are uniformly bounded. Applying Corollary 3.3 again, now with \(\mathbf A_n=\mathbf G_n^{\top}\mathbf G_n\), yields \[ \begin{aligned} \frac{1}{n} \boldsymbol\varepsilon_n^{\top} \mathbf G_n^{\top}\mathbf G_n \boldsymbol\varepsilon_n \overset{p}{\longrightarrow} \sigma_0^2 \lim_{n\to\infty} \frac{1}{n} \operatorname{tr} (\mathbf G_n^{\top}\mathbf G_n) = \sigma_0^2q_G. \end{aligned} \tag{4.14}\] The assumption \(q_G>0\) and the condition \(\sigma_0^2>0\) imply that the limit is strictly positive.

Step 3: Probability limit of the ratio.

Because the denominator in Equation 4.12 converges in probability to the nonzero constant \(\sigma_0^2q_G\), the inversion result in Lemma 3.4 and the product rules in Theorem 3.5 imply \[ \begin{aligned} \widehat\rho_{\mathrm{OLS},n}-\rho_0 &\overset{p}{\longrightarrow} \frac{\sigma_0^2\tau_G}{\sigma_0^2q_G} \\ &= \frac{\tau_G}{q_G}. \end{aligned} \] Therefore, \[ \widehat\rho_{\mathrm{OLS},n} \overset{p}{\longrightarrow} \rho_0+\frac{\tau_G}{q_G}. \] This proves the result.

The proposition clarifies both the usual result and an important qualification. OLS is generally inconsistent because spatial feedback makes \(\tau_G\) nonzero. However, inconsistency is not universal. If the sequence of spatial weights is such that \[ \frac{1}{n}\operatorname{tr}(\mathbf G_n) \longrightarrow0, \] then OLS may be consistent, provided the denominator remains asymptotically nondegenerate. This qualification is central to the analysis in Lee (2002), where the behavior of OLS depends on how the spatial weights change with the sample size.

For the regular spatial designs most commonly used in cross-sectional applications, OLS should not be relied upon to estimate the coefficient on \(\mathbf W_n\mathbf y_n\). Maximum likelihood and instrumental-variable or GMM methods address the simultaneity explicitly. The latter methods are studied in the following chapter.

4.2 Monte Carlo Illustration of OLS Bias

We now illustrate the finite-sample behavior of OLS using a simple Monte Carlo experiment. The data-generating process is \[ \mathbf y_n = \rho_0\mathbf W_n\mathbf y_n + \boldsymbol\varepsilon_n, \qquad \boldsymbol\varepsilon_n \sim \mathcal N(\mathbf{0},\mathbf I_n), \] with \(\rho_0=0.7\). The spatial weights matrix is based on rook contiguity for a \(15\times15\) regular lattice, so \(n=225\), and is row-standardized.

The design holds \(\mathbf W_n\) fixed across replications. This isolates sampling variation arising from the innovation vector. In each replication, we solve the reduced form and estimate the zero-intercept regression of \(\mathbf y_n\) on \(\mathbf W_n\mathbf y_n\), matching the model in Equation 4.1.

# Fix the Monte Carlo design.
set.seed(123)
replications <- 500L
n_side       <- 15L
n            <- n_side^2
rho_0        <- 0.7
sigma_0      <- 1

# Construct a fixed row-standardized rook-contiguity weights matrix for a
# 15 x 15 lattice. The same spatial structure is used in every replication.
rook_neighbors <- spdep::cell2nb(nrow = n_side, ncol = n_side, type = "rook")
rook_weights   <- spdep::nb2listw(rook_neighbors, style = "W")
W              <- spdep::listw2mat(rook_weights)

# Compute the spatial multiplier once because W and rho_0 are fixed across
# Monte Carlo replications.
spatial_multiplier <- solve(diag(n) - rho_0 * W)

# In each replication, generate the innovations, solve the reduced form, and
# estimate the pure SLM by OLS without an intercept.
rho_hat_ols <- replicate(replications, {
  epsilon <- stats::rnorm(n, mean = 0, sd = sigma_0)
  y       <- drop(spatial_multiplier %*% epsilon)
  Wy      <- drop(W %*% y)
  drop(crossprod(Wy, y) / crossprod(Wy))
})

The following table reports the Monte Carlo mean, bias, standard deviation, and root mean squared error of the OLS estimator.

Table 4.1: Monte Carlo performance of OLS in the pure Spatial Lag Model

True_value	Mean_estimate	Bias	Standard_deviation	RMSE
0.7	1.035	0.335	0.054	0.339

Figure 4.1 shows the empirical sampling distribution. The solid vertical line marks the true value \(\rho_0=0.7\), while the dashed line marks the Monte Carlo mean.

Figure 4.1: Sampling distribution of the OLS estimator in a pure Spatial Lag Model

The simulation shows a substantial upward displacement of the OLS sampling distribution. This direction is consistent with Equation 4.11. In the present design, the spatial weights are nonnegative and \(\rho_0=0.7\). The feedback terms described in Equation 4.10 imply a positive trace contribution, so the asymptotic bias term \(\tau_G/q_G\) is positive. The Monte Carlo distribution is therefore centered above the true value.

Some replications also produce OLS estimates outside the admissible parameter space of the Spatial Lag Model. This occurs because OLS treats \(\mathbf W_n\mathbf y_n\) as an ordinary regressor and estimates an unrestricted slope coefficient. It does not impose the requirement that \(\mathbf B_n(\rho)=\mathbf I_n-\rho\mathbf W_n\) be nonsingular. Such values are valid outputs of the unconstrained OLS regression, but they do not necessarily define a valid reduced form for the spatial model.

4.3 Maximum Likelihood Estimation of the Spatial Lag Model

The simultaneity problem documented in Section 4.1 does not make estimation impossible. It means that the joint determination of \(\mathbf y_n\) and \(\mathbf W_n\mathbf y_n\) must be incorporated into the estimation criterion. Maximum likelihood does this by working with the joint distribution of the entire outcome vector rather than treating the spatial lag as an ordinary exogenous regressor.

Maximum likelihood estimation of spatial autoregressive models was developed by Ord (1975) and systematized by Anselin (1988). The likelihood used below is also the starting point of the Gaussian QMLE analysis in Lee (2004).

4.3.1 The Model and Its Spatial Transformation

Consider the Spatial Lag Model \[ \mathbf y_n = \rho\mathbf W_n\mathbf y_n + \mathbf X_n\boldsymbol\beta + \boldsymbol\varepsilon_n, \tag{4.15}\] where \(\mathbf y_n\) is \(n\times1\), \(\mathbf X_n\) is \(n\times k\), \(\boldsymbol\beta\) is \(k\times1\), and \(\mathbf W_n\) is a known \(n\times n\) spatial weights matrix. For the exact Gaussian likelihood, assume \[ \boldsymbol\varepsilon_n \mid (\mathbf X_n,\mathbf W_n) \sim \mathcal N \left( \mathbf 0, \sigma^2\mathbf I_n \right). \tag{4.16}\]

Define \[ \mathbf A_n(\rho) = \mathbf I_n-\rho\mathbf W_n \tag{4.17}\] and the innovation implied by a candidate parameter vector as \[ \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho) = \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta. \tag{4.18}\]

The parameter vector is \[ \boldsymbol\theta = \left( \boldsymbol\beta^{\top}, \rho, \sigma^2 \right)^{\top}. \]

The admissible parameter space for \(\rho\) must contain only values for which \(\mathbf A_n(\rho)\) is nonsingular. This is not merely a numerical restriction: without nonsingularity, the structural system does not determine a unique reduced form (see Section 2.2).

When \(\mathbf A_n(\rho)\) is nonsingular, \[ \mathbf y_n = \mathbf A_n(\rho)^{-1} \left( \mathbf X_n\boldsymbol\beta + \boldsymbol\varepsilon_n \right). \tag{4.19}\]

Consequently, \[ \begin{aligned} \mathbb E \left( \mathbf y_n \mid \mathbf X_n,\mathbf W_n \right) &= \mathbf A_n(\rho)^{-1} \mathbf X_n\boldsymbol\beta, \\ \operatorname{Var} \left( \mathbf y_n \mid \mathbf X_n,\mathbf W_n \right) &= \sigma^2 \mathbf A_n(\rho)^{-1} \mathbf A_n(\rho)^{-\top}. \end{aligned} \tag{4.20}\]

Thus, even though the innovations are conditionally independent, the elements of \(\mathbf y_n\) are generally dependent. The likelihood must account for this joint dependence.

4.3.2 Deriving the Joint Density

The observed random vector is \(\mathbf y_n\), whereas the simple Gaussian density in Equation 4.16 is written for \(\boldsymbol\varepsilon_n\). We therefore use the change-of-variables theorem. For fixed \(\mathbf X_n\) and \(\mathbf W_n\), the transformation from \(\mathbf y_n\) to \(\boldsymbol\varepsilon_n\) is \[ \boldsymbol\varepsilon_n = \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta. \]

Its Jacobian matrix is \[ \frac{ \partial\boldsymbol\varepsilon_n }{ \partial\mathbf y_n^{\top} } = \mathbf A_n(\rho), \tag{4.21}\] so the absolute Jacobian determinant is \[ \left| \det\mathbf A_n(\rho) \right|. \tag{4.22}\]

The density of the innovation vector is \[ f_{\boldsymbol\varepsilon} \left( \boldsymbol\varepsilon_n \mid \mathbf X_n,\mathbf W_n; \sigma^2 \right) = (2\pi\sigma^2)^{-n/2} \exp \left\{ - \frac{ \boldsymbol\varepsilon_n^{\top} \boldsymbol\varepsilon_n }{2\sigma^2} \right\}. \]

Applying the change-of-variables formula gives \[ \begin{aligned} f_{\mathbf y} \left( \mathbf y_n \mid \mathbf X_n,\mathbf W_n; \boldsymbol\theta \right) &= (2\pi\sigma^2)^{-n/2} \left| \det\mathbf A_n(\rho) \right| \\ &\quad\times \exp \left\{ - \frac{1}{2\sigma^2} \left[ \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta \right]^{\top} \left[ \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta \right] \right\}. \end{aligned} \tag{4.23}\]

Proposition 4.2 (Gaussian log-likelihood for the Spatial Lag Model) Under Equation 4.15 and Equation 4.16, the conditional log-likelihood is \[ \begin{aligned} \ell_n(\boldsymbol\theta) &= - \frac n2\log(2\pi) - \frac n2\log(\sigma^2) + \log \left| \det\mathbf A_n(\rho) \right| \\ &\quad- \frac{1}{2\sigma^2} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho)^{\top} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho). \end{aligned} \tag{4.24}\]

The first two terms are the familiar normal-regression terms. The quadratic term measures the fit of the transformed model. The additional log-determinant term adjusts the density for the transformation from \(\boldsymbol\varepsilon_n\) to \(\mathbf y_n\).

The determinant term is necessary because \(\rho\) enters not only the transformed innovations but also the transformation that maps \(\mathbf y_n\) into \(\boldsymbol\varepsilon_n\). When different candidate values of \(\rho\) are considered, they generate different transformations \(\mathbf A_n(\rho)\). Therefore, their quadratic terms cannot be compared as if all transformations had the same Jacobian. If the log-determinant were omitted, estimation would select \(\rho\) only by making \(\boldsymbol\varepsilon_n(\boldsymbol\beta,\rho)^{\top}\boldsymbol\varepsilon_n(\boldsymbol\beta,\rho)\) small. This would define a nonlinear least-squares criterion, not the Gaussian likelihood of the observed vector \(\mathbf y_n\). The factor \(\left|\det\mathbf A_n(\rho)\right|\) is precisely the adjustment required by the change-of-variables formula derived above.

The same likelihood can be recovered directly from the multivariate normal distribution in Equation 4.20. Indeed, \[ \det \left[ \sigma^2 \mathbf A_n(\rho)^{-1} \mathbf A_n(\rho)^{-\top} \right] = (\sigma^2)^n \det\mathbf A_n(\rho)^{-2}, \] which produces the same log-determinant adjustment.

4.3.3 Concentrating Out \(\boldsymbol\beta\) and \(\sigma^2\)

Direct maximization of Equation 4.24 is possible, but it is not necessary (see Chapter 5). For every admissible value of \(\rho\), the maximizing values of \(\boldsymbol\beta\) and \(\sigma^2\) are available in closed form. Substituting them into the likelihood reduces the numerical optimization to one dimension.

For fixed \(\rho\), differentiate Equation 4.24 with respect to \(\boldsymbol\beta\): \[ \frac{ \partial\ell_n(\boldsymbol\theta) }{ \partial\boldsymbol\beta } = \frac{1}{\sigma^2} \mathbf X_n^{\top} \left[ \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta \right]. \tag{4.25}\]

Setting this expression equal to zero gives \[ \widehat{\boldsymbol\beta}_n(\rho) = \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top} \mathbf A_n(\rho)\mathbf y_n. \tag{4.26}\]

Thus, conditional on \(\rho\), one regresses the transformed dependent variable \(\mathbf A_n(\rho)\mathbf y_n\) on the untransformed matrix \(\mathbf X_n\). This point distinguishes the SLM from the SEM considered later.

Differentiating with respect to \(\sigma^2\) yields \[ \frac{ \partial\ell_n(\boldsymbol\theta) }{ \partial\sigma^2 } = - \frac{n}{2\sigma^2} + \frac{ \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho)^{\top} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho) }{2\sigma^4}. \tag{4.27}\]

After substituting \(\widehat{\boldsymbol\beta}_n(\rho)\) and setting the derivative equal to zero, \[ \widehat\sigma_n^2(\rho) = \frac{1}{n} \widehat{\boldsymbol\varepsilon}_n(\rho)^{\top} \widehat{\boldsymbol\varepsilon}_n(\rho), \tag{4.28}\] where \[ \widehat{\boldsymbol\varepsilon}_n(\rho) = \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\widehat{\boldsymbol\beta}_n(\rho). \]

The divisor is \(n\), not \(n-k\), because Equation 4.28 is the Gaussian maximum likelihood estimator conditional on \(\rho\) rather than the usual unbiased OLS estimator.

Let \[ \mathbf M_n = \mathbf I_n - \mathbf X_n \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top} \tag{4.29}\] be the residual-maker matrix. Then \[ \widehat{\boldsymbol\varepsilon}_n(\rho) = \mathbf M_n\mathbf A_n(\rho)\mathbf y_n \] and \[ \widehat\sigma_n^2(\rho) = \frac{1}{n} \mathbf y_n^{\top} \mathbf A_n(\rho)^{\top} \mathbf M_n \mathbf A_n(\rho) \mathbf y_n. \tag{4.30}\]

A computationally useful representation follows from two auxiliary regressions. Define \[ \begin{aligned} \widehat{\boldsymbol\beta}_{O,n} &= \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top}\mathbf y_n, \\ \widehat{\boldsymbol\beta}_{L,n} &= \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top}\mathbf W_n\mathbf y_n. \end{aligned} \tag{4.31}\]

Because \(\mathbf A_n(\rho)\mathbf y_n=\mathbf y_n-\rho\mathbf W_n\mathbf y_n\), \[ \widehat{\boldsymbol\beta}_n(\rho) = \widehat{\boldsymbol\beta}_{O,n} - \rho\widehat{\boldsymbol\beta}_{L,n}. \tag{4.32}\]

Define the corresponding auxiliary residuals \[ \begin{aligned} \mathbf e_{O,n} &= \mathbf y_n - \mathbf X_n\widehat{\boldsymbol\beta}_{O,n}, \\ \mathbf e_{L,n} &= \mathbf W_n\mathbf y_n - \mathbf X_n\widehat{\boldsymbol\beta}_{L,n}. \end{aligned} \tag{4.33}\]

Then \[ \widehat{\boldsymbol\varepsilon}_n(\rho) = \mathbf e_{O,n} - \rho\mathbf e_{L,n} \tag{4.34}\] and \[ \widehat\sigma_n^2(\rho) = \frac{1}{n} \left( \mathbf e_{O,n}-\rho\mathbf e_{L,n} \right)^{\top} \left( \mathbf e_{O,n}-\rho\mathbf e_{L,n} \right). \tag{4.35}\]

Substituting the conditional maximizers into Equation 4.24 gives the concentrated log-likelihood \[ \begin{aligned} \ell_{c,n}(\rho) = - \frac n2 \left[ \log(2\pi)+1 \right] - \frac n2 \log \widehat\sigma_n^2(\rho) + \log \left| \det\mathbf A_n(\rho) \right|. \end{aligned} \tag{4.36}\]

The estimate \(\widehat\rho_n\) maximizes Equation 4.36 over the admissible parameter interval. The remaining estimates are then recovered as \[ \widehat{\boldsymbol\beta}_n = \widehat{\boldsymbol\beta}_n(\widehat\rho_n) \qquad\text{and}\qquad \widehat\sigma_n^2 = \widehat\sigma_n^2(\widehat\rho_n). \]

Profile-likelihood algorithm for the SLM

Construct \(\mathbf W_n\mathbf y_n\).
Regress \(\mathbf y_n\) on \(\mathbf X_n\) and retain \(\mathbf e_{O,n}\).
Regress \(\mathbf W_n\mathbf y_n\) on \(\mathbf X_n\) and retain \(\mathbf e_{L,n}\).
For each admissible \(\rho\), evaluate \(\widehat\sigma_n^2(\rho)\) from Equation 4.35 and then evaluate \(\ell_{c,n}(\rho)\).
Numerically maximize the one-dimensional profile criterion.
Recover \(\widehat{\boldsymbol\beta}_n\) and \(\widehat\sigma_n^2\) at \(\widehat\rho_n\).

The two auxiliary regressions are performed only once. During the numerical search, the residual sum of squares is a quadratic function of \(\rho\); the main spatially specific computational task is repeatedly evaluating the log-determinant.

4.3.4 The Score Vector

The score is needed both for numerical optimization and for the asymptotic analysis developed later. The following identities are useful: \[ \frac{\partial\mathbf A_n(\rho)}{\partial\rho} = -\mathbf W_n, \]

\[ \frac{\partial}{\partial\rho} \log \left| \det\mathbf A_n(\rho) \right| = -\operatorname{tr} \left[ \mathbf A_n(\rho)^{-1}\mathbf W_n \right], \tag{4.37}\] and \[ \frac{ \partial\boldsymbol\varepsilon_n(\boldsymbol\beta,\rho) }{ \partial\rho } = -\mathbf W_n\mathbf y_n. \tag{4.38}\]

Because the matrices \(\mathbf A_n(\rho)\) and \(\mathbf W_n\) commute, \(\mathbf A_n(\rho)^{-1}\mathbf W_n = \mathbf W_n\mathbf A_n(\rho)^{-1}\) whenever the inverse exists. Define \[ \mathbf G_n(\rho) = \mathbf W_n\mathbf A_n(\rho)^{-1}. \tag{4.39}\]

With the parameter order \(\boldsymbol\theta=(\boldsymbol\beta^{\top},\rho,\sigma^2)^{\top}\), the score vector is \[ \mathbf s_n(\boldsymbol\theta) = \frac{ \partial\ell_n(\boldsymbol\theta) }{ \partial\boldsymbol\theta } = \begin{pmatrix} \mathbf s_{\beta,n}(\boldsymbol\theta) \\[4pt] s_{\rho,n}(\boldsymbol\theta) \\[4pt] s_{\sigma^2,n}(\boldsymbol\theta) \end{pmatrix}, \] where \[ \mathbf s_{\beta,n}(\boldsymbol\theta) = \frac{1}{\sigma^2} \mathbf X_n^{\top} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho), \tag{4.40}\]

\[ s_{\rho,n}(\boldsymbol\theta) = - \operatorname{tr}\mathbf G_n(\rho) + \frac{1}{\sigma^2} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho)^{\top} \mathbf W_n\mathbf y_n, \tag{4.41}\] and \[ s_{\sigma^2,n}(\boldsymbol\theta) = - \frac{n}{2\sigma^2} + \frac{1}{2\sigma^4} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho)^{\top} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho). \tag{4.42}\]

To verify Equation 4.41, differentiate the quadratic term: \[ \begin{aligned} \frac{\partial}{\partial\rho} \left[ \boldsymbol\varepsilon_n^{\top} \boldsymbol\varepsilon_n \right] &= 2 \boldsymbol\varepsilon_n^{\top} \frac{\partial\boldsymbol\varepsilon_n}{\partial\rho} \\ &= -2 \boldsymbol\varepsilon_n^{\top} \mathbf W_n\mathbf y_n. \end{aligned} \]

Multiplication by \(-1/(2\sigma^2)\) therefore produces the positive second term in Equation 4.41.

4.3.4.1 The Score at the True Parameter

Let \[ \boldsymbol\theta_0 = \left( \boldsymbol\beta_0^{\top}, \rho_0, \sigma_0^2 \right)^{\top}, \qquad \mathbf A_{n0} = \mathbf A_n(\rho_0), \qquad \mathbf G_{n0} = \mathbf W_n\mathbf A_{n0}^{-1}. \]

At the true parameter, \[ \boldsymbol\varepsilon_n \left( \boldsymbol\beta_0, \rho_0 \right) = \boldsymbol\varepsilon_n. \]

Moreover, the reduced form implies \[ \mathbf W_n\mathbf y_n = \mathbf G_{n0}\mathbf X_n\boldsymbol\beta_0 + \mathbf G_{n0}\boldsymbol\varepsilon_n. \tag{4.43}\]

Substituting Equation 4.43 into Equation 4.41 gives \[ \begin{aligned} s_{\rho,n}(\boldsymbol\theta_0) = \frac{1}{\sigma_0^2} \left( \mathbf G_{n0}\mathbf X_n\boldsymbol\beta_0 \right)^{\top} \boldsymbol\varepsilon_n + \frac{1}{\sigma_0^2} \left[ \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n - \sigma_0^2 \operatorname{tr}\mathbf G_{n0} \right]. \end{aligned} \tag{4.44}\]

This decomposition is central to the asymptotic theory. The first component is a linear form in the innovations. The second is a centered quadratic form. Their moments are governed by the results in Section 3.9, and their joint limiting distribution is covered by Theorem 3.26 under the corresponding regularity conditions.

The score has conditional mean zero at the true parameter. For example, \[ \mathbb E \left[ \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n \mid \mathbf X_n,\mathbf W_n \right] = \sigma_0^2 \operatorname{tr}\mathbf G_{n0} \] by Lemma 3.10. Hence the trace term in Equation 4.44 exactly centers the quadratic form.

4.3.5 The Hessian and Observed Information

The Hessian records the curvature of the log-likelihood. It is used by Newton-type optimization algorithms and provides the observed-information estimator of the covariance matrix.

Continue to order the parameters as \((\boldsymbol\beta^{\top},\rho,\sigma^2)^{\top}\). Differentiating the score components gives \[ \frac{ \partial^2\ell_n }{ \partial\boldsymbol\beta \partial\boldsymbol\beta^{\top} } = - \frac{1}{\sigma^2} \mathbf X_n^{\top}\mathbf X_n, \tag{4.45}\]

\[ \frac{ \partial^2\ell_n }{ \partial\boldsymbol\beta \partial\rho } = - \frac{1}{\sigma^2} \mathbf X_n^{\top}\mathbf W_n\mathbf y_n, \tag{4.46}\]

\[ \frac{ \partial^2\ell_n }{ \partial\boldsymbol\beta \partial\sigma^2 } = - \frac{1}{\sigma^4} \mathbf X_n^{\top} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho), \tag{4.47}\]

\[ \frac{ \partial^2\ell_n }{ \partial\rho^2 } = - \operatorname{tr} \left[ \mathbf G_n(\rho)^2 \right] - \frac{1}{\sigma^2} \mathbf y_n^{\top} \mathbf W_n^{\top}\mathbf W_n \mathbf y_n, \tag{4.48}\]

\[ \frac{ \partial^2\ell_n }{ \partial\rho\partial\sigma^2 } = - \frac{1}{\sigma^4} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho)^{\top} \mathbf W_n\mathbf y_n, \tag{4.49}\] and \[ \frac{ \partial^2\ell_n }{ \partial(\sigma^2)^2 } = \frac{n}{2\sigma^4} - \frac{1}{\sigma^6} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho)^{\top} \boldsymbol\varepsilon_n(\boldsymbol\beta,\rho). \tag{4.50}\]

For completeness, the derivative of the trace term in Equation 4.48 is \[ \begin{aligned} \frac{\partial}{\partial\rho} \operatorname{tr}\mathbf G_n(\rho) &= \operatorname{tr} \left[ \frac{\partial\mathbf A_n(\rho)^{-1}}{\partial\rho} \mathbf W_n \right] \\ &= \operatorname{tr} \left[ \mathbf A_n(\rho)^{-1} \mathbf W_n \mathbf A_n(\rho)^{-1} \mathbf W_n \right] \\ &= \operatorname{tr} \left[ \mathbf G_n(\rho)^2 \right]. \end{aligned} \]

The full Hessian is therefore \[ \mathbf H_n(\boldsymbol\theta) = \begin{pmatrix} - \dfrac{\mathbf X_n^{\top}\mathbf X_n}{\sigma^2} & - \dfrac{\mathbf X_n^{\top}\mathbf W_n\mathbf y_n}{\sigma^2} & - \dfrac{\mathbf X_n^{\top}\boldsymbol\varepsilon_n}{\sigma^4} \\[12pt] - \dfrac{\mathbf y_n^{\top}\mathbf W_n^{\top}\mathbf X_n}{\sigma^2} & - \operatorname{tr}\left[\mathbf G_n(\rho)^2\right] - \dfrac{\mathbf y_n^{\top}\mathbf W_n^{\top}\mathbf W_n\mathbf y_n}{\sigma^2} & - \dfrac{\boldsymbol\varepsilon_n^{\top}\mathbf W_n\mathbf y_n}{\sigma^4} \\[12pt] - \dfrac{\boldsymbol\varepsilon_n^{\top}\mathbf X_n}{\sigma^4} & - \dfrac{\boldsymbol\varepsilon_n^{\top}\mathbf W_n\mathbf y_n}{\sigma^4} & \dfrac{n}{2\sigma^4} - \dfrac{\boldsymbol\varepsilon_n^{\top}\boldsymbol\varepsilon_n}{\sigma^6} \end{pmatrix}, \tag{4.51}\] where \(\boldsymbol\varepsilon_n=\boldsymbol\varepsilon_n(\boldsymbol\beta,\rho)\) in this display.

The observed information matrix is \[ \mathbf J_n(\widehat{\boldsymbol\theta}_n) = - \mathbf H_n(\widehat{\boldsymbol\theta}_n). \]

Under correct Gaussian specification and the regularity conditions developed later, a conventional covariance estimator is \[ \widehat{\operatorname{Var}} \left( \widehat{\boldsymbol\theta}_n \right) = \mathbf J_n(\widehat{\boldsymbol\theta}_n)^{-1}. \tag{4.52}\]

For Gaussian QMLE under nonnormal innovations, the inverse Hessian alone is not in general the correct covariance estimator. The asymptotic section will derive the corresponding sandwich form.

4.3.6 Computing the Spatial Log-Determinant

The term \[ \log \left| \det \left( \mathbf I_n-\rho\mathbf W_n \right) \right| \] must be evaluated repeatedly during optimization. Directly recomputing a dense \(n\times n\) determinant at every candidate value of \(\rho\) can be expensive.

Let \(\omega_{n1},\ldots,\omega_{nn}\) denote the eigenvalues of \(\mathbf W_n\). The determinant identity used by Ord (1975) is \[ \det \left( \mathbf I_n-\rho\mathbf W_n \right) = \prod_{i=1}^{n} \left( 1-\rho\omega_{ni} \right). \tag{4.53}\]

When the eigenvalues are real, this gives \[ \log \left| \det \left( \mathbf I_n-\rho\mathbf W_n \right) \right| = \sum_{i=1}^{n} \log \left| 1-\rho\omega_{ni} \right|. \tag{4.54}\]

This identity is exact; it is not an approximation. Its computational advantage is that the eigenvalues are calculated once and then reused at every evaluation of the profile likelihood.

If the spectrum is real and \(\omega_{n,\min}<0<\omega_{n,\max}\), a standard connected admissible interval containing zero is \[ \frac{1}{\omega_{n,\min}} < \rho < \frac{1}{\omega_{n,\max}}. \tag{4.55}\]

For a nonnegative row-standardized matrix, \(\omega_{n,\max}=1\). The lower endpoint nevertheless depends on the smallest eigenvalue and need not equal \(-1\). Restricting \(\rho\) to a compact subset of the admissible interval keeps the spatial multiplier away from singularity.

For a nonsymmetric spatial weights matrix, some eigenvalues may be complex. The product identity remains valid, but the term-by-term real logarithm in Equation 4.54 requires appropriate handling of conjugate pairs. Reliable software therefore uses methods designed for the structure of the weights matrix rather than applying a real logarithm mechanically to complex eigenvalues (see Section 2.2).

For large sparse systems, alternatives include sparse matrix factorizations, trace approximations, and the characteristic-polynomial approach of Smirnov and Anselin (2001). A broader computational discussion is provided by LeSage and Pace (2010). The important conceptual point is that these methods change how the log-determinant is evaluated, not the likelihood being maximized.

4.4 Maximum Likelihood Estimation of the Spatial Error Model

The Spatial Error Model (SEM) places spatial dependence in the disturbance rather than in the conditional mean (see Section 2.1.3). This changes both the consequences of OLS and the form of the transformation used for maximum likelihood.

4.4.1 Model and Consequences for OLS

Consider \[ \begin{aligned} \mathbf y_n &= \mathbf X_n\boldsymbol\beta + \mathbf u_n, \\ \mathbf u_n &= \lambda\mathbf W_n\mathbf u_n + \boldsymbol\varepsilon_n, \\ \boldsymbol\varepsilon_n \mid (\mathbf X_n,\mathbf W_n) &\sim \mathcal N \left( \mathbf 0, \sigma^2\mathbf I_n \right). \end{aligned} \tag{4.56}\]

Define \[ \mathbf B_n(\lambda) = \mathbf I_n-\lambda\mathbf W_n. \tag{4.57}\]

When \(\mathbf B_n(\lambda)\) is nonsingular, \[ \mathbf u_n = \mathbf B_n(\lambda)^{-1} \boldsymbol\varepsilon_n \] and therefore \[ \mathbf y_n = \mathbf X_n\boldsymbol\beta + \mathbf B_n(\lambda)^{-1} \boldsymbol\varepsilon_n. \tag{4.58}\]

Conditional on \(\mathbf X_n\) and \(\mathbf W_n\),

\[ \mathbb E(\mathbf u_n)=\mathbf 0 \] and \[ \operatorname{Var}(\mathbf u_n) = \sigma^2 \mathbf B_n(\lambda)^{-1} \mathbf B_n(\lambda)^{-\top}. \tag{4.59}\]

The covariance matrix is generally full. Spatially connected units can have correlated composite disturbances even though the primitive innovations are independent. Its diagonal elements can also differ, so the composite error \(\mathbf u_n\) may be heteroskedastic even when \(\boldsymbol\varepsilon_n\) is homoskedastic.

Suppose that \(\mathbf X_n\) is nonstochastic or that strict exogeneity holds. The OLS estimator is \[ \widehat{\boldsymbol\beta}_{\mathrm{OLS}} = \boldsymbol\beta + \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top}\mathbf u_n. \]

Because \(\mathbb E(\mathbf u_n\mid\mathbf X_n,\mathbf W_n)=\mathbf 0\), \[ \mathbb E \left( \widehat{\boldsymbol\beta}_{\mathrm{OLS}} \mid \mathbf X_n,\mathbf W_n \right) = \boldsymbol\beta. \tag{4.60}\]

Thus, unlike in the SLM, the problem is not simultaneity between a spatial lag of the dependent variable and the innovation. Under suitable large-sample conditions, OLS can also remain consistent. However, its conditional covariance matrix is \[ \begin{aligned} \operatorname{Var} \left( \widehat{\boldsymbol\beta}_{\mathrm{OLS}} \mid \mathbf X_n,\mathbf W_n \right) &= \sigma^2 \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top} \mathbf B_n(\lambda)^{-1} \mathbf B_n(\lambda)^{-\top} \mathbf X_n \\ &\quad\times \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1}. \end{aligned} \tag{4.61}\]

This is not the spherical-error formula \(\sigma^2(\mathbf X_n^{\top}\mathbf X_n)^{-1}\). Consequently:

OLS does not exploit the known covariance structure and is generally inefficient relative to GLS;
conventional homoskedastic OLS standard errors are invalid;
spatial dependence must be incorporated into the covariance estimator or into the estimation method.

The phrase “OLS is unbiased but inefficient” therefore requires the exogeneity condition and refers to the coefficient estimator. It does not imply that the usual OLS standard errors remain valid.

4.4.2 Gaussian Likelihood for the SEM

Let \[ \mathbf r_n(\boldsymbol\beta) = \mathbf y_n-\mathbf X_n\boldsymbol\beta. \]

The primitive innovation implied by the SEM is \[ \boldsymbol\varepsilon_n(\boldsymbol\beta,\lambda) = \mathbf B_n(\lambda) \mathbf r_n(\boldsymbol\beta). \tag{4.62}\]

The Jacobian matrix of the transformation from \(\mathbf y_n\) to \(\boldsymbol\varepsilon_n\) is \[ \frac{ \partial\boldsymbol\varepsilon_n }{ \partial\mathbf y_n^{\top} } = \mathbf B_n(\lambda). \]

Therefore, the conditional joint density is \[ \begin{aligned} f_{\mathbf y} \left( \mathbf y_n \mid \mathbf X_n,\mathbf W_n; \boldsymbol\beta,\lambda,\sigma^2 \right) &= (2\pi\sigma^2)^{-n/2} \left| \det\mathbf B_n(\lambda) \right| \\ &\quad\times \exp \left\{ - \frac{1}{2\sigma^2} \mathbf r_n(\boldsymbol\beta)^{\top} \mathbf B_n(\lambda)^{\top} \mathbf B_n(\lambda) \mathbf r_n(\boldsymbol\beta) \right\}. \end{aligned} \tag{4.63}\]

The log-likelihood is \[ \begin{aligned} \ell_n^{\mathrm{SEM}} \left( \boldsymbol\beta,\lambda,\sigma^2 \right) &= - \frac n2\log(2\pi) - \frac n2\log(\sigma^2) + \log \left| \det\mathbf B_n(\lambda) \right| \\ &\quad- \frac{1}{2\sigma^2} \mathbf r_n(\boldsymbol\beta)^{\top} \mathbf B_n(\lambda)^{\top} \mathbf B_n(\lambda) \mathbf r_n(\boldsymbol\beta). \end{aligned} \tag{4.64}\]

Define the precision factor \[ \boldsymbol\Omega_n(\lambda) = \mathbf B_n(\lambda)^{\top} \mathbf B_n(\lambda). \tag{4.65}\]

Apart from the factor \(1/\sigma^2\), \(\boldsymbol\Omega_n(\lambda)\) is the inverse of the covariance matrix in Equation 4.59. This explains why the quadratic term in Equation 4.64 has the Generalized Least Squared (GLS) form.

4.4.3 Concentrated Likelihood and GLS Interpretation

For fixed \(\lambda\), differentiate Equation 4.64 with respect to \(\boldsymbol\beta\). The first-order condition is \[ \mathbf X_n^{\top} \boldsymbol\Omega_n(\lambda) \left( \mathbf y_n- \mathbf X_n\boldsymbol\beta \right) = \mathbf 0. \]

Hence \[ \widehat{\boldsymbol\beta}_n(\lambda) = \left[ \mathbf X_n^{\top} \boldsymbol\Omega_n(\lambda) \mathbf X_n \right]^{-1} \mathbf X_n^{\top} \boldsymbol\Omega_n(\lambda) \mathbf y_n. \tag{4.66}\]

Because \(\boldsymbol\Omega_n(\lambda)=\mathbf B_n(\lambda)^{\top}\mathbf B_n(\lambda)\), this can also be written as \[ \widehat{\boldsymbol\beta}_n(\lambda) = \left[ \mathbf X_n(\lambda)^{\top} \mathbf X_n(\lambda) \right]^{-1} \mathbf X_n(\lambda)^{\top} \mathbf y_n(\lambda), \tag{4.67}\] where \[ \mathbf y_n(\lambda) = \mathbf B_n(\lambda)\mathbf y_n = \mathbf y_n- \lambda\mathbf W_n\mathbf y_n \] and \[ \mathbf X_n(\lambda) = \mathbf B_n(\lambda)\mathbf X_n = \mathbf X_n- \lambda\mathbf W_n\mathbf X_n. \]

For known \(\lambda\), the Gaussian ML estimator of \(\boldsymbol\beta\) is therefore the GLS estimator and can be computed as OLS using the transformed dependent variable and transformed regressors. This is often called the spatial Cochrane–Orcutt transformation.

The profile estimator of the innovation variance is \[ \widehat\sigma_n^2(\lambda) = \frac{1}{n} \widehat{\boldsymbol\varepsilon}_n(\lambda)^{\top} \widehat{\boldsymbol\varepsilon}_n(\lambda), \tag{4.68}\] where \[ \widehat{\boldsymbol\varepsilon}_n(\lambda) = \mathbf B_n(\lambda) \left[ \mathbf y_n- \mathbf X_n \widehat{\boldsymbol\beta}_n(\lambda) \right]. \]

Substitution into Equation 4.64 gives \[ \begin{aligned} \ell_{c,n}^{\mathrm{SEM}}(\lambda) = - \frac n2 \left[ \log(2\pi)+1 \right] - \frac n2 \log \widehat\sigma_n^2(\lambda) + \log \left| \det\mathbf B_n(\lambda) \right|. \end{aligned} \tag{4.69}\]

The minus sign before \(\frac n2\log\widehat\sigma_n^2(\lambda)\) is essential. A larger transformed residual variance lowers, rather than raises, the concentrated likelihood.

Profile-likelihood algorithm for the SEM

For each admissible candidate value of \(\lambda\):

Form \(\mathbf B_n(\lambda)\).
Transform both the outcome and regressors: \(\mathbf y_n(\lambda)=\mathbf B_n(\lambda)\mathbf y_n\) and \(\mathbf X_n(\lambda)=\mathbf B_n(\lambda)\mathbf X_n\).
Regress \(\mathbf y_n(\lambda)\) on \(\mathbf X_n(\lambda)\) to obtain \(\widehat{\boldsymbol\beta}_n(\lambda)\).
Compute \(\widehat\sigma_n^2(\lambda)\) from the transformed residuals.
Evaluate Equation 4.69.

A one-dimensional optimizer selects \(\widehat\lambda_n\). The estimates \(\widehat{\boldsymbol\beta}_n\) and \(\widehat\sigma_n^2\) are then evaluated at \(\widehat\lambda_n\).

This profile procedure does not require a separate manual alternating algorithm. The numerical optimizer supplies candidate values of \(\lambda\), and the closed-form conditional estimates are recalculated at each evaluation.

4.4.4 Score and Information Matrix for the SEM

Define \[ \mathbf G_{B,n}(\lambda) = \mathbf W_n\mathbf B_n(\lambda)^{-1}. \tag{4.70}\]

With parameter order \((\boldsymbol\beta^{\top},\lambda,\sigma^2)^{\top}\), the SEM score is \[ \mathbf s_n^{\mathrm{SEM}} = \begin{pmatrix} \dfrac{1}{\sigma^2} \mathbf X_n^{\top} \mathbf B_n(\lambda)^{\top} \boldsymbol\varepsilon_n \\[10pt] - \operatorname{tr}\mathbf G_{B,n}(\lambda) + \dfrac{1}{\sigma^2} \boldsymbol\varepsilon_n^{\top} \mathbf W_n \mathbf r_n(\boldsymbol\beta) \\[10pt] - \dfrac{n}{2\sigma^2} + \dfrac{1}{2\sigma^4} \boldsymbol\varepsilon_n^{\top} \boldsymbol\varepsilon_n \end{pmatrix}, \tag{4.71}\] where \(\boldsymbol\varepsilon_n=\mathbf B_n(\lambda)\mathbf r_n(\boldsymbol\beta)\).

Under correct Gaussian specification, the Fisher information matrix is \[ \mathcal I_n^{\mathrm{SEM}}(\boldsymbol\theta) = \begin{pmatrix} \dfrac{1}{\sigma^2} \mathbf X_n^{\top} \mathbf B_n(\lambda)^{\top} \mathbf B_n(\lambda) \mathbf X_n & \mathbf 0 & \mathbf 0 \\[12pt] \mathbf 0^{\top} & \operatorname{tr} \left[ \mathbf G_{B,n}(\lambda)^{\top} \mathbf G_{B,n}(\lambda) \right] + \operatorname{tr} \left[ \mathbf G_{B,n}(\lambda)^2 \right] & \dfrac{1}{\sigma^2} \operatorname{tr} \mathbf G_{B,n}(\lambda) \\[12pt] \mathbf 0^{\top} & \dfrac{1}{\sigma^2} \operatorname{tr} \mathbf G_{B,n}(\lambda) & \dfrac{n}{2\sigma^4} \end{pmatrix}. \tag{4.72}\]

Two features are worth emphasizing.

First, the \(\boldsymbol\beta\) block is orthogonal to the covariance parameters under Gaussian normality. This reflects the separation between the conditional mean and covariance structure of the SEM.

Second, the \(\lambda\lambda\) element contains \(\operatorname{tr}[\mathbf G_{B,n}(\lambda)^2]\), not the square of \(\operatorname{tr}\mathbf G_{B,n}(\lambda)\). These are different objects: \[ \operatorname{tr} \left[ \mathbf G_{B,n}^2 \right] \neq \left[ \operatorname{tr} \mathbf G_{B,n} \right]^2 \] in general.

Under correct Gaussian specification and regular identification, \[ \widehat{\operatorname{Var}} \left( \widehat{\boldsymbol\theta}_n \right) = \mathcal I_n^{\mathrm{SEM}} \left( \widehat{\boldsymbol\theta}_n \right)^{-1} \] is the information-matrix covariance estimator. As with the SLM, a QMLE based on a misspecified innovation distribution generally requires a sandwich covariance matrix.

4.5 Comparing the SLM and SEM Likelihoods

The two profile-likelihood procedures look similar because both contain a Gaussian fit term and a spatial log-determinant. Their economic and statistical roles are nevertheless different.

Feature	Spatial Lag Model	Spatial Error Model
Location of spatial dependence	Conditional mean of \(\mathbf y_n\)	Covariance of the composite disturbance
Structural equation	\(\mathbf A_n(\rho)\mathbf y_n=\mathbf X_n\boldsymbol\beta+\boldsymbol\varepsilon_n\)	\(\mathbf B_n(\lambda)(\mathbf y_n-\mathbf X_n\boldsymbol\beta)=\boldsymbol\varepsilon_n\)
OLS problem	Endogenous spatial lag; generally inconsistent	Generally inefficient; conventional standard errors invalid
Transformation for fixed spatial parameter	Transform \(\mathbf y_n\), not \(\mathbf X_n\)	Transform both \(\mathbf y_n\) and \(\mathbf X_n\)
Profile coefficient estimate	OLS of \(\mathbf A_n(\rho)\mathbf y_n\) on \(\mathbf X_n\)	OLS of \(\mathbf B_n(\lambda)\mathbf y_n\) on \(\mathbf B_n(\lambda)\mathbf X_n\)
Spatial determinant	\(\log\|\det\mathbf A_n(\rho)\|\)	\(\log\|\det\mathbf B_n(\lambda)\|\)

In both models, the determinant term is not a penalty added by convention. It is the Jacobian required by the transformation from the observed outcome vector to the primitive innovation vector.

The next section turns from construction of the likelihood to its large-sample behavior. The key tasks will be to establish identification, uniform convergence of the concentrated criterion, convergence of the Hessian, and a central limit theorem for the score. Those arguments will repeatedly invoke the matrix and linear–quadratic-form results developed in Chapter 3.

4.6 Asymptotic Properties of the Spatial Lag Model

The previous sections constructed the Gaussian likelihood for the Spatial Lag Model (SLM) and Spatial Error Model (SEM). This section studies the asymptotic properties of the SLM. The analysis follows Lee (2004). The estimator is called a maximum likelihood estimator when the innovations are Gaussian. It is called a Gaussian quasi-maximum likelihood estimator when the same Gaussian criterion is maximized even though the innovation distribution need not be normal.

The distinction matters for inference. Under correct Gaussian specification, the information equality holds and the inverse information matrix is the asymptotic covariance matrix. Under non-Gaussian innovations, the estimator can remain consistent and asymptotically normal, but its covariance matrix generally has the sandwich form.

Scope of the results developed in this section

The main results below cover the regularly identified SLM. In that case, the spatial parameter is identified partly through variation in the conditional mean, and the estimator has the usual root-\(n\) rate. This regular result may hold whether the sequence \(h_n\) introduced below is bounded or divergent. A divergent \(h_n\) does not, by itself, imply a slower rate. The slower rate found by Lee (2004) arises when the mean-based source of identification becomes asymptotically weak and the information about the spatial parameter is supplied mainly by quadratic forms.

4.6.1 Regularity Conditions

As discussed in Section 3.7 spatial asymptotics require a sequence of models indexed by \(n\). We therefore write \(\mathbf y_n = \rho_0\mathbf W_n\mathbf y_n + \mathbf X_n\boldsymbol\beta_0 +\boldsymbol\varepsilon_n\), where the dimension of every object changes with \(n\). Define \[ \mathbf A_n(\rho) = \mathbf I_n-\rho\mathbf W_n, \qquad \mathbf A_{n0} = \mathbf A_n(\rho_0), \] and \[ \mathbf G_{n0} = \mathbf W_n\mathbf A_{n0}^{-1}. \tag{4.73}\]

The following assumptions provide a transparent version of the regular case in Lee (2004). The number of regressors \(k\) is fixed as \(n\) increases.

Assumption 1 (ML1: Innovations)

For every \(n\), the elements of \(\boldsymbol\varepsilon_n\) are independent and identically distributed. Their common distribution does not depend on \(n\) and satisfies \[ \mathbb E(\varepsilon_{ni})=0, \qquad \mathbb E(\varepsilon_{ni}^2)=\sigma_0^2>0. \]

For some \(\gamma>0\), \[ \mathbb E \left( |\varepsilon_{ni}|^{4+\gamma} \right) <\infty. \]

The fourth moment is needed for variances of quadratic forms. The additional positive moment is used by the linear–quadratic central limit theorem.

Assumption 2 (ML2: Individual spatial weights)

The diagonal of \(\mathbf W_n\) is zero, and its entries satisfy \[ w_{n,ii}=0, \qquad \max_{i,j}|w_{n,ij}| = O\!\left(\frac{1}{h_n}\right), \] where \(h_n\) may be bounded or divergent.

Assumption 2 requires a deeper discussion. The sequence \(h_n\) controls the largest individual spatial weight. The condition \[ \max_{1\leq i,j\leq n} |w_{n,ij}| = O\!\left(\frac{1}{h_n}\right) \] in Lee (2004) means that there is a finite constant \(C\) such that \[ |w_{n,ij}| \leq \frac{C}{h_n} \] uniformly over \(i\), \(j\), and sufficiently large \(n\). This is the deterministic big-\(O\) notation defined in Definition 3.5, applied to a triangular array as discussed in Section 3.7.

Two cases are especially useful.

If \(h_n\) remains bounded, some individual weights need not shrink as the sample expands. This is typical when every unit has only a fixed number of neighbors.
If \(h_n\to\infty\), each individual weight becomes small. This is typical when the number of neighbors increases and the weights are normalized over an expanding neighborhood.

The interpretation is particularly transparent for equal-weight neighbors. If unit \(i\) has \(q_n\) neighbors and each receives weight \(1/q_n\), then \[ \max_j |w_{n,ij}| = \frac{1}{q_n}, \] so one can take \(h_n=q_n\). A fixed number of neighbors corresponds to bounded \(h_n\), whereas an increasing number of neighbors corresponds to divergent \(h_n\).

Inverse-distance weights without a cutoff behave differently because every spatial unit may receive a positive weight from every other unit. To see this, consider \(n\) units located at equally spaced points on a line. Let the distance between two adjacent units be \(d\). From the perspective of an interior unit, there are two units at distance \(d\), two units at distance \(2d\), two units at distance \(3d\), and so on. Let the unnormalized inverse-distance weights be \[ a_{n,ij} = \begin{cases} \dfrac{1}{d_{ij}}, & i\neq j,\\[4pt] 0, & i=j. \end{cases} \]

For an interior unit \(i\), the sum of the unnormalized weights therefore has the approximate form \[ s_{n,i} = \sum_{j\neq i}a_{n,ij} \approx \frac{2}{d} \left( 1+\frac{1}{2}+\frac{1}{3}+\cdots+\frac{1}{m_n} \right), \] where \(m_n\) is the number of distance levels available on either side of the unit. The expression in parentheses is a harmonic sum. It increases without bound and is of order \(\log n\). Units near the boundaries have neighbors on only one side for some distances, but their corresponding sums are still of order \(\log n\).

Consequently, the row sums of the unnormalized inverse-distance matrix do not remain uniformly bounded as the sample size increases. This explains why empirical applications frequently impose a cutoff distance \(d^*\), setting \(a_{n,ij}=0\) whenever \(d_{ij}>d^*\). On a regular lattice with fixed spacing, a fixed cutoff limits the number of neighbors of each unit and prevents the unnormalized row sums from increasing with \(n\).

Now consider the row-standardized inverse-distance matrix, \[ w_{n,ij} = \frac{a_{n,ij}}{s_{n,i}}. \]

Because \(s_{n,i}\) is of order \(\log n\), the largest normalized weight, which is assigned to a nearest neighbor, is of order \(\frac{1}{\log n}\). All other normalized weights are no larger. Hence, \[ \max_{i,j}|w_{n,ij}| = O\left(\frac{1}{\log n}\right), \] so the sequence in the regularity condition may be chosen as \[ h_n\asymp\log n. \]

Assumption 3 (ML3: Relative neighborhood growth)

The rate sequence satisfies \[ \frac{h_n}{n} \longrightarrow0. \]

The restriction \(\frac{h_n}{n}\longrightarrow0\) rules out interactions whose effective neighborhood size grows as fast as the sample itself. Lee shows that when \(h_n/n\) does not vanish, the Gaussian QMLE can fail to be consistent in some spatial designs.

In our previous example, the condition in Assumption 2 holds with \(h_n\asymp\log n\). Moreover, \[ \frac{h_n}{n} \asymp \frac{\log n}{n} \longrightarrow0, \] so Assumption 3 is also satisfied.This calculation verifies the rate of the individual weights and the relative growth condition. Row standardization also implies that the absolute row sums of \(\mathbf W_n\) equal one. It does not, however, automatically establish the uniform column-sum bounds required by Assumption 5. Those bounds, and the corresponding bounds for the spatial multiplier, require a separate argument.

Example 4.1 (Equal-weight neighborhoods) Suppose every unit has \(q_n\) neighbors and the matrix is row-standardized: \[ w_{n,ij} = \begin{cases} 1/q_n, & j\in\mathcal N_{n,i},\\ 0, & j\notin\mathcal N_{n,i}. \end{cases} \]

Then \(h_n=q_n\) is a natural rate sequence. If \(q_n=4\) for every \(n\), then \(h_n\) is bounded. If \(q_n=\lfloor n^{1/3}\rfloor\), then \[ h_n\to\infty, \qquad \frac{h_n}{n} \to0. \]

Assumption 4 (ML4: Existence of the true reduced form)

The matrix \(\mathbf A_{n0}=\mathbf I_n-\rho_0\mathbf W_n\) is nonsingular for every \(n\).

Assumption 5 (ML5: Uniform control of spatial dependence)

The sequences \(\{\mathbf W_n\}\) and \(\{\mathbf A_{n0}^{-1}\}\) are uniformly bounded in both absolute row and column sums.

Note that \[ \|\mathbf W_n\|_2^2 \leq \|\mathbf W_n\|_1 \|\mathbf W_n\|_{\infty} \] and \[ \|\mathbf A_{n0}^{-1}\|_2^2 \leq \|\mathbf A_{n0}^{-1}\|_1 \|\mathbf A_{n0}^{-1}\|_{\infty}. \]

Therefore, the uniform bounds on the absolute row and column sums imply that both sequences are uniformly bounded in spectral norm. These conditions prevent the spatial multiplier from producing uncontrolled cross-sectional amplification as \(n\) increases.

Assumption 6 (ML6: Regressors)

The elements of \(\mathbf X_n\) are nonrandom and uniformly bounded. In addition, \[ \frac{1}{n} \mathbf X_n^{\top}\mathbf X_n \longrightarrow \mathbf Q_{xx}, \] where \(\mathbf Q_{xx}\) is finite and positive definite.

This excludes asymptotic multicollinearity among the columns of \(\mathbf X_n\). A stochastic-regressor version is possible, but it requires conditional moment and probability-limit conditions that are not pursued here.

Assumption 7 (ML7: Uniform stability over the parameter space)

The parameter space \(\mathcal R\) for \(\rho\) is compact, \(\rho_0\) lies in its interior, and \(\mathbf A_n(\rho)\) is nonsingular for every \(\rho\in\mathcal R\). Moreover, \(\mathbf A_n(\rho)^{-1}\) is uniformly bounded in either absolute row sums or absolute column sums, uniformly over \(\rho\in\mathcal R\).

This condition keeps the admissible interval away from values at which the spatial multiplier becomes unstable. For a row-standardized nonnegative matrix, a compact subset of \((-1,1)\) is a common sufficient choice. For a matrix with real eigenvalues, the admissible set must remain inside the interval determined by the extreme eigenvalues, as discussed in Section 2.2.

To state the final condition, define the spatially generated regressor \[ \mathbf d_n = \mathbf G_{n0} \mathbf X_n \boldsymbol\beta_0. \tag{4.74}\]

Assumption 8 (ML8: Regular identification)

The limit \[ \frac{1}{n} \begin{pmatrix} \mathbf X_n & \mathbf d_n \end{pmatrix}^{\top} \begin{pmatrix} \mathbf X_n & \mathbf d_n \end{pmatrix} \longrightarrow \mathbf Q_z \tag{4.75}\] exists and \(\mathbf Q_z\) is positive definite.

Assumption 8 is the defining feature of the regular case. It requires the spatially generated regressor \(\mathbf d_n\) to contain asymptotically nonredundant variation after controlling for \(\mathbf X_n\).

4.6.2 Two Preliminary Matrix Results

Define the projection and residual-maker matrices \[ \mathbf P_n = \mathbf X_n \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top}, \qquad \mathbf M_n = \mathbf I_n-\mathbf P_n. \tag{4.76}\]

The next lemma verifies a property used repeatedly in the consistency proof.

Lemma 4.1 (Uniform boundedness of the regression projectors) Under Assumption 6, the sequences \(\{\mathbf P_n\}\) and \(\{\mathbf M_n\}\) are uniformly bounded in both absolute row and column sums. Moreover, every element of \(\mathbf P_n\) is \(O(1/n)\) uniformly in its row and column indices.

Proof. This result corresponds to Lemma A.5 in Lee (2004). We prove both claims separately. First, we show that every element of \(\mathbf P_n\) is uniformly of order \(O(1/n)\). We then use this elementwise bound to establish uniform bounds for the absolute row and column sums of \(\mathbf P_n\) and \(\mathbf M_n\).

By Assumption 6, \(\mathbf Q_{xx,n} = \frac{1}{n} \mathbf X_n^{\top}\mathbf X_n \longrightarrow\mathbf Q_{xx}\), where \(\mathbf Q_{xx}\) is finite and positive definite. A positive-definite matrix is nonsingular. To see this, suppose instead that \(\mathbf Q_{xx}\) were singular. There would then exist a nonzero vector \(\mathbf a\) such that \(\mathbf Q_{xx}\mathbf a =\mathbf 0\). Premultiplying by \(\mathbf a^{\top}\) would give \(\mathbf a^{\top} \mathbf Q_{xx} \mathbf a =0\), which contradicts positive definiteness because positive definiteness requires \(\mathbf a^{\top}\mathbf Q_{xx}\mathbf a>0\) for every nonzero vector \(\mathbf a\). Therefore, \(\mathbf Q_{xx}\) is nonsingular.

The convergence of a fixed-dimensional matrix sequence is understood element by element, as stated in Definition 3.2. Moreover, matrix inversion is continuous at every nonsingular matrix, as discussed after Definition 3.3. It follows that \(\mathbf Q_{xx,n}\) is nonsingular for all sufficiently large \(n\) and \[ \mathbf Q_{xx,n}^{-1} \longrightarrow \mathbf Q_{xx}^{-1}. \]

Define \[ \mathbf B_n = \mathbf Q_{xx,n}^{-1}. \]

Because \(\mathbf B_n\) is a \(k\times k\) matrix and \(k\) is fixed, the preceding matrix convergence means that \([\mathbf B_n]_{rs}\longrightarrow[\mathbf Q_{xx}^{-1}]_{rs}\) for every \(r,s=1,\ldots,k\). Each scalar sequence \(\{[\mathbf B_n]_{rs}\}_{n\geq1}\) therefore converges to a finite limit. By Definition 3.4, every such sequence is bounded. Hence, for each pair \((r,s)\), there exists a finite constant \(C_{B,rs}\) such that \(\left|[\mathbf B_n]_{rs}\right|\leq C_{B,rs}\) for all sufficiently large \(n\).

There are only \(k^2\) pairs \((r,s)\) because \(k\) is fixed. We may therefore define \[ C_B = \max_{\substack{1\leq r\leq k\\1\leq s\leq k}} C_{B,rs}. \]

The maximum is finite because it is taken over a finite collection of finite constants. Consequently, \(\left|[\mathbf B_n]_{rs}\right|\leq C_B\) uniformly over \(r,s=1,\ldots,k\) and all sufficiently large \(n\).

The elements of \(\mathbf X_n\) are uniformly bounded by Assumption 6. Thus, there exists a finite constant \(C_X\), independent of \(n\), \(i\), and \(r\), such that \(|x_{n,ir}|\leq C_X\) for every \(i=1,\ldots,n\), every \(r=1,\ldots,k\), and every \(n\). Because \[ \mathbf X_n^{\top}\mathbf X_n = n\mathbf Q_{xx,n}, \] we have \[ \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} = \left( n\mathbf Q_{xx,n} \right)^{-1} = \frac{1}{n} \mathbf Q_{xx,n}^{-1} = \frac{1}{n} \mathbf B_n. \]

The second equality follows because \[ \left( n\mathbf Q_{xx,n} \right) \left( \frac{1}{n}\mathbf Q_{xx,n}^{-1} \right) = \mathbf Q_{xx,n} \mathbf Q_{xx,n}^{-1} = \mathbf I_k. \]

Using the definition of \(\mathbf P_n\) in Equation 4.76, we obtain \(\mathbf P_n=\frac{1}{n}\mathbf X_n\mathbf B_n\mathbf X_n^{\top}\). We now derive the \((i,j)\)th element of this matrix product. First, \([\mathbf X_n\mathbf B_n]_{is}=\sum_{r=1}^{k}x_{n,ir}[\mathbf B_n]_{rs}\). Therefore, \[ \begin{aligned} p_{n,ij} &= \frac{1}{n} \left[ \mathbf X_n \mathbf B_n \mathbf X_n^{\top} \right]_{ij} \\ &= \frac{1}{n} \sum_{s=1}^{k} [\mathbf X_n\mathbf B_n]_{is} [\mathbf X_n^{\top}]_{sj} \\ &= \frac{1}{n} \sum_{s=1}^{k} \left( \sum_{r=1}^{k} x_{n,ir} [\mathbf B_n]_{rs} \right) x_{n,js} \\ &= \frac{1}{n} \sum_{r=1}^{k} \sum_{s=1}^{k} x_{n,ir} [\mathbf B_n]_{rs} x_{n,js}. \end{aligned} \]

The third equality uses \([\mathbf X_n^{\top}]_{sj}=x_{n,js}\). We next bound the absolute value of \(p_{n,ij}\): \[ \begin{aligned} |p_{n,ij}| &= \frac{1}{n} \left| \sum_{r=1}^{k} \sum_{s=1}^{k} x_{n,ir} [\mathbf B_n]_{rs} x_{n,js} \right| \\[4pt] &\leq \frac{1}{n} \sum_{r=1}^{k} \sum_{s=1}^{k} \left| x_{n,ir} [\mathbf B_n]_{rs} x_{n,js} \right| && \text{by the triangle inequality} \\[4pt] &= \frac{1}{n} \sum_{r=1}^{k} \sum_{s=1}^{k} |x_{n,ir}| \left| [\mathbf B_n]_{rs} \right| |x_{n,js}| && \text{because absolute values are multiplicative} \\[4pt] &\leq \frac{1}{n} \sum_{r=1}^{k} \sum_{s=1}^{k} C_X C_B C_X && \text{by the uniform bounds on }\mathbf X_n\text{ and }\mathbf B_n \\[4pt] &= \frac{1}{n} \sum_{r=1}^{k} \sum_{s=1}^{k} C_X^2C_B && \text{because }C_XC_BC_X=C_X^2C_B \\[4pt] &= \frac{k^2C_X^2C_B}{n} && \text{because the double sum contains }k^2\text{ terms}. \end{aligned} \]

Define \(C_P =k^2C_X^2C_B.\) Because \(k\), \(C_X\), and \(C_B\) are finite constants that do not depend on \(n\), \(i\), or \(j\), \(C_P\) is also finite and independent of \(n\), \(i\), and \(j\). We have therefore shown that \[ |p_{n,ij}| \leq \frac{C_P}{n} \] uniformly over \(i\) and \(j\). This is precisely the statement that \[ p_{n,ij} = O\left( \frac{1}{n} \right) \] uniformly in the row and column indices.

We now use this elementwise bound to control the absolute row sums of \(\mathbf P_n\). For any row \(i\), \[ \begin{aligned} \sum_{j=1}^{n} |p_{n,ij}| &\leq \sum_{j=1}^{n} \frac{C_P}{n} && \text{because }|p_{n,ij}|\leq\frac{C_P}{n} \\[4pt] &= n\frac{C_P}{n} && \text{because the sum contains }n\text{ terms} \\[4pt] &= C_P. \end{aligned} \]

The bound does not depend on the selected row \(i\). Taking the maximum over all rows therefore gives \[ \max_{1\leq i\leq n} \sum_{j=1}^{n} |p_{n,ij}| \leq C_P. \]

By Definition 3.20, the sequence \(\{\mathbf P_n\}\) is uniformly bounded in absolute row sums.

The column sums are handled separately. For any column \(j\), \[ \begin{aligned} \sum_{i=1}^{n} |p_{n,ij}| &\leq \sum_{i=1}^{n} \frac{C_P}{n} && \text{because }|p_{n,ij}|\leq\frac{C_P}{n} \\[4pt] &= n\frac{C_P}{n} && \text{because the sum contains }n\text{ terms} \\[4pt] &= C_P. \end{aligned} \]

The bound does not depend on the selected column \(j\). Hence, \[ \max_{1\leq j\leq n} \sum_{i=1}^{n} |p_{n,ij}| \leq C_P. \]

By Definition 3.20, the sequence \(\{\mathbf P_n\}\) is uniformly bounded in absolute column sums. It is therefore uniformly bounded in both absolute row and column sums.

It remains to establish the corresponding result for \(\mathbf M_n\). From Equation 4.76, \(\mathbf M_n = \mathbf I_n-\mathbf P_n.\) Let \(\delta_{ij}\) denote the \((i,j)\)th element of \(\mathbf I_n\). Thus, \[ \delta_{ij} = \begin{cases} 1, & i=j,\\ 0, & i\neq j. \end{cases} \]

The \((i,j)\)th element of \(\mathbf M_n\) is therefore \(m_{n,ij} = \delta_{ij}-p_{n,ij}\). For any row \(i\), the triangle inequality gives \[ \begin{aligned} \sum_{j=1}^{n} |m_{n,ij}| &= \sum_{j=1}^{n} |\delta_{ij}-p_{n,ij}| \\[4pt] &\leq \sum_{j=1}^{n} \left( |\delta_{ij}|+|p_{n,ij}| \right) && \text{by }|a-b|\leq|a|+|b| \\[4pt] &= \sum_{j=1}^{n} |\delta_{ij}| + \sum_{j=1}^{n} |p_{n,ij}| && \text{by separating the two finite sums} \\[4pt] &\leq 1+C_P. \end{aligned} \]

The last inequality follows because row \(i\) of \(\mathbf I_n\) contains exactly one element equal to one and all remaining elements equal to zero. Hence, \(\sum_{j=1}^{n} |\delta_{ij}| = 1\).

We already proved that \(\sum_{j=1}^{n}|p_{n,ij}|\leq C_P\).Therefore, \[ \max_{1\leq i\leq n} \sum_{j=1}^{n} |m_{n,ij}| \leq 1+C_P. \]

Thus, \(\{\mathbf M_n\}\) is uniformly bounded in absolute row sums.

For any column \(j\), the same reasoning is written explicitly as \[ \begin{aligned} \sum_{i=1}^{n} |m_{n,ij}| &= \sum_{i=1}^{n} |\delta_{ij}-p_{n,ij}| \\[4pt] &\leq \sum_{i=1}^{n} \left( |\delta_{ij}|+|p_{n,ij}| \right) && \text{by }|a-b|\leq|a|+|b| \\[4pt] &= \sum_{i=1}^{n} |\delta_{ij}| + \sum_{i=1}^{n} |p_{n,ij}| \\[4pt] &\leq 1+C_P. \end{aligned} \]

Column \(j\) of \(\mathbf I_n\) also contains exactly one element equal to one, so \(\sum_{i=1}^{n}|\delta_{ij}|=1\), while the previously established column bound gives \(\sum_{i=1}^{n}|p_{n,ij}|\leq C_P\). Consequently, \[ \max_{1\leq j\leq n} \sum_{i=1}^{n} |m_{n,ij}| \leq 1+C_P. \]

By Definition 3.20, \(\{\mathbf M_n\}\) is uniformly bounded in both absolute row and column sums.

The next result gives a more interpretable characterization of the positive- definiteness requirement in Assumption 8. Define the augmented regressor matrix \[ \mathbf Z_n = \begin{pmatrix} \mathbf X_n & \mathbf d_n \end{pmatrix}. \]

The matrix \(n^{-1}\mathbf Z_n^{\top}\mathbf Z_n\) is the normalized Gram matrix of the columns of \(\mathbf Z_n\). Its entries are the average cross-products between those columns. Moreover, for every \(\mathbf a\in\mathbb R^{k+1}\), \[ \mathbf a^{\top} \left( \frac{1}{n} \mathbf Z_n^{\top}\mathbf Z_n \right) \mathbf a = \frac{1}{n} \left\| \mathbf Z_n\mathbf a \right\|^2. \]

Thus, positive definiteness of its limit rules out asymptotic linear dependence among the columns of \(\mathbf X_n\) and \(\mathbf d_n\).

Proposition 4.3 (Equivalent residualized identification condition) Suppose Assumption 6 holds and the limit in Equation 4.75 exists. Then the limiting augmented Gram matrix \(\mathbf Q_z\) is positive definite if and only if \[ \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n \longrightarrow q_d >0. \tag{4.77}\]

Proof. Write the normalized augmented Gram matrix as \[ \frac{1}{n} \mathbf Z_n^{\top} \mathbf Z_n = \begin{pmatrix} \mathbf Q_{xx,n} & \mathbf q_{xd,n} \\[4pt] \mathbf q_{xd,n}^{\top} & q_{dd,n} \end{pmatrix}, \] where \[ \mathbf Q_{xx,n} = \frac{1}{n} \mathbf X_n^{\top}\mathbf X_n, \qquad \mathbf q_{xd,n} = \frac{1}{n} \mathbf X_n^{\top}\mathbf d_n, \] and \[ q_{dd,n} = \frac{1}{n} \mathbf d_n^{\top}\mathbf d_n. \]

By Assumption 6, \(\mathbf Q_{xx,n}\longrightarrow \mathbf Q_{xx}\), where \(\mathbf Q_{xx}\) is positive definite. It follows that \(\mathbf Q_{xx,n}\) is positive definite and nonsingular for all sufficiently large \(n\). For such \(n\), define \[ s_n = q_{dd,n} - \mathbf q_{xd,n}^{\top} \mathbf Q_{xx,n}^{-1} \mathbf q_{xd,n}. \]

This scalar is called the Schur complement of \(\mathbf Q_{xx,n}\) in the partitioned Gram matrix. We now show that it is the average variation in \(\mathbf d_n\) that remains after projecting \(\mathbf d_n\) onto the columns of \(\mathbf X_n\).

Since \(\mathbf Q_{xx,n}=\frac{1}{n}\mathbf X_n^{\top}\mathbf X_n\), we have \(\mathbf Q_{xx,n}^{-1} =n\left(\mathbf X_n^{\top}\mathbf X_n\right)^{-1}\). Therefore, \[ \begin{aligned} \mathbf q_{xd,n}^{\top} \mathbf Q_{xx,n}^{-1} \mathbf q_{xd,n} &= \left( \frac{1}{n} \mathbf d_n^{\top}\mathbf X_n \right) \left[ n \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \right] \left( \frac{1}{n} \mathbf X_n^{\top}\mathbf d_n \right) \\ &= \frac{1}{n} \mathbf d_n^{\top} \mathbf X_n \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top} \mathbf d_n \\ &= \frac{1}{n} \mathbf d_n^{\top} \mathbf P_n \mathbf d_n. \end{aligned} \]

Consequently, \[ \begin{aligned} s_n &= \frac{1}{n} \mathbf d_n^{\top}\mathbf d_n - \frac{1}{n} \mathbf d_n^{\top} \mathbf P_n \mathbf d_n \\ &= \frac{1}{n} \mathbf d_n^{\top} \left( \mathbf I_n-\mathbf P_n \right) \mathbf d_n \\ &= \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n. \end{aligned} \]

Because \(\mathbf M_n\) is symmetric and idempotent, \[ \mathbf d_n^{\top} \mathbf M_n \mathbf d_n = \mathbf d_n^{\top} \mathbf M_n^{\top} \mathbf M_n \mathbf d_n = \left\| \mathbf M_n\mathbf d_n \right\|^2. \]

Hence \(s_n\) is the average squared residual from projecting \(\mathbf d_n\) onto \(\mathbf X_n\).

By assumption, the augmented Gram matrix converges. Therefore, its blocks satisfy \[ \mathbf q_{xd,n} \longrightarrow \mathbf q_{xd} \qquad\text{and}\qquad q_{dd,n} \longrightarrow q_{dd} \] for some finite \(\mathbf q_{xd}\) and \(q_{dd}\). Since inversion is continuous at the positive-definite matrix \(\mathbf Q_{xx}\), \[ \mathbf Q_{xx,n}^{-1} \longrightarrow \mathbf Q_{xx}^{-1}. \]

It follows that \[ s_n \longrightarrow s = q_{dd} - \mathbf q_{xd}^{\top} \mathbf Q_{xx}^{-1} \mathbf q_{xd}. \]

The limiting augmented Gram matrix is \[ \mathbf Q_z = \begin{pmatrix} \mathbf Q_{xx} & \mathbf q_{xd} \\[4pt] \mathbf q_{xd}^{\top} & q_{dd} \end{pmatrix}. \] Because \(\mathbf Q_{xx}\) is positive definite, the block-matrix positive-definiteness criterion implies that \(\mathbf Q_z\) is positive definite if and only if its Schur complement is strictly positive: \[ s = q_{dd} - \mathbf q_{xd}^{\top} \mathbf Q_{xx}^{-1} \mathbf q_{xd} >0. \]

Using the identity established above, this condition is equivalent to \[ \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n \longrightarrow q_d >0, \] where \(q_d=s\). This proves the result.

Econometric interpretation

The vector \(\mathbf M_n\mathbf d_n\) contains the residuals obtained by projecting the spatially generated regressor \(\mathbf d_n\) onto the columns of \(\mathbf X_n\). Therefore, \[ \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n = \frac{1}{n} \left\| \mathbf M_n\mathbf d_n \right\|^2 \] measures the average variation in \(\mathbf d_n\) that cannot be explained by \(\mathbf X_n\). The condition \(q_d>0\) requires this residual variation to remain asymptotically nonzero. Economically, changes in the spatial parameter must generate variation in the conditional mean that cannot be reproduced by changing only the coefficients on the original regressors.

The quantity in Equation 4.77 measures the variation in \(\mathbf G_{n0}\mathbf X_n\boldsymbol\beta_0\) that cannot be explained by the original regressors. It is this residual variation that supplies regular mean-based information about \(\rho_0\).

4.6.3 The Population Concentrated Criterion

We now construct the deterministic population criterion used to identify the true spatial parameter. All expectations in this subsection are taken under the true data-generating process.

For each candidate value \(\rho\), define \[ Q_n(\rho) = \max_{\boldsymbol\beta,\;\sigma^2>0} \mathbb E \left[ \ell_n \left( \boldsymbol\beta,\rho,\sigma^2 \right) \right], \tag{4.78}\] where \(\ell_n(\boldsymbol\beta,\rho,\sigma^2)\) is the Gaussian log-likelihood defined in Equation 4.36.

The order of the operations in Equation 4.78 is important. We first take the expectation of the log-likelihood under the true data-generating process and then maximize that expected criterion with respect to \(\boldsymbol \beta\) and \(\sigma^2\). Thus, \(Q_n(\rho) = \max_{\boldsymbol\beta,\;\sigma^2>0}\mathbb E\left[\ell_n\left(\boldsymbol\beta,\rho,\sigma^2\right)\right]\), which should not be confused with \(\mathbb E \left[ \max_{\boldsymbol\beta,\;\sigma^2>0} \ell_n \left( \boldsymbol\beta,\rho,\sigma^2 \right) \right]\). In general, maximization and expectation cannot be interchanged.

For a candidate \(\rho\), define \(\delta(\rho) = \rho_0-\rho\) and \[ \mathbf H_n(\rho) = \mathbf A_n(\rho) \mathbf A_{n0}^{-1}. \tag{4.79}\]

Because \(\mathbf A_n(\rho) = \mathbf I_n-\rho\mathbf W_n\) and \(\mathbf A_{n0} = \mathbf I_n-\rho_0\mathbf W_n\), we have \[ \begin{aligned} \mathbf A_n(\rho) &= \mathbf I_n-\rho\mathbf W_n \\ &= \mathbf I_n-\rho_0\mathbf W_n + (\rho_0-\rho)\mathbf W_n \\ &= \mathbf A_{n0} + \delta(\rho)\mathbf W_n. \end{aligned} \]

Premultiplying by \(\mathbf A_{n0}^{-1}\) on the right gives \[ \begin{aligned} \mathbf H_n(\rho) &= \mathbf A_n(\rho) \mathbf A_{n0}^{-1} \\ &= \left[ \mathbf A_{n0} + \delta(\rho)\mathbf W_n \right] \mathbf A_{n0}^{-1} \\ &= \mathbf A_{n0}\mathbf A_{n0}^{-1} + \delta(\rho) \mathbf W_n\mathbf A_{n0}^{-1} \\ &= \mathbf I_n + \delta(\rho)\mathbf G_{n0}, \end{aligned} \tag{4.80}\] where \(\mathbf G_{n0} = \mathbf W_n\mathbf A_{n0}^{-1}\). Using the reduced form, \(\mathbf y_n = \mathbf A_{n0}^{-1} \left(\mathbf X_n\boldsymbol\beta_0 +\boldsymbol\varepsilon_n\right)\), we obtain \[ \begin{aligned} \mathbf A_n(\rho)\mathbf y_n &= \mathbf A_n(\rho) \mathbf A_{n0}^{-1} \left( \mathbf X_n\boldsymbol\beta_0 + \boldsymbol\varepsilon_n \right) \\ &= \mathbf H_n(\rho) \left( \mathbf X_n\boldsymbol\beta_0 + \boldsymbol\varepsilon_n \right) \\ &= \mathbf H_n(\rho) \mathbf X_n\boldsymbol\beta_0 + \mathbf H_n(\rho) \boldsymbol\varepsilon_n. \end{aligned} \]

Using Equation 4.80 in the first term gives \[ \begin{aligned} \mathbf H_n(\rho) \mathbf X_n\boldsymbol\beta_0 &= \left[ \mathbf I_n + \delta(\rho)\mathbf G_{n0} \right] \mathbf X_n\boldsymbol\beta_0 \\ &= \mathbf X_n\boldsymbol\beta_0 + \delta(\rho) \mathbf G_{n0} \mathbf X_n\boldsymbol\beta_0 \\ &= \mathbf X_n\boldsymbol\beta_0 + \delta(\rho)\mathbf d_n, \end{aligned} \] where \(\mathbf d_n = \mathbf G_{n0}\mathbf X_n\boldsymbol\beta_0\). Consequently, \[ \mathbf A_n(\rho)\mathbf y_n = \mathbf X_n\boldsymbol\beta_0 + \delta(\rho)\mathbf d_n + \mathbf H_n(\rho) \boldsymbol\varepsilon_n. \tag{4.81}\]

Taking expectations and using \(\mathbb E\left(\boldsymbol\varepsilon_n\right)=\mathbf 0\) yields \[ \mathbb E \left[ \mathbf A_n(\rho)\mathbf y_n \right] = \mathbf X_n\boldsymbol\beta_0 + \delta(\rho)\mathbf d_n. \tag{4.82}\]

4.6.3.1 Profiling the regression coefficients

For fixed \(\rho\) and \(\sigma^2\), the only part of the expected Gaussian log-likelihood that depends on \(\boldsymbol \beta\) is \[ - \frac{1}{2\sigma^2} \mathbb E \left[ \left\| \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta \right\|^2 \right]. \]

Therefore, maximizing the expected log-likelihood with respect to \(\boldsymbol \beta\) is equivalent to minimizing \[ R_n \left( \boldsymbol\beta;\rho \right) = \mathbb E \left[ \left\| \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta \right\|^2 \right]. \]

Expanding the squared norm gives \[ \begin{aligned} R_n \left( \boldsymbol\beta;\rho \right) &= \mathbb E \left[ \left( \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta \right)^{\top} \left( \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta \right) \right] \\ &= \mathbb E \left[ \mathbf y_n^{\top} \mathbf A_n(\rho)^{\top} \mathbf A_n(\rho) \mathbf y_n \right] - 2\boldsymbol\beta^{\top} \mathbf X_n^{\top} \mathbb E \left[ \mathbf A_n(\rho)\mathbf y_n \right] + \boldsymbol\beta^{\top} \mathbf X_n^{\top} \mathbf X_n \boldsymbol\beta. \end{aligned} \]

The first term does not depend on \(\boldsymbol \beta\). Differentiating the remaining terms with respect to \(\boldsymbol \beta\) gives \[ \frac{\partial R_n \left( \boldsymbol\beta;\rho \right)} {\partial\boldsymbol\beta} = -2 \mathbf X_n^{\top} \mathbb E \left[ \mathbf A_n(\rho)\mathbf y_n \right] + 2 \mathbf X_n^{\top}\mathbf X_n \boldsymbol\beta. \]

The first-order condition is therefore \[ \mathbf X_n^{\top}\mathbf X_n \boldsymbol\beta = \mathbf X_n^{\top} \mathbb E \left[ \mathbf A_n(\rho)\mathbf y_n \right]. \]

Under Assumption 6, \(\mathbf X_n^{\top}\mathbf X_n\) is nonsingular for all sufficiently large \(n\). Hence, the unique minimizer is \[ \boldsymbol\beta_n^{*}(\rho) = \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top} \mathbb E \left[ \mathbf A_n(\rho)\mathbf y_n \right]. \]

Substituting Equation 4.82 gives \[ \begin{aligned} \boldsymbol\beta_n^{*}(\rho) &= \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top} \left[ \mathbf X_n\boldsymbol\beta_0 + \delta(\rho)\mathbf d_n \right] \\ &= \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top}\mathbf X_n \boldsymbol\beta_0 + \delta(\rho) \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top}\mathbf d_n \\ &= \boldsymbol\beta_0 + \delta(\rho) \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top}\mathbf d_n. \end{aligned} \tag{4.83}\]

The Hessian of \(R_n(\boldsymbol\beta;\rho)\) with respect to \(\boldsymbol \beta\) is \[ \frac{\partial^2 R_n \left( \boldsymbol\beta;\rho \right)} {\partial\boldsymbol\beta \partial\boldsymbol\beta^{\top}} = 2\mathbf X_n^{\top}\mathbf X_n. \]

This matrix is positive definite for all sufficiently large \(n\) under Assumption 6. Thus, \(\boldsymbol \beta_n^{*}(\rho)\) is the unique global minimizer of the expected sum of squared residuals and, equivalently, the unique population maximizer of the expected log-likelihood with respect to \(\boldsymbol \beta\).

4.6.3.2 The population residual

Multiplying Equation 4.83 by \(\mathbf X_n\) gives \[ \begin{aligned} \mathbf X_n \boldsymbol\beta_n^{*}(\rho) &= \mathbf X_n\boldsymbol\beta_0 + \delta(\rho) \mathbf X_n \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top} \mathbf d_n \\ &= \mathbf X_n\boldsymbol\beta_0 + \delta(\rho) \mathbf P_n\mathbf d_n, \end{aligned} \] where \(\mathbf P_n\) is defined in Equation 4.76. Subtracting this expression from Equation 4.81 yields \[ \begin{aligned} \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta_n^{*}(\rho) &= \mathbf X_n\boldsymbol\beta_0 + \delta(\rho)\mathbf d_n + \mathbf H_n(\rho)\boldsymbol\varepsilon_n - \mathbf X_n\boldsymbol\beta_0 - \delta(\rho)\mathbf P_n\mathbf d_n \\ &= \delta(\rho) \left( \mathbf I_n-\mathbf P_n \right) \mathbf d_n + \mathbf H_n(\rho)\boldsymbol\varepsilon_n \\ &= \delta(\rho) \mathbf M_n\mathbf d_n + \mathbf H_n(\rho)\boldsymbol\varepsilon_n. \end{aligned} \tag{4.84}\]

4.6.3.3 Profiling the variance

Define the population profile variance associated with the candidate \(\rho\) as \[ \sigma_n^{*2}(\rho) = \frac{1}{n} \mathbb E \left[ \left\| \mathbf A_n(\rho)\mathbf y_n - \mathbf X_n\boldsymbol\beta_n^{*}(\rho) \right\|^2 \right]. \]

The term “profile variance” is used because \(\sigma_n^{*2}(\rho)\) is the population value of \(\sigma^2\) that maximizes the expected Gaussian criterion for a fixed candidate \(\rho\). Except at \(\rho=\rho_0\), it is not the variance of the primitive innovations alone.

Using Equation 4.84, \[ \begin{aligned} \sigma_n^{*2}(\rho) &= \frac{1}{n} \mathbb E \left[ \left\| \delta(\rho)\mathbf M_n\mathbf d_n + \mathbf H_n(\rho)\boldsymbol\varepsilon_n \right\|^2 \right]. \end{aligned} \]

Expanding the squared norm gives \[ \begin{aligned} \sigma_n^{*2}(\rho) = \frac{\delta(\rho)^2}{n} \mathbf d_n^{\top} \mathbf M_n^{\top} \mathbf M_n \mathbf d_n + \frac{2\delta(\rho)}{n} \mathbf d_n^{\top} \mathbf M_n^{\top} \mathbf H_n(\rho) \mathbb E \left( \boldsymbol\varepsilon_n \right) + \frac{1}{n} \mathbb E \left[ \boldsymbol\varepsilon_n^{\top} \mathbf H_n(\rho)^{\top} \mathbf H_n(\rho) \boldsymbol\varepsilon_n \right]. \end{aligned} \]

The middle term is zero because \(\mathbb E\left(\boldsymbol\varepsilon_n\right)=\mathbf 0\). Moreover, \(\mathbf M_n\) is symmetric and idempotent, as used in Proposition 4.3. Therefore, \(\mathbf M_n^{\top}\mathbf M_n = \mathbf M_n^2 =\mathbf M_n,\) and hence \(\mathbf d_n^{\top} \mathbf M_n^{\top}\mathbf M_n \mathbf d_n = \mathbf d_n^{\top}\mathbf M_n\mathbf d_n\).

For the final expectation, apply Equation 3.17 from Lemma 3.10 in Chapter 3 with \(\mathbf A_n = \mathbf H_n(\rho)^{\top} \mathbf H_n(\rho)\), \(\mathbb E\left(\boldsymbol\varepsilon_n\right)=\mathbf 0\), and \(\operatorname{Var}\left(\boldsymbol\varepsilon_n\right)=\sigma_0^2\mathbf I_n\). It follows that \[ \begin{aligned} \mathbb E \left[ \boldsymbol\varepsilon_n^{\top} \mathbf H_n(\rho)^{\top} \mathbf H_n(\rho) \boldsymbol\varepsilon_n \right] &= \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf H_n(\rho) \sigma_0^2\mathbf I_n \right] \\ &= \sigma_0^2 \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf H_n(\rho) \right]. \end{aligned} \]

Combining the three terms gives \[ \begin{aligned} \sigma_n^{*2}(\rho) = \delta(\rho)^2 \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n + \frac{\sigma_0^2}{n} \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf H_n(\rho) \right]. \end{aligned} \tag{4.85}\]

The first term measures the average squared error in the conditional mean created by using \(\rho\) instead of \(\rho_0\). The second term measures the average variance of the transformed innovations.

Define \[ \sigma_{p,n}^2(\rho) = \frac{\sigma_0^2}{n} \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf H_n(\rho) \right]. \tag{4.86}\]

Then \[ \sigma_n^{*2}(\rho) = \sigma_{p,n}^2(\rho) + \delta(\rho)^2 \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n. \tag{4.87}\]

We now verify that \(\sigma_n^{*2}(\rho)\) is the population maximizer with respect to \(\sigma^2\). After substituting \(\boldsymbol \beta_n^{*}(\rho)\) into the expected Gaussian log-likelihood, the terms that depend on \(\sigma^2\) are \[ -\frac{n}{2}\log\sigma^2 - \frac{n}{2\sigma^2} \sigma_n^{*2}(\rho). \]

Let \(v=\sigma^2>0\). Define \[ f_n(v) = - \frac{n}{2}\log v - \frac{n}{2v} \sigma_n^{*2}(\rho). \]

Its derivative is \[ \begin{aligned} f_n'(v) &= - \frac{n}{2v} + \frac{n}{2v^2} \sigma_n^{*2}(\rho) \\ &= \frac{n}{ 2v^2 } \left[ \sigma_n^{*2}(\rho)-v \right]. \end{aligned} \]

Because \(n/(2v^2)>0\), the sign of \(f_n'(v)\) is determined by \(\sigma_n^{*2}(\rho)-v\). Therefore, \[ f_n'(v)>0 \qquad \text{when } 0<v<\sigma_n^{*2}(\rho), \] and \[ f_n'(v)<0 \qquad \text{when } v>\sigma_n^{*2}(\rho). \]

Thus, \(f_n(v)\) increases until \(v=\sigma_n^{*2}(\rho)\) and decreases afterward. Its unique global maximizer is consequently \[ v = \sigma_n^{*2}(\rho). \]

Since \(v=\sigma^2\), the population profile value of the variance parameter is \[ \sigma^2 = \sigma_n^{*2}(\rho). \]

Substituting \(\boldsymbol\beta_n^{*}(\rho)\) and \(\sigma_n^{*2}(\rho)\) into the expected log-likelihood gives \[ \begin{aligned} Q_n(\rho) &= - \frac{n}{2} \log(2\pi) - \frac{n}{2} \log \sigma_n^{*2}(\rho) + \log \left| \det\mathbf A_n(\rho) \right| - \frac{ n\sigma_n^{*2}(\rho) }{ 2\sigma_n^{*2}(\rho) } \\ &= - \frac{n}{2} \left[ \log(2\pi)+1 \right] - \frac{n}{2} \log \sigma_n^{*2}(\rho) + \log \left| \det\mathbf A_n(\rho) \right|. \end{aligned} \tag{4.88}\]

At the true parameter, \(\delta(\rho_0) = \rho_0-\rho_0 =0\). Consequently, \(\mathbf H_n(\rho_0) = \mathbf I_n\), \(\boldsymbol\beta_n^{*}(\rho_0) = \boldsymbol\beta_0\), and \[ \begin{aligned} \sigma_n^{*2}(\rho_0) &= \frac{\sigma_0^2}{n} \operatorname{tr} \left( \mathbf I_n^{\top}\mathbf I_n \right) \\ &= \frac{\sigma_0^2}{n} \operatorname{tr} \left( \mathbf I_n \right) \\ &= \frac{\sigma_0^2}{n}n \\ &= \sigma_0^2. \end{aligned} \]

4.6.4 Why the Population Criterion Is Maximized at \(\rho_0\)

The identification argument can be separated into two parts. First, the covariance structure alone cannot make the population criterion exceed its value at \(\rho_0\). Second, regular mean identification makes the inequality strict away from \(\rho_0\).

Define the auxiliary criterion \[ \begin{aligned} Q_{p,n}(\rho) = - \frac{n}{2} \left[ \log(2\pi)+1 \right] - \frac{n}{2} \log \sigma_{p,n}^{2}(\rho) + \log \left| \det\mathbf A_n(\rho) \right|. \end{aligned} \tag{4.89}\]

Lemma 4.2 (The auxiliary criterion cannot exceed its true value) For every admissible \(\rho\), \[ Q_{p,n}(\rho) \leq Q_{p,n}(\rho_0). \tag{4.90}\]

Proof. Let \(s_{n1}(\rho),\ldots,s_{nn}(\rho)\) be the singular values of \(\mathbf H_n(\rho)\). Since \[ \mathbf A_n(\rho) = \mathbf H_n(\rho) \mathbf A_{n0}, \] we have \[ \log \left| \det\mathbf A_n(\rho) \right| {}- \log \left| \det\mathbf A_{n0} \right| = \log \left| \det\mathbf H_n(\rho) \right|. \]

The determinant equals the product of the singular values in absolute value, so \[ \frac{1}{n} \log \left| \det\mathbf H_n(\rho) \right| = \frac{1}{n} \sum_{i=1}^{n} \log s_{ni}(\rho). \]

Also, \[ \frac{ \sigma_{p,n}^{2}(\rho) }{ \sigma_0^2 } = \frac{1}{n} \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf H_n(\rho) \right] = \frac{1}{n} \sum_{i=1}^{n} s_{ni}(\rho)^2. \]

Therefore, \[ \begin{aligned} \frac{1}{n} \left[ Q_{p,n}(\rho) {}- Q_{p,n}(\rho_0) \right] = - \frac{1}{2} \log \left[ \frac{1}{n} \sum_{i=1}^{n} s_{ni}(\rho)^2 \right] + \frac{1}{2n} \sum_{i=1}^{n} \log s_{ni}(\rho)^2. \end{aligned} \]

Because the logarithm is concave, Jensen’s inequality from Theorem 3.29 implies \[ \frac{1}{n} \sum_{i=1}^{n} \log s_{ni}(\rho)^2 \leq \log \left[ \frac{1}{n} \sum_{i=1}^{n} s_{ni}(\rho)^2 \right]. \]

Substitution into the preceding display yields Equation 4.90.

Now use Equation 4.87: \[ \begin{aligned} \frac{1}{n} \left[ Q_n(\rho)-Q_n(\rho_0) \right] = \frac{1}{n} \left[ Q_{p,n}(\rho)-Q_{p,n}(\rho_0) \right] - \frac{1}{2} \log \left[ \frac{ \sigma_n^{*2}(\rho) }{ \sigma_{p,n}^{2}(\rho) } \right]. \end{aligned} \tag{4.91}\]

Both terms on the right-hand side are nonpositive. The second term equals \[ {}- \frac{1}{2} \log \left[ 1 + \delta(\rho)^2 \frac{ n^{-1}\mathbf d_n^{\top}\mathbf M_n\mathbf d_n }{ \sigma_{p,n}^{2}(\rho) } \right]. \tag{4.92}\]

Under Assumption 5–Assumption 8, the denominator in Equation 4.92 is uniformly bounded above and bounded away from zero. The upper bound follows from the uniform row- and column-sum bounds for \(\mathbf H_n(\rho)\). To see the lower bound, note that \[ \mathbf H_n(\rho)^{-1} = \mathbf A_{n0} \mathbf A_n(\rho)^{-1}. \]

Assumption 5 and Assumption 7 imply that this inverse is uniformly bounded in at least one absolute row- or column-sum norm. Suppose, for example, its row sums are bounded by \(C_H\). Since \[ \mathbf H_n(\rho)^{-1} \mathbf H_n(\rho) = \mathbf I_n, \] for every column \(j\), \[ 1 \leq C_H \max_i |h_{n,ij}(\rho)|. \]

Thus every column of \(\mathbf H_n(\rho)\) contains at least one entry whose absolute value is no smaller than \(1/C_H\). Consequently, \[ \frac{1}{n} \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf H_n(\rho) \right] \geq \frac{1}{C_H^2}. \]

The same argument applies with rows and columns interchanged when the inverse is bounded in column sums.

For any fixed \(\eta>0\), Proposition 4.3 implies that, for sufficiently large \(n\), \[ \inf_{ \rho\in\mathcal R: |\rho-\rho_0|\geq\eta } \delta(\rho)^2 \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n >0. \]

Therefore, there exists \(c_\eta>0\) such that \[ \limsup_{n\to\infty} \sup_{ \rho\in\mathcal R: |\rho-\rho_0|\geq\eta } \frac{1}{n} \left[ Q_n(\rho)-Q_n(\rho_0) \right] \leq -c_\eta. \tag{4.93}\]

Equation Equation 4.93 is the required global identification condition for the concentrated criterion.

4.6.5 Uniform Convergence of the Sample Criterion

The sample concentrated variance from Equation 4.30 is \[ \widehat\sigma_n^2(\rho) = \frac{1}{n} \mathbf y_n^{\top} \mathbf A_n(\rho)^{\top} \mathbf M_n \mathbf A_n(\rho) \mathbf y_n. \]

Using Equation 4.81 and \(\mathbf M_n\mathbf X_n=\mathbf0\), \[ \mathbf M_n \mathbf A_n(\rho) \mathbf y_n = \delta(\rho) \mathbf M_n\mathbf d_n + \mathbf M_n \mathbf H_n(\rho) \boldsymbol\varepsilon_n. \]

Hence \[ \begin{aligned} \widehat\sigma_n^2(\rho) = \delta(\rho)^2 \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n + \frac{2\delta(\rho)}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf H_n(\rho) \boldsymbol\varepsilon_n + \frac{1}{n} \boldsymbol\varepsilon_n^{\top} \mathbf H_n(\rho)^{\top} \mathbf M_n \mathbf H_n(\rho) \boldsymbol\varepsilon_n. \end{aligned} \tag{4.94}\]

Subtracting Equation 4.85 gives \[ \begin{aligned} \widehat\sigma_n^2(\rho) {}- \sigma_n^{*2}(\rho) &= \frac{2\delta(\rho)}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf H_n(\rho) \boldsymbol\varepsilon_n \\ &\quad+ \frac{1}{n} \Big\{ \boldsymbol\varepsilon_n^{\top} \mathbf H_n(\rho)^{\top} \mathbf M_n \mathbf H_n(\rho) \boldsymbol\varepsilon_n \\ &\qquad\qquad- \sigma_0^2 \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf M_n \mathbf H_n(\rho) \right] \Big\} \\ &\quad- \frac{\sigma_0^2}{n} \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf P_n \mathbf H_n(\rho) \right]. \end{aligned} \tag{4.95}\]

Lemma 4.3 (Uniform convergence of the concentrated variance) Under Assumption 1–Assumption 7, \[ \sup_{\rho\in\mathcal R} \left| \widehat\sigma_n^2(\rho) {}- \sigma_n^{*2}(\rho) \right| = o_p(1). \tag{4.96}\]

Proof. We control the three terms in Equation 4.95. First, Assumption 5 and Assumption 6 imply that the elements of \(\mathbf d_n = \mathbf G_{n0} \mathbf X_n \boldsymbol\beta_0\) are uniformly bounded. Therefore, \[ \frac{1}{n} \mathbf d_n^{\top}\mathbf d_n =O(1). \]

The matrices \(\mathbf M_n\) and \(\mathbf H_n(\rho)\) are uniformly bounded in absolute row and column sums. By Equation 3.32, their spectral norms are uniformly bounded. Applying Lemma 3.14 gives \[ \frac{1}{\sqrt n} \mathbf d_n^{\top} \mathbf M_n \mathbf H_n(\rho) \boldsymbol\varepsilon_n =O_p(1) \] for each \(\rho\). Since \[ \mathbf H_n(\rho) = \mathbf I_n + \delta(\rho)\mathbf G_{n0}, \] the expression is affine in \(\delta(\rho)\). The compactness of \(\mathcal R\) therefore gives the uniform result \[ \sup_{\rho\in\mathcal R} \left| \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf H_n(\rho) \boldsymbol\varepsilon_n \right| = o_p(1). \tag{4.97}\]

Second, the matrix \[ \mathbf H_n(\rho)^{\top} \mathbf M_n \mathbf H_n(\rho) \] is uniformly bounded in absolute row and column sums. By Lemma 3.13, \[ \frac{1}{n} \Big\{ \boldsymbol\varepsilon_n^{\top} \mathbf H_n(\rho)^{\top} \mathbf M_n \mathbf H_n(\rho) \boldsymbol\varepsilon_n {}- \sigma_0^2 \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf M_n \mathbf H_n(\rho) \right] \Big\} = o_p(1) \] for each \(\rho\). Expanding \(\mathbf H_n(\rho)=\mathbf I_n+\delta(\rho)\mathbf G_{n0}\) shows that this quadratic form is a polynomial of degree two in \(\delta(\rho)\) with only a finite number of stochastic coefficients. Each coefficient is \(o_p(1)\) after division by \(n\). Since \(\delta(\rho)\) is uniformly bounded on the compact set \(\mathcal R\), the convergence is uniform: \[ \sup_{\rho\in\mathcal R} \left| \frac{1}{n} \Big\{ \boldsymbol\varepsilon_n^{\top} \mathbf H_n(\rho)^{\top} \mathbf M_n \mathbf H_n(\rho) \boldsymbol\varepsilon_n {}- \sigma_0^2 \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf M_n \mathbf H_n(\rho) \right] \Big\} \right| = o_p(1). \tag{4.98}\]

Third, \(\mathbf P_n\) has rank \(k\). Using cyclicity of the trace, \[ \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf P_n \mathbf H_n(\rho) \right] = \operatorname{tr} \left[ \mathbf P_n \mathbf H_n(\rho) \mathbf H_n(\rho)^{\top} \right]. \]

Because \(\mathbf P_n\) is an orthogonal projector, \(\operatorname{tr}(\mathbf P_n)=k\). The spectral norm of \(\mathbf H_n(\rho)\mathbf H_n(\rho)^{\top}\) is uniformly bounded. Hence \[ 0 \leq \operatorname{tr} \left[ \mathbf P_n \mathbf H_n(\rho) \mathbf H_n(\rho)^{\top} \right] \leq k \left\| \mathbf H_n(\rho) \mathbf H_n(\rho)^{\top} \right\|_2 = O(1) \] uniformly over \(\rho\). Therefore, \[ \sup_{\rho\in\mathcal R} \frac{1}{n} \operatorname{tr} \left[ \mathbf H_n(\rho)^{\top} \mathbf P_n \mathbf H_n(\rho) \right] = O\!\left(\frac{1}{n}\right). \tag{4.99}\]

Combining Equation 4.97, Equation 4.98, and Equation 4.99 proves Equation 4.96.

The lower bound established in Section 4.6.4 and Equation 4.96 imply that both \(\sigma_n^{*2}(\rho)\) and \(\widehat\sigma_n^2(\rho)\) are bounded away from zero uniformly over \(\mathcal R\), with probability approaching one. The mean value theorem for the logarithm then gives \[ \sup_{\rho\in\mathcal R} \left| \log\widehat\sigma_n^2(\rho) {}- \log\sigma_n^{*2}(\rho) \right| = o_p(1). \]

Since the determinant terms in the sample and population concentrated criteria are identical, \[ \sup_{\rho\in\mathcal R} \left| \frac{1}{n} \ell_n^c(\rho) {}- \frac{1}{n} Q_n(\rho) \right| = o_p(1). \tag{4.100}\]

4.6.6 Consistency of the Gaussian QMLE

Theorem 4.1 (Consistency of the regular SLM QMLE) Suppose ass-slm-ml1-innovations–Assumption 8 hold. Let \[ \widehat\rho_n \in \arg\max_{\rho\in\mathcal R} \ell_n^c(\rho), \] and define \[ \widehat{\boldsymbol\beta}_n = \widehat{\boldsymbol\beta}_n (\widehat\rho_n), \qquad \widehat\sigma_n^2 = \widehat\sigma_n^2 (\widehat\rho_n). \]

Then \[ \widehat\rho_n \overset{p}{\longrightarrow} \rho_0, \qquad \widehat{\boldsymbol\beta}_n \overset{p}{\longrightarrow} \boldsymbol\beta_0, \qquad \widehat\sigma_n^2 \overset{p}{\longrightarrow} \sigma_0^2. \]

Equivalently, \[ \widehat{\boldsymbol\theta}_n \overset{p}{\longrightarrow} \boldsymbol\theta_0, \] where \[ \boldsymbol\theta_0 = \left( \boldsymbol\beta_0^{\top}, \rho_0, \sigma_0^2 \right)^{\top}. \]

Proof. Fix \(\eta>0\) and define \[ \mathcal R_\eta = \left\{ \rho\in\mathcal R: |\rho-\rho_0|\geq\eta \right\}. \]

By Equation 4.93, there exists \(c_\eta>0\) such that, for sufficiently large \(n\), \[ \sup_{\rho\in\mathcal R_\eta} \frac{1}{n} \left[ Q_n(\rho)-Q_n(\rho_0) \right] \leq -c_\eta. \]

By Equation 4.100, \[ \sup_{\rho\in\mathcal R} \left| \frac{1}{n} \ell_n^c(\rho) {}- \frac{1}{n} Q_n(\rho) \right| = o_p(1). \]

Therefore, with probability approaching one, \[ \sup_{\rho\in\mathcal R} \left| \frac{1}{n} \ell_n^c(\rho) {}- \frac{1}{n} Q_n(\rho) \right| < \frac{c_\eta}{3}. \]

On that event, for every \(\rho\in\mathcal R_\eta\), \[ \begin{aligned} \frac{1}{n} \left[ \ell_n^c(\rho)-\ell_n^c(\rho_0) \right] &\leq \frac{1}{n} \left[ Q_n(\rho)-Q_n(\rho_0) \right] + \frac{2c_\eta}{3} \\ &\leq {}- \frac{c_\eta}{3} <0. \end{aligned} \]

Hence no maximizer of the sample criterion can belong to \(\mathcal R_\eta\) on that event. Thus \[ \Pr \left( |\widehat\rho_n-\rho_0| \geq\eta \right) \longrightarrow0. \]

Since \(\eta>0\) was arbitrary, \[ \widehat\rho_n \overset{p}{\longrightarrow} \rho_0. \]

Next, use the profile coefficient formula Equation 4.26. From Equation 4.81, \[ \begin{aligned} \widehat{\boldsymbol\beta}_n(\rho) {}- \boldsymbol\beta_0 = \delta(\rho) \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top}\mathbf d_n + \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top} \mathbf H_n(\rho) \boldsymbol\varepsilon_n. \end{aligned} \tag{4.101}\]

The first term in Equation 4.101 converges in probability to zero at \(\rho=\widehat\rho_n\) because \(\widehat\rho_n\to_p\rho_0\) and \(\left(\mathbf X_n^{\top}\mathbf X_n\right)^{-1}\mathbf X_n^{\top}\mathbf d_n\) is bounded. For the second term, write \[ \begin{aligned} \left( \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \mathbf X_n^{\top} \mathbf H_n(\rho) \boldsymbol\varepsilon_n = \left( \frac{1}{n} \mathbf X_n^{\top}\mathbf X_n \right)^{-1} \frac{1}{n} \mathbf X_n^{\top} \mathbf H_n(\rho) \boldsymbol\varepsilon_n. \end{aligned} \]

The first factor converges to \(\mathbf Q_{xx}^{-1}\). The second factor is \(o_p(1)\) uniformly over \(\rho\) by the same linear-form argument used in Equation 4.97. Therefore, \[ \widehat{\boldsymbol\beta}_n \overset{p}{\longrightarrow} \boldsymbol\beta_0. \]

Finally, the coefficients in Equation 4.85 are uniformly bounded. Consequently, for a finite constant \(C\) independent of \(n\), \[ \left| \sigma_n^{*2}(\rho_1) {}- \sigma_n^{*2}(\rho_2) \right| \leq C|\rho_1-\rho_2|. \]

This bound holds for all \(\rho_1,\rho_2\in\mathcal R\) and sufficiently large \(n\). Thus, Equation 4.96 and \(\widehat\rho_n\overset{p}{\longrightarrow}\rho_0\) imply \[ \begin{aligned} \widehat\sigma_n^2 &= \widehat\sigma_n^2(\widehat\rho_n) \\ &= \sigma_n^{*2}(\widehat\rho_n) + o_p(1) \\ &\overset{p}{\longrightarrow} \sigma_n^{*2}(\rho_0) = \sigma_0^2. \end{aligned} \]

This proves joint consistency.

The proof shows why the concentrated likelihood is convenient. Consistency of the scalar spatial parameter is established first. Consistency of \(\widehat{\boldsymbol\beta}_n\) and \(\widehat\sigma_n^2\) then follows from their profile formulas.

4.6.7 The Score at the True Parameter

As shown in Equation 4.44, the score for \(\rho\) contains both a linear and a centered quadratic form. Define \[ \mathbf g_n = \operatorname{diag} \left( \mathbf G_{n0} \right), \] where \(\mathbf g_n\) is the \(n\times1\) vector containing the diagonal elements of \(\mathbf G_{n0}\). Also let \(\boldsymbol\iota_n\) denote the \(n\times1\) vector of ones.

At \(\boldsymbol\theta_0\), the normalized score is \[ \frac{1}{\sqrt n} \mathbf s_n(\boldsymbol\theta_0) = \begin{pmatrix} \dfrac{1}{\sigma_0^2\sqrt n} \mathbf X_n^{\top} \boldsymbol\varepsilon_n \\[12pt] \dfrac{1}{\sigma_0^2\sqrt n} \left\{ \mathbf d_n^{\top} \boldsymbol\varepsilon_n + \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n {}- \sigma_0^2 \operatorname{tr} \mathbf G_{n0} \right\} \\[12pt] \dfrac{1}{2\sigma_0^4\sqrt n} \left( \boldsymbol\varepsilon_n^{\top} \boldsymbol\varepsilon_n {}- n\sigma_0^2 \right) \end{pmatrix}. \tag{4.102}\]

The three blocks of the score clarify why the sequence of central limit theorems developed in Section 3.10 is needed here. The \(\boldsymbol\beta\) block is linear in the innovations and, by itself, could be handled by the Kelejian–Prucha linear-form CLT in Theorem 3.22. The \(\rho\) block, however, combines a linear form with a centered quadratic form, while the \(\sigma^2\) block is itself a centered quadratic form. The linear-form theorem is therefore not sufficient for the complete score vector. Its joint limiting distribution must instead be obtained from the Kelejian–Prucha vector linear–quadratic CLT in Theorem 3.26, which is built from the scalar result in Theorem 3.24 and the Cramér–Wold device.

Assume the following limits exist: \[ \boldsymbol\Sigma_{\theta} = \lim_{n\to\infty} \boldsymbol\Sigma_{\theta,n}, \qquad \boldsymbol\Omega_{\theta} = \lim_{n\to\infty} \boldsymbol\Omega_{\theta,n}. \]

The negative expected average Hessian is \[ \boldsymbol\Sigma_{\theta,n} = {}- \mathbb E \left[ \frac{1}{n} \frac{\partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\boldsymbol\theta \partial\boldsymbol\theta^{\top} } \right] \] and equals \[ \boldsymbol\Sigma_{\theta,n} = \begin{pmatrix} \dfrac{1}{n\sigma_0^2} \mathbf X_n^{\top}\mathbf X_n & \dfrac{1}{n\sigma_0^2} \mathbf X_n^{\top}\mathbf d_n & \mathbf0 \\[12pt] \dfrac{1}{n\sigma_0^2} \mathbf d_n^{\top}\mathbf X_n & \dfrac{1}{n\sigma_0^2} \mathbf d_n^{\top}\mathbf d_n + \dfrac{1}{n} \left\{ \operatorname{tr} \left( \mathbf G_{n0}^{\top}\mathbf G_{n0} \right) + \operatorname{tr} \left( \mathbf G_{n0}^2 \right) \right\} & \dfrac{1}{n\sigma_0^2} \operatorname{tr} \mathbf G_{n0} \\[12pt] \mathbf0^{\top} & \dfrac{1}{n\sigma_0^2} \operatorname{tr} \mathbf G_{n0} & \dfrac{1}{2\sigma_0^4} \end{pmatrix}. \tag{4.103}\]

Let \(\mu_3 = \mathbb E\left(\varepsilon_{ni}^3\right)\), and \(\mu_4 = \mathbb E\left(\varepsilon_{ni}^4\right)\). The difference between the score covariance and \(\boldsymbol\Sigma_{\theta,n}\) is \[ \boldsymbol\Omega_{\theta,n} = \begin{pmatrix} \mathbf0 & \dfrac{\mu_3}{n\sigma_0^4} \mathbf X_n^{\top}\mathbf g_n & \dfrac{\mu_3}{2n\sigma_0^6} \mathbf X_n^{\top}\boldsymbol\iota_n \\[12pt] \dfrac{\mu_3}{n\sigma_0^4} \mathbf g_n^{\top}\mathbf X_n & \dfrac{2\mu_3}{n\sigma_0^4} \mathbf d_n^{\top}\mathbf g_n + \dfrac{\mu_4-3\sigma_0^4}{n\sigma_0^4} \mathbf g_n^{\top}\mathbf g_n & \dfrac{ \mu_3\boldsymbol\iota_n^{\top}\mathbf d_n + (\mu_4-3\sigma_0^4) \operatorname{tr}\mathbf G_{n0} }{ 2n\sigma_0^6 } \\[16pt] \dfrac{\mu_3}{2n\sigma_0^6} \boldsymbol\iota_n^{\top}\mathbf X_n & \dfrac{ \mu_3\boldsymbol\iota_n^{\top}\mathbf d_n + (\mu_4-3\sigma_0^4) \operatorname{tr}\mathbf G_{n0} }{ 2n\sigma_0^6 } & \dfrac{\mu_4-3\sigma_0^4}{4\sigma_0^8} \end{pmatrix}. \tag{4.104}\]

The score covariance is therefore \[ \mathbf V_{\theta,n} = \operatorname{Var} \left[ \frac{1}{\sqrt n} \mathbf s_n(\boldsymbol\theta_0) \right] = \boldsymbol\Sigma_{\theta,n} + \boldsymbol\Omega_{\theta,n}. \tag{4.105}\]

The entries in Equation 4.104 follow directly from Lemma 3.12 and Lemma 3.10. For example, \[ \operatorname{Cov} \left( \mathbf X_n^{\top}\boldsymbol\varepsilon_n, \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n \right) = \mu_3 \mathbf X_n^{\top}\mathbf g_n. \]

Similarly, \[ \begin{aligned} \operatorname{Var} \left( \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n \right) = (\mu_4-3\sigma_0^4) \mathbf g_n^{\top}\mathbf g_n + \sigma_0^4 \left\{ \operatorname{tr} \left( \mathbf G_{n0}^{\top}\mathbf G_{n0} \right) + \operatorname{tr} \left( \mathbf G_{n0}^2 \right) \right\}. \end{aligned} \]

If the innovations are Gaussian, then \[ \mu_3=0, \qquad \mu_4=3\sigma_0^4, \] so \[ \boldsymbol\Omega_{\theta,n} = \mathbf0. \tag{4.106}\]

This is the information equality. It generally fails under non-Gaussian innovations.

Lemma 4.4 (Central limit theorem for the regular SLM score) Under Assumption 1 and Assumption 4–Assumption 6, suppose \[ \mathbf V_{\theta,n} \longrightarrow \mathbf V_{\theta}, \] where \(\mathbf V_{\theta}\) is finite and positive definite. Then \[ \frac{1}{\sqrt n} \mathbf s_n(\boldsymbol\theta_0) \overset{d}{\longrightarrow} \mathcal N \left( \mathbf0, \mathbf V_{\theta} \right). \tag{4.107}\]

Proof. We verify directly the primitive conditions of Theorem 3.26. This also makes clear why the linear-form CLT of Theorem 3.22 is not sufficient for the complete likelihood score.

Step 1: Write the score as a centered vector of linear–quadratic forms: Define the symmetric part of \(\mathbf G_{n0}\) by \[ \mathbf G_{n0}^{s} = \frac{ \mathbf G_{n0} + \mathbf G_{n0}^{\top} }{2}. \tag{4.108}\]

By Equation 3.16, \[ \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n = \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0}^{s} \boldsymbol\varepsilon_n. \]

Let \(\mathbf x_{n,r}\) denote column \(r\) of \(\mathbf X_n\). Define the \((k+2)\times1\) vector \[ \widetilde{\mathbf Q}_n = \begin{pmatrix} \dfrac{1}{\sigma_0^2} \mathbf X_n^{\top} \boldsymbol\varepsilon_n \\[12pt] \dfrac{1}{\sigma_0^2} \left( \mathbf d_n^{\top} \boldsymbol\varepsilon_n + \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0}^{s} \boldsymbol\varepsilon_n \right) \\[12pt] \dfrac{1}{2\sigma_0^4} \boldsymbol\varepsilon_n^{\top} \boldsymbol\varepsilon_n \end{pmatrix}. \tag{4.109}\]

Its components have the form required by Theorem 3.26. For \(r=1,\ldots,k\), set \[ \mathbf A_{r,n} = \mathbf0, \qquad \mathbf b_{r,n} = \frac{1}{\sigma_0^2} \mathbf x_{n,r}. \tag{4.110}\]

For the spatial-parameter component, set \[ \mathbf A_{\rho,n} = \frac{1}{\sigma_0^2} \mathbf G_{n0}^{s}, \qquad \mathbf b_{\rho,n} = \frac{1}{\sigma_0^2} \mathbf d_n. \tag{4.111}\]

For the variance-parameter component, set \[ \mathbf A_{\sigma^2,n} = \frac{1}{2\sigma_0^4} \mathbf I_n, \qquad \mathbf b_{\sigma^2,n} = \mathbf0. \tag{4.112}\]

By Equation 3.73 and Assumption 1, \[ \mathbb E \left( \widetilde{\mathbf Q}_n \right) = \begin{pmatrix} \mathbf0_k \\[6pt] \operatorname{tr} \left( \mathbf G_{n0} \right) \\[6pt] \dfrac{n}{2\sigma_0^2} \end{pmatrix}. \tag{4.113}\]

For the spatial component, we used \(\operatorname{tr}\left(\mathbf G_{n0}^{s}\right)=\operatorname{tr}\left(\mathbf G_{n0}\right)\). For the variance component, we used \(\mathbb E\left(\boldsymbol\varepsilon_n^{\top}\boldsymbol\varepsilon_n\right)=n\sigma_0^2\). Comparing Equation 4.109 and Equation 4.113 with Equation 4.102 gives \[ \mathbf s_n(\boldsymbol\theta_0) = \widetilde{\mathbf Q}_n - \mathbb E \left( \widetilde{\mathbf Q}_n \right). \tag{4.114}\]

Thus, the score is exactly the centered vector to which Theorem 3.26 applies.

Step 2: Verify the innovation moment condition: Assumption 1 states that the innovations have a common distribution that does not depend on \(n\) and that, for some \(\gamma>0\), \(\mathbb E\left(|\varepsilon_{ni}|^{4+\gamma}\right)<\infty\). Therefore, \(\sup_{n,i}\mathbb E\left(|\varepsilon_{ni}|^{4+\gamma}\right)<\infty\).

The matrices in Equation 4.111 and Equation 4.112 generally have nonzero diagonal entries. Hence, this is the general-diagonal moment condition in Theorem 3.26.

Step 3: Verify the matrix row- and column-sum conditions: For the \(\boldsymbol\beta\) components, \(\mathbf A_{r,n}=\mathbf0\), so their absolute row and column sums are zero.

For the variance component, every absolute row and column sum of \(\mathbf A_{\sigma^2,n}\) equals \(\frac{1}{2\sigma_0^4}\). For the spatial component, recall that \(\mathbf G_{n0}=\mathbf W_n\mathbf A_{n0}^{-1}\).

By Assumption 5, both factors are uniformly bounded in absolute row and column sums. Therefore, Lemma 3.6 implies that \(\mathbf G_{n0}\) has uniformly bounded absolute row and column sums. The transpose \(\mathbf G_{n0}^{\top}\) has the same property because transposition interchanges row and column sums. The triangle inequality then gives uniform row- and column-sum bounds for \(\mathbf G_{n0}^{s}\) and hence for \(\mathbf A_{\rho,n}\).

Thus, every quadratic-form matrix satisfies the matrix conditions of Theorem 3.26.

Step 4: Verify the condition on the linear coefficients: The elements of \(\mathbf X_n\) are uniformly bounded by Assumption 6. Consequently, the elements of every vector \(\mathbf b_{r,n}\) in Equation 4.110 are uniformly bounded.

We next establish the same property for \(\mathbf d_n\). Recall that \(\mathbf d_n=\mathbf G_{n0}\mathbf X_n\boldsymbol\beta_0\). Because \(k\) is fixed, \(\boldsymbol\beta_0\) is fixed, and the elements of \(\mathbf X_n\) are uniformly bounded, there exists a finite constant \(C_{X\beta}\) such that \[ \left| [\mathbf X_n\boldsymbol\beta_0]_j \right| \leq C_{X\beta} \] uniformly over \(j\) and \(n\). If \(g_{n,ij}\) denotes element \((i,j)\) of \(\mathbf G_{n0}\), then \[ \begin{aligned} |d_{n,i}| &= \left| \sum_{j=1}^{n} g_{n,ij} [\mathbf X_n\boldsymbol\beta_0]_j \right| \\ &\leq C_{X\beta} \sum_{j=1}^{n} |g_{n,ij}|. \end{aligned} \]

The row sums of \(\mathbf G_{n0}\) are uniformly bounded by Step 3. Hence, there exists a finite constant \(C_d\) such that \(|d_{n,i}|\leq C_d\) uniformly over \(i\) and \(n\). Therefore, the elements of \(\mathbf b_{\rho,n}=\mathbf d_n/\sigma_0^2\) are uniformly bounded. The vector \(\mathbf b_{\sigma^2,n}\) is identically zero.

It follows that, for any fixed \(\eta_b>0\) and every score component, \[ \sup_n \frac{1}{n} \sum_{i=1}^{n} |b_{r,ni}|^{2+\eta_b} < \infty. \]

This verifies Equation 3.72.

Step 5: Verify covariance convergence and nondegeneracy: Let \(\boldsymbol\Omega_{Q,n}=\operatorname{Var}\left(\widetilde{\mathbf Q}_n\right)\). Subtracting a deterministic mean does not change a covariance matrix. By Equation 4.114, \[ \boldsymbol\Omega_{Q,n} = \operatorname{Var} \left[ \mathbf s_n(\boldsymbol\theta_0) \right]. \]

Using the definition of \(\mathbf V_{\theta,n}\) in Equation 4.105, \[ \frac{1}{n} \boldsymbol\Omega_{Q,n} = \mathbf V_{\theta,n} \longrightarrow \mathbf V_{\theta}. \tag{4.115}\]

This is the covariance-limit condition Equation 3.77. Because \(\mathbf V_{\theta}\) is positive definite, \(\lambda_{\min}\left(\mathbf V_{\theta}\right)>0\). Eigenvalues are continuous functions of the entries of a fixed-dimensional matrix. Hence, \[ \lambda_{\min} \left( \mathbf V_{\theta,n} \right) \longrightarrow \lambda_{\min} \left( \mathbf V_{\theta} \right). \]

Therefore, for all sufficiently large \(n\), \[ \lambda_{\min} \left( \mathbf V_{\theta,n} \right) \geq \frac{1}{2} \lambda_{\min} \left( \mathbf V_{\theta} \right) > 0. \]

Since \(\boldsymbol\Omega_{Q,n}=n\mathbf V_{\theta,n}\), \[ \frac{1}{n} \lambda_{\min} \left( \boldsymbol\Omega_{Q,n} \right) = \lambda_{\min} \left( \mathbf V_{\theta,n} \right). \]

Thus, the eigenvalue lower bound in Equation 3.75 also holds.

Step 6: Apply the vector linear–quadratic CLT All conditions of Theorem 3.26 have now been verified. Applying Equation 3.78 to \(\widetilde{\mathbf Q}_n\) and using Equation 4.114 gives \[ \frac{1}{\sqrt n} \mathbf s_n(\boldsymbol\theta_0) = \frac{1}{\sqrt n} \left[ \widetilde{\mathbf Q}_n - \mathbb E \left( \widetilde{\mathbf Q}_n \right) \right] \overset{d}{\longrightarrow} \mathcal N \left( \mathbf0, \mathbf V_{\theta} \right). \]

This proves Equation 4.107.

4.6.8 Nonsingularity of the Limiting Information Matrix

Regular identification also prevents the limiting information matrix from becoming singular.

Lemma 4.5 (Positive definiteness of the regular information matrix) Under Assumption 5–Assumption 8, if the limits in Equation 4.103 exist, then \(\boldsymbol\Sigma_{\theta}\) is positive definite.

Proof. Suppose \[ \boldsymbol\Sigma_{\theta} \begin{pmatrix} \mathbf a_{\beta}\\ a_{\rho}\\ a_{\sigma} \end{pmatrix} = \mathbf0. \]

The first block row gives \(\mathbf Q_{xx}\mathbf a_{\beta}+\mathbf q_{xd}a_{\rho}=\mathbf 0\), where \(\mathbf q_{xd} =\lim_{n\to\infty}\frac{1}{n}\mathbf X_n^{\top}\mathbf d_n\). Since \(\mathbf Q_{xx}\) is positive definite, \[ \mathbf a_{\beta} = {}- \mathbf Q_{xx}^{-1} \mathbf q_{xd} a_{\rho}. \tag{4.116}\]

The final block row gives \[ \frac{1}{\sigma_0^2} \left[ \lim_{n\to\infty} \frac{1}{n} \operatorname{tr} \mathbf G_{n0} \right] a_{\rho} + \frac{1}{2\sigma_0^4} a_{\sigma} = 0. \]

Therefore, \[ a_{\sigma} = {}- 2\sigma_0^2 \left[ \lim_{n\to\infty} \frac{1}{n} \operatorname{tr} \mathbf G_{n0} \right] a_{\rho}. \tag{4.117}\]

Substitute Equation 4.116 and Equation 4.117 into the middle block row. After collecting terms, we obtain \[ \begin{aligned} 0 &= \Bigg\{ \frac{1}{\sigma_0^2} \lim_{n\to\infty} \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n \\ &\qquad+ \lim_{n\to\infty} \frac{1}{n} \left[ \operatorname{tr} \left( \mathbf G_{n0}^{\top}\mathbf G_{n0} \right) + \operatorname{tr} \left( \mathbf G_{n0}^2 \right) {}- \frac{2}{n} \left( \operatorname{tr}\mathbf G_{n0} \right)^2 \right] \Bigg\} a_{\rho}. \end{aligned} \tag{4.118}\]

The first term inside braces is strictly positive by Proposition 4.3. The second is nonnegative. To verify the latter, define \[ \mathbf C_n = \mathbf G_{n0} {}- \frac{ \operatorname{tr}\mathbf G_{n0} }{n} \mathbf I_n. \]

Then \[ \begin{aligned} &\operatorname{tr} \left( \mathbf G_{n0}^{\top}\mathbf G_{n0} \right) + \operatorname{tr} \left( \mathbf G_{n0}^2 \right) {}- \frac{2}{n} \left( \operatorname{tr}\mathbf G_{n0} \right)^2 \\ &\qquad= \frac{1}{2} \operatorname{tr} \left[ \left( \mathbf C_n+\mathbf C_n^{\top} \right)^{\top} \left( \mathbf C_n+\mathbf C_n^{\top} \right) \right] \geq0. \end{aligned} \]

Thus the coefficient multiplying \(a_{\rho}\) in Equation 4.118 is strictly positive, so \[ a_{\rho}=0. \]

Equations Equation 4.116 and Equation 4.117 then imply \[ \mathbf a_{\beta}=\mathbf0, \qquad a_{\sigma}=0. \]

Hence the null space of \(\boldsymbol\Sigma_{\theta}\) contains only the zero vector, and the matrix is positive definite.

4.6.9 Convergence of the Hessian

The asymptotic expansion of the estimator requires the Hessian evaluated between \(\widehat{\boldsymbol\theta}_n\) and \(\boldsymbol\theta_0\). Because the parameter is a vector, we use the mean value theorem component by component.

Let \(s_{j,n}(\boldsymbol\theta)\) denote component \(j\) of the score, where \(j=1,\ldots,k+2\). For every \(j\), there is a point \[ \widetilde{\boldsymbol\theta}_{n,j} = \boldsymbol\theta_0 + t_{n,j} \left( \widehat{\boldsymbol\theta}_n {}- \boldsymbol\theta_0 \right), \qquad t_{n,j}\in(0,1), \] such that \[ 0 = s_{j,n} \left( \widehat{\boldsymbol\theta}_n \right) = s_{j,n}(\boldsymbol\theta_0) + \nabla_{\theta} s_{j,n} \left( \widetilde{\boldsymbol\theta}_{n,j} \right)^{\top} \left( \widehat{\boldsymbol\theta}_n {}- \boldsymbol\theta_0 \right). \]

Stack the row gradients into the matrix \[ \overline{\mathbf H}_n = \begin{pmatrix} \nabla_{\theta} s_{1,n} (\widetilde{\boldsymbol\theta}_{n,1})^{\top} \\ \vdots \\ \nabla_{\theta} s_{k+2,n} (\widetilde{\boldsymbol\theta}_{n,k+2})^{\top} \end{pmatrix}. \tag{4.119}\]

Then \[ \sqrt n \left( \widehat{\boldsymbol\theta}_n {}- \boldsymbol\theta_0 \right) = {}- \left( \frac{1}{n} \overline{\mathbf H}_n \right)^{-1} \frac{1}{\sqrt n} \mathbf s_n(\boldsymbol\theta_0), \tag{4.120}\] whenever the matrix is nonsingular.

Lemma 4.6 (Convergence of the intermediate Hessian) Under Assumption 1–Assumption 8, \[ \frac{1}{n} \overline{\mathbf H}_n \overset{p}{\longrightarrow} {}- \boldsymbol\Sigma_{\theta}. \tag{4.121}\]

Proof. We write \[ \begin{aligned} \frac{1}{n} \overline{\mathbf H}_n + \boldsymbol\Sigma_{\theta} = \left[ \frac{1}{n} \overline{\mathbf H}_n {}- \frac{1}{n} \mathbf H_n(\boldsymbol\theta_0) \right] + \left[ \frac{1}{n} \mathbf H_n(\boldsymbol\theta_0) {}- \mathbb E \left( \frac{1}{n} \mathbf H_n(\boldsymbol\theta_0) \right) \right] + \left[ \mathbb E \left( \frac{1}{n} \mathbf H_n(\boldsymbol\theta_0) \right) + \boldsymbol\Sigma_{\theta} \right], \end{aligned} \tag{4.122}\] where \(\mathbf H_n(\boldsymbol\theta)\) denotes the Hessian derived in Section 4.3.5.

The final bracket converges to zero by the definition of \(\boldsymbol\Sigma_{\theta}\). Consider the middle bracket. The random components of the Hessian at the true parameter are sums of terms such as \[ \frac{1}{n} \mathbf X_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n, \]

\[ \frac{1}{n} \left[ \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n {}- \sigma_0^2 \operatorname{tr} \mathbf G_{n0} \right], \] and \[ \frac{1}{n} \left[ \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0}^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n {}- \sigma_0^2 \operatorname{tr} \left( \mathbf G_{n0}^{\top} \mathbf G_{n0} \right) \right]. \]

The first is \(o_p(1)\) by Lemma 3.14. The second and third are \(o_p(1)\) by Lemma 3.13. The same arguments apply to the other blocks in Equation 4.51. Therefore, \[ \frac{1}{n} \mathbf H_n(\boldsymbol\theta_0) {}- \mathbb E \left[ \frac{1}{n} \mathbf H_n(\boldsymbol\theta_0) \right] = o_p(1). \tag{4.123}\]

It remains to control the first bracket in Equation 4.122. Consistency from Theorem 4.1 implies that, for every row \(j\), \(\widetilde{\boldsymbol\theta}_{n,j}\overset{p}{\longrightarrow}\boldsymbol\theta_0\). The Hessian blocks in Section 4.3.5 are continuous functions of \(\boldsymbol\beta\), \(\rho\), and \(\sigma^2\) on a neighborhood of the true parameter. Assumption 5 and Assumption 7 provide uniform bounds for the matrix products and trace terms appearing in these blocks. The random normalized quadratic forms are \(O_p(1)\) by Lemma 3.13, while the normalized linear forms are \(O_p(1)\) by Lemma 3.14.

The most delicate block is the \(\rho\rho\) derivative. From Equation 4.48, \[ \frac{\partial^2\ell_n(\boldsymbol\theta)}{\partial\rho^2} = {}- \operatorname{tr} \left[ \mathbf G_n(\rho)^2 \right] {}- \frac{1}{\sigma^2} \mathbf y_n^{\top} \mathbf W_n^{\top}\mathbf W_n \mathbf y_n. \]

For the trace term, apply the scalar mean value theorem to the function \(\rho\mapsto\operatorname{tr}[\mathbf G_n(\rho)^2]\). Because \[ \frac{\partial\mathbf G_n(\rho)}{\partial\rho} = \mathbf G_n(\rho)^2, \] there is a scalar intermediate value \(\bar\rho_n\) such that \[ \operatorname{tr} \left[ \mathbf G_n(\widetilde\rho_n)^2 \right] {}- \operatorname{tr} \left( \mathbf G_{n0}^2 \right) = 2 \operatorname{tr} \left[ \mathbf G_n(\bar\rho_n)^3 \right] (\widetilde\rho_n-\rho_0). \]

The trace is \(O(n)\) by Lemma 3.9, uniformly in a neighborhood of \(\rho_0\). Therefore, after division by \(n\), this difference is \(o_p(1)\). Moreover, \[ \frac{1}{n} \mathbf y_n^{\top} \mathbf W_n^{\top}\mathbf W_n \mathbf y_n = O_p(1), \] which follows by substituting the reduced form and applying the linear- and quadratic-form bounds from Chapter 3. Since every intermediate variance satisfies \(\widetilde\sigma_n^2\to_p\sigma_0^2\), the change in the second part of the \(\rho\rho\) block is also \(o_p(1)\).

The remaining blocks are simpler. They contain normalized quantities such as \(\mathbf X_n^{\top}\mathbf X_n/n\), \(\mathbf X_n^{\top}\mathbf W_n\mathbf y_n/n\), \(\mathbf X_n^{\top}\boldsymbol\varepsilon_n/n\), and \(\boldsymbol\varepsilon_n^{\top}\boldsymbol\varepsilon_n/n\), multiplied by continuous functions of the parameter. Their stated probability orders and consistency of the intermediate parameters imply an \(o_p(1)\) change in every case. Consequently, changing the parameter from \(\boldsymbol\theta_0\) to the row-specific consistent intermediate points changes every normalized Hessian block by \(o_p(1)\). Hence \[ \frac{1}{n} \overline{\mathbf H}_n {}- \frac{1}{n} \mathbf H_n(\boldsymbol\theta_0) = o_p(1). \tag{4.124}\]

Combining the three terms in Equation 4.122 proves Equation 4.121.

4.6.10 Asymptotic Normality of the Regular QMLE

Theorem 4.2 (Asymptotic normality of the regular SLM QMLE) Suppose Assumption 1–Assumption 8 hold. Assume \[ \boldsymbol\Sigma_{\theta,n} \longrightarrow \boldsymbol\Sigma_{\theta}, \qquad \boldsymbol\Omega_{\theta,n} \longrightarrow \boldsymbol\Omega_{\theta}, \] and define \[ \mathbf V_{\theta} = \boldsymbol\Sigma_{\theta} + \boldsymbol\Omega_{\theta}. \]

Suppose \(\mathbf V_{\theta}\) is finite and positive definite. Then \[ \sqrt n \left( \widehat{\boldsymbol\theta}_n {}- \boldsymbol\theta_0 \right) \overset{d}{\longrightarrow} \mathcal N \left[ \mathbf0, \boldsymbol\Sigma_{\theta}^{-1} \left( \boldsymbol\Sigma_{\theta} + \boldsymbol\Omega_{\theta} \right) \boldsymbol\Sigma_{\theta}^{-1} \right]. \tag{4.125}\]

Equivalently, the asymptotic covariance matrix is \[ \boldsymbol\Sigma_{\theta}^{-1} + \boldsymbol\Sigma_{\theta}^{-1} \boldsymbol\Omega_{\theta} \boldsymbol\Sigma_{\theta}^{-1}. \tag{4.126}\]

If the innovations are Gaussian, then \(\boldsymbol\Omega_{\theta}=\mathbf0\) and \[ \sqrt n \left( \widehat{\boldsymbol\theta}_n {}- \boldsymbol\theta_0 \right) \overset{d}{\longrightarrow} \mathcal N \left( \mathbf0, \boldsymbol\Sigma_{\theta}^{-1} \right). \tag{4.127}\]

Proof. From the componentwise mean-value expansion in Equation 4.120,

\[ \sqrt n \left( \widehat{\boldsymbol\theta}_n {}- \boldsymbol\theta_0 \right) = {}- \left( \frac{1}{n} \overline{\mathbf H}_n \right)^{-1} \frac{1}{\sqrt n} \mathbf s_n(\boldsymbol\theta_0). \]

By Lemma 4.6,

\[ \frac{1}{n} \overline{\mathbf H}_n \overset{p}{\longrightarrow} {}- \boldsymbol\Sigma_{\theta}. \]

By Lemma 4.5, \(\boldsymbol\Sigma_{\theta}\) is nonsingular. The stochastic inversion result Lemma 3.4 therefore gives

\[ {}- \left( \frac{1}{n} \overline{\mathbf H}_n \right)^{-1} \overset{p}{\longrightarrow} \boldsymbol\Sigma_{\theta}^{-1}. \]

By Lemma 4.4,

\[ \frac{1}{\sqrt n} \mathbf s_n(\boldsymbol\theta_0) \overset{d}{\longrightarrow} \mathcal N \left( \mathbf0, \boldsymbol\Sigma_{\theta} + \boldsymbol\Omega_{\theta} \right). \]

Applying the matrix form of Slutsky’s theorem in Corollary 3.2 gives Equation 4.125. Under Gaussian innovations, Equation 4.106 implies \(\boldsymbol\Omega_{\theta}=\mathbf0\), which gives Equation 4.127.

4.6.11 Estimating the Asymptotic Covariance Matrix

Under Gaussian normality, a consistent covariance estimator is obtained from the inverse observed information or the inverse expected information evaluated at the MLE: \[ \widehat{\operatorname{Var}} \left( \widehat{\boldsymbol\theta}_n \right) = \frac{1}{n} \widehat{\boldsymbol\Sigma}_{\theta,n}^{-1}. \]

The factor \(1/n\) is needed because \(\boldsymbol\Sigma_{\theta,n}\) is the information per observation.

For Gaussian QMLE under non-Gaussian innovations, the appropriate estimator is the sandwich matrix

\[ \widehat{\operatorname{Var}} \left( \widehat{\boldsymbol\theta}_n \right) = \frac{1}{n} \widehat{\boldsymbol\Sigma}_{\theta,n}^{-1} \widehat{\mathbf V}_{\theta,n} \widehat{\boldsymbol\Sigma}_{\theta,n}^{-1}, \tag{4.128}\]

where

\[ \widehat{\mathbf V}_{\theta,n} = \widehat{\boldsymbol\Sigma}_{\theta,n} + \widehat{\boldsymbol\Omega}_{\theta,n}. \]

The moments \(\mu_3\) and \(\mu_4\) in Equation 4.104 can be estimated using the third and fourth empirical moments of the fitted innovations. The resulting estimator reflects departures from Gaussian skewness and kurtosis.

The covariance matrix in Equation 4.128 is a model-based spatial QMLE covariance matrix. It should not be confused with an arbitrary cross-sectional heteroskedasticity or dependence correction. The present theory assumes i.i.d. innovations and obtains cross-sectional dependence through the spatial multiplier.

4.6.12 Bounded and Divergent \(h_n\): What Changes?

The regular theorem does not require \(h_n\) to remain bounded. To see why, consider the \(\rho\) score in Equation 4.102:

\[ \frac{1}{\sigma_0^2\sqrt n} \mathbf d_n^{\top} \boldsymbol\varepsilon_n + \frac{1}{\sigma_0^2\sqrt n} \left[ \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n {}- \sigma_0^2 \operatorname{tr}\mathbf G_{n0} \right]. \tag{4.129}\]

Under regular identification,

\[ \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n \longrightarrow q_d>0. \]

The linear term therefore has variance of order \(n\) before division by \(\sqrt n\). It supplies root-\(n\) information about \(\rho_0\) regardless of whether \(h_n\) is bounded or divergent.

The quadratic term behaves differently. Under the elementwise order \(O(1/h_n)\) and the matrix-sum bounds,

\[ \operatorname{Var} \left( \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n \right) = O\!\left(\frac{n}{h_n}\right). \]

Consequently,

\[ \frac{1}{\sqrt n} \left[ \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n {}- \sigma_0^2 \operatorname{tr}\mathbf G_{n0} \right] = O_p \left( \frac{1}{\sqrt{h_n}} \right). \tag{4.130}\]

If \(h_n\) is bounded, the quadratic term can contribute to the root-\(n\) limiting distribution. If \(h_n\to\infty\), it becomes negligible under root-\(n\) normalization. The linear mean component remains, so the regular QMLE still has the root-\(n\) rate.

The slower rate arises when regular mean identification fails. The main regimes in Lee can be summarized as follows.

Identification and spatial design	Main information for \(\rho_0\)	Typical rate for \(\widehat\rho_n\)
ML8 holds; \(h_n\) bounded	Linear and quadratic score components	\(\sqrt n\)
ML8 holds; \(h_n\to\infty\)	Linear score component dominates	\(\sqrt n\)
Mean identification fails; \(h_n\) bounded	Covariance structure can still identify	Often \(\sqrt n\) under an additional variance-identification condition
Mean identification fails; \(h_n\to\infty\) and \(h_n/n\to0\)	Weak quadratic information	Typically \(\sqrt{n/h_n}\)
\(h_n/n\) does not vanish	Information may be too weak	Consistency can fail

The pure SAR process is not covered by ML8

For the pure SAR process,

\[ \mathbf y_n = \rho_0\mathbf W_n\mathbf y_n + \boldsymbol\varepsilon_n, \]

we have \(\mathbf d_n=\mathbf0\). Thus ML8 fails automatically. Identification must come from the covariance structure. With bounded \(h_n\), root-\(n\) estimation can still be possible. With divergent \(h_n\), Lee obtains the slower rate \(\sqrt{n/h_n}\) under additional conditions.

The present chapter proves the regular root-\(n\) case. The irregular theory requires a different normalization of the criterion, score, and Hessian. In particular, when \(h_n\to\infty\) and mean identification is weak, the relevant expansion is of the form

\[ \sqrt{\frac{n}{h_n}} \left( \widehat\rho_n-\rho_0 \right) = {}- \left[ \frac{h_n}{n} \frac{\partial^2 \ell_n^c(\widetilde\rho_n) }{ \partial\rho^2 } \right]^{-1} \sqrt{\frac{h_n}{n}} \frac{\partial \ell_n^c(\rho_0) }{ \partial\rho }. \tag{4.131}\]

This is not merely the regular proof with \(\sqrt n\) replaced mechanically by \(\sqrt{n/h_n}\). The identification condition, variance order, and linear–quadratic CLT must all be reformulated at the new rate.

4.6.13 What the Regular Theory Establishes

The argument can be summarized in five steps.

Assumption 8 ensures that the spatially generated regressor \(\mathbf G_{n0}\mathbf X_n\boldsymbol\beta_0\) contains nonredundant variation after projection on \(\mathbf X_n\).
The population concentrated criterion is separated from its value at \(\rho_0\) outside every neighborhood of the truth.
Linear- and quadratic-form laws of large numbers from Chapter 3 imply uniform convergence of the sample concentrated criterion.
The score is a vector of linear and quadratic forms, so Theorem 3.26 provides its limiting normal distribution.
The componentwise mean-value expansion, Hessian convergence, stochastic inversion, and Slutsky’s theorem produce the limiting distribution of the QMLE.

The key conclusion is

\[ \boxed{ \text{regular identification} \quad\Longrightarrow\quad \sqrt n \left( \widehat{\boldsymbol\theta}_n {}- \boldsymbol\theta_0 \right) =O_p(1) } \] whether \(h_n\) is bounded or divergent. Divergence of \(h_n\) changes the relative importance of the linear and quadratic score components, but it does not reduce the convergence rate while regular mean identification remains available.

The next section applies this asymptotic distribution to nonlinear functions of the SLM parameters, including average direct, indirect, and total effects.

4.7 Statistical Inference in the Regular Spatial Lag Model

The asymptotic results in Theorem 4.2 provide the basis for confidence intervals and hypothesis tests. Throughout this section, the regular root-\(n\) theory is assumed. In particular, Assumptions Assumption 1 through Assumption 8 hold, the true parameter lies in the interior of the parameter space, and the asymptotic covariance matrix is nonsingular.

Let \[ \boldsymbol\theta = \begin{pmatrix} \boldsymbol\beta\\ \rho\\ \sigma^2 \end{pmatrix} \in\mathbb{R}^{k+2}, \] and suppose that \[ \sqrt n \left( \widehat{\boldsymbol\theta}_n - \boldsymbol\theta_0 \right) \overset{d}{\longrightarrow} \mathcal N \left( \mathbf{0}, \mathbf V_{\theta} \right). \tag{4.132}\]

For Gaussian maximum likelihood, \[ \mathbf V_{\theta} = \boldsymbol\Sigma_{\theta}^{-1}, \] whereas Gaussian quasi-maximum likelihood under non-Gaussian innovations gives \[ \mathbf V_{\theta} = \boldsymbol\Sigma_{\theta}^{-1} \left( \boldsymbol\Sigma_{\theta} + \boldsymbol\Omega_{\theta} \right) \boldsymbol\Sigma_{\theta}^{-1}. \]

The distinction matters for inference. The inverse information matrix is valid when the Gaussian likelihood is correctly specified. Under non-Gaussian innovations, the sandwich covariance matrix must be used.

4.7.1 Standard Errors and Confidence Intervals

Let \(\theta_{0,j}\) denote the \(j\)th component of \(\boldsymbol\theta_0\). From Equation 4.132, \[ \frac{ \widehat\theta_{n,j}-\theta_{0,j} }{ \sqrt{ \widehat{\operatorname{Var}} \left( \widehat\theta_{n,j} \right) } } \overset{d}{\longrightarrow} \mathcal N(0,1), \] provided that the estimated variance is consistent.

An asymptotic \(100(1-\alpha)\%\) confidence interval is therefore \[ \widehat\theta_{n,j} \pm z_{1-\alpha/2} \sqrt{ \widehat{\operatorname{Var}} \left( \widehat\theta_{n,j} \right) }, \tag{4.133}\] where \(z_{1-\alpha/2}\) is the corresponding standard-normal quantile.

Two covariance conventions

It is essential to distinguish between the covariance matrix of the estimator and the asymptotic covariance matrix of its root-\(n\) normalization.

If \[ \widehat{\mathbf V}_{\theta,n} \overset{p}{\longrightarrow} \mathbf V_{\theta} \] estimates the covariance of \[ \sqrt n \left( \widehat{\boldsymbol\theta}_n-\boldsymbol\theta_0 \right), \] then \[ \widehat{\operatorname{Var}} \left( \widehat{\boldsymbol\theta}_n \right) = \frac{1}{n} \widehat{\mathbf V}_{\theta,n}. \]

By contrast, software often reports \(\widehat{\operatorname{Var}}(\widehat{\boldsymbol\theta}_n)\) directly. An additional factor \(1/n\) must not be applied to a covariance matrix that is already expressed on the estimator scale.

4.7.2 Wald Tests for Smooth Restrictions

Many hypotheses of interest can be written as \[ H_0: \mathbf r \left( \boldsymbol\theta_0 \right) = \mathbf{0}, \] where \[ \mathbf r: \mathbb{R}^{k+2} \longrightarrow \mathbb{R}^q \] is continuously differentiable in a neighborhood of \(\boldsymbol\theta_0\). Define the \(q\times(k+2)\) Jacobian \[ \mathbf R \left( \boldsymbol\theta \right) = \frac{ \partial \mathbf r \left( \boldsymbol\theta \right) }{ \partial \boldsymbol\theta^{\top} }. \]

Examples include:

\(H_0:\rho_0=0\);
equality of two slope coefficients;
exclusion restrictions on several regressors;
nonlinear restrictions involving \(\rho_0\) and \(\boldsymbol\beta_0\).

Theorem 4.3 (Wald test under regular QMLE) Suppose that Equation 4.132 holds, that \(\mathbf r(\boldsymbol\theta)\) is continuously differentiable in a neighborhood of \(\boldsymbol\theta_0\), and that \[ \mathbf R_0 \mathbf V_{\theta} \mathbf R_0^{\top} \] is nonsingular, where \[ \mathbf R_0 = \mathbf R \left( \boldsymbol\theta_0 \right). \]

Let \[ \widehat{\mathbf C}_{r,n} = \mathbf R \left( \widehat{\boldsymbol\theta}_n \right) \widehat{\operatorname{Var}} \left( \widehat{\boldsymbol\theta}_n \right) \mathbf R \left( \widehat{\boldsymbol\theta}_n \right)^{\top}. \]

Under \(H_0\), \[ W_n = \mathbf r \left( \widehat{\boldsymbol\theta}_n \right)^{\top} \widehat{\mathbf C}_{r,n}^{-1} \mathbf r \left( \widehat{\boldsymbol\theta}_n \right) \overset{d}{\longrightarrow} \chi_q^2. \tag{4.134}\]

Proof. Apply a first-order expansion separately to each component of \(\mathbf r\). For component \(a=1,\ldots,q\), there exists a point \(\widetilde{\boldsymbol\theta}_{n,a}\) on the segment joining \(\widehat{\boldsymbol\theta}_n\) and \(\boldsymbol\theta_0\) such that

\[ r_a \left( \widehat{\boldsymbol\theta}_n \right) - r_a \left( \boldsymbol\theta_0 \right) = \frac{ \partial r_a \left( \widetilde{\boldsymbol\theta}_{n,a} \right) }{ \partial \boldsymbol\theta^{\top} } \left( \widehat{\boldsymbol\theta}_n - \boldsymbol\theta_0 \right). \]

Stacking the \(q\) componentwise expansions gives \[ \sqrt n \left[ \mathbf r \left( \widehat{\boldsymbol\theta}_n \right) - \mathbf r \left( \boldsymbol\theta_0 \right) \right] = \widetilde{\mathbf R}_n \sqrt n \left( \widehat{\boldsymbol\theta}_n - \boldsymbol\theta_0 \right), \] where row \(a\) of \(\widetilde{\mathbf R}_n\) is evaluated at \(\widetilde{\boldsymbol\theta}_{n,a}\).

Consistency of \(\widehat{\boldsymbol\theta}_n\) and continuity of the derivatives imply \[ \widetilde{\mathbf R}_n \overset{p}{\longrightarrow} \mathbf R_0. \]

Therefore, under \(H_0\), \[ \sqrt n \mathbf r \left( \widehat{\boldsymbol\theta}_n \right) \overset{d}{\longrightarrow} \mathcal N \left( \mathbf{0}, \mathbf R_0 \mathbf V_{\theta} \mathbf R_0^{\top} \right). \]

The estimated covariance in Equation 4.134 is consistent for the covariance of \(\mathbf r(\widehat{\boldsymbol\theta}_n)\). Standardization by this covariance and the continuous mapping theorem give the chi-square limit.

For a scalar restriction \(q=1\), the Wald statistic is the square of the corresponding asymptotic \(z\)-statistic.

4.7.3 Likelihood-Ratio and Score Tests

When the innovations are Gaussian and the likelihood is correctly specified, the usual likelihood-ratio and score tests are also available.

Let \(\widehat{\boldsymbol\theta}_n\) be the unrestricted MLE and let \(\widetilde{\boldsymbol\theta}_n\) be the restricted MLE under \[ \mathbf r \left( \boldsymbol\theta \right) = \mathbf{0}. \]

The likelihood-ratio statistic is \[ LR_n = 2 \left[ \ell_n \left( \widehat{\boldsymbol\theta}_n \right) - \ell_n \left( \widetilde{\boldsymbol\theta}_n \right) \right]. \tag{4.135}\]

Under the regular Gaussian ML conditions, \[ LR_n \overset{d}{\longrightarrow} \chi_q^2. \]

Let \[ \mathbf s_n \left( \boldsymbol\theta \right) = \frac{ \partial \ell_n \left( \boldsymbol\theta \right) }{ \partial \boldsymbol\theta } \] and let \[ \mathbf I_n \left( \boldsymbol\theta \right) = - \frac{ \partial^2 \ell_n \left( \boldsymbol\theta \right) }{ \partial \boldsymbol\theta \partial \boldsymbol\theta^{\top} } \] denote the observed information. The unrestricted form of the score statistic is \[ LM_n = \mathbf s_n \left( \widetilde{\boldsymbol\theta}_n \right)^{\top} \mathbf I_n \left( \widetilde{\boldsymbol\theta}_n \right)^{-1} \mathbf s_n \left( \widetilde{\boldsymbol\theta}_n \right), \tag{4.136}\]

with the usual adjustment for nuisance parameters implicit in the inverse full information matrix. Under the null and correct Gaussian specification, \[ LM_n \overset{d}{\longrightarrow} \chi_q^2. \]

Quasi-likelihood and likelihood-ratio tests

Under non-Gaussian innovations, the Gaussian criterion is a quasi-likelihood. The sandwich Wald test remains directly valid under the QMLE theory developed above. The unadjusted likelihood-ratio and score statistics do not generally retain their standard chi-square limits under arbitrary likelihood misspecification. They should therefore not be used as if the Gaussian density were correctly specified unless the additional conditions required for their validity have been established.

Standard chi-square inference also requires that the true parameter be in the interior of the parameter space. Tests at a boundary require a different limiting theory.

4.8 Inference for Spatial Impacts

Chapter 2 derived the impact matrix and the average direct, indirect, and total impacts. The purpose of this section is to derive their sampling uncertainty from the asymptotic distribution of the maximum likelihood or quasi-maximum likelihood estimator.

Two points require care. First, the regular asymptotic theory above was stated for the SLM. Before using it for the SDM, we must show how the SDM fits the same likelihood framework.

Second, the notation must remain consistent with Chapter 2. In that chapter, \(\mathbf D_r\) denotes the Jacobian collecting the \(n^2\) partial derivatives. It must not be reused for the coefficient matrix \(\beta_r\mathbf I_n+\gamma_r\mathbf W_n\).

4.8.1 The SDM as an SLM with an Augmented Regressor Matrix

Consider the SDM \[ \mathbf y_n = \rho \mathbf W_n \mathbf y_n + \alpha \boldsymbol\iota_n + \mathbf X_n \boldsymbol\beta + \mathbf W_n \mathbf X_n \boldsymbol\gamma + \boldsymbol\varepsilon_n, \tag{4.137}\] where \(\mathbf X_n\) does not contain the intercept. Define \[ \mathbf Z_n = \begin{pmatrix} \boldsymbol\iota_n & \mathbf X_n & \mathbf W_n\mathbf X_n \end{pmatrix} \tag{4.138}\] and \[ \boldsymbol\delta = \begin{pmatrix} \alpha\\ \boldsymbol\beta\\ \boldsymbol\gamma \end{pmatrix}. \tag{4.139}\]

Then \[ \mathbf y_n = \rho \mathbf W_n \mathbf y_n + \mathbf Z_n \boldsymbol\delta + \boldsymbol\varepsilon_n. \tag{4.140}\]

Proposition 4.4 (Regular QMLE theory for the SDM) Suppose that the regular SLM assumptions hold after replacing \(\mathbf X_n\) by \(\mathbf Z_n\) and \(\boldsymbol\beta\) by \(\boldsymbol\delta\). In particular, assume that:

the elements of \(\mathbf Z_n\) are uniformly bounded;
\[ \frac{1}{n} \mathbf Z_n^{\top} \mathbf Z_n \longrightarrow \mathbf Q_Z, \] where \(\mathbf Q_Z\) is nonsingular;
the regular identification condition holds for \[ \mathbf d_{Z,n} = \mathbf G_{n0} \mathbf Z_n \boldsymbol\delta_0. \]

Then the consistency and regular root-\(n\) asymptotic-normality results derived for the SLM apply to the SDM after making these substitutions.

Proof. Equation Equation 4.140 has exactly the same algebraic form as the SLM studied above. Conditional on the fixed matrices \(\mathbf W_n\) and \(\mathbf Z_n\), its Gaussian log-likelihood is obtained from the SLM likelihood by replacing \[ \mathbf X_n\boldsymbol\beta \quad\text{with}\quad \mathbf Z_n\boldsymbol\delta. \]

The concentrated estimators, score, Hessian, population criterion, and regular identification argument consequently have the same form with the same replacement. The stated conditions ensure that the augmented design satisfies the regressor and identification assumptions used in the SLM proofs. Therefore, the preceding consistency and asymptotic-normality theorems apply.

If \(\mathbf W_n\) is uniformly bounded in absolute row sums and the elements of \(\mathbf X_n\) are uniformly bounded, then the elements of \(\mathbf W_n\mathbf X_n\) are also uniformly bounded. Nonsingularity of \(\mathbf Z_n^{\top}\mathbf Z_n/n\), however, is an additional rank condition; it does not follow from boundedness alone.

When the weights are row standardized, \[ \mathbf W_n\boldsymbol\iota_n = \boldsymbol\iota_n. \]

Therefore, the spatial lag of the intercept must not be added as a separate column of \(\mathbf Z_n\), because it would be perfectly collinear with the intercept itself.

4.8.2 Impact Definitions from Chapter 2

For compatibility with the likelihood sections of this chapter, define \(\mathbf A_n(\rho)=\mathbf I_n-\rho\mathbf W_n\). Chapter 2 used \(\mathbf B(\rho)\) for this same spatial transformation. Thus, the multiplier denoted there by \(\mathbf B(\rho)^{-1}\) is denoted here by \(\mathbf A_n(\rho)^{-1}\).

For a non-intercept regressor \(r\), define \[ \mathbf K_{r,n} \left( \beta_r,\gamma_r \right) = \beta_r\mathbf I_n + \gamma_r\mathbf W_n. \tag{4.141}\]

The impact matrix is \[ \mathbf S_{r,n} \left( \rho,\beta_r,\gamma_r \right) = \mathbf A_n(\rho)^{-1} \mathbf K_{r,n} \left( \beta_r,\gamma_r \right). \tag{4.142}\]

This is precisely the matrix derived in Proposition 2.4 and Equation 2.43. Its \((i,j)\) element is \[ \frac{ \partial \mathbb E \left( y_{n,i} \mid \mathbf X_n,\mathbf W_n \right) }{ \partial x_{n,jr} }. \]

The orientation established in Chapter 2 remains unchanged:

column \(j\) records impacts from unit \(j\);
row \(i\) records impacts to unit \(i\);
diagonal elements are unit-specific direct effects;
off-diagonal elements are unit-specific indirect effects.

Let \[ \mathcal D_n \left( \mathbf H_n \right) = \frac{1}{n} \operatorname{tr} \left( \mathbf H_n \right) \tag{4.143}\] and \[ \mathcal T_n \left( \mathbf H_n \right) = \frac{1}{n} \boldsymbol\iota_n^{\top} \mathbf H_n \boldsymbol\iota_n. \tag{4.144}\]

Using the definitions in Definition 2.3–Definition 2.5, \[ \operatorname{ADI}_{r,n} = \mathcal D_n \left( \mathbf S_{r,n} \right), \tag{4.145}\]

\[ \operatorname{ATI}_{r,n} = \mathcal T_n \left( \mathbf S_{r,n} \right), \tag{4.146}\] and \[ \operatorname{AII}_{r,n} = \operatorname{ATI}_{r,n} - \operatorname{ADI}_{r,n}. \tag{4.147}\]

As emphasized in Chapter 2, \(\operatorname{ATI}_{r,n}\) is the average of the \(n\) row sums, equivalently the average of the \(n\) column sums. It is not the arithmetic average of all \(n^2\) entries of the impact matrix.

4.8.3 Derivative of the Inverse Spatial Transformation

Lemma 4.7 (Derivative of the inverse spatial transformation) Suppose that \[ \mathbf A_n(\rho) = \mathbf I_n-\rho\mathbf W_n \] is nonsingular on an open interval. Then \[ \frac{ \partial \mathbf A_n(\rho)^{-1} }{ \partial\rho } = \mathbf A_n(\rho)^{-1} \mathbf W_n \mathbf A_n(\rho)^{-1}. \tag{4.148}\]

Proof. The inverse identity is \[ \mathbf A_n(\rho) \mathbf A_n(\rho)^{-1} = \mathbf I_n. \]

Differentiate both sides with respect to \(\rho\): \[ \frac{ \partial\mathbf A_n(\rho) }{ \partial\rho } \mathbf A_n(\rho)^{-1} + \mathbf A_n(\rho) \frac{ \partial\mathbf A_n(\rho)^{-1} }{ \partial\rho } = \mathbf0. \]

Since \[ \frac{ \partial\mathbf A_n(\rho) }{ \partial\rho } = -\mathbf W_n, \] we obtain \[ -\mathbf W_n \mathbf A_n(\rho)^{-1} + \mathbf A_n(\rho) \frac{ \partial\mathbf A_n(\rho)^{-1} }{ \partial\rho } = \mathbf0. \]

Premultiplying by \(\mathbf A_n(\rho)^{-1}\) gives Equation 4.148.

No symmetry of \(\mathbf W_n\) is required for this derivative.

4.8.4 Exact Gradients of the Average Impacts

Define the parameter subvector \[ \boldsymbol\eta_r = \begin{pmatrix} \beta_r\\ \gamma_r\\ \rho \end{pmatrix}. \tag{4.149}\]

From Equation 4.142 and Equation 4.148, \[ \frac{ \partial\mathbf S_{r,n} }{ \partial\beta_r } = \mathbf A_n(\rho)^{-1}, \tag{4.150}\]

\[ \frac{ \partial\mathbf S_{r,n} }{ \partial\gamma_r } = \mathbf A_n(\rho)^{-1} \mathbf W_n, \tag{4.151}\] and \[ \frac{ \partial\mathbf S_{r,n} }{ \partial\rho } = \mathbf A_n(\rho)^{-1} \mathbf W_n \mathbf A_n(\rho)^{-1} \mathbf K_{r,n} \left( \beta_r,\gamma_r \right). \tag{4.152}\]

Proposition 4.5 (Gradients of the average SDM impacts) The gradient of the average direct impact is \[ \mathbf g_{D,r,n}^{\top} = \begin{pmatrix} \mathcal D_n \left[ \mathbf A_n(\rho)^{-1} \right] & \mathcal D_n \left[ \mathbf A_n(\rho)^{-1} \mathbf W_n \right] & \mathcal D_n \left[ \mathbf A_n(\rho)^{-1} \mathbf W_n \mathbf A_n(\rho)^{-1} \mathbf K_{r,n} \right] \end{pmatrix}. \tag{4.153}\]

The gradient of the average total impact is \[ \mathbf g_{T,r,n}^{\top} = \begin{pmatrix} \mathcal T_n \left[ \mathbf A_n(\rho)^{-1} \right] & \mathcal T_n \left[ \mathbf A_n(\rho)^{-1} \mathbf W_n \right] & \mathcal T_n \left[ \mathbf A_n(\rho)^{-1} \mathbf W_n \mathbf A_n(\rho)^{-1} \mathbf K_{r,n} \right] \end{pmatrix}. \tag{4.154}\]

The gradient of the average indirect impact is \[ \mathbf g_{I,r,n} = \mathbf g_{T,r,n} - \mathbf g_{D,r,n}. \tag{4.155}\]

Proof. The functionals \(\mathcal D_n(\cdot)\) and \(\mathcal T_n(\cdot)\) are linear. Therefore, \[ \frac{ \partial\operatorname{ADI}_{r,n} }{ \partial\beta_r } = \mathcal D_n \left( \frac{ \partial\mathbf S_{r,n} }{ \partial\beta_r } \right), \] and the corresponding identities hold for \(\gamma_r\) and \(\rho\). Substitution of Equation 4.150–Equation 4.152 gives Equation 4.153.

Applying the same argument to \(\mathcal T_n(\cdot)\) gives Equation 4.154. Finally, Equation 4.147 implies Equation 4.155.

These are exact analytic gradients. The nonlinearity of the multiplier does not make the delta method unavailable.

4.8.5 SLM Gradients

For the SLM, \(\gamma_r\) is absent and \[ \mathbf S_{r,n}^{\mathrm{SLM}} \left( \rho,\beta_r \right) = \beta_r \mathbf A_n(\rho)^{-1}. \tag{4.156}\]

With \[ \boldsymbol\eta_r^{\mathrm{SLM}} = \begin{pmatrix} \beta_r\\ \rho \end{pmatrix}, \] the direct-impact gradient is \[ \mathbf g_{D,r,n}^{\mathrm{SLM}\top} = \begin{pmatrix} \mathcal D_n \left[ \mathbf A_n(\rho)^{-1} \right] & \beta_r \mathcal D_n \left[ \mathbf A_n(\rho)^{-1} \mathbf W_n \mathbf A_n(\rho)^{-1} \right] \end{pmatrix}, \tag{4.157}\] and the total-impact gradient is \[ \mathbf g_{T,r,n}^{\mathrm{SLM}\top} = \begin{pmatrix} \mathcal T_n \left[ \mathbf A_n(\rho)^{-1} \right] & \beta_r \mathcal T_n \left[ \mathbf A_n(\rho)^{-1} \mathbf W_n \mathbf A_n(\rho)^{-1} \right] \end{pmatrix}. \tag{4.158}\]

The indirect-impact gradient is their difference.

4.8.6 Row-Standardized Weights

Under the row-standardization condition \[ \mathbf W_n \boldsymbol\iota_n = \boldsymbol\iota_n, \]

Proposition 2.6 gives \[ \operatorname{ATI}_{r,n} = \frac{ \beta_r+\gamma_r }{ 1-\rho }. \tag{4.159}\]

Differentiating this result gives \[ \mathbf g_{T,r,n}^{\top} = \begin{pmatrix} \dfrac{1}{1-\rho} & \dfrac{1}{1-\rho} & \dfrac{\beta_r+\gamma_r}{(1-\rho)^2} \end{pmatrix}. \tag{4.160}\]

For the SLM, \[ \operatorname{ATI}_{r,n}^{\mathrm{SLM}} = \frac{\beta_r}{1-\rho}, \] and \[ \mathbf g_{T,r,n}^{\mathrm{SLM}\top} = \begin{pmatrix} \dfrac{1}{1-\rho} & \dfrac{\beta_r}{(1-\rho)^2} \end{pmatrix}. \tag{4.161}\]

Unlike the total impact, the average direct impact generally depends on the complete diagonal structure of \(\mathbf A_n(\rho)^{-1}\). It does not reduce to a function of \((\beta_r,\gamma_r,\rho)\) alone that is independent of \(\mathbf W_n\).

4.9 Delta Method for a Sequence of Spatial Impact Functions

The weights matrix changes with \(n\), so the impact function is itself a sequence. Define \[ \mathbf m_{r,n} \left( \boldsymbol\eta_r \right) = \begin{pmatrix} \operatorname{ADI}_{r,n}\\ \operatorname{AII}_{r,n}\\ \operatorname{ATI}_{r,n} \end{pmatrix}. \tag{4.162}\]

Let \[ \mathbf G_{r,n} \left( \boldsymbol\eta_r \right) = \frac{ \partial \mathbf m_{r,n} \left( \boldsymbol\eta_r \right) }{ \partial \boldsymbol\eta_r^{\top} } \tag{4.163}\] be its Jacobian.

Theorem 4.4 (Delta method for average spatial impacts) Suppose that

\[ \sqrt n \left( \widehat{\boldsymbol\eta}_{r,n} - \boldsymbol\eta_{r,0} \right) \overset{d}{\longrightarrow} \mathcal N \left( \mathbf0, \mathbf V_{\eta_r} \right). \]

Assume that:

every component of \(\mathbf m_{r,n}(\boldsymbol\eta_r)\) is continuously differentiable on a neighborhood of \(\boldsymbol\eta_{r,0}\);
\[ \mathbf G_{r,n} \left( \boldsymbol\eta_{r,0} \right) \longrightarrow \mathbf G_{r,0}; \]
uniformly over a shrinking neighborhood of \(\boldsymbol\eta_{r,0}\), \[ \left\| \mathbf G_{r,n} \left( \boldsymbol\eta_r \right) - \mathbf G_{r,n} \left( \boldsymbol\eta_{r,0} \right) \right\| \longrightarrow 0. \]

Then \[ \sqrt n \left[ \mathbf m_{r,n} \left( \widehat{\boldsymbol\eta}_{r,n} \right) - \mathbf m_{r,n} \left( \boldsymbol\eta_{r,0} \right) \right] \overset{d}{\longrightarrow} \mathcal N \left( \mathbf0, \mathbf G_{r,0} \mathbf V_{\eta_r} \mathbf G_{r,0}^{\top} \right). \tag{4.164}\]

The limiting covariance matrix may be singular.

Proof. Let \(m_{r,n,a}\) be component \(a\) of \(\mathbf m_{r,n}\), for \(a=1,2,3\). Apply the scalar mean value theorem to each component separately: \[ m_{r,n,a} \left( \widehat{\boldsymbol\eta}_{r,n} \right) - m_{r,n,a} \left( \boldsymbol\eta_{r,0} \right) = \frac{ \partial m_{r,n,a} \left( \widetilde{\boldsymbol\eta}_{r,n,a} \right) }{ \partial \boldsymbol\eta_r^{\top} } \left( \widehat{\boldsymbol\eta}_{r,n} - \boldsymbol\eta_{r,0} \right), \] where \(\widetilde{\boldsymbol\eta}_{r,n,a}\) lies on the segment joining \(\widehat{\boldsymbol\eta}_{r,n}\) and \(\boldsymbol\eta_{r,0}\).

Stacking the three componentwise expansions yields \[ \sqrt n \left[ \mathbf m_{r,n} \left( \widehat{\boldsymbol\eta}_{r,n} \right) - \mathbf m_{r,n} \left( \boldsymbol\eta_{r,0} \right) \right] = \widetilde{\mathbf G}_{r,n} \sqrt n \left( \widehat{\boldsymbol\eta}_{r,n} - \boldsymbol\eta_{r,0} \right), \] where row \(a\) of \(\widetilde{\mathbf G}_{r,n}\) is evaluated at \(\widetilde{\boldsymbol\eta}_{r,n,a}\).

Root-\(n\) consistency and the uniform derivative condition imply \[ \widetilde{\mathbf G}_{r,n} - \mathbf G_{r,n} \left( \boldsymbol\eta_{r,0} \right) = o_p(1). \]

Together with the assumed limit of the Jacobian, \[ \widetilde{\mathbf G}_{r,n} \overset{p}{\longrightarrow} \mathbf G_{r,0}. \]

The result follows from Slutsky’s theorem.

The possible singularity is not a technical accident. By definition, \[ \operatorname{ATI}_{r,n} = \operatorname{ADI}_{r,n} + \operatorname{AII}_{r,n}. \]

Therefore, the three elements of \(\mathbf m_{r,n}\) are exactly linearly dependent, and \[ \mathbf g_{T,r,n} = \mathbf g_{D,r,n} + \mathbf g_{I,r,n}. \]

The covariance matrix of \((\operatorname{ADI}_{r,n},\operatorname{AII}_{r,n}, \operatorname{ATI}_{r,n})^{\top}\) consequently has rank at most two. Individual standard errors are well defined, but a joint Wald test should use two nonredundant impacts, such as \((\operatorname{ADI}_{r,n},\operatorname{AII}_{r,n})^{\top}\), or a full-row-rank transformation of the three measures.

4.9.1 Estimated Covariance Matrix

Let \[ \widehat{\mathbf V}_{\eta_r,n} = \widehat{\operatorname{Var}} \left( \widehat{\boldsymbol\eta}_{r,n} \right) \] be the estimated covariance matrix on the estimator scale, and define \[ \widehat{\mathbf G}_{r,n} = \mathbf G_{r,n} \left( \widehat{\boldsymbol\eta}_{r,n} \right). \]

Then \[ \widehat{\operatorname{Var}} \left[ \mathbf m_{r,n} \left( \widehat{\boldsymbol\eta}_{r,n} \right) \right] = \widehat{\mathbf G}_{r,n} \widehat{\mathbf V}_{\eta_r,n} \widehat{\mathbf G}_{r,n}^{\top}. \tag{4.165}\]

The diagonal elements provide the estimated variances of the three reported impact measures.

For example, an asymptotic confidence interval for the average indirect impact is \[ \widehat{\operatorname{AII}}_{r,n} \pm z_{1-\alpha/2} \sqrt{ \left[ \widehat{\operatorname{Var}} \left( \widehat{\mathbf m}_{r,n} \right) \right]_{22} }. \tag{4.166}\]

4.9.2 Extracting the Relevant Parameter Covariance

Let \(\mathbf L_r\) select the entries \((\beta_r,\gamma_r,\rho)^{\top}\) from the complete SDM parameter vector: \[ \widehat{\boldsymbol\eta}_{r,n} = \mathbf L_r \widehat{\boldsymbol\theta}_n. \]

Then \[ \widehat{\mathbf V}_{\eta_r,n} = \mathbf L_r \widehat{\operatorname{Var}} \left( \widehat{\boldsymbol\theta}_n \right) \mathbf L_r^{\top}. \tag{4.167}\]

The covariances between \(\widehat\rho_n\) and the relevant slope estimates must be retained. Treating \(\widehat\rho_n\) as fixed generally understates impact uncertainty.

The analysis conditions on the specified weights matrix. If \(\mathbf W_n\) is itself estimated, the formulas above do not account for that additional source of uncertainty.

4.9.3 Simulation from the Estimated Asymptotic Distribution

An alternative numerical method draws \[ \boldsymbol\eta_{r,n}^{(s)} \sim \mathcal N \left( \widehat{\boldsymbol\eta}_{r,n}, \widehat{\mathbf V}_{\eta_r,n} \right), \qquad s=1,\ldots,S, \] and computes \[ \mathbf m_{r,n}^{(s)} = \mathbf m_{r,n} \left( \boldsymbol\eta_{r,n}^{(s)} \right). \]

The empirical standard deviation of each simulated impact estimates its standard error under the fitted asymptotic distribution.

This simulation does not replace the asymptotic theory. It propagates the same estimated parameter uncertainty through the nonlinear impact function. The delta method uses a local linear approximation; asymptotic simulation evaluates the complete nonlinear function at every draw.

Implementation, treatment of inadmissible draws, and comparison of analytic and simulation-based standard errors belong to the separate R chapter.

Admissible draws of the spatial parameter

A Gaussian draw for \(\rho\) can fall outside the interval on which \(\mathbf A_n(\rho)\) is nonsingular. A computational implementation must use a parameter transformation or a clearly documented rule for invalid draws.

4.9.4 Impacts under Mixed Convergence Rates

The preceding delta method inherits the regular root-\(n\) normalization. It cannot be applied mechanically when different parameter components have different rates.

For example, if \(\widehat\rho_n\) has rate \[ \sqrt{\frac{n}{h_n}} \] while some coefficient estimators have rate \(\sqrt n\), the impact expansion must use a rate matrix rather than a common scalar normalization. Because the impact matrix is nonlinear in \(\rho\), the slower spatial-parameter component may dominate the limiting distribution of the impacts.

A mixed-rate impact theorem therefore requires a separate analysis and is not covered by Theorem 4.4. This is consistent with the irregular regimes discussed earlier in the chapter.

4.10 Marginal Effects in the Spatial Error Model

Chapter 2 already established the relevant special case. For the conventional SEM, \[ \mathbb E \left( \mathbf y_n \mid \mathbf X_n \right) = \mathbf X_n \boldsymbol\beta, \] so \[ \frac{ \partial \mathbb E \left( \mathbf y_n \mid \mathbf X_n \right) }{ \partial \mathbf x_{r,n}^{\top} } = \beta_r\mathbf I_n. \]

Consequently, \[ \operatorname{ADI}_{r,n}^{\mathrm{SEM}} = \operatorname{ATI}_{r,n}^{\mathrm{SEM}} = \beta_r, \qquad \operatorname{AII}_{r,n}^{\mathrm{SEM}} = 0. \tag{4.168}\]

The standard error of the SEM marginal effect is therefore the standard error of \(\widehat\beta_r\). Spatial correlation in the disturbances affects efficient estimation and inference, but it does not create a structural spillover in the conditional mean.

4.11 Main Lessons for Impact Inference

The impact formulas in this section are the same objects derived in Chapter 2. The new contribution is their sampling theory.

The main conclusions are:

the regular SLM QMLE theory extends to the SDM only after the augmented regressor matrix satisfies the corresponding boundedness, rank, and identification conditions;
the exact impact matrix is \[ \mathbf A_n(\rho)^{-1} \left( \beta_r\mathbf I_n+\gamma_r\mathbf W_n \right); \]
the uncertainty of an impact must include the covariance between the spatial-parameter estimate and the relevant slope estimates;
the delta method is available because the required gradients have closed forms;
asymptotic simulation is a complementary numerical method, not a substitute for the underlying asymptotic distribution;
ADI, AII, and ATI are exactly linearly dependent, so their joint three-dimensional covariance matrix is singular;
the conventional SEM has no structural indirect impact in its conditional mean.

Under row standardization, the average total SDM impact remains the result already proved in Chapter 2: \[ \operatorname{ATI}_{r,n} = \frac{ \beta_r+\gamma_r }{ 1-\rho }. \]

The separate R chapter will implement these formulas and compare delta-method and simulation-based standard errors.

4.12 Chapter Summary

This chapter developed maximum likelihood and quasi-maximum likelihood estimation for linear spatial models. The main arguments can be organized around five ideas.

First, the endogenous spatial lag \(\mathbf W_n\mathbf y_n\) prevents ordinary least squares from behaving as it does in the standard linear model. In the pure SLM, \[ \widehat\rho_{\mathrm{OLS},n} - \rho_0 = \frac{ \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n }{ \mathbf y_n^{\top} \mathbf W_n^{\top} \mathbf W_n \mathbf y_n }, \] so consistency depends on the relative limits of two quadratic forms.

Second, the likelihood of a spatial model contains a Jacobian term. For the SLM, \[ \mathbf A_n(\rho) = \mathbf I_n-\rho\mathbf W_n, \] and the log-likelihood contains \[ \log \left| \det \mathbf A_n(\rho) \right|. \]

This term accounts for the transformation from the innovations to the observed outcome. It is not an optional correction.

Third, concentrating the likelihood reduces the optimization problem. For each candidate value of the spatial parameter, the coefficient vector and the innovation variance have closed-form conditional estimators. The remaining optimization is one-dimensional in the conventional SLM and SEM.

Fourth, the asymptotic behavior of the QMLE depends on both identification and the structure of the spatial weights sequence. Under regular identification, \[ \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n \longrightarrow q_d>0, \] where \[ \mathbf d_n = \mathbf G_{n0} \mathbf X_n \boldsymbol\beta_0. \]

The QMLE has the usual root-\(n\) rate even when \(h_n\to\infty\). When this mean-identification condition fails, information about the spatial parameter may come primarily from the covariance structure, and slower rates can arise (Lee 2004).

Fifth, structural coefficients and spatial impacts are different objects. For the SDM, the impact matrix for regressor \(r\) is \[ \mathbf S_{r,n} = \mathbf A_n(\rho)^{-1} \left( \beta_r\mathbf I_n + \gamma_r\mathbf W_n \right). \]

The uncertainty of average direct, indirect, and total impacts must include the sampling uncertainty of the estimated spatial parameter and its covariance with the relevant coefficient estimates.

The theoretical chapter has deliberately separated these arguments from software implementation. Numerical maximization, model estimation in R, diagnostic checks, and computational impact calculations are developed in the separate practical chapter.

4.13 Theoretical Exercises

Exercise 4.1 (Derivatives of the concentrated SLM likelihood) Let \[ \ell_n^c(\rho) = -\frac{n}{2} \left[ \log(2\pi)+1 \right] - \frac{n}{2} \log \widehat\sigma_n^2(\rho) + \log \left| \det \mathbf A_n(\rho) \right|, \] where \[ \widehat\sigma_n^2(\rho) = \frac{1}{n} \mathbf y_n^{\top} \mathbf A_n(\rho)^{\top} \mathbf M_n \mathbf A_n(\rho) \mathbf y_n. \]

Derive the first and second derivatives of \(\ell_n^c(\rho)\) with respect to \(\rho\). Explain why the first-order condition does not generally have a closed-form solution.

Exercise 4.2 (Profile estimator of the innovation variance) For a fixed value of \(\rho\), show that \[ \widehat\sigma_n^2(\rho) = \frac{1}{n} \mathbf z_n(\rho)^{\top} \mathbf M_n \mathbf z_n(\rho), \] where \[ \mathbf z_n(\rho) = \mathbf A_n(\rho) \mathbf y_n \] and \[ \mathbf M_n = \mathbf I_n - \mathbf X_n \left( \mathbf X_n^{\top} \mathbf X_n \right)^{-1} \mathbf X_n^{\top}. \]

Interpret this expression as the residual variance from a transformed regression.

Exercise 4.3 (Derivatives of the Jacobian term) Let \[ \mathbf G_n(\rho) = \mathbf W_n \mathbf A_n(\rho)^{-1}. \]

Show that \[ \frac{\partial}{\partial\rho} \log \left| \det \mathbf A_n(\rho) \right| = -\operatorname{tr} \mathbf G_n(\rho) \] and \[ \frac{\partial^2}{\partial\rho^2} \log \left| \det \mathbf A_n(\rho) \right| = -\operatorname{tr} \left[ \mathbf G_n(\rho)^2 \right]. \]

State the conditions needed for these derivatives to exist.

Exercise 4.4 (OLS in a pure spatial autoregression) Consider \[ \mathbf y_n = \rho_0 \mathbf W_n \mathbf y_n + \boldsymbol\varepsilon_n. \]

Derive the OLS estimator of \(\rho_0\) and show that \[ \widehat\rho_{\mathrm{OLS},n} - \rho_0 = \frac{ \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n }{ \mathbf y_n^{\top} \mathbf W_n^{\top} \mathbf W_n \mathbf y_n }. \]

Using the quadratic-form law of large numbers from Chapter 3, provide sufficient conditions under which the probability limit differs from zero. Discuss why \[ \operatorname{tr} \mathbf G_{n0}=0 \] is not equivalent, in general, to \(\rho_0=0\).

Exercise 4.5 (OLS under a spatial error model) Consider \[ \begin{aligned} \mathbf y_n &= \mathbf X_n \boldsymbol\beta_0 + \mathbf u_n, \\ \mathbf u_n &= \lambda_0 \mathbf W_n \mathbf u_n + \boldsymbol\varepsilon_n. \end{aligned} \]

Assume \[ \mathbb E \left( \boldsymbol\varepsilon_n \mid \mathbf X_n \right) = \mathbf0 \] and \[ \operatorname{Var} \left( \boldsymbol\varepsilon_n \mid \mathbf X_n \right) = \sigma_0^2\mathbf I_n. \]

Show that the OLS estimator of \(\boldsymbol\beta_0\) is conditionally unbiased. Derive its conditional covariance matrix and explain why the conventional homoskedastic OLS covariance estimator is not valid when \(\lambda_0\neq0\).

Exercise 4.6 (Concentrated likelihood of the SEM) For the SEM in Exercise 4.5, define \[ \mathbf B_n(\lambda) = \mathbf I_n-\lambda\mathbf W_n. \]

Derive:

the Gaussian log-likelihood;
the conditional estimator \(\widehat{\boldsymbol\beta}_n(\lambda)\);
the conditional estimator \(\widehat\sigma_n^2(\lambda)\);
the concentrated log-likelihood for \(\lambda\).

Explain why both \(\mathbf y_n\) and \(\mathbf X_n\) must be transformed by \(\mathbf B_n(\lambda)\).

Exercise 4.7 (Regular identification and the information matrix) Let \[ \mathbf d_n = \mathbf G_{n0} \mathbf X_n \boldsymbol\beta_0. \]

Using a partitioned-matrix argument, show that the contribution of \((\boldsymbol\beta,\rho)\) to the limiting information matrix is nonsingular if \[ \frac{1}{n} \mathbf X_n^{\top} \mathbf X_n \longrightarrow \mathbf Q_X>0 \] and \[ \frac{1}{n} \mathbf d_n^{\top} \mathbf M_n \mathbf d_n \longrightarrow q_d>0. \]

Interpret the second condition as a residualized rank condition.

Exercise 4.8 (Linear-quadratic representation of the score) At the true parameter, show that the score for \(\rho\) can be written as \[ s_{\rho,n} \left( \boldsymbol\theta_0 \right) = \frac{1}{\sigma_0^2} \left[ \mathbf d_n^{\top} \boldsymbol\varepsilon_n + \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n - \sigma_0^2 \operatorname{tr} \mathbf G_{n0} \right]. \]

Identify the linear component and the centered quadratic component. Use the results from Chapter 3 to determine the order of each component when:

\(h_n\) is bounded;
\(h_n\to\infty\) and regular identification holds;
\(h_n\to\infty\) and regular mean identification fails.

Exercise 4.9 (Gaussian QMLE under non-Gaussian innovations) Suppose that the first four moments of the innovations exist, but the innovations are not Gaussian. Explain why the inverse expected Hessian alone does not generally estimate the asymptotic covariance matrix of the QMLE.

Starting from \[ \mathbf J_{\theta} = \boldsymbol\Sigma_{\theta} + \boldsymbol\Omega_{\theta}, \] derive the sandwich covariance matrix \[ \boldsymbol\Sigma_{\theta}^{-1} \mathbf J_{\theta} \boldsymbol\Sigma_{\theta}^{-1}. \]

Under what restrictions on the third and fourth moments does the information equality hold?

Exercise 4.10 (The SDM as an augmented SLM) Consider \[ \mathbf y_n = \rho_0 \mathbf W_n \mathbf y_n + \alpha_0 \boldsymbol\iota_n + \mathbf X_n \boldsymbol\beta_0 + \mathbf W_n \mathbf X_n \boldsymbol\gamma_0 + \boldsymbol\varepsilon_n. \]

Construct the augmented design matrix and coefficient vector that write this model in the form of an SLM. State the boundedness, rank, and identification conditions needed to apply the regular QMLE theory.

Explain why \(\mathbf W_n\boldsymbol\iota_n\) must not be included as an additional regressor when \(\mathbf W_n\) is row standardized.

Exercise 4.11 (Delta-method gradients for SDM impacts) For regressor \(r\), let \[ \mathbf S_{r,n} = \mathbf A_n(\rho)^{-1} \left( \beta_r\mathbf I_n + \gamma_r\mathbf W_n \right). \]

Derive the gradients of the average direct and average total impacts with respect to \[ \boldsymbol\eta_r = \begin{pmatrix} \beta_r\\ \gamma_r\\ \rho \end{pmatrix}. \]

Use the result to obtain the gradient of the average indirect impact.

Exercise 4.12 (Linear dependence among average impacts) Prove that \[ \operatorname{ATI}_{r,n} = \operatorname{ADI}_{r,n} + \operatorname{AII}_{r,n}. \]

Show that the delta-method covariance matrix of \[ \begin{pmatrix} \operatorname{ADI}_{r,n}\\ \operatorname{AII}_{r,n}\\ \operatorname{ATI}_{r,n} \end{pmatrix} \] is singular. Propose a nonredundant vector of impacts for a joint Wald test.

Exercise 4.13 (Total impacts under row standardization) Assume \[ \mathbf W_n \boldsymbol\iota_n = \boldsymbol\iota_n. \]

Show that the average total impact in the SDM is \[ \operatorname{ATI}_{r,n} = \frac{ \beta_r+\gamma_r }{ 1-\rho }. \]

Derive its gradient and specialize the result to the SLM.

Exercise 4.14 (A model with two endogenous spatial lags) Consider \[ \mathbf y_n = \rho_1 \mathbf W_{1n} \mathbf y_n + \rho_2 \mathbf W_{2n} \mathbf y_n + \mathbf X_n \boldsymbol\beta + \boldsymbol\varepsilon_n. \]

Define \[ \mathbf A_n \left( \rho_1,\rho_2 \right) = \mathbf I_n - \rho_1\mathbf W_{1n} - \rho_2\mathbf W_{2n}. \]

Derive the Gaussian log-likelihood and the first-order condition for \(\boldsymbol\beta\). Give a sufficient condition for \(\mathbf A_n(\rho_1,\rho_2)\) to be nonsingular. Explain why row normalization of both weights matrices is not, by itself, sufficient to guarantee nonsingularity for every pair \((\rho_1,\rho_2)\).

Exercise 4.15 (Likelihood for a heteroskedastic SAC model) Consider \[ \begin{aligned} \mathbf y_n &= \rho \mathbf W_{1n} \mathbf y_n + \mathbf X_n \boldsymbol\beta + \mathbf u_n, \\ \mathbf u_n &= \lambda \mathbf W_{2n} \mathbf u_n + \boldsymbol\varepsilon_n, \end{aligned} \] with \[ \operatorname{Var} \left( \boldsymbol\varepsilon_n \mid \mathbf Z_n \right) = \boldsymbol\Omega_n(\boldsymbol\alpha), \] where \(\boldsymbol\Omega_n(\boldsymbol\alpha)\) is diagonal and positive definite.

Derive the Gaussian log-likelihood, including both Jacobian terms. Then derive the score with respect to \(\boldsymbol\beta\). State the additional derivative needed to obtain the score for \(\boldsymbol\alpha\).

Technical Appendix: Detailed SLM Derivations

The principal consistency and asymptotic-normality proofs were developed in the main text because they reveal the structure of spatial maximum likelihood. This appendix records the longer algebraic derivations underlying three objects used there:

the derivatives of the concentrated likelihood;
the expected Hessian;
the covariance matrix of the score.

These derivations refine the corresponding calculations in the original notes and use the linear- and quadratic-form identities established in Chapter 3.

4.13.1 Derivatives of the Concentrated Likelihood

For the SLM, define \[ \mathbf A_n(\rho) = \mathbf I_n-\rho\mathbf W_n, \]

\[ \mathbf M_n = \mathbf I_n - \mathbf X_n \left( \mathbf X_n^{\top} \mathbf X_n \right)^{-1} \mathbf X_n^{\top}, \] and \[ \mathbf u_n(\rho) = \mathbf M_n \mathbf A_n(\rho) \mathbf y_n. \]

Also define the vector \[ \mathbf v_n = \mathbf M_n \mathbf W_n \mathbf y_n. \]

The vector \(\mathbf v_n\) does not depend on \(\rho\). Since \[ \frac{\partial\mathbf A_n(\rho)}{\partial\rho} = -\mathbf W_n, \] we have \[ \frac{\partial\mathbf u_n(\rho)}{\partial\rho} = -\mathbf v_n. \tag{4.169}\]

The concentrated innovation variance is \[ \widehat\sigma_n^2(\rho) = \frac{1}{n} \mathbf u_n(\rho)^{\top} \mathbf u_n(\rho). \tag{4.170}\]

Differentiating Equation 4.170 and using Equation 4.169 gives \[ \begin{aligned} \frac{\partial\widehat\sigma_n^2(\rho)}{\partial\rho} &= \frac{1}{n} \left[ \frac{\partial\mathbf u_n(\rho)^{\top}}{\partial\rho} \mathbf u_n(\rho) + \mathbf u_n(\rho)^{\top} \frac{\partial\mathbf u_n(\rho)}{\partial\rho} \right] \\ &= -\frac{2}{n} \mathbf v_n^{\top} \mathbf u_n(\rho). \end{aligned} \tag{4.171}\]

A second differentiation yields \[ \begin{aligned} \frac{\partial^2\widehat\sigma_n^2(\rho)}{\partial\rho^2} &= -\frac{2}{n} \mathbf v_n^{\top} \frac{\partial\mathbf u_n(\rho)}{\partial\rho} \\ &= \frac{2}{n} \mathbf v_n^{\top} \mathbf v_n. \end{aligned} \tag{4.172}\]

The concentrated log-likelihood is \[ \ell_n^c(\rho) = c_n - \frac{n}{2} \log \widehat\sigma_n^2(\rho) + \log \left| \det \mathbf A_n(\rho) \right|, \tag{4.173}\] where \[ c_n = -\frac{n}{2} \left[ \log(2\pi)+1 \right]. \]

Let \[ \mathbf G_n(\rho) = \mathbf W_n \mathbf A_n(\rho)^{-1}. \]

The matrix derivative identity gives \[ \frac{\partial}{\partial\rho} \log \left| \det \mathbf A_n(\rho) \right| = \operatorname{tr} \left[ \mathbf A_n(\rho)^{-1} \frac{\partial\mathbf A_n(\rho)}{\partial\rho} \right] = -\operatorname{tr} \mathbf G_n(\rho). \tag{4.174}\]

Because \[ \frac{\partial\mathbf A_n(\rho)^{-1}}{\partial\rho} = \mathbf A_n(\rho)^{-1} \mathbf W_n \mathbf A_n(\rho)^{-1}, \] we obtain \[ \frac{\partial\mathbf G_n(\rho)}{\partial\rho} = \mathbf G_n(\rho)^2. \tag{4.175}\]

Therefore, \[ \frac{\partial^2}{\partial\rho^2} \log \left| \det \mathbf A_n(\rho) \right| = -\operatorname{tr} \left[ \mathbf G_n(\rho)^2 \right]. \tag{4.176}\]

Substituting Equation 4.171 and Equation 4.174 into the derivative of Equation 4.173 gives \[ \frac{\partial\ell_n^c(\rho)}{\partial\rho} = \frac{ \mathbf v_n^{\top} \mathbf u_n(\rho) }{ \widehat\sigma_n^2(\rho) } - \operatorname{tr} \mathbf G_n(\rho). \tag{4.177}\]

Differentiate the first term in Equation 4.177: \[ \begin{aligned} \frac{\partial}{\partial\rho} \left[ \frac{ \mathbf v_n^{\top} \mathbf u_n(\rho) }{ \widehat\sigma_n^2(\rho) } \right] &= \frac{ -\mathbf v_n^{\top}\mathbf v_n }{ \widehat\sigma_n^2(\rho) } - \frac{ \mathbf v_n^{\top}\mathbf u_n(\rho) }{ \widehat\sigma_n^4(\rho) } \frac{ \partial\widehat\sigma_n^2(\rho) }{ \partial\rho } \\ &= -\frac{ \mathbf v_n^{\top}\mathbf v_n }{ \widehat\sigma_n^2(\rho) } + \frac{2}{n} \frac{ \left[ \mathbf v_n^{\top}\mathbf u_n(\rho) \right]^2 }{ \widehat\sigma_n^4(\rho) }. \end{aligned} \]

Combining this result with Equation 4.176 yields \[ \boxed{ \begin{aligned} \frac{\partial^2\ell_n^c(\rho)}{\partial\rho^2} &= -\operatorname{tr} \left[ \mathbf G_n(\rho)^2 \right] - \frac{ \mathbf v_n^{\top}\mathbf v_n }{ \widehat\sigma_n^2(\rho) } \\ &\quad+ \frac{2}{n} \frac{ \left[ \mathbf v_n^{\top}\mathbf u_n(\rho) \right]^2 }{ \widehat\sigma_n^4(\rho) }. \end{aligned} } \tag{4.178}\]

The score equation in Equation 4.177 combines a rational function of \(\rho\) with the trace of a matrix inverse. It therefore has no general closed-form solution.

4.13.2 Expected Hessian of the SLM Log-Likelihood

Use the parameter ordering \[ \boldsymbol\theta = \begin{pmatrix} \boldsymbol\beta\\ \rho\\ \sigma^2 \end{pmatrix}. \]

At the true parameter, define \[ \mathbf A_{n0} = \mathbf A_n(\rho_0), \]

\[ \mathbf G_{n0} = \mathbf W_n \mathbf A_{n0}^{-1}, \] and \[ \mathbf d_n = \mathbf G_{n0} \mathbf X_n \boldsymbol\beta_0. \]

The reduced form implies \[ \mathbf y_n = \mathbf A_{n0}^{-1} \mathbf X_n \boldsymbol\beta_0 + \mathbf A_{n0}^{-1} \boldsymbol\varepsilon_n. \]

Premultiplication by \(\mathbf W_n\) gives \[ \mathbf W_n \mathbf y_n = \mathbf d_n + \mathbf G_{n0} \boldsymbol\varepsilon_n. \tag{4.179}\]

The identities \[ \mathbb E \left( \boldsymbol\varepsilon_n \right) = \mathbf0 \] and \[ \mathbb E \left( \boldsymbol\varepsilon_n \boldsymbol\varepsilon_n^{\top} \right) = \sigma_0^2 \mathbf I_n \] imply \[ \mathbb E \left( \mathbf W_n\mathbf y_n \right) = \mathbf d_n. \tag{4.180}\]

They also imply \[ \begin{aligned} \mathbb E \left( \mathbf y_n^{\top} \mathbf W_n^{\top} \mathbf W_n \mathbf y_n \right) &= \mathbb E \left[ \left( \mathbf d_n + \mathbf G_{n0}\boldsymbol\varepsilon_n \right)^{\top} \left( \mathbf d_n + \mathbf G_{n0}\boldsymbol\varepsilon_n \right) \right] \\ &= \mathbf d_n^{\top}\mathbf d_n + \sigma_0^2 \operatorname{tr} \left( \mathbf G_{n0}^{\top} \mathbf G_{n0} \right). \end{aligned} \tag{4.181}\]

Finally, \[ \begin{aligned} \mathbb E \left( \boldsymbol\varepsilon_n^{\top} \mathbf W_n \mathbf y_n \right) &= \mathbb E \left[ \boldsymbol\varepsilon_n^{\top} \mathbf d_n + \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n \right] \\ &= \sigma_0^2 \operatorname{tr} \mathbf G_{n0}. \end{aligned} \tag{4.182}\]

The Hessian blocks evaluated at the truth are \[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\boldsymbol\beta \partial\boldsymbol\beta^{\top} } = -\frac{1}{\sigma_0^2} \mathbf X_n^{\top} \mathbf X_n, \tag{4.183}\]

\[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\boldsymbol\beta \partial\rho } = -\frac{1}{\sigma_0^2} \mathbf X_n^{\top} \mathbf W_n \mathbf y_n, \tag{4.184}\]

\[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\boldsymbol\beta \partial\sigma^2 } = -\frac{1}{\sigma_0^4} \mathbf X_n^{\top} \boldsymbol\varepsilon_n, \tag{4.185}\]

\[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\rho^2 } = -\operatorname{tr} \left( \mathbf G_{n0}^2 \right) - \frac{1}{\sigma_0^2} \mathbf y_n^{\top} \mathbf W_n^{\top} \mathbf W_n \mathbf y_n, \tag{4.186}\]

\[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\rho \partial\sigma^2 } = -\frac{1}{\sigma_0^4} \boldsymbol\varepsilon_n^{\top} \mathbf W_n \mathbf y_n, \tag{4.187}\] and \[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial(\sigma^2)^2 } = \frac{n}{2\sigma_0^4} - \frac{1}{\sigma_0^6} \boldsymbol\varepsilon_n^{\top} \boldsymbol\varepsilon_n. \tag{4.188}\]

Taking expectations and using Equation 4.180–Equation 4.182 gives \[ -\mathbb E \left[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\boldsymbol\beta \partial\boldsymbol\beta^{\top} } \right] = \frac{1}{\sigma_0^2} \mathbf X_n^{\top} \mathbf X_n, \]

\[ -\mathbb E \left[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\boldsymbol\beta \partial\rho } \right] = \frac{1}{\sigma_0^2} \mathbf X_n^{\top} \mathbf d_n, \]

\[ -\mathbb E \left[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\boldsymbol\beta \partial\sigma^2 } \right] = \mathbf0, \]

\[ \begin{aligned} -\mathbb E \left[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\rho^2 } \right] &= \frac{1}{\sigma_0^2} \mathbf d_n^{\top} \mathbf d_n \\ &\quad+ \operatorname{tr} \left( \mathbf G_{n0}^{\top} \mathbf G_{n0} \right) + \operatorname{tr} \left( \mathbf G_{n0}^2 \right), \end{aligned} \]

\[ -\mathbb E \left[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial\rho \partial\sigma^2 } \right] = \frac{1}{\sigma_0^2} \operatorname{tr} \mathbf G_{n0}, \] and \[ -\mathbb E \left[ \frac{ \partial^2 \ell_n(\boldsymbol\theta_0) }{ \partial(\sigma^2)^2 } \right] = \frac{n}{2\sigma_0^4}. \]

Therefore, the negative expected average Hessian is \[ \boxed{ \boldsymbol\Sigma_{\theta,n} = \begin{pmatrix} \dfrac{ \mathbf X_n^{\top}\mathbf X_n }{ n\sigma_0^2 } & \dfrac{ \mathbf X_n^{\top}\mathbf d_n }{ n\sigma_0^2 } & \mathbf0 \\[12pt] \dfrac{ \mathbf d_n^{\top}\mathbf X_n }{ n\sigma_0^2 } & \dfrac{ \mathbf d_n^{\top}\mathbf d_n }{ n\sigma_0^2 } + \dfrac{ \operatorname{tr} \left( \mathbf G_{n0}^{\top}\mathbf G_{n0} \right) + \operatorname{tr} \left( \mathbf G_{n0}^2 \right) }{ n } & \dfrac{ \operatorname{tr} \mathbf G_{n0} }{ n\sigma_0^2 } \\[14pt] \mathbf0^{\top} & \dfrac{ \operatorname{tr} \mathbf G_{n0} }{ n\sigma_0^2 } & \dfrac{1}{2\sigma_0^4} \end{pmatrix}. } \tag{4.189}\]

The derivation corrects two common mistakes. First, \[ \mathbb E \left( \mathbf y_n\mathbf y_n^{\top} \right) \] is a second-moment matrix, not the covariance matrix alone. Second, the \(\rho\rho\) entry contains \[ \operatorname{tr} \left( \mathbf G_{n0}^{\top}\mathbf G_{n0} \right) + \operatorname{tr} \left( \mathbf G_{n0}^2 \right), \] not the square of \(\operatorname{tr}\mathbf G_{n0}\).

4.13.3 Covariance Matrix of the Score

At the true parameter, the score blocks are \[ \mathbf s_{\beta,n} = \frac{1}{\sigma_0^2} \mathbf X_n^{\top} \boldsymbol\varepsilon_n, \tag{4.190}\]

\[ s_{\rho,n} = \frac{1}{\sigma_0^2} \left[ \mathbf d_n^{\top} \boldsymbol\varepsilon_n + \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n - \sigma_0^2 \operatorname{tr} \mathbf G_{n0} \right], \tag{4.191}\] and \[ s_{\sigma^2,n} = \frac{1}{2\sigma_0^4} \left[ \boldsymbol\varepsilon_n^{\top} \boldsymbol\varepsilon_n - n\sigma_0^2 \right]. \tag{4.192}\]

Let \[ \mu_3 = \mathbb E \left( \varepsilon_{ni}^3 \right), \qquad \mu_4 = \mathbb E \left( \varepsilon_{ni}^4 \right), \] and define \[ \mathbf g_n = \operatorname{diag} \left( \mathbf G_{n0} \right). \]

Here \(\operatorname{diag}(\mathbf G_{n0})\) denotes the vector formed by the diagonal elements of \(\mathbf G_{n0}\).

The calculations use the identities from Chapter 3: \[ \operatorname{Cov} \left( \mathbf a_n^{\top}\boldsymbol\varepsilon_n, \mathbf b_n^{\top}\boldsymbol\varepsilon_n \right) = \sigma_0^2 \mathbf a_n^{\top} \mathbf b_n, \tag{4.193}\]

\[ \operatorname{Cov} \left( \mathbf a_n^{\top}\boldsymbol\varepsilon_n, \boldsymbol\varepsilon_n^{\top} \mathbf H_n \boldsymbol\varepsilon_n \right) = \mu_3 \mathbf a_n^{\top} \operatorname{diag} \left( \mathbf H_n \right), \tag{4.194}\] and \[ \begin{aligned} \operatorname{Cov} \left( \boldsymbol\varepsilon_n^{\top} \mathbf H_n \boldsymbol\varepsilon_n, \boldsymbol\varepsilon_n^{\top} \mathbf K_n \boldsymbol\varepsilon_n \right) &= \left( \mu_4-3\sigma_0^4 \right) \operatorname{diag} \left( \mathbf H_n \right)^{\top} \operatorname{diag} \left( \mathbf K_n \right) \\ &\quad+ \sigma_0^4 \left[ \operatorname{tr} \left( \mathbf H_n \mathbf K_n^{\top} \right) + \operatorname{tr} \left( \mathbf H_n \mathbf K_n \right) \right]. \end{aligned} \tag{4.195}\]

4.13.3.1 The \(\boldsymbol\beta\boldsymbol\beta\) block

From Equation 4.190, \[ \begin{aligned} \operatorname{Var} \left( \frac{1}{\sqrt n} \mathbf s_{\beta,n} \right) &= \frac{1}{n\sigma_0^4} \mathbf X_n^{\top} \mathbb E \left( \boldsymbol\varepsilon_n \boldsymbol\varepsilon_n^{\top} \right) \mathbf X_n \\ &= \frac{1}{n\sigma_0^2} \mathbf X_n^{\top} \mathbf X_n. \end{aligned} \tag{4.196}\]

4.13.3.2 The \(\boldsymbol\beta\rho\) block

Using Equation 4.193 and Equation 4.194, \[ \begin{aligned} \operatorname{Cov} \left( \frac{1}{\sqrt n} \mathbf s_{\beta,n}, \frac{1}{\sqrt n} s_{\rho,n} \right) &= \frac{1}{n\sigma_0^4} \left[ \sigma_0^2 \mathbf X_n^{\top} \mathbf d_n + \mu_3 \mathbf X_n^{\top} \mathbf g_n \right] \\ &= \frac{ \mathbf X_n^{\top} \mathbf d_n }{ n\sigma_0^2 } + \frac{ \mu_3 \mathbf X_n^{\top} \mathbf g_n }{ n\sigma_0^4 }. \end{aligned} \tag{4.197}\]

4.13.3.3 The \(\boldsymbol\beta\sigma^2\) block

Because \[ \operatorname{diag} \left( \mathbf I_n \right) = \boldsymbol\iota_n, \]

Equation 4.194 gives \[ \operatorname{Cov} \left( \frac{1}{\sqrt n} \mathbf s_{\beta,n}, \frac{1}{\sqrt n} s_{\sigma^2,n} \right) = \frac{ \mu_3 }{ 2n\sigma_0^6 } \mathbf X_n^{\top} \boldsymbol\iota_n. \tag{4.198}\]

4.13.3.4 The \(\rho\rho\) block

Write \[ L_{\rho,n} = \mathbf d_n^{\top} \boldsymbol\varepsilon_n \]

and

\[ Q_{\rho,n} = \boldsymbol\varepsilon_n^{\top} \mathbf G_{n0} \boldsymbol\varepsilon_n - \sigma_0^2 \operatorname{tr} \mathbf G_{n0}. \]

Then \[ s_{\rho,n} = \frac{1}{\sigma_0^2} \left( L_{\rho,n}+Q_{\rho,n} \right). \]

The linear variance is \[ \operatorname{Var} \left( L_{\rho,n} \right) = \sigma_0^2 \mathbf d_n^{\top} \mathbf d_n. \]

The linear-quadratic covariance is \[ \operatorname{Cov} \left( L_{\rho,n}, Q_{\rho,n} \right) = \mu_3 \mathbf d_n^{\top} \mathbf g_n. \]

Finally, Equation 4.195 with \(\mathbf H_n=\mathbf K_n=\mathbf G_{n0}\) gives \[ \begin{aligned} \operatorname{Var} \left( Q_{\rho,n} \right) &= \left( \mu_4-3\sigma_0^4 \right) \mathbf g_n^{\top}\mathbf g_n \\ &\quad+ \sigma_0^4 \left[ \operatorname{tr} \left( \mathbf G_{n0}^{\top}\mathbf G_{n0} \right) + \operatorname{tr} \left( \mathbf G_{n0}^2 \right) \right]. \end{aligned} \]

Therefore, \[ \begin{aligned} \operatorname{Var} \left( \frac{1}{\sqrt n} s_{\rho,n} \right) &= \frac{ \mathbf d_n^{\top}\mathbf d_n }{ n\sigma_0^2 } + \frac{ \operatorname{tr} \left( \mathbf G_{n0}^{\top}\mathbf G_{n0} \right) + \operatorname{tr} \left( \mathbf G_{n0}^2 \right) }{ n } \\ &\quad+ \frac{ 2\mu_3 \mathbf d_n^{\top} \mathbf g_n }{ n\sigma_0^4 } + \frac{ \left( \mu_4-3\sigma_0^4 \right) \mathbf g_n^{\top}\mathbf g_n }{ n\sigma_0^4 }. \end{aligned} \tag{4.199}\]

4.13.3.5 The \(\rho\sigma^2\) block

The covariance between the two linear-quadratic scores is \[ \begin{aligned} \operatorname{Cov} \left( \frac{1}{\sqrt n} s_{\rho,n}, \frac{1}{\sqrt n} s_{\sigma^2,n} \right) &= \frac{ \operatorname{tr} \mathbf G_{n0} }{ n\sigma_0^2 } \\ &\quad+ \frac{ \mu_3 \boldsymbol\iota_n^{\top}\mathbf d_n + \left( \mu_4-3\sigma_0^4 \right) \operatorname{tr} \mathbf G_{n0} }{ 2n\sigma_0^6 }. \end{aligned} \tag{4.200}\]

The first term in Equation 4.200 is the corresponding element of the expected information matrix. The second term is the non-Gaussian correction.

4.13.3.6 The \(\sigma^2\sigma^2\) block

Since \[ \operatorname{Var} \left( \varepsilon_{ni}^2 \right) = \mu_4-\sigma_0^4, \] independence gives \[ \operatorname{Var} \left( \boldsymbol\varepsilon_n^{\top} \boldsymbol\varepsilon_n \right) = n \left( \mu_4-\sigma_0^4 \right). \]

Consequently, \[ \begin{aligned} \operatorname{Var} \left( \frac{1}{\sqrt n} s_{\sigma^2,n} \right) &= \frac{ \mu_4-\sigma_0^4 }{ 4\sigma_0^8 } \\ &= \frac{1}{2\sigma_0^4} + \frac{ \mu_4-3\sigma_0^4 }{ 4\sigma_0^8 }. \end{aligned} \tag{4.201}\]

Collecting the blocks gives \[ \mathbf J_{\theta,n} = \operatorname{Var} \left[ \frac{1}{\sqrt n} \mathbf s_n \left( \boldsymbol\theta_0 \right) \right] = \boldsymbol\Sigma_{\theta,n} + \boldsymbol\Omega_{\theta,n}, \tag{4.202}\] where \(\boldsymbol\Sigma_{\theta,n}\) is given in Equation 4.189 and \[ \boxed{ \boldsymbol\Omega_{\theta,n} = \begin{pmatrix} \mathbf0 & \dfrac{ \mu_3 \mathbf X_n^{\top}\mathbf g_n }{ n\sigma_0^4 } & \dfrac{ \mu_3 \mathbf X_n^{\top}\boldsymbol\iota_n }{ 2n\sigma_0^6 } \\[14pt] \dfrac{ \mu_3 \mathbf g_n^{\top}\mathbf X_n }{ n\sigma_0^4 } & \dfrac{ 2\mu_3 \mathbf d_n^{\top}\mathbf g_n }{ n\sigma_0^4 } + \dfrac{ \left( \mu_4-3\sigma_0^4 \right) \mathbf g_n^{\top}\mathbf g_n }{ n\sigma_0^4 } & \dfrac{ \mu_3 \boldsymbol\iota_n^{\top}\mathbf d_n + \left( \mu_4-3\sigma_0^4 \right) \operatorname{tr} \mathbf G_{n0} }{ 2n\sigma_0^6 } \\[18pt] \dfrac{ \mu_3 \boldsymbol\iota_n^{\top}\mathbf X_n }{ 2n\sigma_0^6 } & \dfrac{ \mu_3 \boldsymbol\iota_n^{\top}\mathbf d_n + \left( \mu_4-3\sigma_0^4 \right) \operatorname{tr} \mathbf G_{n0} }{ 2n\sigma_0^6 } & \dfrac{ \mu_4-3\sigma_0^4 }{ 4\sigma_0^8 } \end{pmatrix}. } \tag{4.203}\]

When the innovations are Gaussian, \[ \mu_3=0 \] and \[ \mu_4=3\sigma_0^4. \]

Therefore, \[ \boldsymbol\Omega_{\theta,n} = \mathbf0 \] and \[ \mathbf J_{\theta,n} = \boldsymbol\Sigma_{\theta,n}. \]

This establishes the information equality for correctly specified Gaussian maximum likelihood.

4.13.4 How the Appendix Enters the Asymptotic Proof

The derivations above are algebraic. The probability limits used in the main proof come from Chapter 3.

The consistency proof uses:

uniform boundedness of the spatial multiplier;
orders for products of spatial matrices;
the law of large numbers for linear forms;
the law of large numbers for centered quadratic forms;
uniform control over the compact spatial-parameter set.

The asymptotic-normality proof uses:

the linear-quadratic representation in Equation 4.190–Equation 4.192;
the vector linear-quadratic central limit theorem;
convergence of the intermediate Hessian;
stochastic inversion;
Slutsky’s theorem.

The expansion is performed component by component. For each component \(\theta_j\) of the score equation, the scalar mean value theorem gives an intermediate point on the segment between \(\widehat{\boldsymbol\theta}_n\) and \(\boldsymbol\theta_0\). Stacking these componentwise equations produces the intermediate Hessian used in the main text. No vector-valued integral form of the mean value theorem is required.

The original appendix supplied the essential structure of these proofs. The main text reorganized that structure so that each stochastic order and each application of a Chapter 3 result is stated at the point where it is used.