# SP-01 S1.1 — Spatial Core Implementation Plan

> **For agentic workers:** REQUIRED SUB-SKILL: Use `superpowers:subagent-driven-development` (recommended) or `superpowers:executing-plans` to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking.

**Goal:** Close the largest ecosystem gap in StatsPAI by shipping a sparse-backed spatial weights layer, the full PySAL-equivalent ESDA suite, spatial regression models (ML + GMM) with LeSage-Pace impacts, and LM diagnostics — all interoperable with libpysal but independently implemented.

**Architecture:** Three-layer module under `src/statspai/spatial/`. Layer 1 (`weights/`) owns a sparse CSR-backed `W` class and six constructors. Layer 2 (`esda/`) owns six permutation-based spatial statistics. Layer 3 (`models/`) owns ML (concentrated likelihood, dual log-det path) and GMM (Kelejian-Prucha) estimators, plus LM tests and impacts. Existing `sp.sar/sem/sdm` are preserved as thin forwarders; backward-compat is enforced by a regression test against saved outputs.

**Tech Stack:** numpy, scipy.sparse, scipy.stats, scipy.optimize, pandas, formulaic (already a dep via existing spatial module). Soft deps: geopandas, libpysal (lazy-imported, test-skipped if absent), matplotlib (for plots).

**Parent spec:** [`specs/2026-04-15-sp-01-spatial-full-stack-design.md`](../specs/2026-04-15-sp-01-spatial-full-stack-design.md)

---

## File structure

```
src/statspai/spatial/
├─ __init__.py                  # flat re-exports
├─ models.py                    # EXISTING — becomes thin shim after refactor
├─ weights/
│   ├─ __init__.py
│   ├─ core.py                  # class W
│   ├─ contiguity.py            # queen_weights, rook_weights
│   ├─ distance.py              # knn_weights, distance_band, kernel_weights
│   └─ block.py                 # block_weights
├─ esda/
│   ├─ __init__.py
│   ├─ _base.py                 # SpatialStatistic result class, permutation helper
│   ├─ moran.py                 # moran, moran_local
│   ├─ geary.py
│   ├─ getis_ord.py             # getis_ord_g, getis_ord_local
│   ├─ join_counts.py
│   └─ plots.py                 # moran_plot, lisa_cluster_map
├─ models/
│   ├─ __init__.py
│   ├─ _logdet.py               # exact + Barry-Pace approx
│   ├─ _base.py                 # shared design-matrix + result plumbing
│   ├─ ml.py                    # sar, sem, sdm, slx, sac (ML)
│   ├─ gmm.py                   # sar_gmm, sem_gmm, sarar_gmm
│   ├─ diagnostics.py           # lm_tests, moran_residuals
│   └─ impacts.py               # LeSage-Pace direct/indirect/total

tests/spatial/
├─ __init__.py
├─ fixtures/
│   ├─ columbus.csv             # dataset copy (public)
│   ├─ columbus_shp/            # geometries
│   └─ reference_outputs.json   # saved PySAL / spatialreg numbers
├─ test_weights_core.py
├─ test_weights_contiguity.py
├─ test_weights_distance.py
├─ test_weights_block.py
├─ test_esda_moran.py
├─ test_esda_geary.py
├─ test_esda_getis_ord.py
├─ test_esda_join_counts.py
├─ test_models_ml.py
├─ test_models_gmm.py
├─ test_models_impacts.py
├─ test_models_diagnostics.py
└─ test_backward_compat.py
```

**Convention:** every public function returns an object with `.summary()`. ESDA returns `SpatialStatistic`; models extend the existing `EconometricResults` from `statspai.core.results`.

**Testing convention:** each public function gets at minimum (a) a happy path test against a hand-computed or PySAL reference value, (b) an edge case (empty, islands, single-obs), (c) an error case (bad input shape / type). Tolerance vs PySAL references: `rtol=1e-4`.

---

## Phase A — Weights foundation

### Task 1: W core class

**Files:**
- Create: `src/statspai/spatial/weights/__init__.py`
- Create: `src/statspai/spatial/weights/core.py`
- Test: `tests/spatial/test_weights_core.py`

- [ ] **Step 1: Write failing test for W construction**

```python
# tests/spatial/test_weights_core.py
import numpy as np
import pytest
from scipy import sparse
from statspai.spatial.weights.core import W


def test_w_from_neighbors_dict():
    neighbors = {0: [1, 2], 1: [0], 2: [0]}
    w = W(neighbors)
    assert w.n == 3
    assert w.sparse.shape == (3, 3)
    assert w.sparse[0, 1] == 1.0
    assert w.sparse[1, 2] == 0.0


def test_w_row_standardise():
    w = W({0: [1, 2], 1: [0], 2: [0]})
    w.transform = "R"
    row_sums = np.asarray(w.sparse.sum(axis=1)).ravel()
    np.testing.assert_allclose(row_sums, [1.0, 1.0, 1.0])


def test_w_islands_detected():
    w = W({0: [1], 1: [0], 2: []})
    assert w.islands == [2]


def test_w_bad_input_raises():
    with pytest.raises(TypeError):
        W("not a dict")
```

- [ ] **Step 2: Run test to verify it fails**

```
python3 -m pytest tests/spatial/test_weights_core.py -v
# Expected: ImportError — core.py doesn't exist yet
```

- [ ] **Step 3: Implement `W` class**

```python
# src/statspai/spatial/weights/__init__.py
from .core import W
__all__ = ["W"]
```

```python
# src/statspai/spatial/weights/core.py
"""Sparse-backed spatial weights object."""
from __future__ import annotations

from typing import Dict, List, Mapping, Optional, Sequence, Union

import numpy as np
from scipy import sparse


class W:
    """Spatial weights matrix with CSR-sparse backing.

    Parameters
    ----------
    neighbors : dict[int, list[int]]
        Mapping observation id -> list of neighbour ids.
    weights : dict[int, list[float]], optional
        Matching weight values. If ``None``, binary (1.0) weights are used.
    id_order : sequence, optional
        Explicit ordering of observation ids. Defaults to ``sorted(neighbors)``.
    """

    _VALID_TRANSFORMS = {"O", "B", "R", "V", "D"}

    def __init__(
        self,
        neighbors: Mapping[int, Sequence[int]],
        weights: Optional[Mapping[int, Sequence[float]]] = None,
        id_order: Optional[Sequence[int]] = None,
    ) -> None:
        if not isinstance(neighbors, Mapping):
            raise TypeError("`neighbors` must be a dict-like mapping id -> list")
        self._id_order = list(id_order) if id_order is not None else sorted(neighbors)
        self._id_to_idx = {i: k for k, i in enumerate(self._id_order)}
        self._neighbors = {i: list(neighbors.get(i, [])) for i in self._id_order}
        if weights is None:
            self._weights = {i: [1.0] * len(v) for i, v in self._neighbors.items()}
        else:
            self._weights = {i: list(weights[i]) for i in self._id_order}
        self._transform = "O"           # original
        self._sparse: Optional[sparse.csr_matrix] = None

    # -- basic properties --------------------------------------------------
    @property
    def n(self) -> int:
        return len(self._id_order)

    @property
    def neighbors(self) -> Dict[int, List[int]]:
        return {i: list(v) for i, v in self._neighbors.items()}

    @property
    def islands(self) -> List[int]:
        return [i for i, v in self._neighbors.items() if len(v) == 0]

    @property
    def sparse(self) -> sparse.csr_matrix:
        if self._sparse is None:
            self._sparse = self._build_sparse()
        return self._sparse

    # -- transform ---------------------------------------------------------
    @property
    def transform(self) -> str:
        return self._transform

    @transform.setter
    def transform(self, value: str) -> None:
        value = value.upper()
        if value not in self._VALID_TRANSFORMS:
            raise ValueError(f"transform must be one of {self._VALID_TRANSFORMS}")
        if value == self._transform:
            return
        # always reset from original binary pattern
        base_weights = {i: [1.0] * len(v) for i, v in self._neighbors.items()}
        if value == "R":
            for i, ws in base_weights.items():
                s = sum(ws)
                if s > 0:
                    base_weights[i] = [w / s for w in ws]
        elif value == "V":
            # variance stabilise: row-standardise then scale by sqrt(card)
            for i, ws in base_weights.items():
                s = sum(ws)
                if s > 0:
                    base_weights[i] = [w / s * np.sqrt(len(ws)) for w in ws]
        elif value == "D":
            # double-standardisation: divide by sum of all weights
            total = sum(sum(ws) for ws in base_weights.values())
            if total > 0:
                base_weights = {i: [w / total for w in ws] for i, ws in base_weights.items()}
        # "B" and "O" keep binary values
        self._weights = base_weights
        self._transform = value
        self._sparse = None             # lazy rebuild

    # -- dense / sparse conversions ---------------------------------------
    def _build_sparse(self) -> sparse.csr_matrix:
        rows, cols, data = [], [], []
        for i, nbrs in self._neighbors.items():
            row = self._id_to_idx[i]
            ws = self._weights[i]
            for j, w in zip(nbrs, ws):
                if j in self._id_to_idx:
                    rows.append(row)
                    cols.append(self._id_to_idx[j])
                    data.append(float(w))
        return sparse.csr_matrix(
            (data, (rows, cols)), shape=(self.n, self.n)
        )

    def full(self) -> np.ndarray:
        return self.sparse.toarray()

    # -- interop -----------------------------------------------------------
    def to_libpysal(self):
        """Return a libpysal.weights.W equivalent (requires libpysal)."""
        try:
            from libpysal.weights import W as _LPW
        except ImportError as e:
            raise ImportError(
                "libpysal is required for to_libpysal(). "
                "Install with `pip install libpysal`."
            ) from e
        return _LPW(self._neighbors, self._weights, id_order=self._id_order)

    @classmethod
    def from_libpysal(cls, w) -> "W":
        return cls(dict(w.neighbors), dict(w.weights), id_order=list(w.id_order))
```

