---
name: "boe-historical-inflation-autocorrelation"
description: "Can 18th-century inflation be reconstructed from the Bank of England's historical price indices, and does it exhibit the same autocorrelation structure as modern CPI? Downloads BoE Millennium Dataset (prices since 1209), computes rolling inflation rates, tests for structural breaks using permutation-based Chow test (2000 shuffles), and compares autocorrelation functions pre-1900 vs post-1900."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget"
tags: ["economic-history", "time-series", "autocorrelation", "structural-breaks", "permutation-test", "claw4s-2026"]
python_version: ">=3.8"
dependencies: []
estimated_runtime: "< 60 seconds on cached data (includes ~15s xlsx parsing, ~2s permutations, ~2s bootstrap)"
data_source: "Bank of England Millennium Macroeconomic Dataset v3.1 (prices since 1209)"
---

# Can Pre-Industrial Inflation Be Modeled with Modern Econometric Tools?

## Motivation

Modern econometric models of inflation (AR, ARIMA, Phillips curve regressions) are
calibrated on post-WWII CPI data. But central banks and economic historians routinely
apply these same models to centuries-old price data. If pre-industrial inflation had
fundamentally different time-series properties — lower persistence, different volatility
clustering, distinct structural breaks — then modern models may produce misleading
conclusions when applied to historical data.

The **methodological hook**: we use a permutation-based Chow test (not the classical
F-test which assumes normality and homoscedasticity) to detect structural breaks in
autocorrelation structure, and bootstrap confidence intervals to quantify the difference
in persistence between historical and modern inflation regimes.

## Prerequisites

- Python 3.8 or later (standard library only — zero pip dependencies)
- Network access for initial data download from Bank of England (cached after first run)
- ~15 MB free disk space

## Step 1: Create workspace

```bash
mkdir -p boe_inflation/cache
```

**Expected:** Exit code 0. Directory `boe_inflation/cache/` exists.

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > boe_inflation/analyze.py
#!/usr/bin/env python3
"""
Bank of England Historical Inflation Autocorrelation Analysis

Reconstructs inflation from the BoE Millennium Dataset (prices since 1209)
and tests whether pre-industrial inflation exhibits the same autocorrelation
structure as modern CPI.

Statistical methods:
- Permutation-based Chow test (2000 shuffles) for structural breaks in AR(1)
- Bootstrap confidence intervals (2000 resamples) for autocorrelation differences
- Fisher z-transform for formal correlation comparison
- Sensitivity analysis across multiple break points and rolling windows

Standard library only. No external dependencies.
"""

import os
import sys
import json
import math
import random
import hashlib
import time
import re
import zipfile
import xml.etree.ElementTree as ET
import urllib.request
import urllib.error
import ssl
import subprocess

###############################################################################
# Configuration
###############################################################################
SEED = 42
N_PERMUTATIONS = 2000
N_BOOTSTRAP = 2000
CACHE_DIR = "cache"
DATA_URL = "https://www.bankofengland.co.uk/-/media/boe/files/statistics/research-datasets/a-millennium-of-macroeconomic-data-for-the-uk.xlsx"
DATA_FILE = os.path.join(CACHE_DIR, "boe_millennium.xlsx")
BREAK_YEAR = 1900
MAX_ACF_LAG = 15
ROLLING_WINDOWS = [10, 20, 30]
SENSITIVITY_BREAK_YEARS = [1750, 1800, 1850, 1900, 1914, 1945, 1971]

###############################################################################
# Math helpers (stdlib only)
###############################################################################
def mean(xs):
    if not xs:
        return 0.0
    return sum(xs) / len(xs)

def variance(xs, ddof=1):
    if len(xs) <= ddof:
        return 0.0
    m = mean(xs)
    return sum((x - m) ** 2 for x in xs) / (len(xs) - ddof)

def std(xs, ddof=1):
    v = variance(xs, ddof)
    return math.sqrt(v) if v > 0 else 0.0

def covariance(xs, ys, ddof=1):
    n = min(len(xs), len(ys))
    if n <= ddof:
        return 0.0
    mx, my = mean(xs[:n]), mean(ys[:n])
    return sum((xs[i] - mx) * (ys[i] - my) for i in range(n)) / (n - ddof)

def correlation(xs, ys):
    sx, sy = std(xs), std(ys)
    if sx == 0 or sy == 0:
        return 0.0
    return covariance(xs, ys) / (sx * sy)

def median(xs):
    s = sorted(xs)
    n = len(s)
    if n % 2 == 1:
        return s[n // 2]
    return (s[n // 2 - 1] + s[n // 2]) / 2

def normal_cdf(x):
    """Standard normal CDF via error function."""
    return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0)))

def percentile(xs, p):
    """Compute the p-th percentile (0-100) of a sorted list."""
    s = sorted(xs)
    k = (p / 100.0) * (len(s) - 1)
    f = math.floor(k)
    c = math.ceil(k)
    if f == c:
        return s[int(k)]
    return s[f] * (c - k) + s[c] * (k - f)

###############################################################################
# Autocorrelation
###############################################################################
def acf(series, max_lag):
    """Compute autocorrelation function up to max_lag."""
    n = len(series)
    m = mean(series)
    denom = sum((x - m) ** 2 for x in series)
    if denom == 0:
        return [1.0] + [0.0] * max_lag
    result = [1.0]
    for lag in range(1, max_lag + 1):
        if lag >= n:
            result.append(0.0)
            continue
        numer = sum((series[t] - m) * (series[t - lag] - m) for t in range(lag, n))
        result.append(numer / denom)
    return result

###############################################################################
# AR(1) Model Fitting
###############################################################################
def ar1_fit(series):
    """Fit AR(1): x_t = a + b*x_{t-1} + e_t. Returns (a, b, rss, residuals)."""
    if len(series) < 4:
        return 0.0, 0.0, float('inf'), []
    x = series[:-1]
    y = series[1:]
    vx = variance(x, ddof=0)
    if vx == 0:
        a = mean(y)
        residuals = [yi - a for yi in y]
        rss = sum(r ** 2 for r in residuals)
        return a, 0.0, rss, residuals
    b = covariance(x, y, ddof=0) / vx
    a = mean(y) - b * mean(x)
    residuals = [y[i] - a - b * x[i] for i in range(len(x))]
    rss = sum(r ** 2 for r in residuals)
    return a, b, rss, residuals

###############################################################################
# Permutation-based Chow Test
###############################################################################
def chow_f_statistic(series, break_idx):
    """Compute Chow F-statistic for AR(1) structural break at break_idx."""
    if break_idx < 5 or break_idx > len(series) - 5:
        return 0.0
    _, _, rss_full, _ = ar1_fit(series)
    _, _, rss_pre, _ = ar1_fit(series[:break_idx])
    _, _, rss_post, _ = ar1_fit(series[break_idx:])
    k = 2  # parameters per sub-model
    n = len(series) - 1  # AR(1) observations
    denom_df = n - 2 * k
    if denom_df <= 0:
        return 0.0
    numer = (rss_full - rss_pre - rss_post) / k
    denom = (rss_pre + rss_post) / denom_df
    if denom <= 0:
        return 0.0
    return numer / denom

