Initial extraction from monorepo

fire_planner/returns/__init__.py

@@ -0,0 +1 @@
+"""Historical-returns loaders and bootstrap samplers."""

fire_planner/returns/bootstrap.py

@@ -0,0 +1,60 @@
+"""Vectorised block bootstrap.
+
+Block bootstrap preserves serial correlation in returns and inflation —
+critical for a FIRE planner where sequence-of-returns risk is the
+dominant failure mode. Drawing IID-shuffled returns understates left
+tails because Depression-era runs of bad years would be impossible.
+
+Default block_size is 5 years per Politis & Romano (1994) guidance for
+asset-return series with multi-year mean-reversion. Block start indices
+are sampled uniformly with replacement from the legal range — i.e.,
+overlapping blocks are allowed, last block extends past series end is
+NOT allowed (we wrap with circular bootstrap).
+"""
+from __future__ import annotations
+
+import numpy as np
+import numpy.typing as npt
+
+from fire_planner.returns.shiller import ReturnsBundle
+
+
+def block_bootstrap(
+    bundle: ReturnsBundle,
+    n_paths: int,
+    n_years: int,
+    block_size: int = 5,
+    rng: np.random.Generator | None = None,
+) -> npt.NDArray[np.float64]:
+    """Return an ndarray of shape (n_paths, n_years, 3).
+
+    The third axis is (stock_nominal, bond_nominal, cpi).
+
+    Implementation: pick `ceil(n_years / block_size)` block starts per
+    path uniformly with replacement from the historical series; index
+    each block circularly; concatenate; truncate to n_years. The whole
+    op is vectorised — no Python-level loops over paths.
+    """
+    if rng is None:
+        rng = np.random.default_rng()
+    if block_size <= 0:
+        raise ValueError("block_size must be positive")
+
+    src = np.stack([bundle.stock_nominal, bundle.bond_nominal, bundle.cpi], axis=-1)
+    src_n = src.shape[0]
+    n_blocks = (n_years + block_size - 1) // block_size
+
+    # (n_paths, n_blocks) of block start indices in [0, src_n)
+    starts = rng.integers(0, src_n, size=(n_paths, n_blocks))
+    # Offsets within each block, broadcast: (1, 1, block_size)
+    offsets = np.arange(block_size).reshape(1, 1, block_size)
+    # (n_paths, n_blocks, block_size) of source indices, mod src_n
+    idx = (starts[:, :, None] + offsets) % src_n
+    # Flatten the inner two axes to (n_paths, n_blocks * block_size)
+    flat_idx = idx.reshape(n_paths, -1)
+    # Trim to exactly n_years
+    flat_idx = flat_idx[:, :n_years]
+    # Gather: result is (n_paths, n_years, 3). Explicit cast — numpy's
+    # advanced-indexing stubs return ndarray[Any], which trips strict mypy.
+    out: npt.NDArray[np.float64] = src[flat_idx]
+    return out

fire_planner/returns/shiller.py

@@ -0,0 +1,99 @@
+"""Shiller historical-returns loader.
+
+Robert Shiller's `ie_data.xls` (http://www.econ.yale.edu/~shiller/data.htm)
+provides monthly S&P 500 prices + dividends, 10-year Treasury rates, and
+CPI from 1871. We fold to annual: stock total return (price + reinvested
+dividends), bond total return (rate + price effect approximation), and
+CPI growth.
+
+The data file isn't shipped — call `load_from_csv` with a path the
+operator has fetched, or use `synthetic_returns()` for tests. The CLI's
+`ingest --source=shiller` command fetches the latest XLS, derives the
+CSV, and caches under `fire_planner/returns/_cache/shiller.csv`.
+
+Expected CSV columns:
+    year, stock_nominal_return, bond_nominal_return, cpi_inflation
+
+Each numeric column is a fraction (0.05 = 5% return), not basis points.
+"""
+from __future__ import annotations
+
+import csv
+from dataclasses import dataclass
+from pathlib import Path
+
+import numpy as np
+import numpy.typing as npt
+
+
+@dataclass(frozen=True)
+class ReturnsBundle:
+    """Aligned historical annual series.
+
+    All four arrays have identical length `n` and share the index axis —
+    `years[i]` is the calendar year of element `i` in the other arrays.
+    """
+    years: npt.NDArray[np.int32]
+    stock_nominal: npt.NDArray[np.float64]
+    bond_nominal: npt.NDArray[np.float64]
+    cpi: npt.NDArray[np.float64]
+
+    def __post_init__(self) -> None:
+        n = len(self.years)
+        if not (len(self.stock_nominal) == n == len(self.bond_nominal) == len(self.cpi)):
+            raise ValueError("ReturnsBundle arrays must share length")
+        if n == 0:
+            raise ValueError("ReturnsBundle cannot be empty")
+
+    @property
+    def n_years(self) -> int:
+        return len(self.years)
+
+    def stock_real(self) -> npt.NDArray[np.float64]:
+        """Real (CPI-adjusted) annual stock return."""
+        return (1 + self.stock_nominal) / (1 + self.cpi) - 1
+
+    def bond_real(self) -> npt.NDArray[np.float64]:
+        return (1 + self.bond_nominal) / (1 + self.cpi) - 1
+
+
+def load_from_csv(path: str | Path) -> ReturnsBundle:
+    p = Path(path)
+    years: list[int] = []
+    stocks: list[float] = []
+    bonds: list[float] = []
+    cpis: list[float] = []
+    with p.open() as f:
+        reader = csv.DictReader(f)
+        for row in reader:
+            years.append(int(row["year"]))
+            stocks.append(float(row["stock_nominal_return"]))
+            bonds.append(float(row["bond_nominal_return"]))
+            cpis.append(float(row["cpi_inflation"]))
+    return ReturnsBundle(
+        years=np.array(years, dtype=np.int32),
+        stock_nominal=np.array(stocks, dtype=np.float64),
+        bond_nominal=np.array(bonds, dtype=np.float64),
+        cpi=np.array(cpis, dtype=np.float64),
+    )
+
+
+def synthetic_returns(seed: int = 42, n_years: int = 150) -> ReturnsBundle:
+    """Deterministic synthetic returns for tests and bootstrap-conver-
+    gence experiments. Calibrated to roughly match Shiller's long-run
+    moments: stocks ~9.5% nominal / 17% vol, bonds ~5% / 8% vol, CPI
+    ~3% / 4% vol.
+
+    NOT for production use — `load_from_csv` is.
+    """
+    rng = np.random.default_rng(seed)
+    stock = rng.normal(0.095, 0.17, n_years)
+    bond = rng.normal(0.05, 0.08, n_years)
+    cpi = rng.normal(0.03, 0.04, n_years)
+    years = np.arange(1871, 1871 + n_years, dtype=np.int32)
+    return ReturnsBundle(
+        years=years,
+        stock_nominal=stock.astype(np.float64),
+        bond_nominal=bond.astype(np.float64),
+        cpi=cpi.astype(np.float64),
+    )