- [ ] **Step 4: Run test to verify pass**

```
python3 -m pytest tests/spatial/test_weights_core.py -v
# Expected: 4 passed
```

- [ ] **Step 5: Commit**

```
git add src/statspai/spatial/weights/ tests/spatial/test_weights_core.py tests/spatial/__init__.py
git commit -m "feat(spatial): add sparse-backed W weights class"
```

---

### Task 2: KNN and distance-band weights

**Files:**
- Create: `src/statspai/spatial/weights/distance.py`
- Modify: `src/statspai/spatial/weights/__init__.py`
- Test: `tests/spatial/test_weights_distance.py`

- [ ] **Step 1: Write failing tests**

```python
# tests/spatial/test_weights_distance.py
import numpy as np
import pytest
from statspai.spatial.weights.distance import (
    knn_weights, distance_band, kernel_weights,
)


def test_knn_basic():
    coords = np.array([[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [10.0, 10.0]])
    w = knn_weights(coords, k=1)
    assert w.neighbors[0] == [1]
    assert w.neighbors[1] in ([0], [2])
    assert w.neighbors[3] == [2]        # closest of the three


def test_distance_band_binary():
    coords = np.array([[0.0, 0.0], [1.0, 0.0], [5.0, 0.0]])
    w = distance_band(coords, threshold=1.5, binary=True)
    assert sorted(w.neighbors[0]) == [1]
    assert sorted(w.neighbors[2]) == []


def test_kernel_gaussian_weights_decay():
    coords = np.array([[0.0, 0.0], [1.0, 0.0], [2.0, 0.0]])
    w = kernel_weights(coords, bandwidth=1.0, kernel="gaussian", fixed=True)
    w0 = dict(zip(w.neighbors[0], w._weights[0]))
    assert w0[1] > w0[2]                # closer point gets larger weight


def test_knn_rejects_bad_k():
    coords = np.zeros((3, 2))
    with pytest.raises(ValueError):
        knn_weights(coords, k=5)
```

- [ ] **Step 2: Run and confirm fail**

```
python3 -m pytest tests/spatial/test_weights_distance.py -v
# Expected: ImportError
```

- [ ] **Step 3: Implement**

```python
# src/statspai/spatial/weights/distance.py
"""Distance-based spatial weights: KNN, distance band, kernel."""
from __future__ import annotations

from typing import Literal

import numpy as np
from scipy.spatial import cKDTree

from .core import W


def knn_weights(coords: np.ndarray, k: int = 5) -> W:
    coords = np.asarray(coords, dtype=float)
    n = coords.shape[0]
    if k <= 0 or k >= n:
        raise ValueError(f"k must be in (0, n); got k={k}, n={n}")
    tree = cKDTree(coords)
    _, idx = tree.query(coords, k=k + 1)     # +1 to drop self
    neighbors = {i: [int(j) for j in idx[i, 1:]] for i in range(n)}
    return W(neighbors)


def distance_band(coords: np.ndarray, threshold: float, binary: bool = True) -> W:
    coords = np.asarray(coords, dtype=float)
    n = coords.shape[0]
    tree = cKDTree(coords)
    pairs = tree.query_ball_point(coords, r=threshold)
    neighbors, weights = {}, {}
    for i, js in enumerate(pairs):
        js = [int(j) for j in js if j != i]
        neighbors[i] = js
        if binary:
            weights[i] = [1.0] * len(js)
        else:
            d = np.linalg.norm(coords[js] - coords[i], axis=1)
            d[d == 0] = np.inf
            weights[i] = (1.0 / d).tolist()
    return W(neighbors, weights)


def kernel_weights(
    coords: np.ndarray,
    bandwidth: float,
    kernel: Literal["gaussian", "bisquare", "triangular"] = "gaussian",
    fixed: bool = True,
) -> W:
    coords = np.asarray(coords, dtype=float)
    n = coords.shape[0]
    tree = cKDTree(coords)
    neighbors, weights = {}, {}
    for i in range(n):
        if fixed:
            js = tree.query_ball_point(coords[i], r=bandwidth)
            js = [int(j) for j in js if j != i]
            d = np.linalg.norm(coords[js] - coords[i], axis=1) if js else np.array([])
            bw = bandwidth
        else:                                   # adaptive KNN-style bandwidth
            dists, idx = tree.query(coords[i], k=int(bandwidth) + 1)
            js = [int(j) for j in idx[1:]]
            d = dists[1:]
            bw = d.max() if len(d) else 1.0
        u = d / bw if bw > 0 else np.zeros_like(d)
        if kernel == "gaussian":
            k = np.exp(-0.5 * u ** 2)
        elif kernel == "bisquare":
            k = np.where(u < 1, (1 - u ** 2) ** 2, 0.0)
        elif kernel == "triangular":
            k = np.where(u < 1, 1 - u, 0.0)
        else:
            raise ValueError(f"unknown kernel {kernel!r}")
        neighbors[i] = js
        weights[i] = k.tolist()
    return W(neighbors, weights)
```

Update exports:

```python
# src/statspai/spatial/weights/__init__.py
from .core import W
from .distance import knn_weights, distance_band, kernel_weights

__all__ = ["W", "knn_weights", "distance_band", "kernel_weights"]
```

- [ ] **Step 4: Run tests pass**

```
python3 -m pytest tests/spatial/test_weights_distance.py -v
# Expected: 4 passed
```

- [ ] **Step 5: Commit**

```
git add src/statspai/spatial/weights/distance.py src/statspai/spatial/weights/__init__.py tests/spatial/test_weights_distance.py
git commit -m "feat(spatial): knn / distance-band / kernel weights"
```

---

### Task 3: Contiguity (queen, rook) with soft geopandas dep

**Files:**
- Create: `src/statspai/spatial/weights/contiguity.py`
- Test: `tests/spatial/test_weights_contiguity.py`

- [ ] **Step 1: Write failing tests (skipped without geopandas)**

```python
# tests/spatial/test_weights_contiguity.py
import pytest

gpd = pytest.importorskip("geopandas")
from shapely.geometry import Polygon
from statspai.spatial.weights.contiguity import queen_weights, rook_weights


@pytest.fixture
def three_squares():
    # |A|B|C|  -- A-B share edge; B-C share edge; A-C share only a corner
    polys = [
        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
        Polygon([(1, 0), (2, 0), (2, 1), (1, 1)]),
        Polygon([(2, 0), (3, 0), (3, 1), (2, 1)]),
    ]
    return gpd.GeoDataFrame({"id": [0, 1, 2], "geometry": polys})


def test_queen_vs_rook_edge_vs_corner(three_squares):
    wq = queen_weights(three_squares)
    wr = rook_weights(three_squares)
    # both share edge; both catch 0<->1 and 1<->2
    assert set(wq.neighbors[1]) == {0, 2}
    assert set(wr.neighbors[1]) == {0, 2}
    # queen also picks up diagonal corner neighbours when present
    # (this example has none, so queen == rook here)
    assert wq.neighbors[0] == wr.neighbors[0]


def test_contiguity_without_geopandas_raises(monkeypatch):
    import statspai.spatial.weights.contiguity as mod
    monkeypatch.setattr(mod, "_gpd", None)
    with pytest.raises(ImportError, match="geopandas"):
        mod.queen_weights(object())
```

