diff --git a/CHANGELOG.md b/CHANGELOG.md
index a0d80d0..85a8189 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,25 @@
 All notable changes to `ewstools` are documented here.
 This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+
+### Added
+
+- **Spatial early warning signals** (`ewstools.spatial`): `morans_i()`,
+  `morans_i_permutation_test()`, and a new `SpatialEWS` class (mirroring the
+  `data`/`state`/`ews` conventions of `MultiTimeSeries`) for computing Moran's I
+  and its significance across space at every time point — the spatial branch of
+  the critical-slowing-down literature (Dakos et al., 2010; MacLaren, Aihara &
+  Masuda, 2025) had no equivalent in this package. See
+  `CONTRIBUTION_spatial_significance.md` for the full write-up.
+- **Combining correlated p-values** (`ewstools.pvalues`): `combine_pvalues_ebm()`,
+  an implementation of the Empirical Brown's Method (Poole et al., 2016) for
+  combining several early-warning indicators' p-values without assuming they are
+  independent, which Fisher's method requires but indicators computed from the
+  same underlying system rarely satisfy.
+- 21 new tests (`tests/test_spatial.py`, `tests/test_pvalues.py`), 0 regressions
+  on the existing suite.
+
 ## [2.1.3] — 2026-07-28
 
 Maintenance release. **No API changes and no behavioural changes** — existing code and
