Applying the System Asset Pricing Model to US Civilian Firearms: Measuring the System Welfare Cost of Lethality-by-Design

{
  "paper_id": "firearms",
  "seed": 42,
  "draws": 100000,
  "correlation": 0.3,
  "annual_revenue": 23.2,
  "channels": {
    "C1_vsl_mortality": {
      "dist": "lognormal",
      "low": 329.97,
      "mid": 433.96,
      "high": 559.95,
      "weight": 0.35
    },
    "C2_healthcare_morbidity": {
      "dist": "normal",
      "low": 7,
      "mid": 8.7,
      "high": 11,
      "weight": 0.15
    },
    "C3_criminal_justice": {
      "dist": "normal",
      "low": 7,
      "mid": 8.8,
      "high": 11,
      "weight": 0.15
    },
    "C4_lost_productivity": {
      "dist": "lognormal",
      "low": 35,
      "mid": 45.7,
      "high": 59.99,
      "weight": 0.2
    },
    "C5_community_damage": {
      "dist": "lognormal",
      "low": 8,
      "mid": 11.3,
      "high": 16,
      "weight": 0.1
    },
    "C6_governance": {
      "dist": "lognormal",
      "low": 0.5,
      "mid": 0.75,
      "high": 1.2,
      "weight": 0.05
    }
  }
}

#!/usr/bin/env python3
"""
Monte Carlo simulation engine for SAPM papers.

Takes a per-domain JSON config, runs N draws with correlated channel
uncertainties, and produces reproducible results for paper injection.

Usage:
    python -m pipeline.monte_carlo tobacco
    python -m pipeline.monte_carlo --config pipeline/mc_configs/tobacco.json

Outputs:
    output/<paper_id>/mc_results.json   — all computed statistics
    output/<paper_id>/mc_histogram.json — binned histogram data for figure generation
"""
import json
import os
import sys
import numpy as np
from pathlib import Path
from scipy import stats

PROJECT_ROOT = Path(__file__).parent.parent
CONFIG_DIR = PROJECT_ROOT / "pipeline" / "mc_configs"
OUTPUT_DIR = PROJECT_ROOT / "output"
LEGACY_DIAGNOSTIC_ENV = "SAPM_ALLOW_LEGACY_DIAGNOSTIC"


def load_config(paper_id: str) -> dict:
    """Load the MC config for a domain."""
    config_path = CONFIG_DIR / f"{paper_id}.json"
    if not config_path.exists():
        raise FileNotFoundError(f"No MC config found: {config_path}")
    return json.loads(config_path.read_text(encoding="utf-8"))


def generate_correlated_draws(config: dict, rng: np.random.Generator) -> np.ndarray:
    """Generate correlated channel draws using Cholesky decomposition.

    Returns array of shape (n_draws, n_channels) with channel welfare costs.
    """
    channels = config["channels"]
    n_draws = config.get("draws", 100_000)
    rho = config.get("correlation", 0.3)
    n_ch = len(channels)

    # Build correlation matrix
    corr = np.full((n_ch, n_ch), rho)
    np.fill_diagonal(corr, 1.0)

    # Cholesky decomposition for correlation
    L = np.linalg.cholesky(corr)

    # Generate independent standard normal draws
    Z = rng.standard_normal((n_draws, n_ch))

    # Apply correlation structure
    Z_corr = Z @ L.T

    # Transform to channel-specific distributions
    results = np.zeros((n_draws, n_ch))
    channel_names = list(channels.keys())

    for i, (name, params) in enumerate(channels.items()):
        dist = params.get("dist", "lognormal")
        low = params["low"]
        mid = params["mid"]
        high = params["high"]

        # Convert uniform quantiles from correlated normals
        U = stats.norm.cdf(Z_corr[:, i])

        if dist == "lognormal":
            # Fit log-normal: median = mid, and 5th/95th pctiles ≈ low/high
            mu = np.log(mid)
            # Use (high/mid) ratio to estimate sigma
            # P95 = exp(mu + 1.645*sigma) → sigma = ln(high/mid) / 1.645
            sigma = np.log(high / mid) / 1.645
            results[:, i] = stats.lognorm.ppf(U, s=sigma, scale=mid)

        elif dist == "normal":
            # mid = mean, (high - low) ≈ 3.29 * sigma (5th to 95th)
            sigma = (high - low) / (2 * 1.645)
            results[:, i] = stats.norm.ppf(U, loc=mid, scale=sigma)

        elif dist == "triangular":
            # Triangular with mode at mid
            c = (mid - low) / (high - low)
            results[:, i] = stats.triang.ppf(U, c, loc=low, scale=high - low)

        elif dist == "uniform":
            results[:, i] = stats.uniform.ppf(U, loc=low, scale=high - low)

        else:
            raise ValueError(f"Unknown distribution '{dist}' for channel {name}")

    return results, channel_names


