diff --git a/swarmmetrics.py b/swarmmetrics.py
index 7c585bf..8079247 100644
--- a/swarmmetrics.py
+++ b/swarmmetrics.py
@@ -284,6 +284,71 @@ def score_echo(messages: list[dict], concept_window_seconds: float = 604800,
     return node_echo
 
 
+def shuffle_test(messages: list[dict], n_shuffles: int = 100,
+                 use_idf: bool = True) -> dict:
+    """
+    Shuffle test for causal echo (Anagnostopoulos, Kumar & Mahdian 2008).
+    Permutes timestamps while holding network fixed.
+    If echo scores survive shuffling, they measured correlation (shared
+    environment) not diffusion (causal influence).
+
+    Returns dict with:
+      - observed: {node: echo_coeff} from actual data
+      - mean_shuffled: {node: mean echo across shuffles}
+      - z_scores: {node: (observed - mean_shuffled) / std_shuffled}
+      - significant: {node: bool} where |z| > 2.0
+
+    Cited: ColonistOne (Colony, 2026-07-28) pointed to this method.
+    """
+    import random
+
+    # Observed echo
+    observed = score_echo(messages, use_idf=use_idf)
+
+    # Run shuffles
+    shuffled_scores = defaultdict(list)
+    timestamps = [m.get("timestamp", 0) for m in messages]
+
+    for _ in range(n_shuffles):
+        # Permute timestamps, keep everything else
+        perm = timestamps.copy()
+        random.shuffle(perm)
+        shuffled_msgs = []
+        for i, m in enumerate(messages):
+            sm = m.copy()
+            sm["timestamp"] = perm[i]
+            shuffled_msgs.append(sm)
+
+        shuffled_echo = score_echo(shuffled_msgs, use_idf=use_idf)
+        all_nodes = set(observed) | set(shuffled_echo)
+        for node in all_nodes:
+            shuffled_scores[node].append(shuffled_echo.get(node, 0.0))
+
+    # Compute z-scores
+    result = {
+        "observed": observed,
+        "mean_shuffled": {},
+        "z_scores": {},
+        "significant": {}
+    }
+
+    for node in set(observed) | set(shuffled_scores):
+        scores = shuffled_scores.get(node, [0.0])
+        mean_s = sum(scores) / len(scores)
+        std_s = (sum((s - mean_s)**2 for s in scores) / len(scores)) ** 0.5
+        obs = observed.get(node, 0.0)
+
+        result["mean_shuffled"][node] = round(mean_s, 4)
+        if std_s > 0:
+            z = (obs - mean_s) / std_s
+        else:
+            z = 0.0 if obs == mean_s else float('inf')
+        result["z_scores"][node] = round(z, 4)
+        result["significant"][node] = abs(z) > 2.0
+
+    return result
+
+
 def detect_shadows(messages: list[dict], now: Optional[float] = None,
                    silence_days: float = 7.0) -> dict:
     """