Function bodies 13 total

def load_grid():
    """Load the specification grid results."""
    path = DATA / "v10_specification_grid.csv"
    if not path.exists():
        raise FileNotFoundError(f"Grid not found at {path}. Run estimate_v10_specification_grid.do first.")
    df = pd.read_csv(path)
    print(f"Loaded {len(df)} rows from specification grid")
    print(f"Outcomes: {df.outcome.nunique()}")
    print(f"Specs: {df.groupby(['shock_var','lag_spec','trend_type','sample_name']).ngroups}")
    return df

def score_significance(group, outcome):
    """Score: fraction of key horizons where coefficient is significant at 95%."""
    sig_count = 0
    total = 0
    for h in KEY_HORIZONS:
        row = group[group.h == h]
        if len(row) == 0:
            continue
        total += 1
        r = row.iloc[0]
        # Significant if CI doesn't cross zero
        if r.ci95_lo > 0 or r.ci95_hi < 0:
            sig_count += 1
    return sig_count / max(total, 1)

def score_sign_consistency(group, outcome):
    """Score: fraction of key horizons with the expected sign."""
    expected = EXPECTED_SIGNS.get(outcome)
    if expected is None:
        return 0.5  # no prior
    correct = 0
    total = 0
    for h in KEY_HORIZONS:
        row = group[group.h == h]
        if len(row) == 0:
            continue
        total += 1
        b = row.iloc[0].b
        if expected == "negative" and b < 0:
            correct += 1
        elif expected == "positive" and b > 0:
            correct += 1
    return correct / max(total, 1)

def score_magnitude(group, outcome, shock_var):
    """Score: 1 if TFP magnitudes are plausible (< 5% for 25bp equivalent)."""
    if "tfp" not in outcome:
        return 1.0  # only check TFP

    # For narrative shock: coefficients are per ±1 event
    # For continuous residual: need FFR first-stage normalization
    max_b = group.b.abs().max()

    if shock_var == "rr_narrative":
        # Per narrative event — no normalization needed
        # A contractionary event typically moves FFR by ~100-200bp
        # So per-25bp equivalent ≈ b / 4 to b / 8
        effective = max_b / 4  # conservative
    else:
        # For rr_resid: raw coefficient. Will be normalized by FFR first stage.
        # Typical b0_FFR ≈ 1.0, so scale_25bp ≈ 0.25
        effective = max_b * 0.25

    # Plausible: < 5% (0.05 in log levels)
    return 1.0 if effective < 0.05 else 0.0