def compute_beta_w(channel_draws: np.ndarray, annual_revenue: float,
                   channel_weights: list = None) -> np.ndarray:
    """Compute β_W for each draw.

    β_W = total_welfare_cost / annual industry revenue

    If channel_weights are provided, apply them. Otherwise, sum all channels.
    """
    if channel_weights:
        weights = np.array(channel_weights)
        total_cost = np.sum(channel_draws * weights, axis=1)
    else:
        total_cost = np.sum(channel_draws, axis=1)

    return total_cost / annual_revenue


def build_histogram(beta_draws: np.ndarray, n_bins: int = 80) -> list:
    """Build histogram data for figure generation."""
    counts, bin_edges = np.histogram(beta_draws, bins=n_bins)
    bins = []
    for i in range(len(counts)):
        bins.append({
            "bin_low": round(float(bin_edges[i]), 3),
            "bin_high": round(float(bin_edges[i + 1]), 3),
            "bin_mid": round(float((bin_edges[i] + bin_edges[i + 1]) / 2), 3),
            "count": int(counts[i]),
            "density": round(float(counts[i] / len(beta_draws)), 6),
        })
    return bins


def build_sensitivity_matrix(config: dict, rng: np.random.Generator) -> list:
    """Build VSL × double-counting sensitivity matrix.

    Varies the largest channel (typically mortality/VSL) and applies
    double-counting adjustments to produce a grid of β_W values.
    """
    # Find the largest channel (assume it's the VSL/mortality channel)
    channels = config["channels"]
    ch_names = list(channels.keys())
    ch_mids = [channels[n]["mid"] for n in ch_names]
    largest_idx = ch_mids.index(max(ch_mids))
    largest_name = ch_names[largest_idx]
    largest_mid = ch_mids[largest_idx]

    annual_revenue = config.get("annual_revenue", config.get("revenue", config["private_payoff"]))

    # VSL multipliers (fraction of central estimate)
    vsl_multipliers = [0.5, 0.75, 1.0, 1.25, 1.5]
    # Double-counting adjustments
    dc_adjustments = [0.0, 0.10, 0.20, 0.30, 0.40]

    matrix = []
    for vsl_mult in vsl_multipliers:
        row = {"vsl_multiplier": vsl_mult, "values": {}}
        for dc_adj in dc_adjustments:
            # Compute deterministic β_W under these parameters
            total = 0
            for i, name in enumerate(ch_names):
                val = ch_mids[i]
                if i == largest_idx:
                    val *= vsl_mult
                total += val
            # Apply double-counting reduction
            total *= (1 - dc_adj)
            beta = round(total / annual_revenue, 2)
            row["values"][f"dc_{int(dc_adj * 100)}pct"] = beta
        matrix.append(row)

    return matrix


