{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-0",
   "metadata": {},
   "source": [
    "# Pre-Trends Power Analysis (Roth 2022)\n",
    "\n",
    "A passing pre-trends test doesn't mean parallel trends holds—it may just mean the test has **low power** to detect violations. **Pre-Trends Power Analysis** (Roth 2022) answers a critical question:\n",
    "\n",
    "> \"What violations could my pre-trends test have detected?\"\n",
    "\n",
    "This notebook covers:\n",
    "\n",
    "1. Motivation: Why pre-trends tests can be misleading\n",
    "2. Basic usage with `PreTrendsPower`\n",
    "3. Computing the Minimum Detectable Violation (MDV)\n",
    "4. Power curves across violation magnitudes\n",
    "5. Different violation types (linear, constant, last period, custom)\n",
    "6. Integration with Honest DiD\n",
    "7. Visualization and reporting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from diff_diff import (\n",
    "    MultiPeriodDiD,\n",
    "    PreTrendsPower,\n",
    "    compute_pretrends_power,\n",
    "    compute_mdv,\n",
    "    plot_pretrends_power,\n",
    ")\n",
    "\n",
    "# For plots\n",
    "try:\n",
    "    import matplotlib.pyplot as plt\n",
    "    plt.style.use('seaborn-v0_8-whitegrid')\n",
    "    HAS_MATPLOTLIB = True\n",
    "except ImportError:\n",
    "    HAS_MATPLOTLIB = False\n",
    "    print(\"matplotlib not installed - visualization examples will be skipped\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-2",
   "metadata": {},
   "source": [
    "## 1. Motivation: The Problem with Pre-Trends Tests\n",
    "\n",
    "Standard practice in DiD analysis is to test for parallel trends by checking if pre-treatment coefficients are jointly zero. However, this approach has a fundamental problem:\n",
    "\n",
    "**A non-significant pre-trends test could mean:**\n",
    "1. Parallel trends actually holds ✓\n",
    "2. There's a violation, but the test lacks power to detect it ✗\n",
    "\n",
    "**Why does this matter?**\n",
    "\n",
    "If your pre-trends test has low power, it provides little reassurance about the validity of your DiD. Even a \"passing\" test (p > 0.05) might miss economically meaningful violations.\n",
    "\n",
    "**Roth (2022) introduces two key concepts:**\n",
    "\n",
    "1. **Power of the pre-trends test**: The probability of rejecting the null (parallel trends) when there *is* a violation of a given magnitude\n",
    "\n",
    "2. **Minimum Detectable Violation (MDV)**: The smallest violation your pre-trends test can detect with a target level of power (e.g., 80%)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-3",
   "metadata": {},
   "source": [
    "## 2. Generate Example Data\n",
    "\n",
    "We'll create panel data suitable for event study analysis with multiple pre-treatment periods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate event study data using the library function\n",
    "from diff_diff import generate_event_study_data\n",
    "\n",
    "# Generate panel data for event study analysis:\n",
    "# - 5 pre-treatment periods (0-4)\n",
    "# - 5 post-treatment periods (5-9)\n",
    "# - Half of units are treated starting at period 5\n",
    "df = generate_event_study_data(\n",
    "    n_units=300,\n",
    "    n_pre=5,\n",
    "    n_post=5,\n",
    "    treatment_fraction=0.5,\n",
    "    treatment_effect=5.0,\n",
    "    unit_fe_sd=2.0,\n",
    "    noise_sd=2.0,\n",
    "    seed=42\n",
    ")\n",
    "\n",
    "print(f\"Generated {len(df)} observations\")\n",
    "print(f\"Units: {df['unit'].nunique()} ({df[df['treated']==1]['unit'].nunique()} treated)\")\n",
    "print(f\"Periods: {df['period'].nunique()} (5 pre, 5 post)\")\n",
    "print(f\"True ATT: 5.0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-5",
   "metadata": {},
   "source": [
    "## 3. Fit Event Study\n",
    "\n",
    "First, we estimate a standard event study to get the pre-period coefficients and their variance-covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-6",
   "metadata": {},
   "outputs": [],
   "source": "# Fit event study with ALL periods (pre and post) relative to reference period\n# For pre-trends power analysis, we need coefficients for pre-periods too\nmp_did = MultiPeriodDiD()\n\n# Use period 4 as the reference period (last pre-period, excluded from estimation)\n# Specify post_periods as the actual post-treatment periods; MultiPeriodDiD\n# automatically estimates pre-period coefficients for the event study.\nevent_results = mp_did.fit(\n    df,\n    outcome='outcome',\n    treatment='treated',\n    time='period',\n    post_periods=[5, 6, 7, 8, 9]\n)\n\nprint(event_results.summary())"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the event study\n",
    "if HAS_MATPLOTLIB:\n",
    "    from diff_diff import plot_event_study\n",
    "    fig, ax = plt.subplots(figsize=(10, 6))\n",
    "    plot_event_study(\n",
    "        event_results,\n",
    "        ax=ax,\n",
    "        title='Event Study: Pre-Trends Look Good',\n",
    "        show=False\n",
    "    )\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-8",
   "metadata": {},
   "source": [
    "The pre-period coefficients (periods 0-3, with period 4 as reference) appear close to zero. But how confident should we be that parallel trends holds? Let's assess the **power** of this pre-trends test."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-9",
   "metadata": {},
   "source": [
    "## 4. Basic Pre-Trends Power Analysis\n",
    "\n",
    "The `PreTrendsPower` class computes the power of the pre-trends test to detect violations of different magnitudes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-10",
   "metadata": {},
   "outputs": [],
   "source": "# Create a PreTrendsPower object\npt = PreTrendsPower(\n    alpha=0.05,      # Significance level for pre-trends test\n    power=0.80,      # Target power for MDV calculation\n    violation_type='linear'  # Type of violation to consider\n)\n\n# Define the actual pre-treatment periods (those before treatment starts at period 5)\n# These are the periods we want to analyze for pre-trends power\npre_treatment_periods = [0, 1, 2, 3]\n\n# Fit to the event study results, specifying which periods are pre-treatment\npt_results = pt.fit(event_results, pre_periods=pre_treatment_periods)\n\nprint(pt_results.summary())"
  },
  {
   "cell_type": "markdown",
   "id": "cell-11",
   "metadata": {},
   "source": [
    "### Interpreting the Results\n",
    "\n",
    "**Key metrics:**\n",
    "\n",
    "1. **MDV (Minimum Detectable Violation)**: The smallest violation magnitude your pre-trends test can detect with 80% power\n",
    "   - Smaller MDV = more informative test\n",
    "   - If MDV is large, even big violations could go undetected\n",
    "\n",
    "2. **Power at specific violations**: How likely is the test to reject when there's a violation?\n",
    "   - Low power = uninformative \"passing\" test\n",
    "   - High power = reassuring \"passing\" test\n",
    "\n",
    "3. **Test informativeness**: Is the MDV small enough to be useful?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-12",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Access key results\n",
    "print(f\"Minimum Detectable Violation (MDV): {pt_results.mdv:.4f}\")\n",
    "print(f\"Target power: {pt_results.target_power:.0%}\")\n",
    "print(f\"Test informativeness: {'Informative' if pt_results.is_informative else 'Uninformative'}\")\n",
    "print(\"\")\n",
    "print(\"Interpretation:\")\n",
    "print(f\"  With 80% power, your pre-trends test can detect violations\")\n",
    "print(f\"  of magnitude {pt_results.mdv:.3f} or larger.\")\n",
    "print(f\"\")\n",
    "print(f\"  Violations smaller than {pt_results.mdv:.3f} would likely go undetected.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-13",
   "metadata": {},
   "source": [
    "## 5. Power at Specific Violation Magnitudes\n",
    "\n",
    "You can compute the power to detect a specific violation magnitude:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute power for specific violation magnitudes\n",
    "violations_to_check = [0.5, 1.0, 2.0, 3.0, 5.0]\n",
    "\n",
    "print(f\"{'Violation':>12} {'Power':>10} {'Detectable?':>15}\")\n",
    "print(\"-\" * 40)\n",
    "\n",
    "for v in violations_to_check:\n",
    "    power = pt_results.power_at(v)\n",
    "    detectable = \"Yes\" if power >= 0.80 else \"No\"\n",
    "    print(f\"{v:>12.1f} {power:>10.1%} {detectable:>15}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-15",
   "metadata": {},
   "source": [
    "## 6. Power Curves\n",
    "\n",
    "A **power curve** shows how the power to detect violations changes with violation magnitude. This is the most useful visualization for understanding your test's informativeness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate power curve\n",
    "curve = pt.power_curve(\n",
    "    event_results,\n",
    "    n_points=50,\n",
    "    pre_periods=pre_treatment_periods\n",
    ")\n",
    "\n",
    "# Preview the data\n",
    "print(\"Power curve data (first 10 points):\")\n",
    "print(curve.to_dataframe().head(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-17",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the power curve\n",
    "if HAS_MATPLOTLIB:\n",
    "    fig, ax = plt.subplots(figsize=(10, 6))\n",
    "    plot_pretrends_power(\n",
    "        curve,\n",
    "        ax=ax,\n",
    "        mdv=pt_results.mdv,  # Show MDV line on plot\n",
    "        target_power=0.80,\n",
    "        title='Pre-Trends Test Power Curve',\n",
    "        show=False\n",
    "    )\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-18",
   "metadata": {},
   "source": [
    "### Reading the Power Curve\n",
    "\n",
    "- **X-axis**: Violation magnitude (larger = worse violation of parallel trends)\n",
    "- **Y-axis**: Power (probability of rejecting when violation exists)\n",
    "- **Horizontal line at 0.80**: Conventional target power\n",
    "- **Vertical line at MDV**: Minimum detectable violation\n",
    "\n",
    "**Key insight**: The curve shows the range of violations your test could miss. If your ATT estimate could be biased by a violation smaller than the MDV, your results may be unreliable."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-19",
   "metadata": {},
   "source": [
    "## 7. Different Violation Types\n",
    "\n",
    "Pre-trends violations can take different forms. The `violation_type` parameter specifies the pattern:\n",
    "\n",
    "1. **Linear** (default): Violation grows linearly over time\n",
    "   - E.g., treated group diverges steadily from control\n",
    "   \n",
    "2. **Constant**: Same violation in all pre-periods\n",
    "   - E.g., level shift between groups\n",
    "   \n",
    "3. **Last period**: Violation only in the period just before treatment\n",
    "   - E.g., anticipation effects\n",
    "   \n",
    "4. **Custom**: User-specified violation pattern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare violation types\n",
    "violation_types = ['linear', 'constant', 'last_period']\n",
    "\n",
    "print(f\"{'Violation Type':>15} {'MDV':>10} {'Power at 2.0':>15}\")\n",
    "print(\"-\" * 45)\n",
    "\n",
    "for vtype in violation_types:\n",
    "    pt_v = PreTrendsPower(violation_type=vtype)\n",
    "    results_v = pt_v.fit(event_results, pre_periods=pre_treatment_periods)\n",
    "    power_at_2 = results_v.power_at(2.0)\n",
    "    print(f\"{vtype:>15} {results_v.mdv:>10.3f} {power_at_2:>15.1%}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Custom violation weights\n",
    "# Example: Violation concentrated in periods 2 and 3 (approaching treatment)\n",
    "# We have pre-periods 0, 1, 2, 3 estimated (reference period 4 is excluded)\n",
    "n_pre = 4  # Periods 0, 1, 2, 3\n",
    "custom_weights = np.zeros(n_pre)\n",
    "custom_weights[-2:] = 1.0  # Weight on last two pre-periods (periods 2 and 3)\n",
    "\n",
    "pt_custom = PreTrendsPower(\n",
    "    violation_type='custom',\n",
    "    violation_weights=custom_weights\n",
    ")\n",
    "results_custom = pt_custom.fit(event_results, pre_periods=pre_treatment_periods)\n",
    "\n",
    "print(f\"Custom violation (last 2 pre-periods): MDV = {results_custom.mdv:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-22",
   "metadata": {},
   "source": [
    "### Visualizing Different Violation Types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-23",
   "metadata": {},
   "outputs": [],
   "source": [
    "if HAS_MATPLOTLIB:\n",
    "    fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
    "    \n",
    "    for ax, vtype in zip(axes, ['linear', 'constant', 'last_period']):\n",
    "        pt_v = PreTrendsPower(violation_type=vtype)\n",
    "        results_v = pt_v.fit(event_results, pre_periods=pre_treatment_periods)\n",
    "        curve_v = pt_v.power_curve(event_results, n_points=50, pre_periods=pre_treatment_periods)\n",
    "        \n",
    "        plot_pretrends_power(\n",
    "            curve_v,\n",
    "            ax=ax,\n",
    "            mdv=results_v.mdv,  # Show MDV line on plot\n",
    "            target_power=0.80,\n",
    "            title=f'Violation Type: {vtype.replace(\"_\", \" \").title()}',\n",
    "            show=False\n",
    "        )\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-24",
   "metadata": {},
   "source": [
    "## 8. Integration with Honest DiD\n",
    "\n",
    "Pre-trends power analysis connects naturally with **Honest DiD** (Rambachan & Roth 2023). The workflow:\n",
    "\n",
    "1. Compute MDV from pre-trends power analysis\n",
    "2. Use MDV to calibrate the violation bound (M) in Honest DiD\n",
    "3. Compute robust confidence intervals under this calibrated bound\n",
    "\n",
    "This answers: \"If violations could be as large as what my pre-trends test could have missed, would my conclusions still hold?\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-25",
   "metadata": {},
   "outputs": [],
   "source": [
    "from diff_diff import HonestDiD\n",
    "\n",
    "# First, compute MDV\n",
    "pt = PreTrendsPower(violation_type='linear')\n",
    "pt_results = pt.fit(event_results, pre_periods=pre_treatment_periods)\n",
    "\n",
    "print(f\"MDV from pre-trends power analysis: {pt_results.mdv:.3f}\")\n",
    "print(\"\")\n",
    "\n",
    "# Use MDV to calibrate Honest DiD\n",
    "# The MDV tells us what violations we couldn't have detected\n",
    "# So we should check robustness to violations up to the MDV\n",
    "honest = HonestDiD(method='smoothness', M=pt_results.mdv)\n",
    "honest_results = honest.fit(event_results)\n",
    "\n",
    "print(\"Honest DiD results (M = MDV):\")\n",
    "print(honest_results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-26",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the built-in sensitivity integration\n",
    "sensitivity_results = pt.sensitivity_to_honest_did(\n",
    "    event_results,\n",
    "    pre_periods=pre_treatment_periods\n",
    ")\n",
    "\n",
    "print(\"Joint sensitivity analysis:\")\n",
    "print(f\"  MDV: {sensitivity_results['mdv']:.3f}\")\n",
    "print(f\"  Max pre-period SE: {sensitivity_results['max_pre_se']:.3f}\")\n",
    "print(f\"  MDV / max(SE): {sensitivity_results['mdv_in_ses']:.2f}\")\n",
    "print(\"\")\n",
    "print(\"Interpretation:\")\n",
    "print(sensitivity_results['interpretation'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-27",
   "metadata": {},
   "source": [
    "## 9. Convenience Functions\n",
    "\n",
    "For quick calculations, use the convenience functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-28",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quick MDV calculation\n",
    "mdv = compute_mdv(event_results, target_power=0.80, violation_type='linear', pre_periods=pre_treatment_periods)\n",
    "print(f\"MDV: {mdv:.3f}\")\n",
    "\n",
    "# Quick power calculation at a specific violation\n",
    "power_result = compute_pretrends_power(event_results, M=2.0, pre_periods=pre_treatment_periods)\n",
    "print(f\"Power at violation=2.0: {power_result.power:.1%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-29",
   "metadata": {},
   "source": [
    "## 10. Working with Real Event Studies\n",
    "\n",
    "In practice, you'll apply pre-trends power analysis to your actual event study estimates. Here's the typical workflow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-30",
   "metadata": {},
   "outputs": [],
   "source": "# Typical workflow for pre-trends power analysis\n\n# Step 1: Estimate event study with proper pre/post period classification\nmp_did = MultiPeriodDiD()\n\n# Specify actual post-treatment periods; pre-period coefficients are\n# estimated automatically by MultiPeriodDiD for the event study\npre_treatment_periods = [0, 1, 2, 3]  # Define which are pre-treatment\n\nresults = mp_did.fit(\n    df, \n    outcome='outcome',\n    treatment='treated', \n    time='period',\n    post_periods=[5, 6, 7, 8, 9]\n)\n\n# Step 2: Assess power of the pre-trends test  \nprint(\"Step 2: Pre-Trends Power Analysis\")\npt = PreTrendsPower(alpha=0.05, power=0.80, violation_type='linear')\npt_results = pt.fit(results, pre_periods=pre_treatment_periods)\nprint(f\"MDV (80% power): {pt_results.mdv:.3f}\")\nprint(\"\")\n\n# Step 3: Interpret\nprint(\"Step 3: Interpretation\")\nprint(f\"Your pre-trends test could only detect violations >= {pt_results.mdv:.3f}\")\nprint(f\"Violations smaller than this would likely go undetected.\")\nprint(\"\")\n\n# Step 4: Connect to Honest DiD for robust inference\nprint(\"Step 4: Robust Inference with Honest DiD\")\nhonest = HonestDiD(method='smoothness', M=pt_results.mdv)\nhonest_results = honest.fit(results)\nprint(f\"Robust 95% CI (M=MDV): [{honest_results.ci_lb:.3f}, {honest_results.ci_ub:.3f}]\")\nprint(f\"Conclusion: {'Effect is robust' if honest_results.is_significant else 'Effect may not be robust'}\")"
  },
  {
   "cell_type": "markdown",
   "id": "cell-31",
   "metadata": {},
   "source": [
    "## 11. Exporting Results\n",
    "\n",
    "Results can be exported to DataFrames for further analysis or reporting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-32",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export single result\n",
    "print(\"Single result as DataFrame:\")\n",
    "print(pt_results.to_dataframe())\n",
    "print(\"\")\n",
    "\n",
    "# Export power curve\n",
    "print(\"Power curve as DataFrame (first 10 rows):\")\n",
    "curve = pt.power_curve(event_results, pre_periods=pre_treatment_periods)\n",
    "print(curve.to_dataframe().head(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-33",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export to dict for JSON serialization\n",
    "result_dict = pt_results.to_dict()\n",
    "print(\"Result as dictionary:\")\n",
    "for key, value in result_dict.items():\n",
    "    print(f\"  {key}: {value}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-34",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "**Key Takeaways:**\n",
    "\n",
    "1. **Pre-trends tests can be misleading**: A \"passing\" test (p > 0.05) doesn't mean parallel trends holds—it may mean the test has low power.\n",
    "\n",
    "2. **MDV quantifies test informativeness**: The Minimum Detectable Violation tells you the smallest violation your test could detect with 80% power.\n",
    "\n",
    "3. **Power curves visualize sensitivity**: See how detection power changes with violation magnitude.\n",
    "\n",
    "4. **Different violation types matter**: Linear, constant, and last-period violations have different detectability.\n",
    "\n",
    "5. **Integration with Honest DiD**: Use MDV to calibrate sensitivity analysis bounds.\n",
    "\n",
    "**Best Practices:**\n",
    "\n",
    "- Always report pre-trends power analysis alongside standard pre-trends tests\n",
    "- Include power curves in supplementary materials\n",
    "- Use MDV to calibrate Honest DiD sensitivity analysis\n",
    "- Consider multiple violation types\n",
    "- Discuss what violation magnitudes would be economically meaningful in your setting\n",
    "\n",
    "**Reference:**\n",
    "\n",
    "Roth, J. (2022). Pretest with Caution: Event-Study Estimates after Testing for Parallel Trends. *American Economic Review: Insights*, 4(3), 305-322. https://doi.org/10.1257/aeri.20210236"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-35",
   "metadata": {},
   "source": [
    "## Related Tutorials\n",
    "\n",
    "- `04_parallel_trends.ipynb` - Testing and visualizing parallel trends\n",
    "- `05_honest_did.ipynb` - Sensitivity analysis for parallel trends violations\n",
    "- `06_power_analysis.ipynb` - Power analysis for study design (sample size, MDE)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}