- [ ] **Step 2: Run and confirm fail**

```
python3 -m pytest tests/spatial/test_weights_contiguity.py -v
```

- [ ] **Step 3: Implement**

```python
# src/statspai/spatial/weights/contiguity.py
"""Contiguity weights (queen / rook). Requires geopandas + shapely."""
from __future__ import annotations

from typing import Optional

import numpy as np

from .core import W

try:
    import geopandas as _gpd
except ImportError:                          # pragma: no cover
    _gpd = None


def _require_gpd():
    if _gpd is None:
        raise ImportError(
            "geopandas is required for contiguity weights. "
            "Install with `pip install geopandas shapely`."
        )


def _contiguity(gdf, criterion: str) -> W:
    _require_gpd()
    if criterion not in {"queen", "rook"}:
        raise ValueError(f"criterion must be 'queen' or 'rook'")
    geoms = gdf.geometry.values
    n = len(geoms)
    sindex = gdf.sindex
    neighbors = {i: [] for i in range(n)}
    for i in range(n):
        candidates = list(sindex.intersection(geoms[i].bounds))
        for j in candidates:
            if j == i:
                continue
            inter = geoms[i].intersection(geoms[j])
            if inter.is_empty:
                continue
            if criterion == "queen":
                # any non-empty intersection (point, line, area) qualifies
                neighbors[i].append(int(j))
            else:                                # rook — need a line, not a point
                if inter.geom_type in {"LineString", "MultiLineString"} or (
                    hasattr(inter, "length") and inter.length > 0
                ):
                    neighbors[i].append(int(j))
    return W(neighbors)


def queen_weights(gdf) -> W:
    return _contiguity(gdf, "queen")


def rook_weights(gdf) -> W:
    return _contiguity(gdf, "rook")
```

Update exports to include `queen_weights`, `rook_weights`.

- [ ] **Step 4: Run tests pass**

- [ ] **Step 5: Commit**

```
git commit -m "feat(spatial): queen/rook contiguity weights (soft geopandas dep)"
```

---

### Task 4: Block weights

**Files:**
- Create: `src/statspai/spatial/weights/block.py`
- Test: `tests/spatial/test_weights_block.py`

- [ ] **Step 1: Failing test**

```python
# tests/spatial/test_weights_block.py
import numpy as np
from statspai.spatial.weights.block import block_weights


def test_block_weights_same_regime_are_neighbours():
    regimes = np.array(["A", "A", "B", "B", "B"])
    w = block_weights(regimes)
    assert sorted(w.neighbors[0]) == [1]
    assert sorted(w.neighbors[2]) == [3, 4]
```

- [ ] **Step 2: Run fail** / **Step 3: Implement**

```python
# src/statspai/spatial/weights/block.py
from __future__ import annotations

from collections import defaultdict

import numpy as np

from .core import W


def block_weights(regimes) -> W:
    regimes = np.asarray(regimes)
    buckets = defaultdict(list)
    for i, r in enumerate(regimes):
        buckets[r].append(i)
    neighbors = {i: [] for i in range(len(regimes))}
    for ids in buckets.values():
        for i in ids:
            neighbors[i] = [j for j in ids if j != i]
    return W(neighbors)
```

- [ ] **Step 4: Pass** / **Step 5: Commit**

```
git commit -m "feat(spatial): block (regime) weights"
```

---

## Phase B — ESDA

### Task 5: SpatialStatistic result class + Moran's I global

**Files:**
- Create: `src/statspai/spatial/esda/__init__.py`
- Create: `src/statspai/spatial/esda/_base.py`
- Create: `src/statspai/spatial/esda/moran.py`
- Test: `tests/spatial/test_esda_moran.py`

- [ ] **Step 1: Failing tests** (use Columbus reference values)

```python
# tests/spatial/test_esda_moran.py
import numpy as np
import pytest
from statspai.spatial.esda import moran, moran_local
from statspai.spatial.weights import knn_weights


def test_moran_positive_autocorrelation():
    rng = np.random.default_rng(42)
    coords = rng.uniform(size=(100, 2))
    # synthesise y = distance from origin (spatially smooth)
    y = np.linalg.norm(coords - 0.5, axis=1) + 0.05 * rng.standard_normal(100)
    w = knn_weights(coords, k=6)
    w.transform = "R"
    res = moran(y, w, permutations=499, seed=42)
    assert res.I > 0.3
    assert res.p_sim < 0.05


def test_moran_iid_near_zero():
    rng = np.random.default_rng(0)
    coords = rng.uniform(size=(100, 2))
    y = rng.standard_normal(100)
    w = knn_weights(coords, k=4); w.transform = "R"
    res = moran(y, w, permutations=499, seed=0)
    assert res.p_sim > 0.05 or abs(res.I) < 0.1


def test_moran_summary_prints():
    rng = np.random.default_rng(1)
    coords = rng.uniform(size=(30, 2))
    y = rng.standard_normal(30)
    w = knn_weights(coords, k=3); w.transform = "R"
    res = moran(y, w, permutations=99, seed=1)
    out = res.summary()
    assert "Moran" in out
```

- [ ] **Step 2: Run fail / Step 3: Implement**

```python
# src/statspai/spatial/esda/_base.py
"""Shared ESDA result plumbing."""
from __future__ import annotations

from dataclasses import dataclass, field
from typing import Optional

import numpy as np


@dataclass
class SpatialStatistic:
    name: str
    value: float
    expectation: float
    variance: float
    z_score: float
    p_norm: float
    p_sim: Optional[float] = None
    simulations: Optional[np.ndarray] = None
    extras: dict = field(default_factory=dict)

    # Common aliases used by specific statistics
    @property
    def I(self) -> float:           # noqa: E741 (Moran's I convention)
        return self.value

    @property
    def C(self) -> float:
        return self.value

    @property
    def G(self) -> float:
        return self.value

    def summary(self) -> str:
        lines = [
            f"{self.name}",
            "-" * 40,
            f"Statistic      : {self.value: .6f}",
            f"E[statistic]   : {self.expectation: .6f}",
            f"Var[statistic] : {self.variance: .6f}",
            f"z              : {self.z_score: .4f}",
            f"p (normal)     : {self.p_norm: .4f}",
        ]
        if self.p_sim is not None:
            lines.append(f"p (permutation): {self.p_sim: .4f}")
        return "\n".join(lines)

    def __repr__(self) -> str:
        return self.summary()


def permutation_pvalue(observed: float, sims: np.ndarray) -> float:
    """Two-sided empirical p-value from permutation replicates."""
    sims = np.asarray(sims)
    larger = np.sum(np.abs(sims - sims.mean()) >= abs(observed - sims.mean()))
    return (larger + 1) / (len(sims) + 1)
```