def evaluate():
    """Main evaluation."""
    df = load_grid()

    # Define spec groups
    spec_cols = ["shock_var", "lag_spec", "trend_type", "sample_name"]

    # Score each spec × outcome combination
    results = []
    for (spec_key, spec_group) in df.groupby(spec_cols):
        shock_var, lag_spec, trend_type, sample_name = spec_key

        for outcome in df.outcome.unique():
            og = spec_group[spec_group.outcome == outcome]
            if len(og) == 0:
                continue

            sig = score_significance(og, outcome)
            sign = score_sign_consistency(og, outcome)
            mag = score_magnitude(og, outcome, shock_var)

            # Composite score (weighted)
            composite = 0.35 * sig + 0.35 * sign + 0.30 * mag

            results.append({
                "shock_var": shock_var,
                "lag_spec": lag_spec,
                "trend_type": trend_type,
                "sample_name": sample_name,
                "outcome": outcome,

def estimate_ffr_first_stage(panel, shock_var="rr_resid", p_s=4, p_y=2):
    """
    Estimate the impact of a unit RR shock on the federal funds rate.
    Uses the same LP specification as the paper.
    Returns b0_FFR (impact coefficient at h=0).
    """
    # FFR LP: FFR_{t+h} - FFR_{t-1} = α + γ*t + β*ε_t + Σ C_l*ε_{t-l} + Σ D_l*FFR_{t-l} + e
    dv = "fedfunds"
    dep = panel[dv] - panel[dv].shift(1)  # h=0: FFR_t - FFR_{t-1}

    # Build regressors
    X_cols = {}
    X_cols[shock_var] = panel[shock_var]
    for l in range(1, p_s + 1):
        X_cols[f"shock_lag{l}"] = panel[shock_var].shift(l)
    for l in range(1, p_y + 1):
        X_cols[f"dv_lag{l}"] = panel[dv].shift(l)
    X_cols["trend"] = np.arange(len(panel))
    X_cols["const"] = 1.0

    X = pd.DataFrame(X_cols)
    valid = dep.notna() & X.notna().all(axis=1)

    y = dep[valid].values
    Xm = X[valid].values

    beta = np.linalg.lstsq(Xm, y, rcond=None)[0]
    b0_ffr = beta[0]  # coefficient on shock_var (impact)

 

def normalize_irfs():
    """Main function: normalize empirical IRFs to 25bp FFR increase."""
    # Load panel
    panel = pd.read_stata(ROOT / "data" / "processed" / "v8_analysis_panel.dta")
    print(f"Panel: {len(panel)} rows")

    # Estimate FFR first stage
    b0_ffr, ffr_path = estimate_ffr_first_stage(panel)
    print(f"\nFFR First Stage (b0_FFR): {b0_ffr:.4f}")
    print(f"FFR path (h=0..4): {[f'{x:.4f}' for x in ffr_path]}")
    print(f"Units: percentage points of FFR per unit RR residual shock")

    # The RR residual and FFR are both in percentage points (quarterly level)
    # A 25bp annual increase = 0.25 percentage points annual ≈ 0.0625 pp quarterly
    # But actually, FFR is in annual percentage points in FRED
    # So b0_FFR is: pp annual FFR change per unit RR shock
    # We want: 25bp annual = 0.25 pp annual
    scale_25bp = 0.25 / b0_ffr
    print(f"\nScaling factor to 25bp annual: {scale_25bp:.4f}")
    print(f"(multiply all empirical IRFs by this factor)")

    #

def get_ffr_first_stage(grid_df, shock_var, lag_spec, trend_type, sample_name):
    """Extract FFR impact coefficient (h=0) for a given specification."""
    mask = (
        (grid_df.outcome == "fedfunds") &
        (grid_df.h == 0) &
        (grid_df.shock_var == shock_var) &
        (grid_df.lag_spec == lag_spec) &
        (grid_df.trend_type == trend_type) &
        (grid_df.sample_name == sample_name)
    )
    rows = grid_df[mask]
    if len(rows) == 0:
        raise ValueError(f"No FFR first-stage found for spec: {shock_var}/{lag_spec}/{trend_type}/{sample_name}")
    return rows.iloc[0].b

def normalize_spec(grid_df, shock_var, lag_spec, trend_type, sample_name):
    """Normalize all IRFs for a given specification to 25bp FFR increase."""
    # Get first-stage
    b0_ffr = get_ffr_first_stage(grid_df, shock_var, lag_spec, trend_type, sample_name)
    scale_25bp = 0.25 / b0_ffr

    print(f"\nSpec: {shock_var} / {lag_spec} / {trend_type} / {sample_name}")
    print(f"  FFR first-stage (h=0): {b0_ffr:.4f}")
    print(f"  Scale to 25bp: {scale_25bp:.4f}")

    # Filter to this spec (excluding fedfunds itself)
    mask = (
        (grid_df.shock_var == shock_var) &
        (grid_df.lag_spec == lag_spec) &
        (grid_df.trend_type == trend_type) &
        (grid_df.sample_name == sample_name) &
        (grid_df.outcome != "fedfunds")
    )
    spec_df = grid_df[mask].copy()

    # Scale
    spec_df["b"] = spec_df["b"] * scale_25bp
    spec_df["se"] = spec_df["se"] * abs(scale_25bp)
    spec_df["ci95_lo"] = spec_df["ci95_lo"] * scale_25bp
    spec_df["ci95_hi"] = spec_df["ci

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--grid", default=str(DATA / "v10_specification_grid.csv"))
    parser.add_argument("--shock", default="rr_resid")
    parser.add_argument("--lags", default="4_2")
    parser.add_argument("--trend", default="none")
    parser.add_argument("--sample", default="full")
    args = parser.parse_args()

    grid = pd.read_csv(args.grid)
    print(f"Loaded {len(grid)} rows from specification grid")

    norm_df, b0_ffr, scale = normalize_spec(
        grid, args.shock, args.lags, args.trend, args.sample
    )

    # Save
    OUT.mkdir(parents=True, exist_ok=True)
    out_path = OUT / "v10_paper_irfs_normalized_25bp.csv"
    norm_df.to_csv(out_path, index=False)
    print(f"\nSaved normalized IRFs to {out_path}")

    # Save metadata
    meta = pd.DataFrame([{
        "shock_var": args.shock,
        "lag_spec": args.lags,
        "trend_type": args.trend,
        "sample_name": args.sample,
        "b0_ffr": b0_ffr,
  

def load_data():
    df = pd.read_csv(DATA)
    print(f"Loaded {len(df)} rows")
    print(f"Outcomes: {df.outcome.unique()}")
    print(f"Sectors: {df.sector.unique()}")
    return df

def plot_irf_panel(df, outcome, ax, title, ylabel, show_legend=False):
    """Plot progressive vs stagnant IRF with confidence bands."""
    for sector, color, fill, label in [
        ("progressive", PROG_COLOR, PROG_FILL, "Progressive (durable mfg)"),
        ("stagnant", STAG_COLOR, STAG_FILL, "Stagnant (services + nondur mfg)"),
    ]:
        sub = df[(df.outcome == outcome) & (df.sector == sector)].sort_values("h")
        if len(sub) == 0:
            continue
        ax.plot(sub.h, sub.b, color=color, linewidth=1.8, label=label)
        ax.fill_between(sub.h, sub.ci95_lo, sub.ci95_hi, color=fill, alpha=0.3)

    ax.axhline(0, color="gray", linewidth=0.5)
    ax.set_xlabel("Years after shock")
    ax.set_ylabel(ylabel)
    ax.set_title(title)
    ax.set_xlim(0, 12)
    if show_legend:
        ax.legend(loc="best", frameon=True, framealpha=0.9)

def main():
    df = load_data()
    OUT.mkdir(parents=True, exist_ok=True)

    # Figure 4: 1x3 panel
    fig, axes = plt.subplots(1, 3, figsize=(14, 4.5))

    plot_irf_panel(df, "Domar", axes[0],
                   "(a) Domar Weights", "Change (×100)", show_legend=True)
    plot_irf_panel(df, "K_re", axes[1],
                   "(b) Capital Share", "Change")
    plot_irf_panel(df, "L_re", axes[2],
                   "(c) Labor Share", "Change")

    fig.suptitle("Long-Run Sectoral Reallocation: Annual KLEMS Panel Evidence",
                 fontsize=12, y=1.02)
    fig.tight_layout()

    outpath = OUT / "paper_v10_klems_reallocation.pdf"
    fig.savefig(outpath, bbox_inches="tight", dpi=300)
    print(f"Saved: {outpath}")
    plt.close(fig)

    # Also make a TFP panel (supplementary)
    fig2, axes2 = plt.subplots(1, 2, figsize=(10, 4.5))
    plot_irf_panel(df, "sec_TFP", axes2[0],
                   "(a) Industry TFP: Progressive", "Percent", show_legend=True)
    # For TFP, also