def permutation_chow_test(series, break_idx, n_perms, rng):
    """Permutation-based Chow test. Returns (observed_F, p_value, n_extreme)."""
    observed_f = chow_f_statistic(series, break_idx)
    count_extreme = 0
    for i in range(n_perms):
        perm = series[:]
        rng.shuffle(perm)
        perm_f = chow_f_statistic(perm, break_idx)
        if perm_f >= observed_f:
            count_extreme += 1
    p_value = (count_extreme + 1) / (n_perms + 1)
    return observed_f, p_value, count_extreme

###############################################################################
# Bootstrap Confidence Intervals
###############################################################################
def bootstrap_stat(data, stat_func, n_boot, rng, ci_level=0.95):
    """Bootstrap CI for stat_func(data). Returns (point_est, lo, hi, dist)."""
    n = len(data)
    point_est = stat_func(data)
    boot_stats = []
    for _ in range(n_boot):
        sample = [data[rng.randint(0, n - 1)] for _ in range(n)]
        boot_stats.append(stat_func(sample))
    boot_stats.sort()
    alpha = 1 - ci_level
    lo = boot_stats[max(0, int(math.floor(alpha / 2 * n_boot)))]
    hi = boot_stats[min(n_boot - 1, int(math.ceil((1 - alpha / 2) * n_boot)) - 1)]
    return point_est, lo, hi, boot_stats

def acf1_func(series):
    """Extract lag-1 autocorrelation."""
    if len(series) < 4:
        return 0.0
    return acf(series, 1)[1]

###############################################################################
# Fisher z-transform
###############################################################################
def fisher_z_test(r1, n1, r2, n2):
    """Compare two correlations using Fisher z-transform."""
    def arctanh_safe(r):
        r = max(-0.9999, min(0.9999, r))
        return 0.5 * math.log((1 + r) / (1 - r))
    z1 = arctanh_safe(r1)
    z2 = arctanh_safe(r2)
    se = math.sqrt(1.0 / max(1, n1 - 3) + 1.0 / max(1, n2 - 3))
    if se == 0:
        return 0.0, 1.0
    z_stat = (z1 - z2) / se
    p_value = 2.0 * (1.0 - normal_cdf(abs(z_stat)))
    return z_stat, p_value

###############################################################################
# Effect Sizes
###############################################################################
def cohens_d(xs, ys):
    """Cohen's d for two independent samples."""
    nx, ny = len(xs), len(ys)
    mx, my = mean(xs), mean(ys)
    vx, vy = variance(xs), variance(ys)
    pooled_var = ((nx - 1) * vx + (ny - 1) * vy) / (nx + ny - 2) if (nx + ny) > 2 else 1.0
    if pooled_var <= 0:
        return 0.0
    return (mx - my) / math.sqrt(pooled_var)

###############################################################################
# Data Download and Caching
###############################################################################
def download_with_retry(url, filepath, max_retries=3):
    """Download file with retry logic and local caching."""
    os.makedirs(os.path.dirname(filepath) or '.', exist_ok=True)
    if os.path.exists(filepath):
        size = os.path.getsize(filepath)
        print(f"  Using cached file: {filepath} ({size:,} bytes)")
        return
    # Browser-like User-Agent required by some CDNs (e.g., Bank of England)
    UA = ('Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 '
          '(KHTML, like Gecko) Chrome/120.0.0.0 Safari/537.36')
    # Try curl first (more reliable with CDN redirects and TLS)
    try:
        print(f"  Downloading via curl...")
        proc = subprocess.run(
            ['curl', '-sL', '-H', f'User-Agent: {UA}', '-o', filepath, url],
            timeout=300, capture_output=True)
        if proc.returncode == 0 and os.path.exists(filepath) and os.path.getsize(filepath) > 1000:
            print(f"  Downloaded {os.path.getsize(filepath):,} bytes via curl")
            return
        else:
            print(f"  curl failed (rc={proc.returncode}), falling back to urllib...")
            if os.path.exists(filepath):
                os.remove(filepath)
    except (FileNotFoundError, subprocess.TimeoutExpired):
        print(f"  curl not available, using urllib...")

    # Fallback to urllib with browser UA
    ctx = ssl.create_default_context()
    for attempt in range(max_retries):
        try:
            print(f"  Downloading via urllib (attempt {attempt + 1}/{max_retries})...")
            req = urllib.request.Request(url, headers={
                'User-Agent': UA,
                'Accept': '*/*',
                'Accept-Language': 'en-US,en;q=0.9'
            })
            with urllib.request.urlopen(req, timeout=180, context=ctx) as resp:
                data = resp.read()
            with open(filepath, 'wb') as f:
                f.write(data)
            print(f"  Downloaded {len(data):,} bytes via urllib")
            return
        except (urllib.error.URLError, OSError, ssl.SSLError) as e:
            print(f"  Attempt {attempt + 1} failed: {e}")
            if attempt < max_retries - 1:
                wait = 2 ** (attempt + 1)
                print(f"  Retrying in {wait}s...")
                time.sleep(wait)
    raise RuntimeError(f"Failed to download {url} after {max_retries} attempts")

def sha256_file(filepath):
    """Compute SHA256 hash of a file."""
    h = hashlib.sha256()
    with open(filepath, 'rb') as f:
        for chunk in iter(lambda: f.read(65536), b''):
            h.update(chunk)
    return h.hexdigest()

###############################################################################
# XLSX Parser (standard library only: zipfile + xml.etree.ElementTree)
###############################################################################
def parse_xlsx_sheet(zf, sheet_path, shared_strings, ns):
    """Parse a single worksheet from an xlsx zip file."""
    if sheet_path not in zf.namelist():
        return []
    tree = ET.parse(zf.open(sheet_path))
    rows_out = []
    for row_el in tree.iter(f'{ns}row'):
        row_num = int(row_el.get('r', '0'))
        cells = {}
        for cell_el in row_el.iter(f'{ns}c'):
            ref = cell_el.get('r', '')
            col = ''.join(c for c in ref if c.isalpha())
            cell_type = cell_el.get('t', '')
            val_el = cell_el.find(f'{ns}v')
            if val_el is not None and val_el.text is not None:
                raw = val_el.text
                if cell_type == 's':
                    idx = int(raw)
                    value = shared_strings[idx] if idx < len(shared_strings) else ''
                elif cell_type == 'b':
                    value = bool(int(raw))
                else:
                    try:
                        value = float(raw)
                        if value == int(value) and '.' not in raw and 'E' not in raw.upper():
                            value = int(value)
                    except ValueError:
                        value = raw
            else:
                # Check for inline string
                is_el = cell_el.find(f'{ns}is')
                if is_el is not None:
                    t_el = is_el.find(f'{ns}t')
                    value = t_el.text if t_el is not None and t_el.text else ''
                else:
                    value = None
            cells[col] = value
        if cells:
            rows_out.append((row_num, cells))
    rows_out.sort(key=lambda x: x[0])
    return rows_out