def run_monte_carlo(paper_id: str):
    """Run the full Monte Carlo simulation for a domain paper."""
    if os.environ.get(LEGACY_DIAGNOSTIC_ENV) != "1":
        raise RuntimeError(
            "Legacy synthetic MC runner is disabled. Use ppt-companion/scripts/run_admitted_mc_20260715.py "
            "with an admitted-mc-packet/1 packet. Set SAPM_ALLOW_LEGACY_DIAGNOSTIC=1 only for "
            "explicit historical diagnostics."
        )
    config = load_config(paper_id)
    rng = np.random.default_rng(seed=config.get("seed", 42))

    n_draws = config.get("draws", 100_000)
    annual_revenue = config.get("annual_revenue", config.get("revenue", config["private_payoff"]))
    channels = config["channels"]

    print(f"Monte Carlo simulation: {paper_id}")
    print(f"  Channels: {len(channels)}")
    print(f"  Draws: {n_draws:,}")
    print(f"  Correlation ρ: {config.get('correlation', 0.3)}")
    print(f"  Annual revenue denominator: ${annual_revenue:,}B")

    # Generate draws
    channel_draws, channel_names = generate_correlated_draws(config, rng)

    # Compute β_W
    weights = config.get("channel_weights")
    beta_draws = compute_beta_w(channel_draws, annual_revenue, weights)

    # Core statistics
    median_beta = float(np.median(beta_draws))
    mean_beta = float(np.mean(beta_draws))
    p5 = float(np.percentile(beta_draws, 5))
    p95 = float(np.percentile(beta_draws, 95))
    p1 = float(np.percentile(beta_draws, 1))
    p99 = float(np.percentile(beta_draws, 99))
    p_below_1 = float(np.mean(beta_draws < 1.0))
    std_beta = float(np.std(beta_draws))

    # Impossibility floor verification (if specified)
    floor = config.get("impossibility_floor")
    p_below_floor = float(np.mean(beta_draws < floor)) if floor else None

    # Total welfare cost statistics
    total_cost = np.sum(channel_draws, axis=1)
    pi_sa = total_cost - annual_revenue

    # Per-channel statistics
    channel_stats = {}
    for i, name in enumerate(channel_names):
        ch_draws = channel_draws[:, i]
        channel_stats[name] = {
            "median": round(float(np.median(ch_draws)), 1),
            "mean": round(float(np.mean(ch_draws)), 1),
            "p5": round(float(np.percentile(ch_draws, 5)), 1),
            "p95": round(float(np.percentile(ch_draws, 95)), 1),
            "std": round(float(np.std(ch_draws)), 1),
        }

    # Histogram data
    histogram = build_histogram(beta_draws)

    # Sensitivity matrix
    sensitivity = build_sensitivity_matrix(config, rng)

    # Assemble results
    results = {
        "paper_id": paper_id,
        "config_file": f"pipeline/mc_configs/{paper_id}.json",
        "seed": config.get("seed", 42),
        "n_draws": n_draws,
        "n_channels": len(channels),
        "correlation_rho": config.get("correlation", 0.3),
        "annual_revenue_B": annual_revenue,
        "private_payoff_B": annual_revenue,
        "beta_w": {
            "median": round(median_beta, 2),
            "mean": round(mean_beta, 2),
            "std": round(std_beta, 2),
            "p1": round(p1, 2),
            "p5": round(p5, 2),
            "p95": round(p95, 2),
            "p99": round(p99, 2),
            "ci_90": [round(p5, 1), round(p95, 1)],
            "ci_99": [round(p1, 1), round(p99, 1)],
            "p_below_1": round(p_below_1, 6),
            "p_below_1_pct": f"{p_below_1 * 100:.4f}%",
        },
        "impossibility_floor": {
            "floor_value": floor,
            "p_below_floor": round(p_below_floor, 6) if p_below_floor is not None else None,
            "p_below_floor_pct": f"{p_below_floor * 100:.4f}%" if p_below_floor is not None else None,
        } if floor else None,
        "welfare_cost": {
            "total_median_B": round(float(np.median(total_cost)), 1),
            "total_mean_B": round(float(np.mean(total_cost)), 1),
            "total_p5_B": round(float(np.percentile(total_cost, 5)), 1),
            "total_p95_B": round(float(np.percentile(total_cost, 95)), 1),
        },
        "pi_sa": {
            "median_B": round(float(np.median(pi_sa)), 1),
            "mean_B": round(float(np.mean(pi_sa)), 1),
            "p5_B": round(float(np.percentile(pi_sa, 5)), 1),
            "p95_B": round(float(np.percentile(pi_sa, 95)), 1),
        },
        "channel_statistics": channel_stats,
        "sensitivity_matrix": sensitivity,
    }

    # Save results
    out_dir = OUTPUT_DIR / paper_id
    out_dir.mkdir(parents=True, exist_ok=True)

    results_path = out_dir / "mc_results.json"
    results_path.write_text(json.dumps(results, indent=2), encoding="utf-8")
    print(f"\n  Results: {results_path}")

    hist_path = out_dir / "mc_histogram.json"
    hist_path.write_text(json.dumps(histogram, indent=2), encoding="utf-8")
    print(f"  Histogram: {hist_path}")

    # Print summary
    print(f"\n  === RESULTS ===")
    print(f"  β_W median:  {results['beta_w']['median']}")
    print(f"  β_W mean:    {results['beta_w']['mean']}")
    print(f"  90% CI:      [{results['beta_w']['ci_90'][0]}, {results['beta_w']['ci_90'][1]}]")
    print(f"  P(β_W < 1):  {results['beta_w']['p_below_1_pct']}")
    if floor:
        print(f"  P(β_W < {floor}): {results['impossibility_floor']['p_below_floor_pct']}")
    print(f"  Π^SA median: −${abs(results['pi_sa']['median_B']):,.0f}B/yr")

    return results


def main():
    if len(sys.argv) < 2:
        print("Usage: python -m pipeline.monte_carlo <paper_id>")
        sys.exit(1)

    paper_id = sys.argv[1]
    run_monte_carlo(paper_id)


if __name__ == "__main__":
    main()