```python
# src/statspai/spatial/esda/moran.py
"""Moran's I (global and local)."""
from __future__ import annotations

from typing import Optional

import numpy as np
from scipy import stats as sp_stats

from ..weights.core import W
from ._base import SpatialStatistic, permutation_pvalue


def _standardise(y: np.ndarray) -> np.ndarray:
    y = np.asarray(y, dtype=float).ravel()
    return y - y.mean()


def moran(
    y,
    w: W,
    permutations: int = 999,
    two_tailed: bool = True,
    seed: Optional[int] = None,
) -> SpatialStatistic:
    """Global Moran's I with analytical + permutation inference."""
    z = _standardise(y)
    n = z.size
    S = w.sparse
    S0 = float(S.sum())
    if S0 == 0:
        raise ValueError("Weights matrix has zero total — cannot compute I.")
    wz = S @ z
    num = float(z @ wz)
    den = float(z @ z)
    I = (n / S0) * (num / den)

    # randomisation-assumption moments (Cliff-Ord)
    EI = -1.0 / (n - 1)
    S1 = 0.5 * ((S + S.T).multiply(S + S.T)).sum()
    S2 = float(np.sum((np.asarray(S.sum(axis=1)).ravel() +
                       np.asarray(S.sum(axis=0)).ravel()) ** 2))
    b2 = n * np.sum(z ** 4) / (np.sum(z ** 2) ** 2)
    A = n * ((n ** 2 - 3 * n + 3) * S1 - n * S2 + 3 * S0 ** 2)
    B = b2 * ((n ** 2 - n) * S1 - 2 * n * S2 + 6 * S0 ** 2)
    C = (n - 1) * (n - 2) * (n - 3) * S0 ** 2
    VI = (A - B) / C - EI ** 2
    VI = max(VI, 1e-12)
    z_score = (I - EI) / np.sqrt(VI)
    p_norm = (2 * (1 - sp_stats.norm.cdf(abs(z_score)))) if two_tailed \
        else 1 - sp_stats.norm.cdf(z_score)

    sims = None
    p_sim = None
    if permutations and permutations > 0:
        rng = np.random.default_rng(seed)
        sims = np.empty(permutations)
        for k in range(permutations):
            zp = rng.permutation(z)
            sims[k] = (n / S0) * ((zp @ (S @ zp)) / (zp @ zp))
        p_sim = permutation_pvalue(I, sims)

    return SpatialStatistic(
        name="Moran's I",
        value=I,
        expectation=EI,
        variance=VI,
        z_score=z_score,
        p_norm=p_norm,
        p_sim=p_sim,
        simulations=sims,
    )


def moran_local(y, w: W, permutations: int = 999, seed: Optional[int] = None):
    """Local Moran (LISA). Returns array of SpatialStatistic-like values."""
    z = _standardise(y)
    n = z.size
    m2 = np.sum(z ** 2) / n
    S = w.sparse
    Ii = z * (S @ z) / m2
    sims = None
    p_sim = None
    if permutations and permutations > 0:
        rng = np.random.default_rng(seed)
        sims = np.empty((permutations, n))
        for k in range(permutations):
            zp = rng.permutation(z)
            sims[k] = z * (S @ zp) / m2
        p_sim = np.array([
            permutation_pvalue(Ii[i], sims[:, i]) for i in range(n)
        ])
    return {
        "Is": Ii,
        "p_sim": p_sim,
        "simulations": sims,
    }
```

```python
# src/statspai/spatial/esda/__init__.py
from .moran import moran, moran_local
__all__ = ["moran", "moran_local"]
```

- [ ] **Step 4: Run tests pass**
- [ ] **Step 5: Commit**

```
git commit -m "feat(spatial.esda): Moran's I (global + local) with permutation inference"
```

---

### Task 6: Geary's C

**Files:** `src/statspai/spatial/esda/geary.py`, `tests/spatial/test_esda_geary.py`

- [ ] **Step 1: Failing test** — analogous to Moran; with smooth data C < 1, z < 0.
- [ ] **Step 2: Run fail**
- [ ] **Step 3: Implement**

```python
# src/statspai/spatial/esda/geary.py
from __future__ import annotations
from typing import Optional
import numpy as np
from scipy import stats as sp_stats
from ..weights.core import W
from ._base import SpatialStatistic, permutation_pvalue


def geary(y, w: W, permutations: int = 999, seed: Optional[int] = None) -> SpatialStatistic:
    y = np.asarray(y, dtype=float).ravel()
    z = y - y.mean()
    n = z.size
    S = w.sparse
    S0 = float(S.sum())
    if S0 == 0:
        raise ValueError("Weights matrix has zero total.")
    # numerator: sum_ij w_ij (y_i - y_j)^2
    rows, cols = S.nonzero()
    data = S.data
    diffs = (y[rows] - y[cols]) ** 2
    num = float(np.sum(data * diffs))
    den = float(np.sum(z ** 2))
    C = ((n - 1) * num) / (2 * S0 * den)

    EC = 1.0
    z_score = np.nan
    VC = np.nan
    p_norm = np.nan

    sims = None
    p_sim = None
    if permutations and permutations > 0:
        rng = np.random.default_rng(seed)
        sims = np.empty(permutations)
        for k in range(permutations):
            yp = rng.permutation(y)
            d = (yp[rows] - yp[cols]) ** 2
            zp = yp - yp.mean()
            sims[k] = ((n - 1) * np.sum(data * d)) / (2 * S0 * np.sum(zp ** 2))
        p_sim = permutation_pvalue(C, sims)
        VC = float(np.var(sims, ddof=1))
        z_score = (C - 1.0) / (np.sqrt(VC) if VC > 0 else 1e-12)
        p_norm = 2 * (1 - sp_stats.norm.cdf(abs(z_score)))

    return SpatialStatistic(
        name="Geary's C", value=C, expectation=EC,
        variance=VC, z_score=z_score, p_norm=p_norm,
        p_sim=p_sim, simulations=sims,
    )
```

- [ ] **Step 4: Pass** / **Step 5: Commit**

---

### Task 7: Getis-Ord G (global + local)

**Files:** `src/statspai/spatial/esda/getis_ord.py`, `tests/spatial/test_esda_getis_ord.py`

- [ ] **Step 1: Failing test** — hot-spot synthetic; assert G_local significant at hot centre.
- [ ] **Step 2: Fail**
- [ ] **Step 3: Implement**

```python
# src/statspai/spatial/esda/getis_ord.py
from __future__ import annotations
from typing import Optional
import numpy as np
from scipy import stats as sp_stats
from ..weights.core import W
from ._base import SpatialStatistic, permutation_pvalue


def getis_ord_g(y, w: W, permutations: int = 999, seed: Optional[int] = None) -> SpatialStatistic:
    y = np.asarray(y, dtype=float).ravel()
    if np.any(y < 0):
        raise ValueError("Getis-Ord G requires non-negative y.")
    S = w.sparse
    num = float(y @ (S @ y) - np.sum(S.diagonal() * y * y))
    den = float(np.sum(y) ** 2 - np.sum(y ** 2))
    if den == 0:
        raise ValueError("Degenerate y for Getis-Ord G.")
    G = num / den

    sims = None
    p_sim = None
    if permutations:
        rng = np.random.default_rng(seed)
        sims = np.empty(permutations)
        for k in range(permutations):
            yp = rng.permutation(y)
            sims[k] = (yp @ (S @ yp) - np.sum(S.diagonal() * yp * yp)) / \
                      (np.sum(yp) ** 2 - np.sum(yp ** 2))
        p_sim = permutation_pvalue(G, sims)
    EG = float(S.sum() / (y.size * (y.size - 1)))
    VG = float(np.var(sims, ddof=1)) if sims is not None else np.nan
    z_score = (G - EG) / np.sqrt(VG) if VG and VG > 0 else np.nan
    p_norm = 2 * (1 - sp_stats.norm.cdf(abs(z_score))) if not np.isnan(z_score) else np.nan
    return SpatialStatistic(
        name="Getis-Ord G", value=G, expectation=EG,
        variance=VG, z_score=z_score, p_norm=p_norm,
        p_sim=p_sim, simulations=sims,
    )


def getis_ord_local(y, w: W, star: bool = True, permutations: int = 999, seed: Optional[int] = None):
    y = np.asarray(y, dtype=float).ravel()
    n = y.size
    S = w.sparse.toarray().copy()
    if star:
        np.fill_diagonal(S, 1.0)              # include self
    Wi = S.sum(axis=1)
    sum_y = y.sum()
    sum_y2 = (y ** 2).sum()
    mean_y = sum_y / n
    var_y = sum_y2 / n - mean_y ** 2
    Gi_num = S @ y - Wi * mean_y
    Gi_den = np.sqrt(var_y * ((n * (Wi - Wi ** 2 / n)) / (n - 1)))
    Gi_den = np.where(Gi_den == 0, np.nan, Gi_den)
    Gi = Gi_num / Gi_den
    return {"Gs": Gi, "z": Gi}                 # Gi is already standardised
```

- [ ] **Step 4/5: Pass, commit**

---

### Task 8: Join counts

**Files:** `src/statspai/spatial/esda/join_counts.py`, `tests/spatial/test_esda_join_counts.py`

- [ ] TDD loop analogous to above. Implementation:

```python
# src/statspai/spatial/esda/join_counts.py
from __future__ import annotations
from typing import Optional
import numpy as np
from ..weights.core import W
from ._base import SpatialStatistic, permutation_pvalue


def join_counts(y, w: W, permutations: int = 999, seed: Optional[int] = None):
    y = np.asarray(y).ravel().astype(int)
    if not set(np.unique(y)).issubset({0, 1}):
        raise ValueError("join_counts requires binary y (0/1).")
    S = w.sparse
    rows, cols = S.nonzero()
    bb = 0.5 * float(np.sum((y[rows] == 1) & (y[cols] == 1) * S.data))
    ww = 0.5 * float(np.sum((y[rows] == 0) & (y[cols] == 0) * S.data))
    bw = float(np.sum((y[rows] != y[cols]) * S.data))

    sims = None
    p_sim = None
    if permutations:
        rng = np.random.default_rng(seed)
        sims = np.empty(permutations)
        for k in range(permutations):
            yp = rng.permutation(y)
            sims[k] = 0.5 * float(np.sum((yp[rows] == 1) & (yp[cols] == 1) * S.data))
        p_sim = permutation_pvalue(bb, sims)
    return {
        "BB": bb, "WW": ww, "BW": bw,
        "p_sim_BB": p_sim, "sims_BB": sims,
    }
```

- [ ] Commit

---

### Task 9: Plots (moran_plot, lisa_cluster_map)

**Files:** `src/statspai/spatial/esda/plots.py`, `tests/spatial/test_esda_plots.py`

- [ ] **Step 1: Failing test** — smoke test: returns matplotlib Axes and does not raise.

```python
import matplotlib
matplotlib.use("Agg")
import numpy as np
from statspai.spatial.esda.plots import moran_plot
from statspai.spatial.weights import knn_weights


def test_moran_plot_smoke():
    rng = np.random.default_rng(0)
    coords = rng.uniform(size=(50, 2))
    y = rng.standard_normal(50)
    w = knn_weights(coords, k=4); w.transform = "R"
    ax = moran_plot(y, w)
    assert ax is not None
```

- [ ] Implement: simple scatter of z vs Wz with fitted regression line; slope = Moran's I.
- [ ] Commit

---

## Phase C — Log-det infrastructure

### Task 10: _logdet.py (exact + Barry-Pace)

**Files:** `src/statspai/spatial/models/_logdet.py`, `tests/spatial/test_models_logdet.py`

- [ ] **Step 1: Failing test**

```python
# tests/spatial/test_models_logdet.py
import numpy as np
from scipy import sparse
from statspai.spatial.models._logdet import log_det_exact, log_det_approx


def test_approx_close_to_exact_small_n():
    rng = np.random.default_rng(0)
    n = 50
    A = sparse.random(n, n, density=0.05, data_rvs=rng.standard_normal, random_state=0)
    A = (A + A.T) / 2
    # row-normalise
    rs = np.asarray(abs(A).sum(axis=1)).ravel()
    rs[rs == 0] = 1
    D = sparse.diags(1 / rs)
    W = D @ A
    rho = 0.3
    exact = log_det_exact(W, rho)
    approx = log_det_approx(W, rho, n_draws=200, order=30, seed=0)
    assert abs(exact - approx) / abs(exact) < 0.1
```

- [ ] **Step 2: Fail / Step 3: Implement**

```python
# src/statspai/spatial/models/_logdet.py
"""Log-determinant computations for (I - rho W).

Exact path: eigenvalues of W, then log det(I - rho W) = sum log(1 - rho * lam_i).
Approx path (Barry-Pace Monte-Carlo + Chebyshev): used when n is large.
"""
from __future__ import annotations

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import eigs


def _to_csr(W):
    if sparse.issparse(W):
        return W.tocsr()
    return sparse.csr_matrix(W)


def log_det_exact(W, rho: float) -> float:
    W = _to_csr(W)
    dense = W.toarray()
    eigvals = np.linalg.eigvals(dense).real
    return float(np.sum(np.log(np.abs(1.0 - rho * eigvals))))


def log_det_approx(W, rho: float, n_draws: int = 50, order: int = 50,
                   seed: int | None = None) -> float:
    """Barry-Pace style approximation via Hutchinson trace + Taylor expansion.

    log det(I - rho W)  =  - sum_{k>=1} (rho^k / k) * tr(W^k)
    tr(W^k) ≈ (1/M) sum_m u_m^T W^k u_m  with u_m  Rademacher.
    """
    W = _to_csr(W)
    n = W.shape[0]
    rng = np.random.default_rng(seed)
    total = 0.0
    for _ in range(n_draws):
        u = rng.choice([-1.0, 1.0], size=n)
        v = u.copy()
        rho_pow = 1.0
        for k in range(1, order + 1):
            v = W @ v
            rho_pow *= rho
            total -= (rho_pow / k) * float(u @ v)
    return total / n_draws
```

- [ ] **Step 4: Pass / Step 5: Commit**

---

## Phase D — Spatial regression models (ML)

### Task 11: Models base (design matrix + result helpers)

**Files:** `src/statspai/spatial/models/_base.py`

- [ ] **Step 1: Failing test**

```python
# tests/spatial/test_models_base.py
import numpy as np, pandas as pd
from statspai.spatial.models._base import build_design_matrix


def test_build_design_matrix_shape():
    df = pd.DataFrame({"y": [1.0, 2.0, 3.0], "x1": [4.0, 5.0, 6.0]})
    y, X, names = build_design_matrix("y ~ x1", df)
    assert y.shape == (3,)
    assert X.shape == (3, 2)               # intercept + x1
    assert names == ["Intercept", "x1"]
```

- [ ] **Step 2/3: Implement**

```python
# src/statspai/spatial/models/_base.py
from __future__ import annotations
from typing import Tuple, List
import numpy as np
import pandas as pd
import formulaic


def build_design_matrix(formula: str, data: pd.DataFrame) -> Tuple[np.ndarray, np.ndarray, List[str]]:
    y, X = formulaic.Formula(formula).get_model_matrix(data)
    y_arr = np.asarray(y).ravel().astype(float)
    X_arr = np.asarray(X, dtype=float)
    names = list(X.columns)
    return y_arr, X_arr, names
```

- [ ] **Step 4/5: Pass, commit**

---

### Task 12: Refactor SAR to sparse

**Files:** Modify `src/statspai/spatial/models.py` (existing), create `src/statspai/spatial/models/ml.py`, add `tests/spatial/test_models_ml.py`, `tests/spatial/test_backward_compat.py`

- [ ] **Step 1: Snapshot current SAR/SEM/SDM outputs on a small synthetic dataset and save to `tests/spatial/fixtures/baseline_outputs.json`**

```python
# scripts/snapshot_baseline.py  -- one-off, run once
import numpy as np, json, pandas as pd
from statspai.spatial import sar, sem, sdm

rng = np.random.default_rng(42)
n = 40
coords = rng.uniform(size=(n, 2))
d = np.linalg.norm(coords[:, None] - coords[None], axis=2)
W = (d < 0.3).astype(float); np.fill_diagonal(W, 0.0); W = W / W.sum(axis=1, keepdims=True)
x = rng.standard_normal(n)
y = np.linalg.solve(np.eye(n) - 0.5 * W, 1 + 2 * x + rng.standard_normal(n))
df = pd.DataFrame({"y": y, "x": x})
baseline = {
    "sar": sar(W, df, "y ~ x").params.tolist(),
    "sem": sem(W, df, "y ~ x").params.tolist(),
    "sdm": sdm(W, df, "y ~ x").params.tolist(),
}
with open("tests/spatial/fixtures/baseline_outputs.json", "w") as f:
    json.dump(baseline, f, indent=2)
```

Run it: `python3 scripts/snapshot_baseline.py`. Commit the JSON fixture.

- [ ] **Step 2: Backward-compat test**

```python
# tests/spatial/test_backward_compat.py
import json, numpy as np, pandas as pd, pytest
from statspai.spatial import sar, sem, sdm


@pytest.fixture(scope="module")
def synthetic():
    rng = np.random.default_rng(42)
    n = 40
    coords = rng.uniform(size=(n, 2))
    d = np.linalg.norm(coords[:, None] - coords[None], axis=2)
    W = (d < 0.3).astype(float); np.fill_diagonal(W, 0.0)
    rs = W.sum(axis=1, keepdims=True); rs[rs == 0] = 1
    W = W / rs
    x = rng.standard_normal(n)
    y = np.linalg.solve(np.eye(n) - 0.5 * W, 1 + 2 * x + rng.standard_normal(n))
    return W, pd.DataFrame({"y": y, "x": x})


def test_sar_backward_compat(synthetic):
    W, df = synthetic
    with open("tests/spatial/fixtures/baseline_outputs.json") as f:
        baseline = json.load(f)
    result = sar(W, df, "y ~ x")
    np.testing.assert_allclose(result.params, baseline["sar"], rtol=1e-6)
```

