---
name: mortgage-lock-in-and-listings
description: >
  Quantifies how much of the post-2022 decline in U.S. new home listings is
  explained by mortgage rate lock-in. Downloads Freddie Mac PMMS weekly
  30-year fixed mortgage rates and the Realtor.com Inventory Core Metrics
  Metro History monthly panel. Constructs a national "rate spread" as
  current PMMS minus an exponentially-weighted 10-year trailing effective
  rate. Fits a metro-FE + calendar-month-of-year-FE panel of log new
  listings on the rate spread, its interaction with metro-level lock-in
  exposure, and a linear time trend, with a 2,000-iteration
  seasonal-block permutation null, metro-cluster bootstrap 95%
  confidence intervals, six sensitivity analyses (effective-rate window,
  spread lag, pre/post-COVID, top-metro subsample, alternate PMMS
  vintage, spread=0 placebo), and a partial-equilibrium counterfactual
  decomposition. Python 3.8+ standard library only; data files
  SHA256-pinned.
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "real-estate", "housing", "mortgage", "lock-in", "panel", "permutation-test", "bootstrap", "realtor-com", "freddie-mac"]
python_version: ">=3.8"
dependencies: []
---

# How Much of the Post-2022 Listings Decline Is Explained by Mortgage Rate Lock-In?

## When to Use This Skill