def parse_xlsx(filepath):
    """Parse xlsx into {sheet_name: [(row_num, {col: value}), ...]}."""
    NS = '{http://schemas.openxmlformats.org/spreadsheetml/2006/main}'
    RELS_NS = '{http://schemas.openxmlformats.org/package/2006/relationships}'

    with zipfile.ZipFile(filepath, 'r') as zf:
        # Shared strings
        shared_strings = []
        if 'xl/sharedStrings.xml' in zf.namelist():
            ss_tree = ET.parse(zf.open('xl/sharedStrings.xml'))
            for si in ss_tree.iter(f'{NS}si'):
                parts = []
                for t in si.iter(f'{NS}t'):
                    if t.text:
                        parts.append(t.text)
                shared_strings.append(''.join(parts))

        # Workbook: sheet names
        wb_tree = ET.parse(zf.open('xl/workbook.xml'))
        sheets_info = []
        for sheet_el in wb_tree.iter(f'{NS}sheet'):
            name = sheet_el.get('name', '')
            rid = sheet_el.get('{http://schemas.openxmlformats.org/officeDocument/2006/relationships}id', '')
            sheets_info.append((name, rid))

        # Relationships
        rid_map = {}
        rels_path = 'xl/_rels/workbook.xml.rels'
        if rels_path in zf.namelist():
            rels_tree = ET.parse(zf.open(rels_path))
            for rel in rels_tree.iter(f'{RELS_NS}Relationship'):
                target = rel.get('Target', '')
                if not target.startswith('/'):
                    target = 'xl/' + target
                rid_map[rel.get('Id', '')] = target

        workbook = {}
        for name, rid in sheets_info:
            path = rid_map.get(rid)
            if not path:
                continue
            rows = parse_xlsx_sheet(zf, path, shared_strings, NS)
            workbook[name] = rows
    return workbook

###############################################################################
# Data Extraction
###############################################################################
def extract_price_series(workbook):
    """Find and extract annual price-level data from the BoE Millennium Dataset.

    Strategy: Target sheet 'A47. Wages and prices' which contains the CPI
    preferred measure (column D) from 1209-2016. Column A has years as integers.
    Falls back to pattern-based detection on all sheets if A47 is not found.
    """
    all_sheets = list(workbook.keys())
    print(f"  Available sheets ({len(all_sheets)}): {all_sheets[:10]}{'...' if len(all_sheets) > 10 else ''}")

    # Priority order: A47 (has CPI 1209-2016), then A1, then others with price keywords
    priority_sheets = []
    other_sheets = []
    for name in all_sheets:
        nl = name.lower()
        if 'a47' in nl or 'wages and prices' in nl:
            priority_sheets.insert(0, name)
        elif 'a1' in nl and 'headline' in nl:
            priority_sheets.append(name)
        elif any(k in nl for k in ['price', 'cpi', 'rpi', 'inflation']):
            priority_sheets.append(name)
        else:
            other_sheets.append(name)
    candidate_sheets = priority_sheets + other_sheets

    # Step 1: Header-based detection (look for CPI/price columns)
    price_keywords = [
        'consumer price index', 'cpi', 'price index', 'composite price',
        'cost of living', 'retail price', 'rpi', 'gdp deflator',
        'prices', 'price level'
    ]

    for sheet_name in candidate_sheets:
        rows = workbook[sheet_name]
        if len(rows) < 50:
            continue

        # Scan first 10 rows for ALL price-keyword column matches
        # Rank: "prefer" > "consumer price index" > "cpi" > other keywords
        best_match = None  # (priority, h_idx, col, label)
        for h_idx in range(min(10, len(rows))):
            _, header_cells = rows[h_idx]
            for col, val in header_cells.items():
                if val is None:
                    continue
                vs = str(val).lower().strip()
                for pk in price_keywords:
                    if pk in vs:
                        # Assign priority: lower = better
                        if 'prefer' in vs:
                            prio = 0
                        elif 'consumer price index' in vs:
                            prio = 1
                        elif pk == 'cpi':
                            prio = 2
                        else:
                            prio = 3
                        if best_match is None or prio < best_match[0]:
                            best_match = (prio, h_idx, col, str(val).strip())
                        break

        if best_match is None:
            continue

        _, best_h_idx, price_col, price_label = best_match

        # Year column: try column A (most common in BoE dataset)
        year_col = 'A'
        years = []
        prices = []
        # Start extracting from the row after the best header match
        for r_idx in range(best_h_idx + 1, len(rows)):
            _, cells = rows[r_idx]
            y_raw = cells.get(year_col)
            p_raw = cells.get(price_col)
            if y_raw is None or p_raw is None:
                continue
            try:
                yr = int(float(y_raw)) if isinstance(y_raw, (int, float)) else int(str(y_raw).strip())
                pr = float(p_raw)
                if 1200 <= yr <= 2025 and pr > 0:
                    years.append(yr)
                    prices.append(pr)
            except (ValueError, TypeError):
                continue

        if len(years) >= 200:
            print(f"  Found series in sheet '{sheet_name}': column '{price_label}'")
            print(f"  Year range: {min(years)}-{max(years)}, N={len(years)}")
            return years, prices, sheet_name, price_label

    # Step 2: Pattern-based fallback
    print("  Header search failed; trying pattern-based detection...")
    best_result = None
    best_span = 0
    for sheet_name in candidate_sheets:
        rows = workbook[sheet_name]
        if len(rows) < 50:
            continue
        all_cols = set()
        for _, cells in rows:
            all_cols.update(cells.keys())
        all_cols = sorted(all_cols)

        for yc in all_cols:
            y_vals = []
            for _, cells in rows:
                v = cells.get(yc)
                if v is not None:
                    try:
                        iv = int(float(v))
                        if 1200 <= iv <= 2025:
                            y_vals.append(iv)
                    except (ValueError, TypeError):
                        pass
            if len(y_vals) < 100:
                continue
            for pc in all_cols:
                if pc == yc:
                    continue
                pairs = []
                for _, cells in rows:
                    yv = cells.get(yc)
                    pv = cells.get(pc)
                    if yv is None or pv is None:
                        continue
                    try:
                        yr = int(float(yv))
                        pr = float(pv)
                        if 1200 <= yr <= 2025 and pr > 0:
                            pairs.append((yr, pr))
                    except (ValueError, TypeError):
                        continue
                if len(pairs) >= 100:
                    span = max(p[0] for p in pairs) - min(p[0] for p in pairs)
                    if span > best_span:
                        best_span = span
                        pairs.sort()
                        best_result = (
                            [p[0] for p in pairs], [p[1] for p in pairs],
                            sheet_name, f"Column {pc}"
                        )
    if best_result:
        years, prices, sname, plabel = best_result
        print(f"  Pattern-detected series in '{sname}': {plabel}")
        print(f"  Year range: {min(years)}-{max(years)}, N={len(years)}")
        return best_result

    raise RuntimeError("Could not find suitable price-level series in the dataset")