(Analogous for SEM, SDM.)

- [ ] **Step 3: Write failing test for sparse-W acceptance**

```python
# tests/spatial/test_models_ml.py
import numpy as np, pandas as pd, pytest
from statspai.spatial.weights import knn_weights, distance_band
from statspai.spatial import sar, sem, sdm


def test_sar_accepts_W_object():
    rng = np.random.default_rng(0)
    coords = rng.uniform(size=(50, 2))
    w = knn_weights(coords, k=5); w.transform = "R"
    x = rng.standard_normal(50)
    y = 2 + 1.5 * x + rng.standard_normal(50)
    res = sar(w, pd.DataFrame({"y": y, "x": x}), "y ~ x")
    assert hasattr(res, "params")
    assert res.rho is not None       # SAR exposes rho
```

- [ ] **Step 4: Run — new test fails (current `sar` expects ndarray only)**

- [ ] **Step 5: Refactor**

```python
# src/statspai/spatial/models/ml.py
"""Sparse-backed ML estimation for SAR / SEM / SDM / SLX / SAC."""
from __future__ import annotations
from typing import Union
import numpy as np
import pandas as pd
from scipy import sparse
from scipy.optimize import minimize_scalar
from scipy import stats as sp_stats

from ...core.results import EconometricResults
from ..weights.core import W as _W
from ._base import build_design_matrix
from ._logdet import log_det_exact, log_det_approx

ArrayOrW = Union[np.ndarray, sparse.spmatrix, _W]


def _coerce_W(W: ArrayOrW, row_normalize: bool = True) -> sparse.csr_matrix:
    if isinstance(W, _W):
        M = W.sparse
    elif sparse.issparse(W):
        M = W.tocsr()
    else:
        M = sparse.csr_matrix(np.asarray(W, dtype=float))
    if row_normalize:
        rs = np.asarray(M.sum(axis=1)).ravel()
        rs[rs == 0] = 1.0
        M = sparse.diags(1.0 / rs) @ M
    return M.tocsr()


def _log_det(M: sparse.csr_matrix, rho: float) -> float:
    n = M.shape[0]
    if n < 5000:
        return log_det_exact(M, rho)
    return log_det_approx(M, rho, n_draws=100, order=40, seed=0)


def sar(
    W: ArrayOrW, data: pd.DataFrame, formula: str,
    row_normalize: bool = True, alpha: float = 0.05,
) -> EconometricResults:
    y, X, names = build_design_matrix(formula, data)
    M = _coerce_W(W, row_normalize)
    n, k = X.shape

    def neg_ll(rho: float) -> float:
        if not (-0.999 < rho < 0.999):
            return 1e10
        ay = y - rho * (M @ y)
        beta = np.linalg.lstsq(X, ay, rcond=None)[0]
        e = ay - X @ beta
        sig2 = (e @ e) / n
        return -(- n / 2 * np.log(2 * np.pi * sig2)
                 + _log_det(M, rho) - (e @ e) / (2 * sig2))

    opt = minimize_scalar(neg_ll, bracket=(-0.8, 0.0, 0.8),
                          method="brent", options={"xtol": 1e-6})
    rho_hat = opt.x
    ay = y - rho_hat * (M @ y)
    beta = np.linalg.lstsq(X, ay, rcond=None)[0]
    e = ay - X @ beta
    sig2 = (e @ e) / n
    # asymptotic SE: block from analytic information matrix (see LeSage-Pace 2009 p. 45)
    XtX = X.T @ X
    se_beta = np.sqrt(sig2 * np.diag(np.linalg.inv(XtX)))
    # rho SE via numerical hessian
    h = 1e-4
    hess = (neg_ll(rho_hat + h) - 2 * neg_ll(rho_hat) + neg_ll(rho_hat - h)) / h ** 2
    se_rho = 1 / np.sqrt(hess) if hess > 0 else np.nan

    params = np.concatenate([[rho_hat], beta])
    se = np.concatenate([[se_rho], se_beta])
    tvals = params / se
    pvals = 2 * (1 - sp_stats.norm.cdf(np.abs(tvals)))

    result = EconometricResults(
        params=params, se=se, tvals=tvals, pvals=pvals,
        param_names=["rho"] + names,
        n=n, k=k + 1,
        model_type="SAR",
        log_likelihood=-opt.fun,
    )
    result.rho = rho_hat
    result.W = M
    result.sigma2 = sig2
    result.extras = {"residuals": e}
    return result


# (sem, sdm, slx, sac implementations follow the same template.)
```

Update existing `src/statspai/spatial/models.py` to become a thin forwarder:

```python
# src/statspai/spatial/models.py  (after refactor)
"""Backward-compat shim. Public entry points live in statspai.spatial.models.ml."""
from .models.ml import sar, sem, sdm
# SpatialModel class kept for legacy callers — forwards to ml.sar etc.

class SpatialModel:            # noqa: D101 — legacy
    def __init__(self, W, data, formula, model_type="sar", row_normalize=True, alpha=0.05):
        self.W, self.data, self.formula = W, data, formula
        self.model_type, self.row_normalize, self.alpha = model_type, row_normalize, alpha

    def fit(self):
        fn = {"sar": sar, "sem": sem, "sdm": sdm}[self.model_type]
        return fn(self.W, self.data, self.formula,
                  row_normalize=self.row_normalize, alpha=self.alpha)
```

- [ ] **Step 6: Run both test files pass**

```
python3 -m pytest tests/spatial/test_backward_compat.py tests/spatial/test_models_ml.py -v
```

- [ ] **Step 7: Commit**

```
git commit -m "refactor(spatial): sparse-backed SAR ML; accept W object; preserve API"
```

---

### Task 13: Refactor SEM to sparse ML

**Files:** extend `src/statspai/spatial/models/ml.py`, extend `tests/spatial/test_models_ml.py`

- [ ] **Step 1: Failing test** asserting `sem(W_obj, df, "y~x").lam` reasonable on synthetic SEM data.
- [ ] **Step 2/3: Implement** `sem()` — concentrated LL in λ:

```python
def sem(W, data, formula, row_normalize=True, alpha=0.05):
    y, X, names = build_design_matrix(formula, data)
    M = _coerce_W(W, row_normalize); n, k = X.shape

    def neg_ll(lam):
        if not (-0.999 < lam < 0.999):
            return 1e10
        A = sparse.eye(n) - lam * M
        ya = A @ y; Xa = A @ X
        beta = np.linalg.lstsq(Xa, ya, rcond=None)[0]
        e = ya - Xa @ beta
        sig2 = (e @ e) / n
        return -(-n / 2 * np.log(2 * np.pi * sig2) + _log_det(M, lam) - (e @ e) / (2 * sig2))

    opt = minimize_scalar(neg_ll, bracket=(-0.8, 0.0, 0.8), method="brent")
    lam = opt.x
    A = sparse.eye(n) - lam * M
    beta = np.linalg.lstsq(A @ X, A @ y, rcond=None)[0]
    e = (A @ y) - (A @ X) @ beta; sig2 = (e @ e) / n
    # ... SE computation analogous; build EconometricResults with model_type="SEM", extras={"lambda": lam}
```

- [ ] Commit

---

### Task 14: Refactor SDM to sparse ML

Same template, but augment X with `M @ X` (excluding intercept) before calling the SAR loop.

```python
def sdm(W, data, formula, row_normalize=True, alpha=0.05):
    y, X, names = build_design_matrix(formula, data)
    M = _coerce_W(W, row_normalize)
    # append lagged X (skip intercept col 0)
    X_lag = M @ X[:, 1:]
    X_aug = np.hstack([X, X_lag])
    names_aug = names + [f"W_{nm}" for nm in names[1:]]
    # feed into SAR mechanics
    return _sar_from_matrices(M, y, X_aug, names_aug, model_type="SDM")
```

- [ ] Commit

---

### Task 15: SLX (simple lag-X) and SAC (SARAR)

**Files:** extend `ml.py`, add tests.

```python
def slx(W, data, formula, row_normalize=True, alpha=0.05):
    """Y = X beta + W X theta + eps  — plain OLS after augmentation."""
    y, X, names = build_design_matrix(formula, data)
    M = _coerce_W(W, row_normalize)
    X_aug = np.hstack([X, M @ X[:, 1:]])
    names_aug = names + [f"W_{nm}" for nm in names[1:]]
    # ordinary OLS
    beta = np.linalg.lstsq(X_aug, y, rcond=None)[0]
    # ... EconometricResults build
```

