Likelihood Estimation of Spatial Autoregressive Models

Why the Spatial Likelihood Has an Extra Term

In an ordinary linear regression, the Gaussian likelihood is built from the residuals and their variance. Spatial autoregressive models require one more ingredient: the response or disturbance is transformed by a matrix that depends on the spatial coefficient. That transformation changes the density’s scale, so its Jacobian contributes a log-determinant to the likelihood. Omitting it, or using its sign incorrectly, changes the parameter being optimized.

The objective here is likelihood estimation for two familiar specifications: the spatial-lag model, in which the outcome depends on neighboring outcomes, and the spatial-error model, in which the disturbances have spatial dependence. The likelihood calculations condition on the observed spatial weights matrix and the regressors.

Models and Assumptions

Let \(y\) be an \(n\times 1\) outcome vector, \(X\) an \(n\times k\) matrix of regressors with full column rank, and \(W\) a fixed \(n\times n\) spatial weights matrix. Write \(I\) for the \(n\times n\) identity matrix. The entries of \(W\) and its normalization are treated as given during estimation.

For the spatial-lag model define \(A_\rho=I-\rho W\). Rearranging gives \(A_\rho y=X\beta+\varepsilon\). For the spatial-error model define \(A_\lambda=I-\lambda W\); rearrangement gives \(A_\lambda(y-X\beta)=\varepsilon\). In both cases the innovation vector is assumed to have independent, equal-variance Gaussian components conditional on \(X,W\). The transformations must be invertible so that the density of the observed outcome is well defined.

Deriving the Gaussian Log-Likelihood

If a random vector \(z\) has density \(f_z\), and \(z=g(y)\) is a one-to-one differentiable transformation, then the density of \(y\) is \(f_y(y)=f_z(g(y))|\det(\partial g(y)/\partial y')|\). Taking logarithms adds the log absolute determinant of the transformation’s Jacobian. In the spatial-lag model, the transformation from \(y\) to innovations is \(\varepsilon=A_\rho y-X\beta\), whose derivative with respect to \(y\) is \(A_\rho\). In the spatial-error model it is \(\varepsilon=A_\lambda(y-X\beta)\), with derivative \(A_\lambda\). Thus each model contributes its own log-determinant term.

Here \(|A|\) denotes the absolute value of the determinant of a square matrix \(A\), so \(\log|A|\) is defined whenever \(A\) is nonsingular. The sign is important: the log-likelihood contains \(+\log|A|\), not its negative. If a transformation compresses volume, its determinant has magnitude below one and this term is negative; if it expands volume, the term is positive.

Concentrating Out Regression Coefficients and Variance

For any fixed admissible spatial coefficient, if the minimized transformed sum of squared residuals is positive, the likelihood can be maximized over the ordinary regression coefficients and innovation variance. If that sum is zero, the likelihood is unbounded as \(\sigma^2\to0\) (with \(\sigma^2>0\)), so there is no finite variance MLE or attained maximum. In the spatial-lag model, the transformed outcome is \(A_\rho y\) and the regression matrix remains \(X\). The conditional maximum-likelihood coefficient is the ordinary least-squares coefficient from regressing \(A_\rho y\) on \(X\). In the spatial-error model, the transformed residual is \(A_\lambda(y-X\beta)\); minimizing its squared length gives a transformed least-squares coefficient.

These expressions require the indicated matrices to be invertible; full column rank of \(X\) suffices for the spatial-lag expression, while the spatial-error expression also uses nonsingularity of \(A_\lambda\). When the relevant \(SSR\) is positive, the variance estimate at a fixed spatial coefficient is the maximum-likelihood estimate \(SSR/n\), not a degrees-of-freedom-adjusted unbiased estimate. Substituting it into the likelihood gives the concentrated log-likelihood.

The formula follows by differentiating the Gaussian log-likelihood with respect to \(\sigma^2\). If \(SSR\) is the relevant transformed sum of squared residuals, the variance-dependent part is \(-\frac{n}{2}\log(\sigma^2)-\frac{SSR}{2\sigma^2}\). Its derivative is \(-n/(2\sigma^2)+SSR/(2(\sigma^2)^2)\), which is zero at \(\widehat{\sigma}^2=SSR/n\). The second derivative there is negative, so this is a maximum. Substitution yields the displayed concentrated expression. Estimation then searches over the admissible spatial coefficient, recomputing both the determinant and the transformed residual sum of squares at each candidate.

Numerical Evaluation Over the Admissible Region

A spatial coefficient is excluded whenever \(I-\theta W\) is singular, where \(\theta\) denotes \(\rho\) or \(\lambda\). At such a value, the transformation has no inverse and the Gaussian density transformation used above fails. In applied work, the optimization region is typically a specified admissible interval containing zero and bounded away from determinant singularities. Additional stability restrictions may also be imposed. For example, a condition ensuring a convergent spatial multiplier is not identical to the minimal requirement that the likelihood transformation be nonsingular; the model’s stated parameter restrictions should determine the region used.

For a real eigenvalue \(\omega\) of \(W\), a singularity occurs when \(1-\theta\omega=0\). This gives a practical way to locate possible boundaries. The determinant can be evaluated directly using a numerically stable factorization, such as pivoted LU. Alternatively, when eigenvalues are available, \(\log|I-\theta W|\) can be computed from the factors \(1-\theta\omega_j\), combining complex conjugate pairs to obtain the real log absolute determinant. Near a singularity, both the determinant term and numerical conditioning require particular care.

Practical Optimization and Common Pitfalls

A reliable numerical procedure treats the determinant and residual calculations as parts of one objective, not as optional corrections. For each admissible coefficient, form the transformation matrix, compute its log absolute determinant, estimate the conditional regression coefficients, calculate the transformed residual sum of squares, and evaluate the concentrated likelihood. A grid can help locate promising regions, followed by a bounded one-dimensional optimizer. If the admissible set has multiple components, optimize within each component and compare the resulting maxima rather than allowing an algorithm to cross a singularity.

The spatial-error model uses the same overall strategy but a different transformed residual and a different conditional coefficient estimate. Its residual is \(A_\lambda(y-X\beta)\), not \(A_\lambda y-X\beta\). Applying the spatial-lag coefficient formula to the error model would therefore optimize the wrong sum of squares. Likewise, a determinant evaluated for one model or parameter value cannot be reused for another.

• Wrong: “The determinant is below one, so I add a positive adjustment.” Correction: Use \(\log|A|\) exactly as calculated; if \(|A|<1\), its log is negative.
• Wrong: “A small residual sum of squares must give the largest likelihood.” Correction: The likelihood also includes the log-determinant and the variance term, so compare the complete concentrated objective.
• Wrong: “Any coefficient with a finite numerical output is admissible.” Correction: Check that the transformation is nonsingular and that the coefficient satisfies the model’s stated admissibility restrictions.
• Wrong: “The variance estimate is \(SSR/(n-k)\) in the concentrated Gaussian likelihood.” Correction: Maximum likelihood uses \(SSR/n\); degrees-of-freedom adjustments serve a different purpose.