###############################################################################
# Inflation Computation
###############################################################################
def compute_inflation(years, prices):
    """Compute year-over-year inflation rates for consecutive years."""
    inf_years = []
    inf_rates = []
    for i in range(1, len(years)):
        if years[i] == years[i - 1] + 1 and prices[i - 1] > 0:
            rate = (prices[i] - prices[i - 1]) / prices[i - 1]
            inf_years.append(years[i])
            inf_rates.append(rate)
    return inf_years, inf_rates

def rolling_mean_std(values, window):
    """Compute rolling mean and std with given window size."""
    results = []
    for i in range(len(values) - window + 1):
        chunk = values[i:i + window]
        results.append((mean(chunk), std(chunk)))
    return results

###############################################################################
# Main Analysis
###############################################################################
def run_analysis():
    """Execute the full analysis pipeline."""
    results = {}
    rng = random.Random(SEED)

    # =========================================================================
    print("=" * 72)
    print("[1/8] Downloading and caching BoE Millennium Dataset")
    print("=" * 72)
    download_with_retry(DATA_URL, DATA_FILE)
    file_hash = sha256_file(DATA_FILE)
    file_size = os.path.getsize(DATA_FILE)
    print(f"  SHA256: {file_hash}")
    print(f"  Size:   {file_size:,} bytes")
    results['data'] = {
        'url': DATA_URL,
        'sha256': file_hash,
        'size_bytes': file_size,
        'cache_path': DATA_FILE
    }

    # Store hash for verification
    hash_file = os.path.join(CACHE_DIR, "boe_millennium.sha256")
    if os.path.exists(hash_file):
        with open(hash_file) as f:
            expected_hash = f.read().strip()
        if expected_hash != file_hash:
            print(f"  WARNING: Hash mismatch! Expected {expected_hash}, got {file_hash}")
            print(f"  Dataset may have been updated. Results may differ from original run.")
        else:
            print(f"  Hash verified: matches cached checksum")
    else:
        with open(hash_file, 'w') as f:
            f.write(file_hash)
        print(f"  Stored hash for future verification")

    # =========================================================================
    print()
    print("=" * 72)
    print("[2/8] Parsing XLSX and extracting price-level series")
    print("=" * 72)
    t0 = time.time()
    workbook = parse_xlsx(DATA_FILE)
    print(f"  Parsed {len(workbook)} sheets in {time.time() - t0:.1f}s")

    years, prices, sheet_name, price_label = extract_price_series(workbook)

    # Sort by year and remove duplicates (keep first)
    paired = sorted(zip(years, prices))
    seen = set()
    clean_years, clean_prices = [], []
    for y, p in paired:
        if y not in seen:
            seen.add(y)
            clean_years.append(y)
            clean_prices.append(p)
    years, prices = clean_years, clean_prices

    results['series'] = {
        'sheet': sheet_name,
        'column': price_label,
        'n_years': len(years),
        'year_min': min(years),
        'year_max': max(years),
        'price_min': round(min(prices), 4),
        'price_max': round(max(prices), 4)
    }
    print(f"  Clean series: {len(years)} years, {min(years)}-{max(years)}")

    # =========================================================================
    print()
    print("=" * 72)
    print("[3/8] Computing annual inflation rates")
    print("=" * 72)
    inf_years, inf_rates = compute_inflation(years, prices)
    print(f"  Inflation series: {len(inf_rates)} annual observations")
    print(f"  Year range: {inf_years[0]}-{inf_years[-1]}")
    print(f"  Mean inflation: {mean(inf_rates)*100:.2f}%")
    print(f"  Median inflation: {median(inf_rates)*100:.2f}%")
    print(f"  Std deviation: {std(inf_rates)*100:.2f}%")
    print(f"  Min: {min(inf_rates)*100:.2f}%, Max: {max(inf_rates)*100:.2f}%")

    results['inflation_summary'] = {
        'n_obs': len(inf_rates),
        'year_range': [inf_years[0], inf_years[-1]],
        'mean_pct': round(mean(inf_rates) * 100, 4),
        'median_pct': round(median(inf_rates) * 100, 4),
        'std_pct': round(std(inf_rates) * 100, 4),
        'min_pct': round(min(inf_rates) * 100, 4),
        'max_pct': round(max(inf_rates) * 100, 4)
    }

    # =========================================================================
    print()
    print("=" * 72)
    print("[4/8] Computing rolling inflation statistics")
    print("=" * 72)
    rolling_results = {}
    for w in ROLLING_WINDOWS:
        roll = rolling_mean_std(inf_rates, w)
        start_idx = w - 1
        roll_years = inf_years[start_idx:start_idx + len(roll)]
        roll_means = [r[0] for r in roll]
        roll_stds = [r[1] for r in roll]
        rolling_results[f'window_{w}'] = {
            'n_windows': len(roll),
            'mean_of_means_pct': round(mean(roll_means) * 100, 4),
            'mean_of_stds_pct': round(mean(roll_stds) * 100, 4),
            'max_rolling_mean_pct': round(max(roll_means) * 100, 4),
            'min_rolling_mean_pct': round(min(roll_means) * 100, 4)
        }
        print(f"  Window={w}yr: {len(roll)} windows, "
              f"avg_mean={mean(roll_means)*100:.2f}%, avg_vol={mean(roll_stds)*100:.2f}%")
    results['rolling'] = rolling_results

    # =========================================================================
    print()
    print("=" * 72)
    print("[5/8] Permutation-based Chow test for structural breaks (2000 shuffles)")
    print("=" * 72)

    # Find break index in inflation series
    break_idx = None
    for i, y in enumerate(inf_years):
        if y >= BREAK_YEAR:
            break_idx = i
            break
    if break_idx is None:
        break_idx = len(inf_rates) // 2

    print(f"  Primary break year: {BREAK_YEAR} (index {break_idx}/{len(inf_rates)})")
    print(f"  Pre-break: {break_idx} obs, Post-break: {len(inf_rates) - break_idx} obs")
    print(f"  Running {N_PERMUTATIONS} permutations... (this may take several minutes)")

    t0 = time.time()
    obs_f, p_val, n_extreme = permutation_chow_test(inf_rates, break_idx, N_PERMUTATIONS, rng)
    elapsed = time.time() - t0

    print(f"  Observed F-statistic: {obs_f:.4f}")
    print(f"  Permutation p-value: {p_val:.6f} ({n_extreme}/{N_PERMUTATIONS} as extreme)")
    print(f"  Elapsed: {elapsed:.1f}s")

    results['chow_test'] = {
        'break_year': BREAK_YEAR,
        'break_index': break_idx,
        'pre_n': break_idx,
        'post_n': len(inf_rates) - break_idx,
        'observed_F': round(obs_f, 6),
        'p_value': round(p_val, 6),
        'n_permutations': N_PERMUTATIONS,
        'n_extreme': n_extreme,
        'significant_005': p_val < 0.05,
        'elapsed_seconds': round(elapsed, 1)
    }

    # =========================================================================
    print()
    print("=" * 72)
    print("[6/8] Autocorrelation analysis: pre-1900 vs post-1900")
    print("=" * 72)

    pre_inflation = [inf_rates[i] for i in range(len(inf_rates)) if inf_years[i] < BREAK_YEAR]
    post_inflation = [inf_rates[i] for i in range(len(inf_rates)) if inf_years[i] >= BREAK_YEAR]

    print(f"  Pre-{BREAK_YEAR}: {len(pre_inflation)} observations")
    print(f"  Post-{BREAK_YEAR}: {len(post_inflation)} observations")

    acf_pre = acf(pre_inflation, MAX_ACF_LAG)
    acf_post = acf(post_inflation, MAX_ACF_LAG)

    print(f"
  {'Lag':<6} {'Pre-1900':>10} {'Post-1900':>10} {'Diff':>10}")
    print(f"  {'-'*38}")
    for lag in range(MAX_ACF_LAG + 1):
        diff = acf_post[lag] - acf_pre[lag]
        print(f"  {lag:<6} {acf_pre[lag]:>10.4f} {acf_post[lag]:>10.4f} {diff:>10.4f}")

    # Fisher z-test for lag-1 autocorrelation
    z_stat, z_pval = fisher_z_test(
        acf_pre[1], len(pre_inflation),
        acf_post[1], len(post_inflation)
    )
    print(f"
  Fisher z-test for lag-1 ACF difference:")
    print(f"    z-statistic: {z_stat:.4f}")
    print(f"    p-value: {z_pval:.6f}")

    # AR(1) coefficients
    _, ar1_pre, _, _ = ar1_fit(pre_inflation)
    _, ar1_post, _, _ = ar1_fit(post_inflation)
    print(f"
  AR(1) coefficient pre-{BREAK_YEAR}:  {ar1_pre:.4f}")
    print(f"  AR(1) coefficient post-{BREAK_YEAR}: {ar1_post:.4f}")

    # Cohen's d for inflation levels
    cd = cohens_d(pre_inflation, post_inflation)
    print(f"
  Cohen's d (mean inflation): {cd:.4f}")

    results['acf_comparison'] = {
        'pre_period': f'<{BREAK_YEAR}',
        'post_period': f'>={BREAK_YEAR}',
        'pre_n': len(pre_inflation),
        'post_n': len(post_inflation),
        'acf_pre': {f'lag_{i}': round(v, 6) for i, v in enumerate(acf_pre)},
        'acf_post': {f'lag_{i}': round(v, 6) for i, v in enumerate(acf_post)},
        'lag1_pre': round(acf_pre[1], 6),
        'lag1_post': round(acf_post[1], 6),
        'lag1_diff': round(acf_post[1] - acf_pre[1], 6),
        'fisher_z_stat': round(z_stat, 6),
        'fisher_z_pval': round(z_pval, 6),
        'ar1_coef_pre': round(ar1_pre, 6),
        'ar1_coef_post': round(ar1_post, 6),
        'cohens_d_inflation': round(cd, 6),
        'pre_mean_pct': round(mean(pre_inflation) * 100, 4),
        'post_mean_pct': round(mean(post_inflation) * 100, 4),
        'pre_std_pct': round(std(pre_inflation) * 100, 4),
        'post_std_pct': round(std(post_inflation) * 100, 4)
    }

    # =========================================================================
    print()
    print("=" * 72)
    print("[7/8] Bootstrap confidence intervals (2000 resamples)")
    print("=" * 72)
    print(f"  Computing bootstrap CIs for ACF(1) in each period...")

    t0 = time.time()

    # Block bootstrap to preserve temporal dependence structure.
    # Block size chosen as ~ sqrt(n) following Politis & Romano (1994).
    def block_bootstrap_acf1(series, n_boot, rng, block_size=None):
        """Block bootstrap for ACF(1). Preserves temporal structure."""
        n = len(series)
        if block_size is None:
            block_size = max(3, int(math.sqrt(n)))
        n_blocks = (n + block_size - 1) // block_size
        boot_acf1s = []
        for _ in range(n_boot):
            sample = []
            while len(sample) < n:
                start = rng.randint(0, n - block_size)
                sample.extend(series[start:start + block_size])
            sample = sample[:n]
            boot_acf1s.append(acf1_func(sample))
        boot_acf1s.sort()
        lo = boot_acf1s[int(0.025 * n_boot)]
        hi = boot_acf1s[int(0.975 * n_boot) - 1]
        return boot_acf1s, lo, hi

    # Block bootstrap CI for pre-period ACF(1)
    pre_boot, pre_acf1_lo, pre_acf1_hi = block_bootstrap_acf1(
        pre_inflation, N_BOOTSTRAP, random.Random(SEED + 1))
    pre_acf1_point = acf_pre[1]
    print(f"  Pre-{BREAK_YEAR} ACF(1): {pre_acf1_point:.4f} [{pre_acf1_lo:.4f}, {pre_acf1_hi:.4f}]")

    # Block bootstrap CI for post-period ACF(1)
    post_boot, post_acf1_lo, post_acf1_hi = block_bootstrap_acf1(
        post_inflation, N_BOOTSTRAP, random.Random(SEED + 2))
    post_acf1_point = acf_post[1]
    print(f"  Post-{BREAK_YEAR} ACF(1): {post_acf1_point:.4f} [{post_acf1_lo:.4f}, {post_acf1_hi:.4f}]")

    # Bootstrap CI for ACF(1) difference via paired block bootstraps
    diff_boot = [post_boot[i] - pre_boot[i] for i in range(N_BOOTSTRAP)]
    diff_boot.sort()
    diff_lo = diff_boot[int(0.025 * N_BOOTSTRAP)]
    diff_hi = diff_boot[int(0.975 * N_BOOTSTRAP) - 1]
    diff_point = acf_post[1] - acf_pre[1]
    elapsed = time.time() - t0

    print(f"  Difference (post - pre): {diff_point:.4f} [{diff_lo:.4f}, {diff_hi:.4f}]")
    ci_excludes_zero = (diff_lo > 0 and diff_hi > 0) or (diff_lo < 0 and diff_hi < 0)
    print(f"  95% CI excludes zero: {ci_excludes_zero}")
    print(f"  Elapsed: {elapsed:.1f}s")

    # i.i.d. bootstrap for mean inflation difference (mean is robust to resampling)
    n_pre, n_post = len(pre_inflation), len(post_inflation)
    mean_diff_boot = []
    rng_mean = random.Random(SEED + 4)
    for _ in range(N_BOOTSTRAP):
        pre_s = [pre_inflation[rng_mean.randint(0, n_pre - 1)] for _ in range(n_pre)]
        post_s = [post_inflation[rng_mean.randint(0, n_post - 1)] for _ in range(n_post)]
        mean_diff_boot.append(mean(post_s) - mean(pre_s))
    mean_diff_boot.sort()
    md_lo = mean_diff_boot[int(0.025 * N_BOOTSTRAP)]
    md_hi = mean_diff_boot[int(0.975 * N_BOOTSTRAP) - 1]
    md_point = mean(post_inflation) - mean(pre_inflation)
    print(f"
  Mean inflation diff: {md_point*100:.2f}% [{md_lo*100:.2f}%, {md_hi*100:.2f}%]")

    # Block bootstrap for volatility (std) difference
    vol_diff_boot = []
    rng_vol = random.Random(SEED + 5)
    pre_blk = max(3, int(math.sqrt(n_pre)))
    post_blk = max(3, int(math.sqrt(n_post)))
    for _ in range(N_BOOTSTRAP):
        # Block-resample pre
        pre_s = []
        while len(pre_s) < n_pre:
            st = rng_vol.randint(0, n_pre - pre_blk)
            pre_s.extend(pre_inflation[st:st + pre_blk])
        pre_s = pre_s[:n_pre]
        # Block-resample post
        post_s = []
        while len(post_s) < n_post:
            st = rng_vol.randint(0, n_post - post_blk)
            post_s.extend(post_inflation[st:st + post_blk])
        post_s = post_s[:n_post]
        vol_diff_boot.append(std(post_s) - std(pre_s))
    vol_diff_boot.sort()
    vd_lo = vol_diff_boot[int(0.025 * N_BOOTSTRAP)]
    vd_hi = vol_diff_boot[int(0.975 * N_BOOTSTRAP) - 1]
    vd_point = std(post_inflation) - std(pre_inflation)
    print(f"  Volatility diff: {vd_point*100:.2f}% [{vd_lo*100:.2f}%, {vd_hi*100:.2f}%]")

    results['bootstrap'] = {
        'n_bootstrap': N_BOOTSTRAP,
        'acf1_pre': {'point': round(pre_acf1_point, 6), 'ci_lo': round(pre_acf1_lo, 6), 'ci_hi': round(pre_acf1_hi, 6)},
        'acf1_post': {'point': round(post_acf1_point, 6), 'ci_lo': round(post_acf1_lo, 6), 'ci_hi': round(post_acf1_hi, 6)},
        'acf1_diff': {'point': round(diff_point, 6), 'ci_lo': round(diff_lo, 6), 'ci_hi': round(diff_hi, 6), 'ci_excludes_zero': ci_excludes_zero},
        'mean_diff_pct': {'point': round(md_point * 100, 4), 'ci_lo': round(md_lo * 100, 4), 'ci_hi': round(md_hi * 100, 4)},
        'vol_diff_pct': {'point': round(vd_point * 100, 4), 'ci_lo': round(vd_lo * 100, 4), 'ci_hi': round(vd_hi * 100, 4)}
    }

    # =========================================================================
    print()
    print("=" * 72)
    print("[8/8] Sensitivity analysis across break points and parameters")
    print("=" * 72)

    sensitivity = {}

    # 8a: Sensitivity to break-year choice
    print("  Testing multiple candidate break years...")
    break_year_results = []
    for by in SENSITIVITY_BREAK_YEARS:
        bi = None
        for i, y in enumerate(inf_years):
            if y >= by:
                bi = i
                break
        if bi is None or bi < 10 or bi > len(inf_rates) - 10:
            continue
        pre = [inf_rates[i] for i in range(len(inf_rates)) if inf_years[i] < by]
        post = [inf_rates[i] for i in range(len(inf_rates)) if inf_years[i] >= by]
        if len(pre) < 10 or len(post) < 10:
            continue
        acf1_p = acf(pre, 1)[1]
        acf1_q = acf(post, 1)[1]
        ar1_p = ar1_fit(pre)[1]
        ar1_q = ar1_fit(post)[1]
        cd_val = cohens_d(pre, post)
        z_s, z_p = fisher_z_test(acf1_p, len(pre), acf1_q, len(post))
        entry = {
            'break_year': by,
            'pre_n': len(pre),
            'post_n': len(post),
            'acf1_pre': round(acf1_p, 6),
            'acf1_post': round(acf1_q, 6),
            'acf1_diff': round(acf1_q - acf1_p, 6),
            'ar1_pre': round(ar1_p, 6),
            'ar1_post': round(ar1_q, 6),
            'fisher_z': round(z_s, 4),
            'fisher_p': round(z_p, 6),
            'cohens_d': round(cd_val, 4)
        }
        break_year_results.append(entry)
        sig = "*" if z_p < 0.05 else ""
        print(f"    {by}: ACF1 pre={acf1_p:.3f} post={acf1_q:.3f} diff={acf1_q-acf1_p:+.3f} "
              f"Fisher p={z_p:.4f}{sig}")
    sensitivity['break_years'] = break_year_results

    # 8b: Sensitivity to ACF lag choice
    print("
  Testing ACF differences across multiple lags...")
    lag_results = []
    for lag in [1, 2, 3, 5, 10]:
        if lag > MAX_ACF_LAG:
            continue
        pre_acf_val = acf(pre_inflation, lag)[lag]
        post_acf_val = acf(post_inflation, lag)[lag]
        diff = post_acf_val - pre_acf_val
        lag_results.append({
            'lag': lag,
            'pre': round(pre_acf_val, 6),
            'post': round(post_acf_val, 6),
            'diff': round(diff, 6)
        })
        print(f"    Lag {lag}: pre={pre_acf_val:.4f} post={post_acf_val:.4f} diff={diff:+.4f}")
    sensitivity['acf_lags'] = lag_results

    # 8c: Scan for optimal break point (max Chow F-statistic without permutation)
    print("
  Scanning for optimal break point (max F-statistic)...")
    best_f, best_year = 0.0, inf_years[0]
    scan_results = []
    step = max(1, len(inf_rates) // 100)
    for idx in range(20, len(inf_rates) - 20, step):
        f_val = chow_f_statistic(inf_rates, idx)
        yr = inf_years[idx]
        scan_results.append((yr, f_val))
        if f_val > best_f:
            best_f = f_val
            best_year = yr
    print(f"  Optimal break year: {best_year} (F={best_f:.4f})")

    # Top 5 break points
    scan_results.sort(key=lambda x: -x[1])
    top5 = scan_results[:5]
    print(f"  Top 5 break points:")
    for yr, fv in top5:
        print(f"    {yr}: F={fv:.4f}")
    sensitivity['optimal_break'] = {
        'best_year': best_year,
        'best_F': round(best_f, 6),
        'top5': [{'year': yr, 'F': round(fv, 6)} for yr, fv in top5]
    }

    results['sensitivity'] = sensitivity

    # =========================================================================
    # Write results
    # =========================================================================
    print()
    print("=" * 72)
    print("Writing results...")
    print("=" * 72)

    with open('results.json', 'w') as f:
        json.dump(results, f, indent=2)
    print(f"  Wrote results.json ({os.path.getsize('results.json'):,} bytes)")

    # Write report
    write_report(results)
    print(f"  Wrote report.md ({os.path.getsize('report.md'):,} bytes)")

    print()
    print("=" * 72)
    print("ANALYSIS COMPLETE")
    print("=" * 72)

    return results

###############################################################################
# Report Generation
###############################################################################
def write_report(r):
    """Generate a human-readable Markdown report."""
    lines = []
    a = lines.append

    a("# Bank of England Historical Inflation: Autocorrelation Analysis")
    a("")
    a("## Data Source")
    a(f"- **URL**: {r['data']['url']}")
    a(f"- **SHA256**: `{r['data']['sha256']}`")
    a(f"- **Size**: {r['data']['size_bytes']:,} bytes")
    a(f"- **Sheet**: {r['series']['sheet']}")
    a(f"- **Column**: {r['series']['column']}")
    a(f"- **Years**: {r['series']['year_min']}-{r['series']['year_max']} (N={r['series']['n_years']})")
    a("")

    a("## Inflation Summary")
    s = r['inflation_summary']
    a(f"- Observations: {s['n_obs']} ({s['year_range'][0]}-{s['year_range'][1]})")
    a(f"- Mean: {s['mean_pct']:.2f}%")
    a(f"- Median: {s['median_pct']:.2f}%")
    a(f"- Std Dev: {s['std_pct']:.2f}%")
    a(f"- Range: [{s['min_pct']:.2f}%, {s['max_pct']:.2f}%]")
    a("")

    a("## Structural Break Test (Permutation-based Chow Test)")
    c = r['chow_test']
    a(f"- Break year: {c['break_year']}")
    a(f"- Pre-break: {c['pre_n']} obs, Post-break: {c['post_n']} obs")
    a(f"- Observed F: {c['observed_F']:.4f}")
    a(f"- Permutation p-value: {c['p_value']:.6f} ({c['n_permutations']} permutations)")
    a(f"- **Significant at 0.05**: {'Yes' if c['significant_005'] else 'No'}")
    a("")

    a("## Autocorrelation Comparison")
    ac = r['acf_comparison']
    a(f"| Metric | Pre-{r['chow_test']['break_year']} | Post-{r['chow_test']['break_year']} |")
    a("|--------|----------|-----------|")
    a(f"| N | {ac['pre_n']} | {ac['post_n']} |")
    a(f"| Mean inflation | {ac['pre_mean_pct']:.2f}% | {ac['post_mean_pct']:.2f}% |")
    a(f"| Std inflation | {ac['pre_std_pct']:.2f}% | {ac['post_std_pct']:.2f}% |")
    a(f"| ACF lag-1 | {ac['lag1_pre']:.4f} | {ac['lag1_post']:.4f} |")
    a(f"| AR(1) coef | {ac['ar1_coef_pre']:.4f} | {ac['ar1_coef_post']:.4f} |")
    a("")
    a(f"- ACF(1) difference: {ac['lag1_diff']:.4f}")
    a(f"- Fisher z-test: z={ac['fisher_z_stat']:.4f}, p={ac['fisher_z_pval']:.6f}")
    a(f"- Cohen's d (inflation levels): {ac['cohens_d_inflation']:.4f}")
    a("")

    a("## Bootstrap Confidence Intervals")
    b = r['bootstrap']
    a(f"- {b['n_bootstrap']} bootstrap resamples")
    a(f"- Pre ACF(1): {b['acf1_pre']['point']:.4f} [{b['acf1_pre']['ci_lo']:.4f}, {b['acf1_pre']['ci_hi']:.4f}]")
    a(f"- Post ACF(1): {b['acf1_post']['point']:.4f} [{b['acf1_post']['ci_lo']:.4f}, {b['acf1_post']['ci_hi']:.4f}]")
    a(f"- ACF(1) diff: {b['acf1_diff']['point']:.4f} [{b['acf1_diff']['ci_lo']:.4f}, {b['acf1_diff']['ci_hi']:.4f}]")
    a(f"- **95% CI excludes zero**: {b['acf1_diff']['ci_excludes_zero']}")
    a(f"- Mean inflation diff: {b['mean_diff_pct']['point']:.2f}% [{b['mean_diff_pct']['ci_lo']:.2f}%, {b['mean_diff_pct']['ci_hi']:.2f}%]")
    a(f"- Volatility diff: {b['vol_diff_pct']['point']:.2f}% [{b['vol_diff_pct']['ci_lo']:.2f}%, {b['vol_diff_pct']['ci_hi']:.2f}%]")
    a("")

    a("## Sensitivity Analysis")
    if 'break_years' in r['sensitivity']:
        a("### Break Year Sensitivity")
        a("| Break Year | Pre N | Post N | ACF1 Pre | ACF1 Post | Diff | Fisher p |")
        a("|-----------|-------|--------|----------|-----------|------|----------|")
        for e in r['sensitivity']['break_years']:
            sig = "*" if e['fisher_p'] < 0.05 else ""
            a(f"| {e['break_year']} | {e['pre_n']} | {e['post_n']} | {e['acf1_pre']:.4f} | {e['acf1_post']:.4f} | {e['acf1_diff']:+.4f} | {e['fisher_p']:.4f}{sig} |")
        a("")

    if 'optimal_break' in r['sensitivity']:
        ob = r['sensitivity']['optimal_break']
        a(f"### Optimal Break Point")
        a(f"- Best year: {ob['best_year']} (F={ob['best_F']:.4f})")
        a(f"- Top 5: {', '.join(str(t['year']) for t in ob['top5'])}")
        a("")

    a("---")
    a(f"*Analysis completed with seed={SEED}, {N_PERMUTATIONS} permutations, {N_BOOTSTRAP} bootstrap resamples.*")

    with open('report.md', 'w') as f:
        f.write('
'.join(lines))

###############################################################################
# Verification Mode
###############################################################################
def verify():
    """Run 8+ machine-checkable assertions on the results."""
    print("=" * 72)
    print("VERIFICATION MODE")
    print("=" * 72)

    assert os.path.exists('results.json'), "results.json must exist"
    with open('results.json') as f:
        r = json.load(f)
    print("  [OK] results.json exists and is valid JSON")

    assert os.path.exists('report.md'), "report.md must exist"
    with open('report.md') as f:
        report_text = f.read()
    assert len(report_text) > 500, "report.md should be substantial"
    print("  [OK] report.md exists and has content")

    # Assertion 1: Sufficient data coverage
    n_years = r['series']['n_years']
    assert n_years >= 200, f"Should have 200+ years of data, got {n_years}"
    print(f"  [OK] Data coverage: {n_years} years (>= 200)")

    # Assertion 2: Year range spans centuries
    yr_min = r['series']['year_min']
    yr_max = r['series']['year_max']
    span = yr_max - yr_min
    assert span >= 300, f"Year span should be 300+, got {span}"
    print(f"  [OK] Year span: {span} years ({yr_min}-{yr_max})")

    # Assertion 3: Inflation series is reasonable
    n_inf = r['inflation_summary']['n_obs']
    assert n_inf >= 150, f"Should have 150+ inflation observations, got {n_inf}"
    print(f"  [OK] Inflation observations: {n_inf}")

    # Assertion 4: Mean inflation is in plausible range
    mean_inf = r['inflation_summary']['mean_pct']
    assert -5 < mean_inf < 20, f"Mean inflation {mean_inf}% seems implausible"
    print(f"  [OK] Mean inflation: {mean_inf:.2f}% (plausible range)")

    # Assertion 5: Chow test was computed with correct N_PERMUTATIONS
    assert r['chow_test']['n_permutations'] == N_PERMUTATIONS
    print(f"  [OK] Chow test used {N_PERMUTATIONS} permutations")

    # Assertion 6: p-value is valid
    p_val = r['chow_test']['p_value']
    assert 0 < p_val <= 1, f"p-value {p_val} out of range"
    print(f"  [OK] Chow test p-value: {p_val:.6f}")

    # Assertion 7: ACF values are in [-1, 1]
    for period in ['acf1_pre', 'acf1_post']:
        v = r['bootstrap'][period]['point']
        assert -1 <= v <= 1, f"ACF(1) {period} = {v} out of [-1,1]"
    print(f"  [OK] ACF values in valid range")

    # Assertion 8: Bootstrap CIs have correct structure
    diff_ci = r['bootstrap']['acf1_diff']
    assert diff_ci['ci_lo'] < diff_ci['ci_hi'], "CI lower must be < upper"
    print(f"  [OK] Bootstrap CI: [{diff_ci['ci_lo']:.4f}, {diff_ci['ci_hi']:.4f}]")

    # Assertion 9: Sensitivity analysis covers multiple break years
    n_breaks = len(r['sensitivity']['break_years'])
    assert n_breaks >= 3, f"Sensitivity should test 3+ break years, got {n_breaks}"
    print(f"  [OK] Sensitivity: {n_breaks} break years tested")

    # Assertion 10: SHA256 hash recorded
    assert len(r['data']['sha256']) == 64, "SHA256 should be 64 hex chars"
    print(f"  [OK] SHA256 hash recorded: {r['data']['sha256'][:16]}...")

    # Assertion 11: Results are deterministic (check key values haven't changed)
    # This checks that the analysis produces the same structure on re-run
    assert 'rolling' in r, "Rolling stats should be present"
    assert 'sensitivity' in r, "Sensitivity analysis should be present"
    print(f"  [OK] All expected result sections present")

    print()
    print(f"  ALL 11 ASSERTIONS PASSED")
    print("  VERIFICATION COMPLETE")
    return True

###############################################################################
# Entry Point
###############################################################################
if __name__ == '__main__':
    if '--verify' in sys.argv:
        success = verify()
        sys.exit(0 if success else 1)
    else:
        run_analysis()

SCRIPT_EOF
```

**Expected:** Exit code 0. File `boe_inflation/analyze.py` exists and is non-empty.

## Step 3: Run analysis

```bash
cd boe_inflation && python3 analyze.py
```

**Expected output format:**
```
======================================================================
[1/8] Downloading and caching BoE Millennium Dataset
======================================================================
  ...
[2/8] Parsing XLSX and extracting price-level series
  ...
[3/8] Computing annual inflation rates
  ...
[4/8] Computing rolling inflation statistics
  ...
[5/8] Permutation-based Chow test for structural breaks (2000 shuffles)
  ...
[6/8] Autocorrelation analysis: pre-1900 vs post-1900
  ...
[7/8] Bootstrap confidence intervals (2000 resamples)
  ...
[8/8] Sensitivity analysis across break points and parameters
  ...
ANALYSIS COMPLETE
```

**Expected files produced:**
- `results.json` — structured results with all statistics
- `report.md` — human-readable Markdown report

**Expected runtime:** Under 60 seconds on cached data. First run includes ~30MB download.

## Step 4: Verify results

```bash
cd boe_inflation && python3 analyze.py --verify
```

**Expected output:**
```
VERIFICATION MODE
  [OK] results.json exists and is valid JSON
  [OK] report.md exists and has content
  [OK] Data coverage: ... years (>= 200)
  [OK] Year span: ... years
  [OK] Inflation observations: ...
  [OK] Mean inflation: ...% (plausible range)
  [OK] Chow test used 2000 permutations
  [OK] Chow test p-value: ...
  [OK] ACF values in valid range
  [OK] Bootstrap CI: [...]
  [OK] Sensitivity: ... break years tested
  [OK] SHA256 hash recorded: ...
  [OK] All expected result sections present

  ALL 11 ASSERTIONS PASSED
  VERIFICATION COMPLETE
```

**Exit code:** 0 on success, non-zero on failure.

## Success Criteria

1. Script downloads and caches the BoE Millennium Dataset with SHA256 verification
2. Extracts a price-level series spanning 200+ years
3. Computes inflation rates and rolling statistics
4. Runs permutation-based Chow test with 2000 shuffles
5. Compares autocorrelation structure pre/post 1900 with Fisher z-test
6. Provides bootstrap confidence intervals for all key statistics
7. Performs sensitivity analysis across 7 candidate break years
8. All 11 verification assertions pass
9. `results.json` and `report.md` are produced

## Failure Conditions

- Network failure: Download fails after 3 retries → script exits with error message
- Data parsing failure: No suitable price series found → RuntimeError with list of sheets searched
- Insufficient data: Fewer than 200 years or 150 inflation observations → assertion failure
- Statistical anomaly: ACF values outside [-1, 1] → assertion failure