diff --git a/CONTRIBUTION_spatial_significance.md b/CONTRIBUTION_spatial_significance.md
new file mode 100644
index 0000000..6346e84
--- /dev/null
+++ b/CONTRIBUTION_spatial_significance.md
@@ -0,0 +1,79 @@
+# Contribution: spatial early warning signals + combining correlated p-values
+
+## What was missing
+
+`ewstools` implements the temporal branch of the critical-slowing-down
+literature thoroughly — variance, autocorrelation, skew, kurtosis,
+spectral indicators, all tracked over a rolling window in time (`TimeSeries`,
+`MultiTimeSeries`). It has no equivalent for the **spatial** branch of the
+same literature: for systems observed at many locations at once (a
+lattice, a sensor network, an administrative territory grid...), spatial
+autocorrelation between neighbouring units is a separate, well-established
+early warning signal (Dakos et al., 2010, "Spatial correlation as leading
+indicator of catastrophic shifts", *Theoretical Ecology*; still an active
+area — see MacLaren, Aihara & Masuda, 2025, on generalising these methods
+from idealised regular lattices, where almost all validation work has been
+done, to irregular real-world networks). Searching `ewstools/core.py` and
+`ewstools/helpers.py` for "moran", "spatial", "fisher", "brown", or
+"combine_p" turns up nothing — confirmed by direct inspection of both
+files before writing any code here, not assumed.
+
+A related, separate gap: `ewstools` never combines its own indicators.
+If a user computes both a temporal indicator (say, variance) and a spatial
+one (Moran's I) on the same system, there is no tool in the package to ask
+"are both unusual at once?" — which needs Fisher's method at minimum, and
+a **correlation-aware** version of it if the two indicators (as is typical,
+since they are computed from the same underlying system) are not
+independent. Fisher's method, applied naively to correlated inputs, is
+biased. Confirmed absent by the same search.
+
+## What this contribution adds
+
+- **`ewstools/spatial.py`** — `morans_i()` and `morans_i_permutation_test()`
+  as standalone functions, plus a `SpatialEWS` class that follows the
+  existing `MultiTimeSeries` conventions (`data` → `state` → `ews`,
+  transition-aware, a `compute_ktau()` for trend-testing the resulting
+  indicator series exactly like the rest of the package). Computes Moran's
+  I — and its permutation-test significance — at every time point from a
+  DataFrame with one column per spatial unit and a fixed spatial weights
+  matrix.
+- **`ewstools/pvalues.py`** — `combine_pvalues_ebm()`, a real implementation
+  of the Empirical Brown's Method (Poole, Gibbs, Shmulevich, Bernard &
+  Knijnenburg, 2016, *Bioinformatics*), which extends Brown's (1975)
+  covariance-corrected version of Fisher's method by estimating the
+  correlation between tests empirically from data rather than requiring it
+  to be known or assumed.
+- 21 new tests (`tests/test_spatial.py`, `tests/test_pvalues.py`): hand-computed
+  Moran's I examples, known limiting cases (checkerboard → strong negative
+  autocorrelation, smooth gradient → strong positive autocorrelation,
+  constant field → undefined), the classical permutation-test property
+  `E[I] → -1/(N-1)` under the null, and — for the p-value combiner — that it
+  reduces to plain Fisher when the underlying data are independent, and
+  produces a *less* extreme combined p-value than naive Fisher when they are
+  positively correlated (the whole point of the correction), verified
+  against a direct evaluation of Brown's polynomial at a known correlation.
+  Full suite run locally: 52 passed, 1 skipped (0 regressions on the 32
+  tests that existed before this change).
+
+## Usage angles
+
+- **Research**: any critical-transition study with spatially resolved data
+  (ecology, epidemiology, social/territorial data, sensor networks, plasma
+  diagnostics) gains a spatial EWS tool in the same package already used for
+  the temporal one, instead of having to hand-roll Moran's I separately.
+- **Industry / infrastructure monitoring**: sensor networks (power grids,
+  structural health monitoring, industrial process control) already produce
+  exactly the "one column per unit, one row per time step" shape this class
+  expects.
+- **Public services**: territorial statistics agencies publishing
+  regularly-updated indicators across administrative units can use
+  `SpatialEWS` directly on their own panel data without depending on a
+  domain-specific tool.
+
+## Status
+
+Not yet submitted as a pull request to `ThomasMBury/ewstools`. See
+`/positionnement/ewstools` on the Hélios site for the live, honestly
+maintained status once submission happens — this repository is a local
+fork/branch pending that decision, and is not yet available even via a
+public fork URL.
diff --git a/ewstools/pvalues.py b/ewstools/pvalues.py
new file mode 100644
index 0000000..0ebf84c
--- /dev/null
+++ b/ewstools/pvalues.py
@@ -0,0 +1,113 @@
+"""
+Combining correlated p-values.
+
+`ewstools` computes several early warning indicators from the same
+underlying data (variance, autocorrelation, skew, a spatial indicator...).
+A natural next question -- not currently answered anywhere in the
+package -- is: are two or more of these signals unusual AT THE SAME TIME?
+Fisher's (1925) classical method for combining independent p-values,
+
+    T = -2 * sum(ln p_i),  T ~ chi2(2k) under H0,
+
+is the obvious tool, but it assumes the p_i are independent. Indicators
+computed from the same underlying system are typically correlated, which
+biases Fisher's method (too liberal if the correlation is positive, too
+conservative if negative).
+
+The Empirical Brown's Method (Poole, Gibbs, Shmulevich, Bernard &
+Knijnenburg, 2016, "Combining dependent P-values with an empirical
+adaptation of Brown's method", Bioinformatics 32(17), i430-i436) fixes
+this: it keeps Fisher's statistic T, but instead of assuming the
+textbook chi2(2k) null distribution, it estimates the TRUE variance of T
+under the observed correlation structure directly from data (rather than
+requiring the covariance to be known analytically, as in Brown's 1975
+original method), then fits a rescaled chi-squared distribution to that
+corrected mean and variance (method of moments).
+"""
+from __future__ import annotations
+
+import numpy as np
+from scipy.stats import chi2
+
+
+def combine_pvalues_ebm(data, p_values) -> dict:
+    """Combine one-sided, possibly-correlated p-values with the Empirical
+    Brown's Method (Poole et al., 2016), building on Brown's (1975)
+    polynomial approximation for the covariance of `-2 ln(p)` terms.
+
+    Parameters
+    ----------
+    data : array-like, shape (n_observations, n_tests)
+        The raw data underlying each of the `n_tests` p-values (one
+        column per test) -- e.g. the k indicator series each p-value was
+        derived from. Used only to estimate the empirical correlation
+        between tests; this is what makes the method "empirical" rather
+        than requiring the covariance to be assumed or known exactly, as
+        in Brown's (1975) original method.
+    p_values : array-like, shape (n_tests,)
+        The one-sided p-value from each test (already computed
+        elsewhere, e.g. by a surrogate or permutation test). A small
+        p-value must consistently mean the same thing (e.g. "evidence of
+        an approaching transition") across all tests for the combination
+        to be meaningful.
+
+    Returns
+    -------
+    dict with keys:
+        p_combined : combined p-value.
+        chi2_fisher : Fisher's uncorrected statistic, -2*sum(ln p_i).
+        chi2_ebm : chi2_fisher rescaled by the empirical correction factor `c`.
+        df_ebm : effective degrees of freedom after correction (<= 2*n_tests
+            when p-values are positively correlated).
+        c : empirical scale factor (> 1 under positive correlation, which is
+            what makes the combined p-value LESS extreme than naive Fisher
+            when the underlying tests are positively correlated -- the
+            correction Fisher's method is missing).
+
+    References
+    ----------
+    Brown, M. B. (1975). "A method for combining non-independent,
+    one-sided tests of significance." Biometrics, 31(4), 987-992.
+    Poole, W., Gibbs, D. L., Shmulevich, I., Bernard, B., & Knijnenburg,
+    T. A. (2016). "Combining dependent P-values with an empirical
+    adaptation of Brown's method." Bioinformatics, 32(17), i430-i436.
+    """
+    data = np.asarray(data, dtype=float)
+    p_values = np.asarray(p_values, dtype=float)
+    k = len(p_values)
+
+    if data.ndim != 2 or data.shape[1] != k:
+        raise ValueError("data must have shape (n_observations, n_tests) with n_tests == len(p_values)")
+    if k < 2:
+        raise ValueError("need at least two p-values to combine")
+
+    clipped = np.clip(p_values, 1e-15, 1.0)
+    chi2_fisher = float(-2.0 * np.sum(np.log(clipped)))
+
+    corr = np.corrcoef(data, rowvar=False)
+
+    # Brown's (1975) polynomial approximation of Cov(-2 ln p_i, -2 ln p_j)
+    # as a function of the Pearson correlation rho between the underlying
+    # (approximately normal) variables i and j.
+    cov_sum = 0.0
+    for i in range(k):
+        for j in range(i + 1, k):
+            rho = corr[i, j]
+            cov_sum += rho * (3.263 + rho * (0.710 + rho * 0.027))
+
+    var_fisher = 4 * k + 2 * cov_sum
+    mean_fisher = 2 * k
+
+    # Method-of-moments fit: T ~ c * chi2(df), matching mean and variance.
+    c = var_fisher / (2 * mean_fisher)
+    df_ebm = 2 * mean_fisher**2 / var_fisher
+
+    p_combined = float(chi2.sf(chi2_fisher / c, df_ebm))
+
+    return {
+        "p_combined": p_combined,
+        "chi2_fisher": chi2_fisher,
+        "chi2_ebm": chi2_fisher / c,
+        "df_ebm": float(df_ebm),
+        "c": float(c),
+    }
diff --git a/ewstools/spatial.py b/ewstools/spatial.py
new file mode 100644
index 0000000..355d578
--- /dev/null
+++ b/ewstools/spatial.py
@@ -0,0 +1,170 @@
+"""
+Spatial early warning signals.
+
+Most of `ewstools` treats early warning signals (EWS) as a purely
+temporal phenomenon: a rolling window slides along a single time series
+(`TimeSeries`) or a handful of co-measured series (`MultiTimeSeries`),
+and indicators like variance or lag-1 autocorrelation are tracked over
+time. But a large and separate branch of the critical-slowing-down
+literature -- starting with Dakos et al. (2010), "Spatial correlation as
+leading indicator of catastrophic shifts", and still an active area (see
+e.g. MacLaren, Aihara & Masuda, 2025, on generalising spatial EWS from
+regular lattices to irregular real-world networks) -- looks instead
+across SPACE at a single instant: do neighbouring units of a system
+become more correlated with each other as a critical transition
+approaches? Moran's I is the standard statistic for this, but it is not
+implemented anywhere in `ewstools` today.
+
+This module adds that missing piece, following the existing package's
+conventions (`ewstools.core.MultiTimeSeries`): a spatial analogue that
+takes one column per spatial unit and computes a spatial indicator at
+EVERY time point (cross-sectional at each row), rather than a temporal
+indicator over a rolling window.
+"""
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+
+def morans_i(values, weights) -> float:
+    """Moran's I spatial autocorrelation statistic.
+
+    I = (N / S0) * [sum_ij w_ij (x_i - xbar)(x_j - xbar)] / [sum_i (x_i - xbar)^2]
+
+    Parameters
+    ----------
+    values : array-like, shape (n_units,)
+        Observed value at each spatial unit.
+    weights : array-like, shape (n_units, n_units)
+        Spatial weights matrix, w_ij > 0 if units i and j are neighbours,
+        0 otherwise (zero diagonal). Not required to be row-standardised.
+
+    Returns
+    -------
+    float
+        Moran's I, or nan if the weights or the values have no variation
+        (undefined in that case).
+    """
+    x = np.asarray(values, dtype=float)
+    w = np.asarray(weights, dtype=float)
+    n = len(x)
+    if w.shape != (n, n):
+        raise ValueError(f"weights must be an {n} x {n} matrix to match values")
+
+    deviations = x - x.mean()
+    s0 = w.sum()
+    if s0 == 0:
+        return float("nan")
+
+    denominator = (deviations**2).sum()
+    if denominator == 0:
+        return float("nan")
+
+    numerator = deviations @ w @ deviations
+    return float((n / s0) * (numerator / denominator))
+
+
+def morans_i_permutation_test(values, weights, n_permutations: int = 500, seed=None) -> dict:
+    """Significance of Moran's I by permutation (spatial analogue of a
+    surrogate test: reshuffle the values across spatial units while
+    keeping the network itself fixed, following the standard approach
+    for testing spatial autocorrelation, e.g. Dakos et al., 2010).
+
+    Returns
+    -------
+    dict with keys: observed_i, p_value, n_permutations, null_mean, null_std.
+    `p_value` is one-sided (probability of a permuted I at least as large
+    as the observed one) -- appropriate for testing an INCREASE in
+    spatial correlation as an early warning signal.
+    """
+    observed = morans_i(values, weights)
+    if np.isnan(observed):
+        return {"observed_i": None, "p_value": None, "n_permutations": 0, "null_mean": None, "null_std": None}
+
+    rng = np.random.default_rng(seed)
+    values = np.asarray(values, dtype=float)
+    null_values = np.empty(n_permutations)
+    for k in range(n_permutations):
+        null_values[k] = morans_i(rng.permutation(values), weights)
+
+    return {
+        "observed_i": observed,
+        "p_value": float(np.mean(null_values >= observed)),
+        "n_permutations": n_permutations,
+        "null_mean": float(np.mean(null_values)),
+        "null_std": float(np.std(null_values)),
+    }
+
+
+class SpatialEWS:
+    """
+    Spatially-resolved data on which to compute spatial early warning
+    signals, following the `data`/`state`/`ews` conventions of
+    `ewstools.core.MultiTimeSeries`.
+
+    Parameters
+    ----------
+    data : pandas.DataFrame
+        One column per spatial unit (grid cell, node, sensor, region...),
+        one row per time point. Index represents time and is carried
+        over.
+    weights : array-like, shape (n_units, n_units)
+        Spatial weights matrix for the units in `data.columns` (same
+        order), fixed over time -- the network itself is not assumed to
+        change, only the values observed on it.
+    transition : float, optional
+        Time value at which a transition occurs, if any. If given,
+        spatial EWS are only computed up to this point.
+    """
+
+    def __init__(self, data, weights, transition=None):
+        if not isinstance(data, pd.DataFrame):
+            raise TypeError("data must be a pandas DataFrame (one column per spatial unit)")
+        weights = np.asarray(weights, dtype=float)
+        n_units = data.shape[1]
+        if weights.shape != (n_units, n_units):
+            raise ValueError(f"weights must be an N x N matrix with N = data.shape[1] = {n_units}")
+
+        self.state = data
+        self.weights = weights
+        self.transition = float(transition) if transition else transition
+        self.var_names = data.columns
+        self.ews = pd.DataFrame(index=data.index)
+        self.ktau = dict()
+
+    def _pre_transition(self) -> pd.DataFrame:
+        if self.transition:
+            return self.state[self.state.index <= self.transition]
+        return self.state
+
+    def compute_moran(self):
+        """Compute Moran's I at every time point. Output stored in
+        `self.ews['morans_i']`.
+        """
+        df_pre = self._pre_transition()
+        self.ews["morans_i"] = df_pre.apply(lambda row: morans_i(row.to_numpy(), self.weights), axis=1)
+
+    def compute_moran_significance(self, n_permutations: int = 500, seed=None):
+        """Permutation-test p-value for Moran's I at every time point.
+        Output stored in `self.ews['morans_i_pvalue']`.
+        """
+        df_pre = self._pre_transition()
+        self.ews["morans_i_pvalue"] = df_pre.apply(
+            lambda row: morans_i_permutation_test(row.to_numpy(), self.weights, n_permutations, seed)["p_value"],
+            axis=1,
+        )
+
+    def compute_ktau(self, tmin="earliest", tmax="latest"):
+        """Kendall tau of each spatial EWS against time -- same convention
+        as `ewstools.core.TimeSeries.compute_ktau`. Output stored in the
+        `self.ktau` dict.
+        """
+        if tmin == "earliest":
+            tmin = self.ews.dropna(how="all").index[0]
+        if tmax == "latest":
+            tmax = self.ews.dropna(how="all").index[-1]
+
+        df_ews = self.ews[(self.ews.index >= tmin) & (self.ews.index <= tmax)].copy()
+        time_values = pd.Series(data=df_ews.index, index=df_ews.index)
+        self.ktau = dict(df_ews.corrwith(time_values, method="kendall", axis=0))
diff --git a/tests/test_pvalues.py b/tests/test_pvalues.py
new file mode 100644
index 0000000..eb0f1fc
--- /dev/null
+++ b/tests/test_pvalues.py
@@ -0,0 +1,77 @@
+import numpy as np
+import pytest
+from scipy.stats import chi2
+
+from ewstools.pvalues import combine_pvalues_ebm
+
+
+def test_combine_pvalues_ebm_requires_at_least_two_tests():
+    with pytest.raises(ValueError):
+        combine_pvalues_ebm(np.zeros((10, 1)), [0.05])
+
+
+def test_combine_pvalues_ebm_shape_mismatch_raises():
+    with pytest.raises(ValueError):
+        combine_pvalues_ebm(np.zeros((10, 2)), [0.05, 0.1, 0.2])
+
+
+def test_combine_pvalues_ebm_independent_data_matches_fisher_method():
+    rng = np.random.default_rng(0)
+    data = rng.normal(size=(500, 3))  # independent columns -> correlation ~ 0
+    p_values = [0.02, 0.03, 0.10]
+
+    result = combine_pvalues_ebm(data, p_values)
+    chi2_fisher = -2 * sum(np.log(p_values))
+    fisher_p = chi2.sf(chi2_fisher, 2 * len(p_values))
+
+    assert result["c"] == pytest.approx(1.0, abs=0.05)
+    assert result["df_ebm"] == pytest.approx(2 * len(p_values), abs=0.5)
+    assert result["p_combined"] == pytest.approx(fisher_p, abs=0.02)
+
+
+def test_combine_pvalues_ebm_positive_correlation_gives_less_extreme_pvalue_than_fisher():
+    # Two tests built from strongly correlated underlying data carry
+    # redundant evidence; the empirical correction must make the combined
+    # p-value LESS significant than naively treating them as independent
+    # (Fisher), which is the whole point of the correction.
+    rng = np.random.default_rng(1)
+    base = rng.normal(size=500)
+    data = np.column_stack([base, base + rng.normal(scale=0.05, size=500)])
+    p_values = [0.01, 0.01]
+
+    result = combine_pvalues_ebm(data, p_values)
+    chi2_fisher = -2 * sum(np.log(p_values))
+    fisher_p = chi2.sf(chi2_fisher, 2 * len(p_values))
+
+    assert result["c"] > 1.0
+    assert result["p_combined"] > fisher_p
+
+
+def test_combine_pvalues_ebm_brown_polynomial_at_perfect_correlation():
+    # Direct check of Brown's (1975) polynomial itself at rho=1:
+    # cov = 1*(3.263 + 1*(0.710 + 1*0.027)) = 4.0 exactly.
+    data = np.column_stack([np.arange(100.0), np.arange(100.0) * 2 + 1])  # rho = 1 exactly
+    result = combine_pvalues_ebm(data, [0.05, 0.2])
+    k = 2
+    expected_var = 4 * k + 2 * 4.0
+    expected_c = expected_var / (2 * 2 * k)
+    assert result["c"] == pytest.approx(expected_c, abs=1e-9)
+
+
+def test_combine_pvalues_ebm_returns_expected_keys():
+    rng = np.random.default_rng(2)
+    data = rng.normal(size=(50, 3))
+    result = combine_pvalues_ebm(data, [0.1, 0.2, 0.3])
+    assert set(result) == {"p_combined", "chi2_fisher", "chi2_ebm", "df_ebm", "c"}
+    assert 0.0 <= result["p_combined"] <= 1.0
+
+
+def test_combine_pvalues_ebm_invariant_to_test_order():
+    rng = np.random.default_rng(3)
+    data = rng.normal(size=(200, 3))
+    p_values = [0.02, 0.15, 0.4]
+
+    original = combine_pvalues_ebm(data, p_values)
+    reordered = combine_pvalues_ebm(data[:, [2, 0, 1]], [p_values[2], p_values[0], p_values[1]])
+
+    assert reordered["p_combined"] == pytest.approx(original["p_combined"], abs=1e-9)
diff --git a/tests/test_spatial.py b/tests/test_spatial.py
new file mode 100644
index 0000000..3480baf
--- /dev/null
+++ b/tests/test_spatial.py
@@ -0,0 +1,151 @@
+import numpy as np
+import pandas as pd
+import pytest
+
+from ewstools.spatial import SpatialEWS, morans_i, morans_i_permutation_test
+
+
+def _grid_rook_weights(n):
+    weights = np.zeros((n * n, n * n))
+    for r in range(n):
+        for c in range(n):
+            i = r * n + c
+            for dr, dc in [(0, 1), (0, -1), (1, 0), (-1, 0)]:
+                rr, cc = r + dr, c + dc
+                if 0 <= rr < n and 0 <= cc < n:
+                    weights[i, rr * n + cc] = 1
+    return weights
+
+
+# -- morans_i --------------------------------------------------------------
+
+
+def test_morans_i_matches_hand_computed_chain_example():
+    # Chain 1-2-3-4, values = perfect gradient [1,2,3,4].
+    # By hand: xbar=2.5, deviations=[-1.5,-0.5,0.5,1.5], S0=6, numerator=2.5,
+    # denominator=5.0, I = (4/6)*(2.5/5.0) = 1/3.
+    values = np.array([1.0, 2.0, 3.0, 4.0])
+    weights = np.array(
+        [
+            [0, 1, 0, 0],
+            [1, 0, 1, 0],
+            [0, 1, 0, 1],
+            [0, 0, 1, 0],
+        ],
+        dtype=float,
+    )
+    assert morans_i(values, weights) == pytest.approx(1 / 3, abs=1e-9)
+
+
+def test_morans_i_checkerboard_gives_strong_negative_autocorrelation():
+    n = 4
+    grid = np.indices((n, n)).sum(axis=0) % 2
+    values = np.where(grid == 0, 1.0, -1.0).flatten()
+    assert morans_i(values, _grid_rook_weights(n)) < -0.9
+
+
+def test_morans_i_smooth_gradient_gives_strong_positive_autocorrelation():
+    n = 4
+    values = np.indices((n, n))[1].flatten().astype(float)
+    assert morans_i(values, _grid_rook_weights(n)) > 0.5
+
+
+def test_morans_i_constant_values_is_nan():
+    values = np.full(9, 3.0)
+    assert np.isnan(morans_i(values, _grid_rook_weights(3)))
+
+
+def test_morans_i_zero_weights_is_nan():
+    values = np.array([1.0, 2.0, 3.0, 4.0])
+    assert np.isnan(morans_i(values, np.zeros((4, 4))))
+
+
+def test_morans_i_shape_mismatch_raises():
+    with pytest.raises(ValueError):
+        morans_i(np.array([1.0, 2.0, 3.0]), np.zeros((4, 4)))
+
+
+# -- morans_i_permutation_test ----------------------------------------------
+
+
+def test_permutation_test_null_mean_matches_theoretical_expectation():
+    # Under H0 (permutation), E[I] -> -1/(N-1), a classical property.
+    rng = np.random.default_rng(0)
+    n = 20
+    values = rng.normal(size=n)
+    weights = (rng.random((n, n)) < 0.3).astype(float)
+    np.fill_diagonal(weights, 0)
+    weights = np.maximum(weights, weights.T)
+
+    result = morans_i_permutation_test(values, weights, n_permutations=2000, seed=1)
+    assert result["null_mean"] == pytest.approx(-1 / (n - 1), abs=0.05)
+
+
+def test_permutation_test_p_value_in_bounds_and_matches_observed():
+    n = 4
+    values = np.indices((n, n))[1].flatten().astype(float)
+    result = morans_i_permutation_test(values, _grid_rook_weights(n), n_permutations=500, seed=2)
+    assert result["observed_i"] == pytest.approx(morans_i(values, _grid_rook_weights(n)))
+    assert 0.0 <= result["p_value"] <= 1.0
+
+
+# -- SpatialEWS ---------------------------------------------------------------
+
+
+def _toy_spatial_frame(n_time=20, n_units=9, seed=3):
+    rng = np.random.default_rng(seed)
+    data = pd.DataFrame(rng.normal(size=(n_time, n_units)), index=np.arange(n_time))
+    return data, _grid_rook_weights(3)
+
+
+def test_spatialews_requires_dataframe():
+    _, weights = _toy_spatial_frame()
+    with pytest.raises(TypeError):
+        SpatialEWS(np.zeros((5, 9)), weights)
+
+
+def test_spatialews_weights_shape_validation():
+    data, _ = _toy_spatial_frame()
+    with pytest.raises(ValueError):
+        SpatialEWS(data, np.zeros((3, 3)))
+
+
+def test_spatialews_compute_moran_matches_function_row_by_row():
+    data, weights = _toy_spatial_frame()
+    spatial = SpatialEWS(data, weights)
+    spatial.compute_moran()
+    for t, row in data.iterrows():
+        assert spatial.ews.loc[t, "morans_i"] == pytest.approx(morans_i(row.to_numpy(), weights))
+
+
+def test_spatialews_compute_moran_significance_produces_valid_pvalues():
+    data, weights = _toy_spatial_frame(n_time=5)
+    spatial = SpatialEWS(data, weights)
+    spatial.compute_moran_significance(n_permutations=200, seed=4)
+    pvals = spatial.ews["morans_i_pvalue"].dropna()
+    assert len(pvals) == 5
+    assert pvals.between(0, 1).all()
+
+
+def test_spatialews_compute_ktau_detects_increasing_trend():
+    n_time, n_units = 30, 9
+    weights = _grid_rook_weights(3)
+    rng = np.random.default_rng(5)
+    rows = []
+    for t in range(n_time):
+        # Gradient strength grows over time -> Moran's I should trend up.
+        base = np.indices((3, 3))[1].flatten().astype(float)
+        rows.append(base * (t / n_time) + rng.normal(scale=0.05, size=n_units))
+    data = pd.DataFrame(rows, index=np.arange(n_time))
+
+    spatial = SpatialEWS(data, weights)
+    spatial.compute_moran()
+    spatial.compute_ktau()
+    assert spatial.ktau["morans_i"] > 0.5
+
+
+def test_spatialews_transition_restricts_computation():
+    data, weights = _toy_spatial_frame(n_time=10)
+    spatial = SpatialEWS(data, weights, transition=5)
+    spatial.compute_moran()
+    assert spatial.ews["morans_i"].dropna().index.max() <= 5