Spatial Econometrics — Lag, Error, and Spillovers

Python · PyMC · R (spatialreg, spdep)  ·  Download spatial-econometrics module

Spillover or Nuisance?

The areal examples put spatial structure in a random effect — a CAR prior that smooths a map. Spatial econometrics builds the neighbour structure into the regression itself, and the difference is substantive rather than stylistic. Under a spatial lag a region's outcome depends on its neighbours' outcomes: a genuine spillover, with ρ\rho measuring its strength. Under a spatial error the regressors are fine and omitted spatially-correlated factors merely leak through the residuals, so λ\lambda is a nuisance to correct, not an effect to interpret. Same WW, opposite meanings.

SAR: y=ρWy+Xβ+εSEM: y=Xβ+u,  u=λWu+ε\text{SAR: } y = \rho W y + X\beta + \varepsilon \qquad\qquad \text{SEM: } y = X\beta + u,\ \ u = \lambda W u + \varepsilon

On Anselin's Columbus data — 49 neighbourhoods, crime against income and housing value — OLS leaves Moran's I of 0.22 in its residuals (p=0.002p=0.002), so it is mis-specified. Both spatial models improve on it, and both engines agree on the estimates: ρ\rho is 0.39 from scratch, 0.41 in PyMC and 0.42 in R's lagsarlm; λ\lambda is 0.56 and 0.55.

Why the coefficient is not the answer

A spatial lag has a consequence ordinary regression lacks. Because y=(I−ρW)−1(Xβ+ε)y=(I-\rho W)^{-1}(X\beta+\varepsilon), changing a covariate in one region moves the outcome everywhere through the spatial multiplier, so each covariate has a direct effect on its own region, an indirect spillover onto others, and a total — the LeSage–Pace decomposition — and the coefficient βk\beta_k is none of them. Income's coefficient is −1.09, but its total effect is −1.82: reporting the coefficient alone would understate the policy impact by the entire indirect column. The from-scratch and R engines land on the same total to two decimals.

Columbus — income on crimeOLSSAR coefdirectindirecttotal
from scratch−1.60−1.09−1.14−0.67−1.82
R spatialreg−1.60−1.05−1.10−0.72−1.82

Choosing between them

Which specification, then? Here the notebook's textbook intuition — that a spatial error leaves the coefficient near OLS while a lag shrinks it — does not survive contact with the data. Income runs −1.60 (OLS) → −1.09 (SAR) → −0.95 (SEM), so SEM moves it further, not less, and R reproduces the same ordering (−1.60, −1.05, −0.96). The intuition describes expectations under a true model; one 49-neighbourhood sample need not obey it.

The formal answer comes from the Lagrange-multiplier tests, which exist precisely for this choice. Their robust versions are the ones to read: adjRSerr p=0.83p=0.83 says no error term is needed once a lag is allowed, while adjRSlag p=0.053p=0.053 says a lag is still needed once an error term is. AIC agrees — SAR 375 < SEM 377 < OLS 383. On Columbus the evidence favours the lag. That is a statistical verdict and not a causal one: the tests choose between two error structures, and only theory can license reading ρ\rho as neighbours actually causing crime.

Notebooks

Python Bayesian SAR and SEM from scratch — conjugate draws conditional on the spatial parameter, Metropolis on it with the exact log-determinant Jacobian, the LeSage–Pace effects decomposition, and a PyMC cross-check (Python) Open → R Independent R implementation via spatialreg — lagsarlm and errorsarlm by maximum likelihood, impacts for the effects, and the Lagrange-multiplier tests that adjudicate lag against error (R) Open →

Downloads

sp_sarsem.py Bayesian SAR and SEM sharing one sampler — conjugate β and σ² conditional on the spatial parameter, Metropolis on that parameter with the exact log-determinant Jacobian, and the LeSage–Pace direct/indirect/total effects (NumPy) columbus.csv Anselin's Columbus, Ohio data — 49 neighbourhoods with crime, household income and housing value col_adj.csv The 49 × 49 binary contiguity matrix columbus.geojson Neighbourhood boundaries for the maps

Spatial-Econometrics Module — Source Code

"""
spatecon.py -- SPATIAL ECONOMETRICS: spatial lag (SAR) and spatial error (SEM) models (from scratch).

Backs the notebooks in  "Spatial Econometrics -- Lag, Error, and Spillovers".

The areal projects put spatial structure in a random EFFECT (a CAR prior that smooths a map). Spatial
ECONOMETRICS instead builds the neighbour structure into the REGRESSION itself, and the distinction
is substantive, not cosmetic:

  SPATIAL LAG (SAR):   y = rho W y + X beta + eps
      the outcome in region i depends on the outcome in its NEIGHBOURS -- a genuine SPILLOVER
      (crime, prices, policy diffuse across borders). rho is the strength of the spillover.

  SPATIAL ERROR (SEM): y = X beta + u ,   u = lambda W u + eps
      the REGRESSORS are fine, but omitted spatially-correlated factors leak through the errors.
      lambda is a nuisance to be corrected, not a spillover to be interpreted.

The two look similar but mean opposite things: under SAR a shock to one region propagates to others
(a substantive multiplier); under SEM the spatial term is only about efficient/unbiased inference.
Choosing between them is the identification problem this literature is built around.

SAR has a second consequence that ordinary regression lacks: because y = (I - rho W)^{-1}(X beta +
eps), a change in x at region i moves y EVERYWHERE through the spatial multiplier (I - rho W)^{-1}.
So each covariate has a DIRECT effect (on its own region), an INDIRECT effect (spillover onto
others) and a TOTAL effect -- the LeSage-Pace decomposition -- and the naive coefficient beta_k is
none of these.

Estimation is Bayesian with a row-standardised W. Conditional on the spatial parameter the model is
an ordinary linear regression on the spatially-filtered data, so beta and sigma^2 are conjugate;
the spatial parameter itself is sampled by Metropolis using the exact log-determinant
log|I - rho W| (the Jacobian of the spatial transformation). Contrast with the areal CAR: same W,
but a simultaneous autoregression of the data rather than a hierarchical prior on an effect.
"""