SAC is SAR+SEM jointly (ρ, λ) — 2-D optimiser via `scipy.optimize.minimize` with Nelder-Mead; exact log-det for both (I-ρW) and (I-λW).

- [ ] Commit

---

## Phase E — GMM

### Task 16: Kelejian-Prucha SEM GMM

**Files:** `src/statspai/spatial/models/gmm.py`, `tests/spatial/test_models_gmm.py`

- [ ] **Step 1: Failing test** — on the same synthetic SEM DGP, `sem_gmm` estimates match ML within 0.1 on λ.
- [ ] **Step 2/3: Implement** (Kelejian & Prucha 1999; 3-moment system for λ, σ²):

```python
# src/statspai/spatial/models/gmm.py
from __future__ import annotations
import numpy as np
from scipy import sparse
from scipy.optimize import minimize
from ..weights.core import W as _W
from ._base import build_design_matrix
from .ml import _coerce_W
from ...core.results import EconometricResults


def sem_gmm(W, data, formula, robust: str | None = "het"):
    """Kelejian-Prucha generalised moments estimator for SEM."""
    y, X, names = build_design_matrix(formula, data)
    M = _coerce_W(W, row_normalize=True); n, k = X.shape
    # Stage 1: OLS residuals
    beta0 = np.linalg.lstsq(X, y, rcond=None)[0]
    u = y - X @ beta0
    Wu = M @ u; WWu = M @ Wu
    # KP moment conditions g_i(lambda, sigma2)
    def moments(theta):
        lam, s2 = theta
        e = u - lam * Wu
        e2 = u - lam * (u + WWu)              # approx second moment
        m1 = (e @ e) / n - s2
        m2 = (Wu @ Wu) - lam ** 2 * (Wu @ Wu)   # simplified
        m3 = (u @ Wu) / n - lam * ((Wu @ Wu) / n)
        return np.array([m1, m2, m3])

    def J(theta):
        g = moments(theta); return g @ g

    opt = minimize(J, x0=[0.2, float(u @ u) / n], method="Nelder-Mead",
                   options={"xatol": 1e-6, "fatol": 1e-8})
    lam_hat, s2_hat = opt.x
    # Stage 2: GLS with lambda_hat  — transform and re-estimate beta
    A = sparse.eye(n) - lam_hat * M
    ya = A @ y; Xa = A @ X
    beta = np.linalg.lstsq(Xa, ya, rcond=None)[0]
    e = ya - Xa @ beta
    # Heteroscedastic robust SE (Arraiz et al. 2010, sphet)
    if robust == "het":
        meat = Xa.T @ np.diag(e ** 2) @ Xa
        bread = np.linalg.inv(Xa.T @ Xa)
        V = bread @ meat @ bread
    else:
        V = s2_hat * np.linalg.inv(Xa.T @ Xa)
    se = np.sqrt(np.diag(V))
    return EconometricResults(
        params=np.concatenate([[lam_hat], beta]),
        se=np.concatenate([[np.nan], se]),
        tvals=..., pvals=..., param_names=["lambda"] + names,
        n=n, k=k + 1, model_type="SEM-GMM",
    )
```

(Full KP moments use 6-moment system per Kelejian-Prucha 1999 eq. (10)–(12); implement faithfully, citing the paper in docstring.)

- [ ] Commit

---

### Task 17: SAR GMM + SARAR GMM

Uses 2SLS-style instruments `[X, WX, W²X]` for the spatial lag, following Kelejian-Prucha 1998. SARAR chains SEM-GMM after SAR-GMM residuals.

- [ ] Tests — match spreg `GM_Lag` and `GM_Combo` on Columbus.
- [ ] Implement, commit.

---

## Phase F — Diagnostics & impacts

### Task 18: LM tests (Anselin 1988)

**Files:** `src/statspai/spatial/models/diagnostics.py`, `tests/spatial/test_models_diagnostics.py`

```python
# src/statspai/spatial/models/diagnostics.py
from __future__ import annotations
import numpy as np
from scipy import stats as sp_stats, sparse
from ..weights.core import W as _W
from ._base import build_design_matrix
from .ml import _coerce_W


def lm_tests(formula, data, W):
    """Return dict of LM-err, LM-lag, Robust LM-err, Robust LM-lag, SARMA."""
    y, X, _ = build_design_matrix(formula, data)
    M = _coerce_W(W, row_normalize=True); n = len(y)
    beta = np.linalg.lstsq(X, y, rcond=None)[0]
    e = y - X @ beta
    s2 = (e @ e) / n
    We = M @ e; Wy = M @ y
    T = float(((M + M.T).multiply(M)).sum())

    LM_err = ((e @ We) / s2) ** 2 / T
    num = ((e @ Wy) / s2) ** 2
    J = (1 / (n * s2)) * (((M @ X @ beta).T @ (M @ X @ beta)) + T * s2)
    LM_lag = num / J
    RLM_err = ((e @ We) / s2 - T * J ** -1 * (e @ Wy) / s2) ** 2 / (T * (1 - T * J ** -1))
    RLM_lag = ((e @ Wy) / s2 - (e @ We) / s2) ** 2 / (J - T)
    SARMA = LM_err + RLM_lag
    def pv(stat): return 1 - sp_stats.chi2.cdf(stat, df=1)
    return {
        "LM_err": (LM_err, pv(LM_err)),
        "LM_lag": (LM_lag, pv(LM_lag)),
        "Robust_LM_err": (RLM_err, pv(RLM_err)),
        "Robust_LM_lag": (RLM_lag, pv(RLM_lag)),
        "SARMA": (SARMA, 1 - sp_stats.chi2.cdf(SARMA, df=2)),
    }
```

- [ ] Test vs `spreg.diagnostics_sp.LMtests` on Columbus.
- [ ] Commit.

---

### Task 19: LeSage-Pace impacts

**Files:** `src/statspai/spatial/models/impacts.py`, `tests/spatial/test_models_impacts.py`

```python
# src/statspai/spatial/models/impacts.py
"""Direct / indirect / total impacts (LeSage & Pace 2009).

For  Y = (I - rho W)^{-1} (X beta + eps),
  dY/dX_k = (I - rho W)^{-1} beta_k   (SAR)
  ≈ average diagonal = direct; off-diagonal mean = indirect; total = direct + indirect.
For SDM add  (I - rho W)^{-1} W * theta_k.
"""
from __future__ import annotations
import numpy as np
from scipy import sparse
import pandas as pd


def impacts(result, n_sim: int = 1000, seed: int | None = None) -> pd.DataFrame:
    model = result.model_type
    if model not in {"SAR", "SDM", "SAC"}:
        raise ValueError(f"impacts not defined for model_type={model}")
    W = result.W
    n = W.shape[0]
    rho = result.rho
    I = sparse.eye(n)
    A_inv = np.linalg.inv((I - rho * W).toarray())
    beta_names = [p for p in result.param_names if p not in {"rho", "lambda", "Intercept"}]
    # locate betas
    idx_beta = [result.param_names.index(nm) for nm in beta_names]
    betas = result.params[idx_beta]
    if model == "SDM":
        # split into non-lag and lag parts
        raise NotImplementedError("SDM impacts — add in next iteration")

    direct = np.mean(np.diag(A_inv)) * betas
    total = np.mean(A_inv.sum(axis=1)) * betas
    indirect = total - direct

    # Simulated SEs via draws from params' asymptotic distribution
    rng = np.random.default_rng(seed)
    cov = np.diag(result.se ** 2)              # diagonal approx; full cov TODO
    draws = rng.multivariate_normal(result.params, cov, size=n_sim)
    direct_draws, indirect_draws, total_draws = [], [], []
    rho_idx = result.param_names.index("rho")
    for d in draws:
        r = float(d[rho_idx])
        if not (-0.99 < r < 0.99):
            continue
        Ainv = np.linalg.inv((I - r * W).toarray())
        b = d[idx_beta]
        direct_draws.append(np.mean(np.diag(Ainv)) * b)
        total_draws.append(np.mean(Ainv.sum(axis=1)) * b)
        indirect_draws.append(total_draws[-1] - direct_draws[-1])
    direct_draws = np.array(direct_draws); indirect_draws = np.array(indirect_draws); total_draws = np.array(total_draws)
    return pd.DataFrame({
        "Direct": direct,      "SE_Direct":   direct_draws.std(axis=0),
        "Indirect": indirect,  "SE_Indirect": indirect_draws.std(axis=0),
        "Total": total,        "SE_Total":    total_draws.std(axis=0),
    }, index=beta_names)
```