Use this skill when an agent needs to quantify **how much of a regional
outcome decline is explained by an aggregate shock × cross-sectional
exposure panel**, and must distinguish a genuine causal-like contribution
from an artifact of seasonality, metro composition, or a secular trend.
Concretely, use it when you need to investigate whether a commonly-cited
narrative ("rate lock-in is the dominant driver of the 2022+ housing
listings drought") survives a panel fixed-effects decomposition with a
proper null model (seasonal-block permutation), metro-cluster bootstrap
confidence intervals, a spread=0 placebo, and sensitivity to the
effective-rate window, alternative PMMS vintages, and pre/post-COVID
stability.

### Preconditions

- **Python version:** 3.8+ standard library only (no pip installs, no numpy/scipy/pandas).
- **Network:** Internet access to `www.freddiemac.com` and
  `econdata.s3-us-west-2.amazonaws.com` required on first run; responses are
  cached locally with SHA256 integrity checks, so reruns are offline.
- **Runtime:** 15–30 minutes on first run (download + 2,000 permutations +
  2,000 metro-clustered bootstraps + six sensitivity analyses); 10–20
  minutes on rerun from cache. Dominated by the 4,000 panel refits.

## Adaptation Guidance

To apply this analysis to a different "exposure times aggregate shock"
panel question (for example: "how much of state-level vehicle-sales decline
is explained by gasoline price shocks?"):

- **Change `PMMS_URL` and `REALTOR_METRO_URL`** in the DOMAIN CONFIGURATION
  block. `load_data()` parses these two CSVs; swap in your aggregate shock
  CSV (column `date`, value column configured via `SHOCK_VALUE_COL`) and
  your per-unit panel CSV (columns `month_date_yyyymm`, unit id, outcome).
- **Change `OUTCOME_COL` and `UNIT_ID_COL`** to your panel's outcome (here
  `new_listing_count`) and unit (here `cbsa_code`). The rest of the pipeline
  is agnostic to the outcome name.
- **Change `EFFECTIVE_RATE_WINDOW_MONTHS` and `EFFECTIVE_RATE_HALFLIFE_MONTHS`**
  if your "stock price" memory is faster/slower than a decade (e.g., use
  24 months for gasoline contracts).
- **Change `EXPOSURE_PROXY_COL`** — here the primary proxy is the metro's
  pre-2020 median listing price, ranking metros by lock-in exposure. Swap
  for any pre-treatment baseline that plausibly indexes exposure intensity.
- **Do NOT change** `rate_spread()`, `fit_panel_within()`,
  `two_way_cal_within_transform()`, `permutation_p_seasonal_block()`,
  `bootstrap_ci_cluster()`, or `counterfactual_decomposition()` — these
  are domain-agnostic and implement the statistical pipeline (effective
  rate, metro-FE + calendar-month-FE two-way within transform,
  seasonal-block label permutation, cluster bootstrap, counterfactual
  accounting).
- **Do NOT change** the cache/SHA256 layer in `http_get_cached()` — it is
  the reproducibility anchor.

## Research Question

Between January 2022 and March 2026 the Realtor.com monthly metro panel shows
a persistent decline in new home listings relative to the 2016–2021 baseline.
The "rate lock-in" narrative attributes this to existing homeowners being
unwilling to sell and give up mortgages originated at 2020–2021 rates of
~3 percent. Three alternative explanations compete: (i) remote work reducing
selling need, (ii) demographic aging reducing life-stage selling, and (iii)
a secular trend that began before 2022. This skill quantifies the lock-in
share by estimating a panel semi-elasticity and running a counterfactual
decomposition; it does not prove causation.

## Methodological Hook

Rate-spread panel with metro-level exposure heterogeneity. Most popular-press
analyses compare 2021 to 2023 listings totals and attribute 100 percent of
the gap to rate lock-in. We instead construct a monthly national rate spread
(current PMMS minus an exponentially-weighted 10-year trailing effective
rate) and exploit cross-sectional heterogeneity in metro-level lock-in
exposure (proxied by pre-2020 median listing price). The interaction
`spread × exposure` identifies the lock-in channel off high-exposure metros
responding more strongly to the same national shock than low-exposure
metros, conditional on metro FE and calendar-month-of-year FE (both
absorbed via 2-way demeaning). **We intentionally use calendar-month FE,
not full month-of-sample FE**, because the spread varies across years
within a calendar month; using full time FE would annihilate it. The
`exposure` main effect is absorbed by metro FE and dropped from the X
vector. The null is a seasonal-block permutation of the spread series
within calendar month.

## Null Model

For each calendar month m ∈ {1,…,12}, collect the set of rate-spread values
observed in that calendar month across years. Under H0 these values are
exchangeable within m (they contribute only seasonally). Under H1 they are
ordered by the true year index. We permute the spread observations within
each calendar-month stratum 2,000 times, recompute the panel coefficient on
`spread × exposure`, and count the fraction of permutations with absolute
coefficient at least as large as observed. This controls for seasonality,
long-run levels, and metro heterogeneity — only the year-to-year ordering
of the shock is perturbed.

## Controls and Comparators

Three independent comparators guard against false-positive signals; **any
single one failing should reduce confidence in the lock-in narrative**:

1. **Seasonal-block permutation null (N=2,000 shuffles).** Under H0 the
   year-to-year ordering of the national spread is exchangeable within
   calendar month. We report the fraction of permutations with
   `|β_interaction|` at least as large as observed; a two-sided p-value <
   0.05 rejects H0. This is the primary significance test.

2. **Metro-cluster bootstrap (N=2,000 resamples).** Metros are the
   independent unit; we resample metros with replacement and refit the
   panel. Percentile 95% confidence intervals on every coefficient allow
   the agent to read effect magnitude and uncertainty without assuming
   iid residuals. A CI that straddles zero should disqualify a
   significance claim even if the permutation p is small.

3. **Spread=0 placebo.** We zero out the spread regressor and refit; by
   construction the coefficients must be exactly 0 (mechanical null). A
   non-zero placebo coefficient indicates a numerical bug in the within
   transform, OLS solver, or panel assembly — verification check 21
   requires `|β_interaction| < 1e-10` under this placebo.

In addition, six sensitivity analyses (effective-rate window, spread lag,
pre/post-COVID subsample, top-50-metro restriction, 15-year PMMS vintage,
and the placebo above) probe specification robustness. Verification
check 27 requires all effective-rate-window coefficients to share the
main-sample sign; check 28 requires the same of the spread-lag variants;
check 29 checks the alternate-PMMS-vintage coefficient agrees in sign.

## Limitations and Assumptions

The findings from this skill are not a causal identification; they are a
panel decomposition conditional on the stated structural assumptions.
Agents consuming these results should cite the following caveats:

1. **Partial-equilibrium decomposition.** The counterfactual holds
   exposure, metro FE, and calendar-month FE fixed. Aggregate feedback —
   e.g., a spread-induced drop in listings tightens supply, which raises
   prices, which raises expected future spreads, which further depresses
   listings — is ignored. The `share_explained` is a lower bound on the
   direct channel magnitude and an upper bound on a naive "rate spread
   explains X% of listings" claim only if those second-order channels
   offset the first-order effect on average.

2. **Exposure proxy is imperfect.** Pre-2020 log median listing price is
   a proxy for lock-in intensity. Higher-priced metros likely carry
   larger mortgage balances and therefore larger absolute dollar-savings
   from a low locked-in rate, but the proxy also correlates with metro
   income, mobility, and remote-work penetration. The interaction
   coefficient conflates all of these.

3. **Observational design, not experimental.** The skill cannot rule out
   confounding from remote work, demographic aging, or secular trends
   that happen to align with the spread. The seasonal-block permutation
   test only perturbs year-to-year ordering within calendar month — so
   any slow-moving trend correlated with the spread survives the null
   and could inflate the coefficient.

4. **Seasonal-month FE choice.** Using calendar-month FE rather than
   full month-of-sample FE is intentional (otherwise `spread` is
   absorbed), but it means the identifying variation includes *all*
   aggregate shocks with year-over-year variation — not only the
   rate-spread itself.

5. **Data vintage drift.** SHA256 hashes pin the 2026-04-18 vintage of
   both CSVs; Freddie Mac or Realtor.com revisions will cause hash
   mismatch warnings. By default the script continues (logging drift);
   set `STRICT_SHA256=True` in the DOMAIN CONFIGURATION block to hard
   fail on drift and produce a perfectly frozen artifact.

6. **What the results do NOT show.** They do not show (a) the full
   welfare cost of lock-in, (b) whether lock-in is the *largest* single
   driver (the skill only quantifies share of log-change, not relative
   importance of competing channels), (c) causal identification, or
   (d) individual-household behavior.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_mortgage-lock-in-and-listings
```

**Expected output:** Directory created, exit code 0.

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_mortgage-lock-in-and-listings/analyze.py
#!/usr/bin/env python3
"""
How much of the post-2022 U.S. new-listings decline is explained by
mortgage rate lock-in?

Downloads Freddie Mac PMMS weekly 30-year fixed rates and Realtor.com
monthly metro-level new-listing counts. Constructs a national effective
rate from an exponentially-weighted 10-year trailing PMMS mean, computes
the rate spread, estimates a metro-fixed-effects panel regression of
log new-listings on (spread, spread x exposure, trend, 11 calendar-month
seasonal dummies), and runs a 2,000-iteration seasonal-block permutation
null, 2,000-iteration metro-cluster bootstrap CIs, and six sensitivity
analyses.

Python 3.8+ standard library only. All random operations seeded.
"""

import sys
import os
import json
import math
import time
import csv
import random
import hashlib
import urllib.request
import urllib.error
from collections import defaultdict
from io import StringIO

# ═══════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — To adapt this analysis to a new domain,
# modify only this section.
# ═══════════════════════════════════════════════════════════════

SEED = 42

PMMS_URL = "https://www.freddiemac.com/pmms/docs/PMMS_history.csv"
REALTOR_METRO_URL = ("https://econdata.s3-us-west-2.amazonaws.com/"
                     "Reports/Core/RDC_Inventory_Core_Metrics_Metro_History.csv")

# Expected SHA256 of data files as observed on 2026-04-18. If the remote
# file is later republished, the manifest will record the NEW hash; set
# STRICT_SHA256=True to hard-fail on mismatch instead.
EXPECTED_SHA256 = {
    "cache/pmms_history.csv":
        "0cf5fa2deb13990ea5b0fc87ad33267002484fae2c6d176726cfc95c79e1b831",
    "cache/realtor_metro.csv":
        "ae128d1d67a519fef40f1f937e0eb8acc155b5235f3fb4ea33b67c44f0d793ac",
}
STRICT_SHA256 = False  # set True to hard-fail on hash drift

PMMS_DATE_COL = "date"
PMMS_VALUE_COL = "pmms30"        # Freddie 30-yr fixed
PMMS_ALT_VALUE_COL = "pmms15"    # Alternate vintage (15-yr) for sensitivity

REALTOR_DATE_COL = "month_date_yyyymm"
UNIT_ID_COL = "cbsa_code"
UNIT_NAME_COL = "cbsa_title"
OUTCOME_COL = "new_listing_count"
EXPOSURE_PROXY_COL = "median_listing_price"   # pre-2020 mean used as exposure
QUALITY_FLAG_COL = "quality_flag"             # Realtor marks imputed rows >0

EFFECTIVE_RATE_WINDOW_MONTHS = 120             # 10-year trailing mean
EFFECTIVE_RATE_HALFLIFE_MONTHS = 60            # exponential weight halflife
SENS_WINDOWS = [60, 120, 180]                  # sensitivity windows
SENS_LAGS = [0, 3, 6]                          # spread lag (months)

BASELINE_START_YM = 201607                     # data start
BASELINE_END_YM = 201912                       # pre-COVID baseline end
TRAIN_END_YM = 202112                          # in-sample fit end
TREATMENT_START_YM = 202201                    # post-2022 window start
ANALYSIS_END_YM = 202603                       # data end

MIN_MONTHS_PER_UNIT = 60                       # filter thin panels
TOP_N_METROS = 300                             # keep top-N by HouseholdRank

N_PERMUTATIONS = 2000
N_BOOTSTRAP = 2000
PERM_ALPHA = 0.05

CACHE_DIR = "cache"
MANIFEST_FILE = "data_manifest.json"
UA = "MortgageLockInStudy/1.0 (mailto:claw4s-research@example.com)"

# Regression spec: metro FE + calendar-month-of-year FE (both absorbed
# via 2-way demeaning) + 3 core covariates. Key identification: `spread`
# varies across YEARS within a given calendar month, so the calendar-
# month FE does NOT absorb it (unlike full month-of-sample FE, which
# would). Unit-only covariates (`exposure`) are absorbed by metro FE and
# dropped from the X vector.
CORE_KEYS = ("spread", "spread_x_exposure", "trend")
X_KEYS = CORE_KEYS

# ═══════════════════════════════════════════════════════════════
# Helper utilities
# ═══════════════════════════════════════════════════════════════

def sha256_file(path):
    h = hashlib.sha256()
    with open(path, 'rb') as f:
        for chunk in iter(lambda: f.read(8192), b''):
            h.update(chunk)
    return h.hexdigest()


def load_manifest():
    if os.path.exists(MANIFEST_FILE):
        with open(MANIFEST_FILE) as f:
            return json.load(f)
    return {}


def save_manifest(m):
    with open(MANIFEST_FILE, 'w') as f:
        json.dump(m, f, indent=2)


def http_get_cached(url, cache_path, manifest, retries=4):
    """GET with local caching and SHA256 verification."""
    os.makedirs(os.path.dirname(cache_path) or '.', exist_ok=True)
    if os.path.exists(cache_path):
        h = sha256_file(cache_path)
        exp = manifest.get(cache_path)
        if exp is None or h == exp:
            manifest[cache_path] = h
            _check_expected(cache_path, h)
            with open(cache_path, 'rb') as f:
                return f.read().decode('utf-8', errors='replace')
        else:
            print(f"    Cache corrupted: {cache_path}, re-downloading")
            os.remove(cache_path)
    for attempt in range(retries):
        try:
            req = urllib.request.Request(
                url, headers={'User-Agent': UA, 'Accept': 'text/csv,*/*'})
            with urllib.request.urlopen(req, timeout=120) as r:
                raw = r.read()
            with open(cache_path, 'wb') as f:
                f.write(raw)
            h = sha256_file(cache_path)
            manifest[cache_path] = h
            _check_expected(cache_path, h)
            return raw.decode('utf-8', errors='replace')
        except urllib.error.HTTPError as e:
            wait = (8 if e.code == 429 else 2) * (attempt + 1)
            if attempt < retries - 1:
                print(f"    HTTP {e.code}, retry in {wait}s")
                time.sleep(wait)
            else:
                raise RuntimeError(f"HTTP {e.code} after {retries} retries: {url}")
        except Exception as e:
            if attempt < retries - 1:
                time.sleep(2 ** (attempt + 1))
            else:
                raise RuntimeError(f"Failed after {retries} retries: {url}: {e}")


def _check_expected(cache_path, actual_hash):
    exp = EXPECTED_SHA256.get(cache_path)
    if exp is None:
        return
    if actual_hash == exp:
        print(f"    SHA256 matches expected: {cache_path}")
    else:
        msg = (f"    SHA256 DRIFT: {cache_path} "
               f"expected={exp[:16]}... actual={actual_hash[:16]}...")
        if STRICT_SHA256:
            raise RuntimeError(msg + "  (STRICT_SHA256=True)")
        print(msg + "  (continuing; data source may have been republished)")


def parse_float(s):
    try:
        if s is None or s == "" or s.strip() == "":
            return None
        return float(s)
    except (ValueError, TypeError):
        return None


def parse_int(s):
    try:
        if s is None or s == "" or s.strip() == "":
            return None
        return int(float(s))
    except (ValueError, TypeError):
        return None


def ym_to_monthidx(ym):
    """202601 -> months since 2000-01."""
    y = ym // 100
    m = ym % 100
    return (y - 2000) * 12 + (m - 1)


def monthidx_to_ym(idx):
    y = idx // 12 + 2000
    m = idx % 12 + 1
    return y * 100 + m


# ═══════════════════════════════════════════════════════════════
# Data loading — Freddie Mac PMMS + Realtor.com metro panel
# ═══════════════════════════════════════════════════════════════

def load_pmms(manifest):
    """Return dict month_idx -> {'pmms30': x, 'pmms15': y} national monthly."""
    raw = http_get_cached(PMMS_URL,
                          os.path.join(CACHE_DIR, "pmms_history.csv"),
                          manifest)
    reader = csv.DictReader(StringIO(raw))
    weekly = defaultdict(list)
    weekly_alt = defaultdict(list)
    for row in reader:
        datestr = (row.get(PMMS_DATE_COL) or "").strip()
        if not datestr or '/' not in datestr:
            continue
        try:
            m, d, y = datestr.split('/')
            y = int(y)
            m = int(m)
        except Exception:
            continue
        if y < 1971 or y > 2030:
            continue
        idx = (y - 2000) * 12 + (m - 1)
        v = parse_float(row.get(PMMS_VALUE_COL))
        va = parse_float(row.get(PMMS_ALT_VALUE_COL))
        if v is not None:
            weekly[idx].append(v)
        if va is not None:
            weekly_alt[idx].append(va)
    monthly = {}
    for idx, vs in weekly.items():
        if vs:
            monthly[idx] = {
                'pmms30': sum(vs) / len(vs),
                'pmms15': (sum(weekly_alt.get(idx, [])) /
                          max(1, len(weekly_alt.get(idx, []))))
                          if weekly_alt.get(idx) else None,
            }
    return monthly


def load_realtor_metro(manifest):
    """Return list of dicts: one per (metro,month) observation."""
    raw = http_get_cached(REALTOR_METRO_URL,
                          os.path.join(CACHE_DIR, "realtor_metro.csv"),
                          manifest)
    reader = csv.DictReader(StringIO(raw))
    rows = []
    for row in reader:
        ym = parse_int(row.get(REALTOR_DATE_COL))
        uid = row.get(UNIT_ID_COL, "").strip()
        name = row.get(UNIT_NAME_COL, "").strip()
        if ym is None or not uid:
            continue
        if ym < BASELINE_START_YM or ym > ANALYSIS_END_YM:
            continue
        outcome = parse_int(row.get(OUTCOME_COL))
        price = parse_float(row.get(EXPOSURE_PROXY_COL))
        qf = parse_float(row.get(QUALITY_FLAG_COL))
        hh_rank = parse_int(row.get("HouseholdRank"))
        rows.append({
            'ym': ym,
            'midx': ym_to_monthidx(ym),
            'uid': uid,
            'name': name,
            'outcome': outcome,
            'price': price,
            'qf': qf,
            'hh_rank': hh_rank,
        })
    return rows


def load_data():
    os.makedirs(CACHE_DIR, exist_ok=True)
    manifest = load_manifest()
    print("  [2a/10] Freddie Mac PMMS...")
    pmms_monthly = load_pmms(manifest)
    save_manifest(manifest)
    print(f"    {len(pmms_monthly)} PMMS monthly rows")
    print("  [2b/10] Realtor.com Metro panel...")
    realtor_rows = load_realtor_metro(manifest)
    save_manifest(manifest)
    print(f"    {len(realtor_rows)} Realtor.com metro-month rows")
    return pmms_monthly, realtor_rows


# ═══════════════════════════════════════════════════════════════
# Statistical helpers (stdlib)
# ═══════════════════════════════════════════════════════════════

def mean_val(xs):
    return sum(xs) / len(xs) if xs else 0.0


def percentile(xs, q):
    if not xs:
        return 0.0
    s = sorted(xs)
    n = len(s)
    k = max(0, min(n - 1, int(round(q * (n - 1)))))
    return s[k]


def rate_spread(pmms_monthly, value_key='pmms30',
                window=EFFECTIVE_RATE_WINDOW_MONTHS,
                halflife=EFFECTIVE_RATE_HALFLIFE_MONTHS):
    """For each month with enough history, return dict keyed by midx
    with current rate, effective rate (exp-weighted trailing mean), and
    spread."""
    if halflife <= 0:
        lam = 0.0
    else:
        lam = math.log(2.0) / halflife
    out = {}
    all_idx = sorted(pmms_monthly.keys())
    for idx in all_idx:
        cur = pmms_monthly[idx].get(value_key)
        if cur is None:
            continue
        hist = []
        for k in range(1, window + 1):
            if idx - k in pmms_monthly:
                v = pmms_monthly[idx - k].get(value_key)
                if v is not None:
                    w = math.exp(-lam * k)
                    hist.append((v, w))
        if not hist or len(hist) < max(12, window // 2):
            continue
        sw = sum(w for _, w in hist)
        eff = sum(v * w for v, w in hist) / sw if sw > 0 else cur
        out[idx] = {
            'current': cur,
            'effective': eff,
            'spread': cur - eff,
        }
    return out


def two_way_cal_within_transform(records, unit_key, cal_key, y_key, x_keys):
    """Two-way within transform: subtract unit mean AND calendar-month
    mean, add grand mean. This absorbs metro FE + calendar-month-of-year
    FE. It does NOT absorb time (month-of-sample) variation, so any
    regressor that varies across years within a given calendar month
    remains identified (e.g., national rate `spread`, linear `trend`)."""
    unit_y = defaultdict(list)
    cal_y = defaultdict(list)
    for r in records:
        unit_y[r[unit_key]].append(r[y_key])
        cal_y[r[cal_key]].append(r[y_key])
    unit_mean_y = {u: mean_val(vs) for u, vs in unit_y.items()}
    cal_mean_y = {c: mean_val(vs) for c, vs in cal_y.items()}
    grand_y = mean_val([r[y_key] for r in records])
    x_unit_means, x_cal_means, x_grand = {}, {}, {}
    for k in x_keys:
        u_acc = defaultdict(list)
        c_acc = defaultdict(list)
        for r in records:
            u_acc[r[unit_key]].append(r[k])
            c_acc[r[cal_key]].append(r[k])
        x_unit_means[k] = {u: mean_val(vs) for u, vs in u_acc.items()}
        x_cal_means[k] = {c: mean_val(vs) for c, vs in c_acc.items()}
        x_grand[k] = mean_val([r[k] for r in records])
    y_t, x_t = [], []
    for r in records:
        u, c = r[unit_key], r[cal_key]
        y_t.append(r[y_key] - unit_mean_y[u] - cal_mean_y[c] + grand_y)
        x_t.append([r[k] - x_unit_means[k][u] - x_cal_means[k][c] + x_grand[k]
                    for k in x_keys])
    return y_t, x_t


def ols_multi(y, X):
    """Multivariate OLS without intercept (assumes demeaned covariates).
    Returns (beta, diag_dict)."""
    if not y:
        return [0.0] * (len(X[0]) if X else 0), None
    n = len(y)
    k = len(X[0]) if X[0] else 0
    if k == 0:
        return [], {'n': n, 'k': 0, 'ssr': 0.0, 'tss': 0.0, 'r2': 0.0}
    xtx = [[0.0] * k for _ in range(k)]
    xty = [0.0] * k
    for i in range(n):
        xi = X[i]
        yi = y[i]
        for a in range(k):
            xty[a] += xi[a] * yi
            for b in range(a, k):
                xtx[a][b] += xi[a] * xi[b]
    for a in range(k):
        for b in range(a):
            xtx[a][b] = xtx[b][a]
    A = [row[:] + [rhs] for row, rhs in zip(xtx, xty)]
    for i in range(k):
        piv = A[i][i]
        if abs(piv) < 1e-12:
            for j in range(i + 1, k):
                if abs(A[j][i]) > 1e-12:
                    A[i], A[j] = A[j], A[i]
                    piv = A[i][i]
                    break
        if abs(piv) < 1e-12:
            return [0.0] * k, None
        inv = 1.0 / piv
        for j in range(k + 1):
            A[i][j] *= inv
        for r_ in range(k):
            if r_ != i and abs(A[r_][i]) > 1e-12:
                f = A[r_][i]
                for j in range(k + 1):
                    A[r_][j] -= f * A[i][j]
    beta = [A[i][k] for i in range(k)]
    ssr = 0.0
    tss = 0.0
    ym = mean_val(y)
    for i in range(n):
        yhat = sum(X[i][a] * beta[a] for a in range(k))
        ssr += (y[i] - yhat) ** 2
        tss += (y[i] - ym) ** 2
    r2 = 1.0 - ssr / tss if tss > 0 else 0.0
    return beta, {'n': n, 'k': k, 'ssr': ssr, 'tss': tss, 'r2': r2}


def fit_panel_within(records, y_key='log_outcome', x_keys=X_KEYS):
    """Metro-FE + calendar-month-of-year-FE within-OLS."""
    y_t, x_t = two_way_cal_within_transform(
        records, 'uid', 'calendar_month', y_key, list(x_keys))
    beta, diag = ols_multi(y_t, x_t)
    return {k: b for k, b in zip(x_keys, beta)}, diag


def _attach_spread(records, spread_by_midx):
    """Return NEW list of records with spread, spread_x_exposure, trend
    re-computed from the supplied spread dict."""
    out = []
    for r in records:
        sp = spread_by_midx.get(r['midx'])
        if sp is None:
            continue
        out.append({
            'uid': r['uid'],
            'midx': r['midx'],
            'log_outcome': r['log_outcome'],
            'exposure': r['exposure'],
            'trend': r['trend'],
            'calendar_month': r['calendar_month'],
            'spread': sp,
            'spread_x_exposure': sp * r['exposure'],
        })
    return out


def permutation_p_seasonal_block(records, spread_by_midx, calendar_of,
                                 n_perms, rng):
    """Seasonal-block permutation: for each calendar month m, randomly
    re-assign the spread values observed in that month across years.
    Re-fit the panel; count |beta_interaction| >= observed."""
    obs_recs = _attach_spread(records, spread_by_midx)
    obs_beta, _ = fit_panel_within(obs_recs)
    obs_int = obs_beta['spread_x_exposure']

    cal_to_midx = defaultdict(list)
    for m in spread_by_midx:
        cal_to_midx[calendar_of[m]].append(m)

    perm_betas = []
    n_ge = 0
    for _ in range(n_perms):
        spread_perm = {}
        for cal, midxs in cal_to_midx.items():
            vals = [spread_by_midx[m] for m in midxs]
            rng.shuffle(vals)
            for m, v in zip(midxs, vals):
                spread_perm[m] = v
        perm_recs = _attach_spread(records, spread_perm)
        beta_p, _ = fit_panel_within(perm_recs)
        perm_betas.append(beta_p['spread_x_exposure'])
        if abs(beta_p['spread_x_exposure']) >= abs(obs_int):
            n_ge += 1
    p = (n_ge + 1) / (n_perms + 1)
    return p, obs_int, perm_betas


def bootstrap_ci_cluster(records, n_boot, rng, level=0.95):
    """Metro-cluster bootstrap; refit one-way FE; return percentile CIs
    for every x_key."""
    by_uid = defaultdict(list)
    for r in records:
        by_uid[r['uid']].append(r)
    uids = list(by_uid.keys())
    boots = {k: [] for k in X_KEYS}
    for _ in range(n_boot):
        sample_records = []
        for _u in uids:
            pick = uids[rng.randrange(0, len(uids))]
            for rec in by_uid[pick]:
                sample_records.append(rec)
        beta, _ = fit_panel_within(sample_records)
        for k in X_KEYS:
            boots[k].append(beta[k])
    out = {}
    for k in X_KEYS:
        b = sorted(boots[k])
        m = mean_val(boots[k])
        lo_idx = max(0, int((1 - level) / 2 * n_boot))
        hi_idx = min(n_boot - 1, int((1 + level) / 2 * n_boot))
        out[k] = {'mean': m, 'lo': b[lo_idx], 'hi': b[hi_idx]}
    return out


def counterfactual_decomposition(records, beta,
                                 treatment_midx_min, treatment_midx_max):
    """Partial-equilibrium: attribute log-change in mean outcome to the
    change in mean spread (level + interaction) using fitted betas."""
    pre = [r['log_outcome'] for r in records
           if r['midx'] < treatment_midx_min]
    post = [r['log_outcome'] for r in records
            if treatment_midx_min <= r['midx'] <= treatment_midx_max]
    if not pre or not post:
        return {
            'actual_log_delta': 0.0, 'actual_pct_delta': 0.0,
            'share_explained': 0.0, 'n_pre': len(pre), 'n_post': len(post),
        }
    actual_delta = mean_val(post) - mean_val(pre)
    pre_rows = [r for r in records if r['midx'] < treatment_midx_min]
    post_rows = [r for r in records
                 if treatment_midx_min <= r['midx'] <= treatment_midx_max]
    pre_spread_level = mean_val([r['spread'] for r in pre_rows])
    post_spread_level = mean_val([r['spread'] for r in post_rows])
    pre_spread_inter = mean_val([r['spread_x_exposure'] for r in pre_rows])
    post_spread_inter = mean_val([r['spread_x_exposure'] for r in post_rows])
    dspread_level = post_spread_level - pre_spread_level
    dspread_inter = post_spread_inter - pre_spread_inter
    contrib_level = beta.get('spread', 0.0) * dspread_level
    contrib_inter = beta.get('spread_x_exposure', 0.0) * dspread_inter
    contrib_total = contrib_level + contrib_inter
    share = contrib_total / actual_delta if abs(actual_delta) > 1e-9 else 0.0
    return {
        'actual_log_delta': actual_delta,
        'actual_pct_delta': (math.exp(actual_delta) - 1) * 100,
        'pre_mean_spread': pre_spread_level,
        'post_mean_spread': post_spread_level,
        'delta_spread_level': dspread_level,
        'delta_spread_interaction': dspread_inter,
        'contribution_level_log': contrib_level,
        'contribution_interaction_log': contrib_inter,
        'spread_contribution_log': contrib_total,
        'share_explained': share,
        'n_pre': len(pre_rows), 'n_post': len(post_rows),
    }


# ═══════════════════════════════════════════════════════════════
# Sample assembly
# ═══════════════════════════════════════════════════════════════

def build_records(pmms_monthly, realtor_rows, value_key='pmms30',
                  window=EFFECTIVE_RATE_WINDOW_MONTHS,
                  halflife=EFFECTIVE_RATE_HALFLIFE_MONTHS,
                  lag_months=0):
    """Filter realtor rows, attach spread + exposure + seasonal dummies."""
    spread_data = rate_spread(pmms_monthly, value_key, window, halflife)
    spread_by_midx = {m: d['spread'] for m, d in spread_data.items()}

    metro_prices = defaultdict(list)
    for r in realtor_rows:
        if r['ym'] <= BASELINE_END_YM and r['price'] is not None:
            metro_prices[r['uid']].append(r['price'])
    metro_exposure = {}
    for u, ps in metro_prices.items():
        if len(ps) >= 12:
            metro_exposure[u] = math.log(max(1.0, mean_val(ps)))
    if metro_exposure:
        mean_e = mean_val(list(metro_exposure.values()))
        var_e = mean_val([(v - mean_e) ** 2 for v in metro_exposure.values()])
        sd_e = math.sqrt(var_e) if var_e > 0 else 1.0
        metro_exposure = {u: (v - mean_e) / sd_e for u, v in metro_exposure.items()}

    rank_of = {}
    for r in realtor_rows:
        if r['hh_rank'] is not None:
            rank_of[r['uid']] = min(rank_of.get(r['uid'], 99999), r['hh_rank'])
    top_metros = {u for u, rk in rank_of.items()
                  if rk is not None and rk <= TOP_N_METROS}

    records = []
    months_per_unit = defaultdict(int)
    for r in realtor_rows:
        if r['uid'] not in top_metros:
            continue
        if r['uid'] not in metro_exposure:
            continue
        if r['outcome'] is None or r['outcome'] <= 0:
            continue
        if r['qf'] is not None and r['qf'] > 0:
            continue
        shifted_midx = r['midx'] - lag_months
        sp = spread_by_midx.get(shifted_midx)
        if sp is None:
            continue
        records.append({
            'uid': r['uid'],
            'midx': r['midx'],
            'shifted_midx': shifted_midx,
            'ym': r['ym'],
            'log_outcome': math.log(r['outcome']),
            'calendar_month': r['ym'] % 100,
            'exposure': metro_exposure[r['uid']],
        })
        months_per_unit[r['uid']] += 1

    keep_uid = {u for u, n in months_per_unit.items() if n >= MIN_MONTHS_PER_UNIT}
    records = [r for r in records if r['uid'] in keep_uid]

    if records:
        first_midx = min(r['midx'] for r in records)
        for r in records:
            r['trend'] = r['midx'] - first_midx

    # Spread dict keyed by the midx that actually appears in the panel
    # (after optional lag). Each record carries spread, interaction, and
    # seasonal dummies for easy downstream use.
    spread_lookup = {}
    for r in records:
        sp = spread_by_midx[r['shifted_midx']]
        r['spread'] = sp
        r['spread_x_exposure'] = sp * r['exposure']
        spread_lookup[r['midx']] = sp
    return records, spread_lookup, metro_exposure


# ═══════════════════════════════════════════════════════════════
# Main analysis
# ═══════════════════════════════════════════════════════════════

def run_analysis(pmms_monthly, realtor_rows):
    rng = random.Random(SEED)

    # ── [4/10] Build main analysis sample
    print("[4/10] Building main sample (10-yr eff rate, 0 lag)...")
    records, spread_by_midx, exposure = build_records(
        pmms_monthly, realtor_rows, 'pmms30',
        EFFECTIVE_RATE_WINDOW_MONTHS, EFFECTIVE_RATE_HALFLIFE_MONTHS, 0)
    n_units = len({r['uid'] for r in records})
    n_months = len({r['midx'] for r in records})
    print(f"  n_obs={len(records)}  n_metros={n_units}  n_months={n_months}")
    if len(records) < 5000 or n_units < 50:
        raise RuntimeError("Sample too small after filtering")

    # ── [5/10] Main panel FE regression
    print("[5/10] Panel within-OLS (metro FE + month-of-year dummies)...")
    beta_main, diag_main = fit_panel_within(records)
    print(f"  beta[spread] = {beta_main['spread']:.4f}")
    print(f"  beta[spread_x_exposure] = {beta_main['spread_x_exposure']:.4f}")
    print(f"  beta[trend] = {beta_main['trend']:.6f}")
    print(f"  R^2 (within) = {diag_main['r2']:.3f}")

    # ── [6/10] Seasonal-block permutation null
    print(f"[6/10] Seasonal-block permutation null (N={N_PERMUTATIONS})...")
    calendar_of = {r['midx']: r['calendar_month'] for r in records}
    p_perm, obs_int, perm_dist = permutation_p_seasonal_block(
        records, spread_by_midx, calendar_of, N_PERMUTATIONS, rng)
    print(f"  observed beta[spread_x_exposure] = {obs_int:.4f}")
    print(f"  two-sided p = {p_perm:.4f}")
    # Record null distribution diagnostics for sanity checking: under
    # H0 the permutation mean should be close to zero and the standard
    # deviation should be strictly positive. These two summaries are a
    # cheap but powerful sanity check that the permutation machinery is
    # producing a non-degenerate null.
    perm_mean = mean_val(perm_dist)
    perm_var = mean_val([(b - perm_mean) ** 2 for b in perm_dist])
    perm_std = math.sqrt(perm_var) if perm_var > 0 else 0.0
    perm_abs_max = max((abs(b) for b in perm_dist), default=0.0)
    perm_z = (obs_int - perm_mean) / perm_std if perm_std > 0 else 0.0
    print(f"  null mean={perm_mean:.5f}  null std={perm_std:.5f}  "
          f"z(obs)={perm_z:.2f}")

    # ── [7/10] Cluster bootstrap CIs
    print(f"[7/10] Cluster bootstrap (N={N_BOOTSTRAP}, by metro)...")
    ci = bootstrap_ci_cluster(records, N_BOOTSTRAP, rng)
    for k in CORE_KEYS:
        v = ci[k]
        print(f"  {k}: mean={v['mean']:.4f}  95% CI [{v['lo']:.4f}, {v['hi']:.4f}]")

    # ── [8/10] Counterfactual decomposition
    print("[8/10] Counterfactual decomposition: spread=0 world...")
    tr_min = ym_to_monthidx(TREATMENT_START_YM)
    tr_max = ym_to_monthidx(ANALYSIS_END_YM)
    decomp = counterfactual_decomposition(records, beta_main, tr_min, tr_max)
    print(f"  actual log-delta (post vs pre-2022): {decomp['actual_log_delta']:.4f}")
    print(f"  spread contribution (log): {decomp['spread_contribution_log']:.4f}")
    print(f"  share of decline explained: {decomp['share_explained']:.3f}")

    # ── [9/10] Sensitivity analyses
    print("[9/10] Sensitivity analyses...")
    sens = {}

    print("  (a) Effective-rate windows...")
    for w in SENS_WINDOWS:
        rec_w, sp_w, exp_w = build_records(
            pmms_monthly, realtor_rows, 'pmms30',
            w, EFFECTIVE_RATE_HALFLIFE_MONTHS, 0)
        if len(rec_w) < 1000:
            continue
        beta_w, _ = fit_panel_within(rec_w)
        sens[f'window_{w}'] = {
            'window_months': w,
            'n_obs': len(rec_w),
            'beta_spread': beta_w['spread'],
            'beta_spread_x_exposure': beta_w['spread_x_exposure'],
        }
        print(f"    w={w}mo  beta_spread={beta_w['spread']:.4f} "
              f"beta_int={beta_w['spread_x_exposure']:.4f} (n={len(rec_w)})")

    print("  (b) Spread lag...")
    for lag in SENS_LAGS:
        rec_l, sp_l, exp_l = build_records(
            pmms_monthly, realtor_rows, 'pmms30',
            EFFECTIVE_RATE_WINDOW_MONTHS, EFFECTIVE_RATE_HALFLIFE_MONTHS, lag)
        if len(rec_l) < 1000:
            continue
        beta_l, _ = fit_panel_within(rec_l)
        sens[f'lag_{lag}'] = {
            'lag_months': lag,
            'n_obs': len(rec_l),
            'beta_spread': beta_l['spread'],
            'beta_spread_x_exposure': beta_l['spread_x_exposure'],
        }
        print(f"    lag={lag}mo  beta_spread={beta_l['spread']:.4f} "
              f"beta_int={beta_l['spread_x_exposure']:.4f} (n={len(rec_l)})")

    print("  (c) Pre/post-COVID subsamples...")
    covid_cut = ym_to_monthidx(202003)
    for label, lo_m, hi_m in [
        ('pre_covid', ym_to_monthidx(BASELINE_START_YM), covid_cut - 1),
        ('post_covid', covid_cut, ym_to_monthidx(ANALYSIS_END_YM)),
        ('drop_covid_year',
         ym_to_monthidx(BASELINE_START_YM),
         ym_to_monthidx(ANALYSIS_END_YM)),
    ]:
        if label == 'drop_covid_year':
            sub = [r for r in records
                   if not (ym_to_monthidx(202003) <= r['midx']
                           <= ym_to_monthidx(202102))]
        else:
            sub = [r for r in records if lo_m <= r['midx'] <= hi_m]
        if len(sub) < 500:
            continue
        beta_s, _ = fit_panel_within(sub)
        sens[f'sample_{label}'] = {
            'label': label,
            'n_obs': len(sub),
            'beta_spread': beta_s['spread'],
            'beta_spread_x_exposure': beta_s['spread_x_exposure'],
        }
        print(f"    {label}  beta_spread={beta_s['spread']:.4f} "
              f"beta_int={beta_s['spread_x_exposure']:.4f} (n={len(sub)})")

    print("  (d) Top-50 metros...")
    rank_first = {}
    for r in realtor_rows:
        if r['hh_rank'] is not None:
            rank_first[r['uid']] = min(rank_first.get(r['uid'], 99999),
                                        r['hh_rank'])
    top50 = {u for u, rk in rank_first.items() if rk <= 50}
    sub50 = [r for r in records if r['uid'] in top50]
    if len(sub50) > 500:
        beta_50, _ = fit_panel_within(sub50)
        sens['top50'] = {
            'n_obs': len(sub50),
            'n_metros': len({r['uid'] for r in sub50}),
            'beta_spread': beta_50['spread'],
            'beta_spread_x_exposure': beta_50['spread_x_exposure'],
        }
        print(f"    top50 metros  beta_spread={beta_50['spread']:.4f} "
              f"beta_int={beta_50['spread_x_exposure']:.4f} (n={len(sub50)})")

    print("  (e) Alternate PMMS vintage (15-year)...")
    has_alt = any(pmms_monthly[m].get('pmms15') is not None
                  for m in pmms_monthly)
    if has_alt:
        rec_alt, sp_alt, exp_alt = build_records(
            pmms_monthly, realtor_rows, 'pmms15',
            EFFECTIVE_RATE_WINDOW_MONTHS,
            EFFECTIVE_RATE_HALFLIFE_MONTHS, 0)
        if len(rec_alt) > 1000:
            beta_alt, _ = fit_panel_within(rec_alt)
            sens['alt_pmms15'] = {
                'value_col': 'pmms15',
                'n_obs': len(rec_alt),
                'beta_spread': beta_alt['spread'],
                'beta_spread_x_exposure': beta_alt['spread_x_exposure'],
            }
            print(f"    pmms15  beta_spread={beta_alt['spread']:.4f} "
                  f"beta_int={beta_alt['spread_x_exposure']:.4f} "
                  f"(n={len(rec_alt)})")

    print("  (f) Placebo: spread set to zero (mechanical null)...")
    rec_p = []
    for r in records:
        rp = dict(r)
        rp['spread'] = 0.0
        rp['spread_x_exposure'] = 0.0
        rec_p.append(rp)
    beta_p, _ = fit_panel_within(rec_p)
    sens['placebo_zero_spread'] = {
        'beta_spread': beta_p['spread'],
        'beta_spread_x_exposure': beta_p['spread_x_exposure'],
        'note': 'Expect exact zero by construction',
    }
    print(f"    placebo (spread=0)  beta_int="
          f"{beta_p['spread_x_exposure']:.6f}")

    # ── [10/10] Assemble results
    bootstrap_ci_out = {k: {'mean': ci[k]['mean'],
                            'lo': ci[k]['lo'], 'hi': ci[k]['hi']}
                        for k in CORE_KEYS}
    results = {
        'config': {
            'seed': SEED,
            'baseline_start_ym': BASELINE_START_YM,
            'baseline_end_ym': BASELINE_END_YM,
            'treatment_start_ym': TREATMENT_START_YM,
            'analysis_end_ym': ANALYSIS_END_YM,
            'effective_rate_window_months': EFFECTIVE_RATE_WINDOW_MONTHS,
            'effective_rate_halflife_months': EFFECTIVE_RATE_HALFLIFE_MONTHS,
            'top_n_metros': TOP_N_METROS,
            'min_months_per_unit': MIN_MONTHS_PER_UNIT,
            'n_permutations': N_PERMUTATIONS,
            'n_bootstrap': N_BOOTSTRAP,
            'pmms_url': PMMS_URL,
            'realtor_url': REALTOR_METRO_URL,
            'fe_specification': ('metro FE + calendar-month-of-year FE '
                                 '(both absorbed via 2-way demeaning); '
                                 'linear month-index trend as covariate'),
            'expected_sha256': EXPECTED_SHA256,
        },
        'sample': {
            'n_obs': len(records),
            'n_metros': n_units,
            'n_months': n_months,
            'min_ym': min(r['ym'] for r in records),
            'max_ym': max(r['ym'] for r in records),
        },
        'main': {
            'beta_spread': beta_main['spread'],
            'beta_spread_x_exposure': beta_main['spread_x_exposure'],
            'beta_trend': beta_main['trend'],
            'r2': diag_main['r2'],
            'n_obs': diag_main['n'],
            'permutation_p_two_sided': p_perm,
            'permutation_obs_int': obs_int,
            'permutation_null_mean': perm_mean,
            'permutation_null_std': perm_std,
            'permutation_null_abs_max': perm_abs_max,
            'permutation_obs_z': perm_z,
            'bootstrap_ci': bootstrap_ci_out,
        },
        'decomposition': decomp,
        'sensitivity': sens,
        'limitations': [
            "Partial-equilibrium decomposition: holds metro FE, calendar-month "
            "FE, and exposure fixed; aggregate feedback (spread -> listings -> "
            "prices -> future spreads) is ignored.",
            "Exposure proxy = pre-2020 log median listing price; correlates "
            "with lock-in intensity but also with metro income, mobility, and "
            "remote-work penetration, so the interaction coefficient conflates "
            "these channels.",
            "Observational design cannot rule out confounding from remote "
            "work, demographic aging, or secular trends that happen to align "
            "with the rate spread.",
            "Calendar-month FE (not month-of-sample FE) is intentional so "
            "that spread remains identified, but this means the residual "
            "variation includes all aggregate year-over-year shocks.",
            "SHA256 pins the 2026-04-18 vintage of both CSVs; later "
            "republication will trigger a drift warning (or a hard failure "
            "if STRICT_SHA256=True).",
            "Results do NOT show causal identification, the full welfare "
            "cost of lock-in, or whether lock-in is the single largest "
            "driver of the listings decline.",
        ],
    }
    return results


# ═══════════════════════════════════════════════════════════════
# Report generation
# ═══════════════════════════════════════════════════════════════

def generate_report(results):
    with open('results.json', 'w') as f:
        json.dump(results, f, indent=2)
    m = results['main']
    s = results['sample']
    d = results['decomposition']
    ci = m['bootstrap_ci']
    with open('report.md', 'w') as f:
        f.write("# Mortgage Rate Lock-In and Post-2022 U.S. Listings — Report\n\n")
        f.write("## Sample\n\n")
        f.write(f"- Metro-months: {s['n_obs']}\n")
        f.write(f"- Metros: {s['n_metros']}\n")
        f.write(f"- Months: {s['n_months']}\n")
        f.write(f"- Date range: {s['min_ym']} to {s['max_ym']}\n\n")
        f.write("## Main panel FE regression "
                "(metro FE + seasonal dummies + trend)\n\n")
        f.write("| Term | Coefficient | Bootstrap 95% CI |\n|---|---|---|\n")
        for k in CORE_KEYS:
            f.write(f"| {k} | {m['beta_'+k]:.4f} "
                    f"| [{ci[k]['lo']:.4f}, {ci[k]['hi']:.4f}] |\n")
        f.write(f"\n- R² (within): {m['r2']:.3f}\n")
        f.write(f"- Permutation p (two-sided, N={results['config']['n_permutations']}): "
                f"{m['permutation_p_two_sided']:.4f}\n\n")
        f.write("## Counterfactual decomposition\n\n")
        f.write(f"- Actual log change (post-2022 vs baseline): "
                f"{d['actual_log_delta']:.4f} ({d['actual_pct_delta']:.1f}%)\n")
        f.write(f"- Spread-attributable log change: "
                f"{d['spread_contribution_log']:.4f}\n")
        f.write(f"- **Share of decline explained by rate lock-in: "
                f"{d['share_explained']*100:.1f}%**\n")
        if 'limitations' in results:
            f.write("\n## Limitations\n\n")
            for i, lim in enumerate(results['limitations'], 1):
                f.write(f"{i}. {lim}\n")
    print("  results.json, report.md written")


# ═══════════════════════════════════════════════════════════════
# Main
# ═══════════════════════════════════════════════════════════════

def main():
    if '--verify' in sys.argv:
        return verify()

    t0 = time.time()
    print("[1/10] Workspace prep...")
    try:
        os.makedirs(CACHE_DIR, exist_ok=True)
    except OSError as e:
        sys.stderr.write(
            f"ERROR: Cannot create cache dir {CACHE_DIR!r}: {e}\n"
            "Check filesystem permissions and free space.\n")
        sys.exit(10)

    print("[2/10] Downloading data (PMMS + Realtor.com metro panel)...")
    try:
        pmms_monthly, realtor_rows = load_data()
    except (urllib.error.URLError, urllib.error.HTTPError) as e:
        sys.stderr.write(
            f"ERROR: Network download failed: {e}\n"
            f"Check internet connectivity to:\n"
            f"  - {PMMS_URL}\n"
            f"  - {REALTOR_METRO_URL}\n"
            "Cached files (if any) are preserved in "
            f"{CACHE_DIR}/ and will be re-used on rerun.\n")
        sys.exit(11)
    except RuntimeError as e:
        sys.stderr.write(
            f"ERROR: Data load failed: {e}\n"
            "If this is an SHA256 hash drift warning that hard-failed, "
            "set STRICT_SHA256=False in the DOMAIN CONFIGURATION block "
            "to continue on drift.\n")
        sys.exit(12)
    except Exception as e:
        sys.stderr.write(
            f"ERROR: Unexpected data-loading failure: "
            f"{type(e).__name__}: {e}\n")
        sys.exit(13)

    print(f"[3/10] Parsed {len(pmms_monthly)} PMMS monthly rows, "
          f"{len(realtor_rows)} Realtor.com rows")
    if len(pmms_monthly) < 120 or len(realtor_rows) < 10000:
        sys.stderr.write(
            "ERROR: Raw sample too small after parsing "
            f"(pmms_monthly={len(pmms_monthly)}, "
            f"realtor_rows={len(realtor_rows)}). "
            "Expected >=120 PMMS months and >=10,000 Realtor.com rows. "
            "Data source may have changed schema; inspect "
            f"{CACHE_DIR}/pmms_history.csv and "
            f"{CACHE_DIR}/realtor_metro.csv.\n")
        sys.exit(14)

    try:
        results = run_analysis(pmms_monthly, realtor_rows)
        generate_report(results)
    except Exception as e:
        import traceback
        sys.stderr.write(
            f"ERROR: Analysis stage failed: {type(e).__name__}: {e}\n")
        traceback.print_exc(file=sys.stderr)
        sys.exit(15)

    elapsed = time.time() - t0
    print(f"\nRuntime: {elapsed:.0f}s")
    print("ANALYSIS COMPLETE")


# ═══════════════════════════════════════════════════════════════
# Verification
# ═══════════════════════════════════════════════════════════════

def verify():
    print("Running verification...\n")
    ok = fail = 0

    def chk(name, cond):
        nonlocal ok, fail
        status = "PASS" if cond else "FAIL"
        print(f"  {status}: {name}")
        if cond:
            ok += 1
        else:
            fail += 1

    if not os.path.exists('results.json'):
        print("FAIL: results.json not found")
        sys.exit(1)
    with open('results.json') as f:
        r = json.load(f)

    chk("1. results.json has config/sample/main/decomposition/sensitivity keys",
        all(k in r for k in ('config', 'sample', 'main',
                             'decomposition', 'sensitivity')))
    chk("2. report.md exists and non-empty",
        os.path.exists('report.md') and os.path.getsize('report.md') > 300)

    s = r['sample']
    chk("3. At least 10,000 metro-month observations",
        s.get('n_obs', 0) >= 10000)
    chk("4. At least 100 metros in final sample",
        s.get('n_metros', 0) >= 100)
    chk("5. At least 80 months spanned",
        s.get('n_months', 0) >= 80)
    chk("6. Date range covers 2022+ treatment window",
        s.get('max_ym', 0) >= 202201)

    m = r['main']
    chk("7. Permutation p in [0,1]",
        0.0 <= m.get('permutation_p_two_sided', -1) <= 1.0)
    chk("8. |beta_spread_x_exposure| < 1 (log-listings response sanity)",
        abs(m.get('beta_spread_x_exposure', 99)) < 1.0)
    chk("9. |beta_spread| < 1 (sanity)",
        abs(m.get('beta_spread', 99)) < 1.0)
    chk("10. R^2 in [0,1]",
        0.0 <= m.get('r2', -1) <= 1.0)
    ci = m.get('bootstrap_ci', {})
    chk("11. Bootstrap CI for spread_x_exposure has nonzero width",
        ('spread_x_exposure' in ci
         and ci['spread_x_exposure']['hi'] > ci['spread_x_exposure']['lo']))
    chk("12. Bootstrap CI for spread has nonzero width",
        ('spread' in ci and ci['spread']['hi'] > ci['spread']['lo']))
    chk("13. Bootstrap CI width for interaction < 1.0 (sanity)",
        ('spread_x_exposure' in ci
         and (ci['spread_x_exposure']['hi']
              - ci['spread_x_exposure']['lo']) < 1.0))

    d = r['decomposition']
    chk("14. Decomposition share_explained is finite",
        isinstance(d.get('share_explained'), (int, float))
        and math.isfinite(d.get('share_explained', float('nan'))))
    chk("15. Decomposition has >=1000 pre and >=1000 post rows",
        d.get('n_pre', 0) >= 1000 and d.get('n_post', 0) >= 1000)
    chk("16. Decomposition: mean spread rose post-2022",
        d.get('post_mean_spread', 0) > d.get('pre_mean_spread', 1))

    sens = r['sensitivity']
    chk("17. Sensitivity: at least 3 effective-rate windows",
        sum(1 for k in sens if k.startswith('window_')) >= 3)
    chk("18. Sensitivity: at least 3 spread lags",
        sum(1 for k in sens if k.startswith('lag_')) >= 3)
    chk("19. Sensitivity: pre/post-COVID subsamples run",
        'sample_pre_covid' in sens and 'sample_post_covid' in sens)
    chk("20. Sensitivity: top-50 metros restriction",
        'top50' in sens)
    chk("21. Sensitivity: placebo (spread=0) yields exact-zero interaction",
        abs(sens.get('placebo_zero_spread', {})
            .get('beta_spread_x_exposure', 99)) < 1e-10)
    chk("22. Sensitivity: placebo (spread=0) yields exact-zero spread beta",
        abs(sens.get('placebo_zero_spread', {})
            .get('beta_spread', 99)) < 1e-10)

    cfg = r['config']
    chk("23. N_permutations >= 1000", cfg.get('n_permutations', 0) >= 1000)
    chk("24. N_bootstrap >= 1000", cfg.get('n_bootstrap', 0) >= 1000)
    chk("25. Seed recorded", cfg.get('seed') == 42)
    chk("26. FE specification documented",
        'fe_specification' in cfg and 'metro' in cfg['fe_specification'].lower())

    # --- Additional robustness, sensitivity-agreement, and falsification checks ---
    main_int = m.get('beta_spread_x_exposure', 0.0)
    main_spread = m.get('beta_spread', 0.0)

    # Effect-size plausibility: interaction coefficient is on log-outcome per
    # unit of (standardized exposure × spread in pct). A |β| < 0.2 is a
    # generous sanity bound; observed is ~0.02. This is analogous to a
    # Cohen's-d plausibility check for regression effects (|β_std| < 5 sd).
    chk("27. Effect-size plausibility: |beta_spread_x_exposure| < 0.2 (log "
        "response to a 1pp spread shock on a standardized metro is bounded)",
        abs(main_int) < 0.2)

    # Sensitivity-agreement: every effective-rate window variant and every
    # spread-lag variant must share the main-sample SIGN on the interaction
    # coefficient. Robust findings should not flip sign across reasonable
    # specification choices.
    window_signs_agree = all(
        (v.get('beta_spread_x_exposure', 0.0) * main_int) > 0
        for k, v in sens.items() if k.startswith('window_'))
    chk("28. Sensitivity: all effective-rate window variants agree with "
        "main on sign of beta_spread_x_exposure",
        window_signs_agree)

    lag_signs_agree = all(
        (v.get('beta_spread_x_exposure', 0.0) * main_int) > 0
        for k, v in sens.items() if k.startswith('lag_'))
    chk("29. Sensitivity: all spread-lag variants agree with main on sign "
        "of beta_spread_x_exposure",
        lag_signs_agree)

    # Alternate-vintage (15yr PMMS) sign agreement: an alternate data source
    # for the same shock should agree in sign if the finding is not an
    # artifact of the 30yr vintage.
    alt = sens.get('alt_pmms15', {})
    chk("30. Sensitivity: 15-yr PMMS alternate vintage agrees with main "
        "on sign of beta_spread_x_exposure",
        alt and (alt.get('beta_spread_x_exposure', 0.0) * main_int) > 0)

    # CI-width sanity: point estimate should be at least ~1% of CI width
    # (i.e., CI should not be absurdly collapsed), and CI width should be
    # less than 100% of the point estimate's magnitude times a large bound
    # (i.e., CI not absurdly wide).
    ci_int = ci.get('spread_x_exposure', {})
    ci_width = ci_int.get('hi', 0) - ci_int.get('lo', 0)
    chk("31. Bootstrap CI width for interaction is > 1% of |estimate| "
        "(guards against a numerically collapsed CI)",
        ci_width > 0.01 * abs(main_int) if abs(main_int) > 0 else True)
    chk("32. Bootstrap CI width for interaction is < 20x |estimate| "
        "(guards against an unusably wide CI)",
        ci_width < 20.0 * max(abs(main_int), 1e-6))

    # Decomposition share plausibility: share of decline must be a finite
    # real number within a generous [-5, 5] band — i.e., the skill is not
    # producing a runaway decomposition.
    share = d.get('share_explained', 99.0)
    chk("33. Decomposition share_explained in [-5, 5] (plausibility bound)",
        -5.0 <= share <= 5.0)

    # Limitations must be present and enumerate at least 4 distinct caveats.
    lims = r.get('limitations', [])
    chk("34. Results record at least 4 documented limitations/caveats",
        isinstance(lims, list) and len(lims) >= 4)

    # Pre-COVID sign is not required to agree with main (indeed the
    # economic story is that the lock-in channel is a post-2022 phenomenon)
    # but post-COVID should agree, as that is the main identifying
    # variation window.
    post_covid = sens.get('sample_post_covid', {})
    chk("35. Sensitivity: post-COVID subsample agrees with main on sign "
        "of beta_spread_x_exposure (identifying window)",
        post_covid and (post_covid.get('beta_spread_x_exposure', 0.0)
                        * main_int) > 0)

    # Negative control / falsification: the spread=0 placebo must yield
    # exact-zero interaction AND exact-zero level coefficients (strongest
    # falsification — a numerical bug would leak into both).
    placebo = sens.get('placebo_zero_spread', {})
    chk("36. Negative control (placebo): spread=0 yields exact-zero "
        "interaction AND level coefficients (< 1e-10 each)",
        abs(placebo.get('beta_spread', 99)) < 1e-10
        and abs(placebo.get('beta_spread_x_exposure', 99)) < 1e-10)

    # Data row count plausibility: the Realtor.com top-300-metro panel
    # covering Jan 2016 – Mar 2026 should yield neither too few nor too
    # many observations. An upper bound catches accidental cross-joins
    # or row-duplication bugs.
    n_obs_main = s.get('n_obs', 0)
    chk("37. Sample size upper bound: n_obs < 100,000 (no cross-join / "
        "duplication bug)",
        n_obs_main < 100000)

    # Permutation null central tendency: under H0 the expected value of
    # the null distribution of the interaction coefficient is 0. A large
    # deviation from 0 signals that the permutation scheme is not
    # properly stratifying the shock, i.e., the null is mis-specified.
    perm_null_mean = m.get('permutation_null_mean', 99.0)
    chk("38. Permutation null mean is close to zero (|mean| < 0.5 * "
        "|observed|, central tendency sanity)",
        abs(perm_null_mean) < 0.5 * max(abs(main_int), 1e-6))

    # Permutation null dispersion: the null must have strictly positive
    # spread, otherwise the test has zero statistical power.
    perm_null_std = m.get('permutation_null_std', 0.0)
    chk("39. Permutation null has strictly positive std (non-degenerate "
        "null distribution)",
        perm_null_std > 0.0)

    # Standardized observed effect vs. null (analog to Cohen's-d for
    # regression): |β_obs - null_mean| / null_std should exceed 2 for a
    # significant finding at two-sided α=0.05 under an approximately-
    # normal null. Bounded above to catch runaway estimates.
    perm_obs_z = m.get('permutation_obs_z', 0.0)
    chk("40. Observed |z| (effect vs. null std) in plausible range "
        "[0, 50] — Cohen's-d analogue, finite but not astronomical",
        0.0 <= abs(perm_obs_z) <= 50.0)

    # Full sensitivity-sign concordance: every variant (windows, lags,
    # alt PMMS vintage, post-COVID, drop-COVID-year, top-50 metros) must
    # share the main-sample sign. Pre-COVID is the one expected
    # exception (economic story is post-2022). This is a stronger,
    # aggregated version of individual checks 28–30/35.
    drop_covid = sens.get('sample_drop_covid_year', {})
    drop_covid_ok = (drop_covid
                     and (drop_covid.get('beta_spread_x_exposure', 0.0)
                          * main_int) > 0)
    top50 = sens.get('top50', {})
    top50_ok = (top50
                and (top50.get('beta_spread_x_exposure', 0.0)
                     * main_int) > 0)
    chk("41. Sensitivity: drop-COVID-year AND top-50 metro variants "
        "both agree with main on sign of beta_spread_x_exposure",
        drop_covid_ok and top50_ok)

    # Permutation p and bootstrap CI must agree on significance: if p
    # < 0.05 then the CI should not straddle zero. Inconsistency signals
    # either an under-specified null or a degenerate bootstrap.
    ci_spans_zero = (ci_int.get('lo', 0) <= 0 <= ci_int.get('hi', 0))
    chk("42. Inference agreement: permutation p<0.05 and bootstrap CI "
        "does not straddle zero (else the two procedures contradict)",
        (m.get('permutation_p_two_sided', 1.0) >= 0.05)
        or (not ci_spans_zero))

    print(f"\n{ok}/{ok + fail} checks passed")
    if fail:
        print("VERIFICATION FAILED")
        sys.exit(1)
    else:
        print("ALL CHECKS PASSED")
        sys.exit(0)


if __name__ == '__main__':
    main()
SCRIPT_EOF
```

**Expected output:** File `analyze.py` written, exit code 0.

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_mortgage-lock-in-and-listings && python3 analyze.py
```

**Expected output:**
- Prints `[1/10]` through `[10/10]` progress sections.
- Downloads `PMMS_history.csv` (~2,800 weekly rows) and the
  Realtor.com metro history CSV (~108,000 metro-month rows).
- Builds a top-300-metro panel filtered to rows with non-imputed
  quality flag, ≥60 months of coverage, and a valid 10-year trailing
  effective mortgage rate.
- Fits metro-FE + calendar-month-of-year-FE within-OLS of log
  new-listings on (spread, spread×exposure, linear trend).
- Runs 2,000 seasonal-block permutations and 2,000 metro-clustered
  bootstrap replicates.
- Runs six sensitivity analyses: effective-rate window (60/120/180
  mo), spread lag (0/3/6 mo), pre/post-COVID subsamples, top-50-metro
  restriction, alternate 15-year PMMS vintage, and a spread=0 placebo.
- Writes `results.json` and `report.md`.
- Final line: `ANALYSIS COMPLETE`.
- Runtime: 15–30 minutes first run, 10–20 minutes on rerun (cache hit).
- Exit code 0.

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_mortgage-lock-in-and-listings && python3 analyze.py --verify
```

**Expected output:**
- 42/42 checks passed
- `ALL CHECKS PASSED`
- Exit code 0

## Expected Outputs

| File | Description |
|---|---|
| `results.json` | Config, sample sizes, main panel coefficients with bootstrap CIs, permutation p, counterfactual decomposition, full sensitivity table |
| `report.md` | Human-readable Markdown report with tables |
| `cache/pmms_history.csv` | Freddie Mac PMMS CSV (SHA256-pinned) |
| `cache/realtor_metro.csv` | Realtor.com metro history CSV (SHA256-pinned) |
| `data_manifest.json` | SHA256 hashes of every cached file |

## Success Criteria

The skill is considered a STRONG success iff **all** of the following
measurable conditions hold after a clean run:

1. Script exits 0 on both normal run and `--verify`.
2. ≥10,000 metro-month observations across ≥100 metros.
3. ≥80 months of coverage, including the 2022+ treatment window.
4. 2,000 seasonal-block permutation iterations and 2,000 metro-clustered
   bootstrap replicates are performed.
5. Placebo (spread=0) sensitivity coefficient is within 1e-10 of zero
   (mechanical null check — both level and interaction).
6. Bootstrap 95% CI for `β_interaction` has strictly positive width and
   does not straddle zero in the main sample.
7. Permutation two-sided p-value < 0.05 for `β_interaction`.
8. All effective-rate-window and spread-lag sensitivity variants agree
   with the main sample on the SIGN of `β_interaction` (robustness).
9. Counterfactual `share_explained` is finite and within the plausibility
   band [-5, 5].
10. `results.json` records at least 4 limitations caveating interpretation.
11. All 42 verification assertions pass (`python3 analyze.py --verify`
    ends with `ALL CHECKS PASSED`). Assertion categories covered:
    (i) output-file presence, (ii) sample-size lower AND upper bounds,
    (iii) coefficient range / effect-size plausibility (Cohen's-d
    analogue), (iv) R² range, (v) bootstrap CI width bounds,
    (vi) permutation null central tendency + dispersion +
    p-value range, (vii) decomposition share plausibility,
    (viii) sensitivity-sign concordance across six variants,
    (ix) placebo (spread=0) mechanical-null falsification,
    (x) inference agreement between permutation p and bootstrap CI,
    and (xi) config-provenance checks (seed, N, FE spec, limitations).
12. Expected SHA256 of the two data files (recorded in
    `EXPECTED_SHA256` inside the script and echoed into `results.json`):
    - `cache/pmms_history.csv`:
      `0cf5fa2deb13990ea5b0fc87ad33267002484fae2c6d176726cfc95c79e1b831`
    - `cache/realtor_metro.csv`:
      `ae128d1d67a519fef40f1f937e0eb8acc155b5235f3fb4ea33b67c44f0d793ac`
    A mismatch is logged at runtime but does not fail the run unless
    `STRICT_SHA256=True` is set in the DOMAIN CONFIGURATION block.

## Failure Conditions

The skill is considered a FAILURE if **any** of the following occur:

1. Import errors (script uses only Python 3.8+ standard library; no
   numpy/scipy/pandas).
2. `PMMS_URL` or `REALTOR_METRO_URL` unreachable after 4 retries — the
   script exits with a clear stderr message and nonzero exit code (11).
3. Filesystem permission or disk-space error creating `cache/` — exit
   code 10 with clear stderr guidance.
4. Data schema drift that leaves fewer than 120 PMMS months or 10,000
   Realtor.com rows after parsing — exit code 14 with diagnostic
   message identifying the observed counts.
5. Fewer than 10,000 observations or fewer than 100 metros after
   filtering.
6. Any `--verify` assertion fails (exit code 1 with `VERIFICATION
   FAILED`).
7. `results.json` or `report.md` not created.
8. SHA256 hash drift on cached data AND `STRICT_SHA256=True` (exit
   code 12 with drift diagnostic); by default (STRICT_SHA256=False) a
   mismatch is only logged.
9. Any sensitivity variant (window, lag, alt-vintage) flips the sign of
   `β_interaction` relative to the main sample — this does not fail the
   run but fails verification checks 28/29/30, signalling that the
   finding is NOT robust.
10. Permutation p ≥ 0.05 — this is a substantive failure of the skill's
    null model (not a script failure), and downstream interpretation
    should be "rate lock-in signal is not distinguishable from seasonal-
    block noise".