import numpy as np


def row_standardise(W):
    W = np.asarray(W, float); rs = W.sum(1, keepdims=True); rs[rs == 0] = 1
    return W / rs


def morans_i(y, W):
    y = np.asarray(y, float); z = y - y.mean(); return (len(y) / W.sum()) * (z @ W @ z) / (z @ z)


def ols(y, X):
    XtXi = np.linalg.inv(X.T @ X); b = XtXi @ (X.T @ y)
    r = y - X @ b; s2 = (r @ r) / (len(y) - X.shape[1])
    return b, np.sqrt(np.diag(s2 * XtXi)), r


def _spatial_gibbs(y, X, W, rng, draws, burn, model, b_sd=100.0):
    """shared Bayesian sampler for SAR (model='lag') and SEM (model='error') on row-standardised W."""
    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape
    Wrs = row_standardise(W); evals = np.linalg.eigvals(Wrs).real
    lo, hi = 1.0 / evals.min() + 1e-3, 1.0 / evals.max() - 1e-3      # propriety bounds for the parameter
    par = 0.0; beta = np.zeros(k); sig2 = np.var(y); step = 0.1
    I = np.eye(n); P0 = np.eye(k) / b_sd ** 2
    B = np.empty((draws, k)); PAR = np.empty(draws); SIG = np.empty(draws); acc = 0
    def logdet(p):
        s, ld = np.linalg.slogdet(I - p * Wrs); return ld
    for it in range(draws + burn):
        A = I - par * Wrs
        if model == "lag":
            ystar, Xstar = A @ y, X                                  # (I-rhoW)y = Xb + eps
        else:
            ystar, Xstar = A @ y, A @ X                              # (I-lamW)(y-Xb)=eps -> By=BXb+eps
        # beta | par, sig2   (conjugate)
        V = np.linalg.inv(Xstar.T @ Xstar / sig2 + P0); m = V @ (Xstar.T @ ystar / sig2)
        beta = m + np.linalg.cholesky(V) @ rng.standard_normal(k)
        # sig2 | .   (inverse-gamma)
        resid = ystar - Xstar @ beta; sig2 = (resid @ resid) / rng.chisquare(n)
        # par | .   (Metropolis with exact log-determinant Jacobian)
        pp = par + step * rng.standard_normal()
        if lo < pp < hi:
            Ap = I - pp * Wrs
            if model == "lag":
                rp = Ap @ y - X @ beta
            else:
                rp = Ap @ (y - X @ beta)
            rc = (A @ y - X @ beta) if model == "lag" else (A @ (y - X @ beta))
            dlp = (logdet(pp) - logdet(par)) - (rp @ rp - rc @ rc) / (2 * sig2)
            if np.log(rng.random()) < dlp:
                par = pp; acc += 1
        if it < burn and it % 200 == 199:
            step *= np.exp((acc / 200 - 0.4) * 0.5); acc = 0
        if it >= burn:
            j = it - burn; B[j] = beta; PAR[j] = par; SIG[j] = np.sqrt(sig2)
    key = "rho" if model == "lag" else "lambda"
    return {"beta": B, key: PAR, "sigma": SIG}


def sar_gibbs(y, X, W, rng, draws=4000, burn=2000):
    return _spatial_gibbs(y, X, W, rng, draws, burn, "lag")


def sem_gibbs(y, X, W, rng, draws=4000, burn=2000):
    return _spatial_gibbs(y, X, W, rng, draws, burn, "error")


def sar_effects(beta_draws, rho_draws, W, col):
    """LeSage-Pace direct / indirect / total effects for covariate `col` from SAR posterior draws.
    Returns (direct, indirect, total) posterior-draw arrays."""
    Wrs = row_standardise(W); n = Wrs.shape[0]; I = np.eye(n)
    D = np.empty(len(rho_draws)); T = np.empty(len(rho_draws))
    for s in range(len(rho_draws)):
        M = np.linalg.inv(I - rho_draws[s] * Wrs)
        D[s] = beta_draws[s, col] * np.trace(M) / n                 # avg own-region effect
        T[s] = beta_draws[s, col] * M.sum() / n                     # avg total (own + spillovers)
    return D, T - D, T

References