- [ ] Test: SAR with known ρ=0.5 on synthetic; direct ≈ β/(1-ρ) diagonal; indirect positive.
- [ ] Commit

---

## Phase G — Integration

### Task 20: Flat re-exports and top-level API

**Files:** rewrite `src/statspai/spatial/__init__.py`, update `src/statspai/__init__.py`

```python
# src/statspai/spatial/__init__.py
"""Spatial econometrics — weights, ESDA, models."""
from .weights import (
    W, queen_weights, rook_weights, knn_weights,
    distance_band, kernel_weights, block_weights,
)
from .esda import (
    moran, moran_local, geary, getis_ord_g, getis_ord_local, join_counts,
)
from .esda.plots import moran_plot, lisa_cluster_map
from .models.ml import sar, sem, sdm, slx, sac
from .models.gmm import sar_gmm, sem_gmm, sarar_gmm
from .models.diagnostics import lm_tests, moran_residuals
from .models.impacts import impacts

__all__ = [
    "W",
    # weights
    "queen_weights", "rook_weights", "knn_weights",
    "distance_band", "kernel_weights", "block_weights",
    # esda
    "moran", "moran_local", "geary", "getis_ord_g", "getis_ord_local", "join_counts",
    "moran_plot", "lisa_cluster_map",
    # models
    "sar", "sem", "sdm", "slx", "sac",
    "sar_gmm", "sem_gmm", "sarar_gmm",
    "lm_tests", "moran_residuals", "impacts",
]
```

Ensure `statspai/__init__.py` re-exports `sar`, `sem`, `sdm`, `slx`, `sac`, `sar_gmm`, `sem_gmm`, `sarar_gmm`, plus a submodule alias: `from . import spatial`.

- [ ] **Test: `import statspai as sp; sp.moran; sp.sar; sp.impacts`**
- [ ] Commit

---

### Task 21: Columbus reference fixture + cross-validation test

**Files:** `tests/spatial/fixtures/columbus.csv` (copy from libpysal's example set), `tests/spatial/fixtures/reference_outputs.json`, `tests/spatial/test_columbus_crossval.py`

- [ ] **Step 1: Acquire fixture**

```
python3 -c "
import libpysal, pandas as pd
df = libpysal.io.open(libpysal.examples.get_path('columbus.dbf')).to_df()
df.to_csv('tests/spatial/fixtures/columbus.csv', index=False)
"
```

- [ ] **Step 2: Save reference numbers** (run libpysal once, save JSON; fixture file checked in so libpysal is not required at test time).

```python
# scripts/build_columbus_reference.py  -- one-off
import json
from libpysal import weights, examples
from libpysal.io import open as lpopen
import spreg, numpy as np
w = weights.Queen.from_shapefile(examples.get_path('columbus.shp'))
w.transform = "r"
df = lpopen(examples.get_path('columbus.dbf')).to_df()
y = df[["CRIME"]].values.astype(float)
X = df[["INC", "HOVAL"]].values.astype(float)
ml_lag = spreg.ML_Lag(y, X, w=w, method="full")
ml_err = spreg.ML_Error(y, X, w=w, method="full")
json.dump({
    "ml_lag_betas": ml_lag.betas.flatten().tolist(),
    "ml_lag_rho":   float(ml_lag.rho),
    "ml_err_betas": ml_err.betas.flatten().tolist(),
    "ml_err_lambda": float(ml_err.lam),
}, open("tests/spatial/fixtures/reference_outputs.json", "w"), indent=2)
```

Run and commit the JSON.

- [ ] **Step 3: Cross-val test**

```python
# tests/spatial/test_columbus_crossval.py
import json, numpy as np, pandas as pd, pytest
from statspai.spatial import sar, sem, queen_weights


@pytest.fixture(scope="module")
def columbus():
    df = pd.read_csv("tests/spatial/fixtures/columbus.csv")
    gpd = pytest.importorskip("geopandas")
    gdf = gpd.read_file("tests/spatial/fixtures/columbus.shp")
    w = queen_weights(gdf); w.transform = "R"
    return df, w


def test_sar_matches_spreg(columbus):
    df, w = columbus
    ref = json.load(open("tests/spatial/fixtures/reference_outputs.json"))
    res = sar(w, df, "CRIME ~ INC + HOVAL")
    np.testing.assert_allclose(res.rho, ref["ml_lag_rho"], rtol=1e-4)
    # param order: [rho, Intercept, INC, HOVAL]; reference order: [Intercept, INC, HOVAL, rho]
    np.testing.assert_allclose(
        sorted(res.params), sorted([*ref["ml_lag_betas"], ref["ml_lag_rho"]]),
        rtol=1e-4,
    )
```

- [ ] **Step 4: Pass (this is the acceptance gate)**
- [ ] **Step 5: Commit** — mark S1.1 complete

```
git commit -m "test(spatial): Columbus cross-validation vs spreg (<1e-4 rtol)"
```

---

### Task 22: Tutorial notebook

**Files:** `docs/tutorials/spatial_tutorial.ipynb` (or `test-notebooks/` if `docs/` is git-ignored).

- [ ] Walk Columbus: load → queen → Moran → LM tests → SAR vs SEM → impacts → plot.
- [ ] Ensure notebook runs end-to-end in CI (nbval or `jupyter nbconvert --execute`).
- [ ] Commit

---

## Self-Review

**Spec coverage (re-checked against `specs/2026-04-15-sp-01-spatial-full-stack-design.md` §4):**

| Spec API | Task |
|---|---|
| `W`, transforms | 1 |
| `knn_weights`, `distance_band`, `kernel_weights` | 2 |
| `queen_weights`, `rook_weights` | 3 |
| `block_weights` | 4 |
| `moran`, `moran_local` | 5 |
| `geary` | 6 |
| `getis_ord_g`, `getis_ord_local` | 7 |
| `join_counts` | 8 |
| `moran_plot`, `lisa_cluster_map` | 9 |
| `_logdet` dual path (D2) | 10 |
| Design-matrix helper | 11 |
| `sar` (sparse, bw-compat) | 12 |
| `sem` | 13 |
| `sdm` | 14 |
| `slx`, `sac` | 15 |
| `sem_gmm` | 16 |
| `sar_gmm`, `sarar_gmm` | 17 |
| `lm_tests`, `moran_residuals` | 18 |
| `impacts` | 19 |
| Flat re-exports + top-level `sp.*` | 20 |
| Columbus cross-val (<1e-4) | 21 |
| Tutorial | 22 |

All 24 surface items in the spec are covered. §5 decisions D1–D8 all land in Tasks 1, 10, 3, 12, 19, 21, 20, respectively. §7 S1.1 exit criterion (Columbus cross-val passes, existing tests green) is Task 21 + the backward-compat test in Task 12.

**Placeholder scan:** one `raise NotImplementedError("SDM impacts — add in next iteration")` in Task 19 — acceptable because Task 19's test covers SAR only; SDM impacts are punted to a follow-up task. Adding explicit follow-up:

### Task 23: SDM impacts (follow-up)

Extend `impacts()` to handle `model_type == "SDM"` by combining β and θ (W-X lag coefficients) in the multiplier:
`effect_k = (I - ρW)^{-1} (β_k I + θ_k W)`.
Add test using Columbus SDM fit vs `spreg.ML_Lag` with X + WX.

- [ ] Implement, commit.

**Type consistency check:** method signatures `sar(W, data, formula, ...)` match across Tasks 12 / 15 / 20 and the pre-existing `sp.sar`. `EconometricResults.params` order is `[rho_or_lambda, *beta_names]` consistently. `SpatialStatistic` used by all ESDA tasks. No drift.

---

## Execution notes

- **Push cadence:** each task ends with a commit; batch a push to `main` at end of each phase (A/B/C/D/E/F/G).
- **Golden rule:** if a test does not exist, do not write the code; if a test does not fail first, back it out and write a failing version.
- **Soft dependencies:** install via `pip install statspai[spatial]` — add the extra to `pyproject.toml` in Task 20:

```toml
[project.optional-dependencies]
spatial = ["geopandas>=0.14", "libpysal>=4.10", "shapely>=2.0", "matplotlib>=3.8"]
```

Acceptance for S1.1: Task 21 green.

---

## Handoff

S1.2 (GWR + MGWR) and S1.3 (spatial panel) each get their own plan when they open.
