{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a86464bc",
   "metadata": {},
   "source": [
    "# Texas Borderlands Enforcement and 2019 Disclosure Reform\n",
    "\n",
    "## Study Scope and Testable Design\n",
    "\n",
    "This notebook estimates **Texas-internal borderlands gaps** and **border-specific post-2019 disclosure shifts** at the district-year level.\n",
    "\n",
    "### In-Scope Research Questions\n",
    "\n",
    "**RQ1 (Border gaps):**\n",
    "Do border-exposed Texas districts differ from non-border districts in enforcement intensity and pipeline outcomes?\n",
    "\n",
    "**RQ2 (Disclosure heterogeneity):**\n",
    "Did the 2019 disclosure reform change enforcement outcomes differently in border districts versus non-border districts (level shift and post-policy trend differential)?\n",
    "\n",
    "### Outcomes (Operationalized)\n",
    "\n",
    "- **Inspection intensity:** inspections per well\n",
    "- **Detection intensity:** violations per inspection\n",
    "- **Enforcement timing:** average days from violation discovery to enforcement action date\n",
    "- **Follow-through proxy:** resolution/compliance-on-reinspection rate\n",
    "\n",
    "### Hypotheses Tested in This Notebook\n",
    "\n",
    "**H1 (Border inspection gap):**\n",
    "Border districts have lower inspection intensity than non-border districts.\n",
    "\n",
    "**H2 (Border pipeline disadvantage):**\n",
    "Border districts show weaker enforcement pipeline outcomes (higher violations per inspection and/or slower timing and/or lower resolution rates).\n",
    "\n",
    "**H3 (Disclosure heterogeneity in levels):**\n",
    "Post-2019 level shifts differ between border and non-border districts (`post_2019:border`).\n",
    "\n",
    "**H4 (Disclosure heterogeneity in trends):**\n",
    "Post-2019 trend shifts differ between border and non-border districts (`post_trend:border`).\n",
    "\n",
    "### Out of Scope for Current Data\n",
    "\n",
    "Direct Woods reaction-function tests are not estimated here because cross-jurisdiction competitor stringency data are unavailable.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4d51d985",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Setup complete\n"
     ]
    }
   ],
   "source": [
    "# Core setup\n",
    "from pathlib import Path\n",
    "import os\n",
    "import warnings\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.formula.api as smf\n",
    "from sqlalchemy import create_engine, text\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "pd.set_option('display.max_columns', 100)\n",
    "pd.set_option('display.float_format', '{:.4f}'.format)\n",
    "\n",
    "sns.set_theme(style='whitegrid', context='notebook')\n",
    "\n",
    "print('Setup complete')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e67af62c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Connected to PostgreSQL\n",
      "  Host: localhost\n",
      "  Database: texas_data\n",
      "  User: postgres\n",
      "\n",
      "Tables (first 15):\n",
      "['census_tract_demographics', 'census_tracts_2021', 'inspections', 'oil_gas_basins', 'rrc_well_api_raw', 'shale_plays', 'spatial_ref_sys', 'texmex_regions', 'violations', 'well_geo_features', 'well_shape_tract', 'well_shapes', 'well_shapes_backup_dedup', 'well_with_demographics_table']\n"
     ]
    }
   ],
   "source": [
    "# Explicit PostgreSQL connection (district notebook pattern)\n",
    "PGHOST = os.getenv('PGHOST', 'localhost')\n",
    "PGPORT = os.getenv('PGPORT', '5432')\n",
    "PGUSER = os.getenv('PGUSER', 'postgres')\n",
    "PGPASSWORD = os.getenv('PGPASSWORD', '')\n",
    "PGDATABASE = os.getenv('PGDATABASE', 'texas_data')\n",
    "\n",
    "pg_url = f'postgresql+psycopg2://{PGUSER}:{PGPASSWORD}@{PGHOST}:{PGPORT}/{PGDATABASE}'\n",
    "engine = create_engine(pg_url)\n",
    "\n",
    "try:\n",
    "    with engine.begin() as conn:\n",
    "        conn.execute(text('SELECT 1'))\n",
    "        tables = pd.read_sql(\n",
    "            text(\n",
    "                \"SELECT table_name \"\n",
    "                \"FROM information_schema.tables \"\n",
    "                \"WHERE table_schema = 'public' \"\n",
    "                \"AND table_type = 'BASE TABLE' \"\n",
    "                \"ORDER BY table_name\"\n",
    "            ),\n",
    "            conn,\n",
    "        )\n",
    "\n",
    "    print('Connected to PostgreSQL')\n",
    "    print(f'  Host: {PGHOST}')\n",
    "    print(f'  Database: {PGDATABASE}')\n",
    "    print(f'  User: {PGUSER}')\n",
    "    print('\\nTables (first 15):')\n",
    "    print(tables['table_name'].head(15).tolist())\n",
    "except Exception as exc:\n",
    "    print('Failed to connect to PostgreSQL with current PG* environment variables.')\n",
    "    raise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f76e1881",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2026-02-20 22:20:34,409 - INFO - Loaded 1010432 wells from public.well_shape_tract\n",
      "2026-02-20 22:22:30,525 - INFO - Recomputed border metrics for 1010432 wells (21112 within 25km, 40339 within 50km)\n",
      "2026-02-20 22:22:37,763 - INFO - Loaded 1878764 inspections from public.inspections\n",
      "2026-02-20 22:22:38,991 - INFO - Loaded 193338 violations from public.violations\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wells: 1,010,432\n",
      "Inspections: 1,878,764\n",
      "Violations: 193,338\n",
      "\n",
      "Well columns: ['api_norm', 'census_tract_geoid', 'distance_to_texmex_km', 'latitude', 'longitude', 'texmex_name', 'within_25km_texmex', 'within_50km_texmex']\n",
      "\n",
      "Inspection columns: ['api_norm', 'compliance', 'county', 'days_since_last_inspection', 'district', 'field_name', 'inspection_date', 'operator_name']\n",
      "\n",
      "Violation columns: ['api_norm', 'compliant_on_reinsp', 'district', 'drilling_permit_no', 'field_name', 'last_enf_action', 'last_enf_action_date', 'lease_fac_name', 'major_viol_ind', 'oil_lease_gas_well_id', 'operator_name', 'p5_operator_no', 'total_violations', 'violated_rule', 'violated_rule_desc', 'violation_disc_date', 'violation_number', 'violation_row_id', 'well_no']\n"
     ]
    }
   ],
   "source": [
    "# Load analyzer from local module and pull source data\n",
    "from well_analyzer import WellAnalyzer\n",
    "\n",
    "try:\n",
    "    analyzer = WellAnalyzer(engine=engine, chunk_size=50_000)\n",
    "except Exception as exc:\n",
    "    print('Failed to initialize WellAnalyzer after DB connection.')\n",
    "    print('If connection above succeeded, verify required source tables exist (well_shape_tract, inspections, violations).')\n",
    "    raise\n",
    "\n",
    "well_data = analyzer.data['well_data'].copy()\n",
    "inspections = analyzer.data['inspections'].copy()\n",
    "violations = analyzer.data['violations'].copy()\n",
    "\n",
    "print(f'Wells: {len(well_data):,}')\n",
    "print(f'Inspections: {len(inspections):,}')\n",
    "print(f'Violations: {len(violations):,}')\n",
    "print('\\nWell columns:', sorted(well_data.columns.tolist()))\n",
    "print('\\nInspection columns:', sorted(inspections.columns.tolist()))\n",
    "print('\\nViolation columns:', sorted(violations.columns.tolist()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a5599a25",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Config:\n",
      "  Years: 2015-2025\n",
      "  Policy year: 2019\n",
      "  Border districts: ['01', '02', '06', '08', '09', '10', '8A']\n",
      "  Output directory: /home/dadams/Repos/texas-borderlands/analysis/output_borderlands\n"
     ]
    }
   ],
   "source": [
    "# Analysis configuration\n",
    "START_YEAR = 2015\n",
    "END_YEAR = 2025\n",
    "POLICY_YEAR = 2019\n",
    "\n",
    "OUTPUT_DIR = Path('output_borderlands')\n",
    "OUTPUT_DIR.mkdir(exist_ok=True)\n",
    "\n",
    "# Border-adjacent district set for baseline binary treatment coding.\n",
    "# You can revise this list if your substantive definition changes.\n",
    "BORDER_DISTRICTS = {'01', '02', '06', '08', '8A', '09', '10'}\n",
    "\n",
    "print('Config:')\n",
    "print(f'  Years: {START_YEAR}-{END_YEAR}')\n",
    "print(f'  Policy year: {POLICY_YEAR}')\n",
    "print(f'  Border districts: {sorted(BORDER_DISTRICTS)}')\n",
    "print(f'  Output directory: {OUTPUT_DIR.resolve()}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d890fb98",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filtered sample:\n",
      "  Inspections: 1,867,859\n",
      "  Violations: 191,762\n"
     ]
    }
   ],
   "source": [
    "# Prepare core identifiers and time fields\n",
    "inspections = inspections.copy()\n",
    "violations = violations.copy()\n",
    "well_data = well_data.copy()\n",
    "\n",
    "inspections['api_norm'] = inspections['api_norm'].astype(str).str.strip()\n",
    "violations['api_norm'] = violations['api_norm'].astype(str).str.strip()\n",
    "well_data['api_norm'] = well_data['api_norm'].astype(str).str.strip()\n",
    "\n",
    "inspections['inspection_date'] = pd.to_datetime(inspections.get('inspection_date'), errors='coerce')\n",
    "violations['violation_disc_date'] = pd.to_datetime(violations.get('violation_disc_date'), errors='coerce')\n",
    "violations['last_enf_action_date'] = pd.to_datetime(violations.get('last_enf_action_date'), errors='coerce')\n",
    "\n",
    "inspections['year'] = inspections['inspection_date'].dt.year\n",
    "violations['year'] = violations['violation_disc_date'].dt.year\n",
    "\n",
    "inspections = inspections[inspections['year'].between(START_YEAR, END_YEAR, inclusive='both')].copy()\n",
    "violations = violations[violations['year'].between(START_YEAR, END_YEAR, inclusive='both')].copy()\n",
    "\n",
    "print('Filtered sample:')\n",
    "print(f'  Inspections: {len(inspections):,}')\n",
    "print(f'  Violations: {len(violations):,}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fa1727b1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Well coordinate ranges: lon [-106.537, -93.535], lat [25.857, 36.499]\n",
      "County CRS: EPSG:4326\n",
      "Well border map built:\n",
      "well_border_exposed\n",
      "non-border    840912\n",
      "border        169520\n",
      "Name: count, dtype: int64\n",
      "\n",
      "Border subtype counts (50km):\n",
      "  TX-MX: 40339\n",
      "  TX-NM: 81567\n",
      "  TX-OK: 19643\n",
      "  TX-LA: 29675\n"
     ]
    }
   ],
   "source": [
    "# Build well-level border exposure map + state-border flags\n",
    "from shapely.geometry import LineString\n",
    "\n",
    "well_border = well_data[['api_norm', 'latitude', 'longitude']].drop_duplicates('api_norm').copy()\n",
    "well_border['latitude'] = pd.to_numeric(well_border['latitude'], errors='coerce')\n",
    "well_border['longitude'] = pd.to_numeric(well_border['longitude'], errors='coerce')\n",
    "\n",
    "# Start with Mexico-border metrics already computed in WellAnalyzer\n",
    "for c in ['distance_to_texmex_km', 'within_25km_texmex', 'within_50km_texmex', 'texmex_name']:\n",
    "    if c in well_data.columns:\n",
    "        well_border = well_border.merge(\n",
    "            well_data[['api_norm', c]].drop_duplicates('api_norm'),\n",
    "            on='api_norm',\n",
    "            how='left'\n",
    "        )\n",
    "\n",
    "for c in ['within_25km_texmex', 'within_50km_texmex']:\n",
    "    if c not in well_border.columns:\n",
    "        well_border[c] = 0\n",
    "    well_border[c] = pd.to_numeric(well_border[c], errors='coerce').fillna(0).astype(int)\n",
    "\n",
    "# Initialize state-border fields\n",
    "for c in [\n",
    "    'distance_to_nm_border_km', 'distance_to_ok_border_km', 'distance_to_la_border_km',\n",
    "    'within_25km_nm_border', 'within_50km_nm_border',\n",
    "    'within_25km_ok_border', 'within_50km_ok_border',\n",
    "    'within_25km_la_border', 'within_50km_la_border',\n",
    "]:\n",
    "    well_border[c] = np.nan if c.startswith('distance_') else 0\n",
    "\n",
    "try:\n",
    "    import geopandas as gpd\n",
    "\n",
    "    def _resolve_path(candidates):\n",
    "        for c in candidates:\n",
    "            if c.exists():\n",
    "                return c\n",
    "        return None\n",
    "\n",
    "    county_shp = _resolve_path([\n",
    "        Path('../data/texas_county_shape/tl_2025_48_cousub/tl_2025_48_cousub.shp'),\n",
    "        Path('data/texas_county_shape/tl_2025_48_cousub/tl_2025_48_cousub.shp'),\n",
    "    ])\n",
    "\n",
    "    if county_shp is None:\n",
    "        print('State-border flags skipped: Texas county shapefile not found.')\n",
    "    else:\n",
    "        counties = gpd.read_file(county_shp)\n",
    "        if counties.crs is None:\n",
    "            counties = counties.set_crs('EPSG:4269', allow_override=True)\n",
    "        counties = counties.to_crs('EPSG:4326')\n",
    "\n",
    "        tx_union = counties.geometry.union_all()\n",
    "        tx_boundary = tx_union.boundary\n",
    "\n",
    "        # Fast coordinate sanity checks to catch ingest issues\n",
    "        lon_min, lon_max = well_border['longitude'].min(), well_border['longitude'].max()\n",
    "        lat_min, lat_max = well_border['latitude'].min(), well_border['latitude'].max()\n",
    "        print(f'Well coordinate ranges: lon [{lon_min:.3f}, {lon_max:.3f}], lat [{lat_min:.3f}, {lat_max:.3f}]')\n",
    "        print(f'County CRS: {counties.crs}')\n",
    "\n",
    "        # Build explicit border seed lines and clip to actual Texas boundary.\n",
    "        # This avoids accidental use of wrong boundary sectors.\n",
    "        nm_seed = LineString([(-103.02, 31.78), (-103.02, 36.50)])\n",
    "        ok_seed = LineString([(-103.02, 36.50), (-100.00, 36.50), (-94.04, 33.57)])\n",
    "        la_seed = LineString([(-94.04, 33.02), (-93.70, 31.00), (-93.75, 29.70)])\n",
    "\n",
    "        nm_geom = tx_boundary.intersection(nm_seed.buffer(0.35))\n",
    "        ok_geom = tx_boundary.intersection(ok_seed.buffer(0.35))\n",
    "        la_geom = tx_boundary.intersection(la_seed.buffer(0.35))\n",
    "\n",
    "        # Fallback to seed lines if an intersection unexpectedly collapses.\n",
    "        if nm_geom.is_empty:\n",
    "            nm_geom = nm_seed\n",
    "        if ok_geom.is_empty:\n",
    "            ok_geom = ok_seed\n",
    "        if la_geom.is_empty:\n",
    "            la_geom = la_seed\n",
    "\n",
    "        valid_pts = well_border.dropna(subset=['longitude', 'latitude']).copy()\n",
    "        if len(valid_pts) > 0:\n",
    "            pts = gpd.GeoDataFrame(\n",
    "                valid_pts[['api_norm', 'longitude', 'latitude']].copy(),\n",
    "                geometry=gpd.points_from_xy(valid_pts['longitude'], valid_pts['latitude']),\n",
    "                crs='EPSG:4326'\n",
    "            ).to_crs('EPSG:5070')\n",
    "\n",
    "            nm_m = gpd.GeoSeries([nm_geom], crs='EPSG:4326').to_crs('EPSG:5070').iloc[0]\n",
    "            ok_m = gpd.GeoSeries([ok_geom], crs='EPSG:4326').to_crs('EPSG:5070').iloc[0]\n",
    "            la_m = gpd.GeoSeries([la_geom], crs='EPSG:4326').to_crs('EPSG:5070').iloc[0]\n",
    "\n",
    "            pts['distance_to_nm_border_km'] = pts.geometry.distance(nm_m) / 1000.0\n",
    "            pts['distance_to_ok_border_km'] = pts.geometry.distance(ok_m) / 1000.0\n",
    "            pts['distance_to_la_border_km'] = pts.geometry.distance(la_m) / 1000.0\n",
    "\n",
    "            pts['within_25km_nm_border'] = (pts['distance_to_nm_border_km'] <= 25).astype(int)\n",
    "            pts['within_50km_nm_border'] = (pts['distance_to_nm_border_km'] <= 50).astype(int)\n",
    "            pts['within_25km_ok_border'] = (pts['distance_to_ok_border_km'] <= 25).astype(int)\n",
    "            pts['within_50km_ok_border'] = (pts['distance_to_ok_border_km'] <= 50).astype(int)\n",
    "            pts['within_25km_la_border'] = (pts['distance_to_la_border_km'] <= 25).astype(int)\n",
    "            pts['within_50km_la_border'] = (pts['distance_to_la_border_km'] <= 50).astype(int)\n",
    "\n",
    "            cols = [\n",
    "                'api_norm',\n",
    "                'distance_to_nm_border_km', 'distance_to_ok_border_km', 'distance_to_la_border_km',\n",
    "                'within_25km_nm_border', 'within_50km_nm_border',\n",
    "                'within_25km_ok_border', 'within_50km_ok_border',\n",
    "                'within_25km_la_border', 'within_50km_la_border',\n",
    "            ]\n",
    "            well_border = well_border.drop(\n",
    "                columns=[c for c in cols if c in well_border.columns and c != 'api_norm'],\n",
    "                errors='ignore'\n",
    "            )\n",
    "            well_border = well_border.merge(pts[cols], on='api_norm', how='left')\n",
    "\n",
    "except ImportError:\n",
    "    print('State-border flags skipped: geopandas not installed.')\n",
    "\n",
    "# Fill remaining flags/distances\n",
    "for c in [\n",
    "    'within_25km_nm_border', 'within_50km_nm_border',\n",
    "    'within_25km_ok_border', 'within_50km_ok_border',\n",
    "    'within_25km_la_border', 'within_50km_la_border',\n",
    "]:\n",
    "    well_border[c] = pd.to_numeric(well_border[c], errors='coerce').fillna(0).astype(int)\n",
    "\n",
    "for c in ['distance_to_nm_border_km', 'distance_to_ok_border_km', 'distance_to_la_border_km']:\n",
    "    well_border[c] = pd.to_numeric(well_border[c], errors='coerce')\n",
    "\n",
    "well_border['within_50km_state_border_any'] = (\n",
    "    (well_border['within_50km_nm_border'] > 0) |\n",
    "    (well_border['within_50km_ok_border'] > 0) |\n",
    "    (well_border['within_50km_la_border'] > 0)\n",
    ").astype(int)\n",
    "\n",
    "well_border['well_border_exposed'] = (\n",
    "    (well_border['within_50km_texmex'] > 0) |\n",
    "    (well_border['within_50km_state_border_any'] > 0)\n",
    ").astype(int)\n",
    "\n",
    "print('Well border map built:')\n",
    "print(well_border['well_border_exposed'].value_counts(dropna=False).rename(index={0:'non-border',1:'border'}))\n",
    "print('\\nBorder subtype counts (50km):')\n",
    "print('  TX-MX:', int((well_border['within_50km_texmex'] > 0).sum()))\n",
    "print('  TX-NM:', int((well_border['within_50km_nm_border'] > 0).sum()))\n",
    "print('  TX-OK:', int((well_border['within_50km_ok_border'] > 0).sum()))\n",
    "print('  TX-LA:', int((well_border['within_50km_la_border'] > 0).sum()))\n",
    "\n",
    "if int((well_border['within_50km_nm_border'] > 0).sum()) < 500:\n",
    "    print('WARNING: TX-NM 50km count is unexpectedly low; inspect border geometry and coordinate ranges above.')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f6c09f4f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 542472 characters]",
      "text/plain": [
       "<Figure size 1000x1000 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved map: /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/well_border_exposure_map.png\n",
      "Points plotted (trimmed to Texas): 50,000\n",
      "Border points: 8,556; Non-border points: 41,444\n"
     ]
    }
   ],
   "source": [
    "# Map: border vs non-border wells\n",
    "try:\n",
    "    import geopandas as gpd\n",
    "    from matplotlib.lines import Line2D\n",
    "\n",
    "    def _resolve_path(candidates):\n",
    "        for c in candidates:\n",
    "            if c.exists():\n",
    "                return c\n",
    "        return None\n",
    "\n",
    "    # Toggle border line display\n",
    "    SHOW_BORDER_LINE = False\n",
    "\n",
    "    # Support running notebook from either analysis/ or repo root\n",
    "    county_shp = _resolve_path([\n",
    "        Path('../data/texas_county_shape/tl_2025_48_cousub/tl_2025_48_cousub.shp'),\n",
    "        Path('data/texas_county_shape/tl_2025_48_cousub/tl_2025_48_cousub.shp'),\n",
    "    ])\n",
    "    texmex_shp = _resolve_path([\n",
    "        Path('../data/texmex_shape/tl_2023_us_internationalboundary.shp'),\n",
    "        Path('data/texmex_shape/tl_2023_us_internationalboundary.shp'),\n",
    "    ])\n",
    "\n",
    "    if county_shp is None or texmex_shp is None:\n",
    "        print('Required shapefiles not found from current working directory.')\n",
    "        print('Checked county paths under ../data and ./data')\n",
    "        print('Checked texmex paths under ../data and ./data')\n",
    "    else:\n",
    "        counties = gpd.read_file(county_shp)\n",
    "        texmex = gpd.read_file(texmex_shp)\n",
    "\n",
    "        # Build well points from available coordinates\n",
    "        coord_cols = {'longitude', 'latitude'}\n",
    "        if not coord_cols.issubset(set(well_data.columns)):\n",
    "            print('Cannot map wells: longitude/latitude columns are missing from well_data.')\n",
    "        else:\n",
    "            points_df = well_data[['api_norm', 'longitude', 'latitude']].dropna().copy()\n",
    "            points_df['longitude'] = pd.to_numeric(points_df['longitude'], errors='coerce')\n",
    "            points_df['latitude'] = pd.to_numeric(points_df['latitude'], errors='coerce')\n",
    "            points_df = points_df.dropna(subset=['longitude', 'latitude'])\n",
    "\n",
    "            points_df = points_df.merge(\n",
    "                well_border[['api_norm', 'well_border_exposed']],\n",
    "                on='api_norm',\n",
    "                how='left'\n",
    "            )\n",
    "            points_df['well_border_exposed'] = points_df['well_border_exposed'].fillna(0).astype(int)\n",
    "\n",
    "            # Fallback if no wells are marked border-exposed: use district border coding\n",
    "            if points_df['well_border_exposed'].sum() == 0:\n",
    "                well_district = inspections[['api_norm', 'district']].dropna().copy()\n",
    "                well_district['api_norm'] = well_district['api_norm'].astype(str).str.strip()\n",
    "                well_district['district'] = well_district['district'].astype(str).str.strip()\n",
    "                well_district = (\n",
    "                    well_district.groupby('api_norm')['district']\n",
    "                    .agg(lambda s: s.mode().iat[0] if not s.mode().empty else s.iloc[0])\n",
    "                    .reset_index()\n",
    "                )\n",
    "                points_df = points_df.merge(well_district, on='api_norm', how='left')\n",
    "                points_df['well_border_exposed'] = points_df['district'].isin(BORDER_DISTRICTS).astype(int)\n",
    "                print('No positive within_25km/50km flags found; using district-based fallback for map coloring.')\n",
    "\n",
    "            wells_gdf = gpd.GeoDataFrame(\n",
    "                points_df,\n",
    "                geometry=gpd.points_from_xy(points_df['longitude'], points_df['latitude']),\n",
    "                crs='EPSG:4326'\n",
    "            )\n",
    "\n",
    "            # Keep only Texas border segments if NAME exists\n",
    "            if 'NAME' in texmex.columns:\n",
    "                texmex_plot = texmex[texmex['NAME'].astype(str).str.contains('Texas', case=False, na=False)].copy()\n",
    "                if texmex_plot.empty:\n",
    "                    texmex_plot = texmex.copy()\n",
    "            else:\n",
    "                texmex_plot = texmex.copy()\n",
    "\n",
    "            # Harmonize CRS\n",
    "            if counties.crs is None:\n",
    "                counties = counties.set_crs('EPSG:4326', allow_override=True)\n",
    "            if texmex_plot.crs is None:\n",
    "                texmex_plot = texmex_plot.set_crs('EPSG:4326', allow_override=True)\n",
    "            if wells_gdf.crs is None:\n",
    "                wells_gdf = wells_gdf.set_crs('EPSG:4326', allow_override=True)\n",
    "\n",
    "            if counties.crs != wells_gdf.crs:\n",
    "                counties = counties.to_crs(wells_gdf.crs)\n",
    "            if texmex_plot.crs != wells_gdf.crs:\n",
    "                texmex_plot = texmex_plot.to_crs(wells_gdf.crs)\n",
    "\n",
    "            # Trim everything to Texas counties footprint\n",
    "            texas_union = counties.geometry.union_all()\n",
    "            texmex_plot = gpd.clip(texmex_plot, texas_union)\n",
    "            wells_gdf = gpd.clip(wells_gdf, texas_union)\n",
    "\n",
    "            # Sample for legibility and speed (after clip)\n",
    "            max_points = 50000\n",
    "            if len(wells_gdf) > max_points:\n",
    "                wells_gdf = wells_gdf.sample(max_points, random_state=42)\n",
    "\n",
    "            fig, ax = plt.subplots(figsize=(10, 10))\n",
    "            counties.plot(ax=ax, color='#f2efe8', edgecolor='#d7d2c8', linewidth=0.35)\n",
    "            if SHOW_BORDER_LINE and not texmex_plot.empty:\n",
    "                texmex_plot.plot(ax=ax, color='#9ca3af', linewidth=0.8, alpha=0.8)\n",
    "\n",
    "            non_border = wells_gdf[wells_gdf['well_border_exposed'] == 0]\n",
    "            border = wells_gdf[wells_gdf['well_border_exposed'] == 1]\n",
    "\n",
    "            if not non_border.empty:\n",
    "                non_border.plot(ax=ax, markersize=3, color='#6b7280', alpha=0.25, zorder=2)\n",
    "            if not border.empty:\n",
    "                border.plot(ax=ax, markersize=6, color='#ef4444', alpha=0.65, zorder=3)\n",
    "\n",
    "            # Force Texas extent so map does not zoom to continental US\n",
    "            minx, miny, maxx, maxy = counties.total_bounds\n",
    "            pad_x = (maxx - minx) * 0.02\n",
    "            pad_y = (maxy - miny) * 0.02\n",
    "            ax.set_xlim(minx - pad_x, maxx + pad_x)\n",
    "            ax.set_ylim(miny - pad_y, maxy + pad_y)\n",
    "\n",
    "            n_non = len(non_border)\n",
    "            n_border = len(border)\n",
    "\n",
    "            legend_items = [\n",
    "                Line2D([0], [0], marker='o', color='none', markerfacecolor='#6b7280', markersize=6,\n",
    "                       label=f'Non-border wells (n={n_non:,})'),\n",
    "                Line2D([0], [0], marker='o', color='none', markerfacecolor='#ef4444', markersize=7,\n",
    "                       label=f'Border wells (n={n_border:,})'),\n",
    "            ]\n",
    "\n",
    "            ax.legend(handles=legend_items, loc='lower left', frameon=True)\n",
    "            ax.set_title('Texas Wells: Border vs Non-Border Exposure', fontsize=13, weight='bold')\n",
    "            ax.set_axis_off()\n",
    "\n",
    "            map_path = OUTPUT_DIR / 'well_border_exposure_map.png'\n",
    "            fig.savefig(map_path, dpi=300, bbox_inches='tight')\n",
    "            plt.show()\n",
    "\n",
    "            print(f'Saved map: {map_path.resolve()}')\n",
    "            print(f'Points plotted (trimmed to Texas): {len(wells_gdf):,}')\n",
    "            print(f'Border points: {n_border:,}; Non-border points: {n_non:,}')\n",
    "\n",
    "except ImportError:\n",
    "    print('geopandas/shapely not installed; skipping map cell.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b33ae546",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Panel built:\n",
      "  Rows: 143\n",
      "  Districts: 13\n",
      "  Years: 2015-2025\n",
      "  district  year  total_inspections  unique_wells  compliant_inspections  \\\n",
      "0       01  2015               3498          2816                   3027   \n",
      "1       01  2016               6499          4055                   5028   \n",
      "2       01  2017               8649          6153                   7613   \n",
      "3       01  2018              10966          9109                   9668   \n",
      "4       01  2019               8097          6447                   6818   \n",
      "\n",
      "   border_exposed_insp  inspection_intensity  compliance_rate  \\\n",
      "0                  319                1.2422           0.8654   \n",
      "1                  625                1.6027           0.7737   \n",
      "2                 1042                1.4057           0.8802   \n",
      "3                 1914                1.2039           0.8816   \n",
      "4                 2096                1.2559           0.8420   \n",
      "\n",
      "   share_border_exposed_insp  total_violations  major_violations  \\\n",
      "0                     0.0912               592                 0   \n",
      "1                     0.0962              1902                 0   \n",
      "2                     0.1205              1439                 0   \n",
      "3                     0.1745              1771                 0   \n",
      "4                     0.2589              1506                 2   \n",
      "\n",
      "   compliant_on_reinsp  avg_days_to_enforcement  med_days_to_enforcement  \\\n",
      "0                  284                 124.9122                  19.0000   \n",
      "1                  472                 230.0910                  20.0000   \n",
      "2                  740                 326.0945                  33.0000   \n",
      "3                 1012                 290.2191                  82.0000   \n",
      "4                  771                 261.2776                  73.0000   \n",
      "\n",
      "   violations_per_inspection  major_share  resolution_rate  post_2019  \\\n",
      "0                     0.1692       0.0000           0.4797          0   \n",
      "1                     0.2927       0.0000           0.2482          0   \n",
      "2                     0.1664       0.0000           0.5142          0   \n",
      "3                     0.1615       0.0000           0.5714          0   \n",
      "4                     0.1860       0.0013           0.5120          1   \n",
      "\n",
      "   year_from_policy  post_trend  border_district  border_exposure_district  \n",
      "0                -4           0                1                         0  \n",
      "1                -3           0                1                         0  \n",
      "2                -2           0                1                         0  \n",
      "3                -1           0                1                         0  \n",
      "4                 0           0                1                         1  \n"
     ]
    }
   ],
   "source": [
    "# District-year outcome panel\n",
    "inspections_merged = inspections.merge(\n",
    "    well_border[['api_norm', 'well_border_exposed']],\n",
    "    on='api_norm',\n",
    "    how='left'\n",
    ")\n",
    "inspections_merged['well_border_exposed'] = inspections_merged['well_border_exposed'].fillna(0).astype(int)\n",
    "\n",
    "violations_merged = violations.merge(\n",
    "    well_border[['api_norm', 'well_border_exposed']],\n",
    "    on='api_norm',\n",
    "    how='left'\n",
    ")\n",
    "violations_merged['well_border_exposed'] = violations_merged['well_border_exposed'].fillna(0).astype(int)\n",
    "\n",
    "# Timing metric used in district analysis\n",
    "days_to_enf = (violations_merged['last_enf_action_date'] - violations_merged['violation_disc_date']).dt.days\n",
    "violations_merged['days_to_enforcement'] = days_to_enf.where((days_to_enf >= 0) & (days_to_enf <= 3650))\n",
    "\n",
    "insp_panel = inspections_merged.groupby(['district', 'year'], as_index=False).agg(\n",
    "    total_inspections=('api_norm', 'size'),\n",
    "    unique_wells=('api_norm', 'nunique'),\n",
    "    compliant_inspections=('compliance', lambda x: (x == 'Yes').sum()),\n",
    "    border_exposed_insp=('well_border_exposed', 'sum')\n",
    ")\n",
    "insp_panel['inspection_intensity'] = insp_panel['total_inspections'] / insp_panel['unique_wells'].replace(0, np.nan)\n",
    "insp_panel['compliance_rate'] = insp_panel['compliant_inspections'] / insp_panel['total_inspections'].replace(0, np.nan)\n",
    "insp_panel['share_border_exposed_insp'] = insp_panel['border_exposed_insp'] / insp_panel['total_inspections'].replace(0, np.nan)\n",
    "\n",
    "viol_panel = violations_merged.groupby(['district', 'year'], as_index=False).agg(\n",
    "    total_violations=('api_norm', 'size'),\n",
    "    major_violations=('major_viol_ind', lambda x: (x == 'Y').sum()),\n",
    "    compliant_on_reinsp=('compliant_on_reinsp', lambda x: (x == 'Y').sum()),\n",
    "    avg_days_to_enforcement=('days_to_enforcement', 'mean'),\n",
    "    med_days_to_enforcement=('days_to_enforcement', 'median')\n",
    ")\n",
    "\n",
    "panel = insp_panel.merge(viol_panel, on=['district', 'year'], how='left')\n",
    "for c in ['total_violations', 'major_violations', 'compliant_on_reinsp']:\n",
    "    panel[c] = panel[c].fillna(0)\n",
    "\n",
    "panel['violations_per_inspection'] = panel['total_violations'] / panel['total_inspections'].replace(0, np.nan)\n",
    "panel['major_share'] = panel['major_violations'] / panel['total_violations'].replace(0, np.nan)\n",
    "panel['resolution_rate'] = panel['compliant_on_reinsp'] / panel['total_violations'].replace(0, np.nan)\n",
    "\n",
    "panel['post_2019'] = (panel['year'] >= POLICY_YEAR).astype(int)\n",
    "panel['year_from_policy'] = panel['year'] - POLICY_YEAR\n",
    "panel['post_trend'] = np.where(panel['year'] >= POLICY_YEAR, panel['year_from_policy'], 0)\n",
    "panel['border_district'] = panel['district'].astype(str).str.strip().isin(BORDER_DISTRICTS).astype(int)\n",
    "\n",
    "# Exposure-based fallback/alternative: district is border if >=25% of inspections are border-exposed.\n",
    "panel['border_exposure_district'] = (panel['share_border_exposed_insp'] >= 0.25).astype(int)\n",
    "\n",
    "print('Panel built:')\n",
    "print(f'  Rows: {len(panel):,}')\n",
    "print(f'  Districts: {panel.district.nunique()}')\n",
    "print(f'  Years: {panel.year.min()}-{panel.year.max()}')\n",
    "print(panel.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "40777bbd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "District competitor links created:\n",
      "  rows: 11\n",
      "  file: /home/dadams/Repos/texas-borderlands/data/district_competitor_links.csv\n",
      "district competitor_jurisdiction border_type  n_border_wells  n_wells_district  raw_weight  weight\n",
      "      01                     MEX       TX-MX            7377             31898      0.2313  1.0000\n",
      "      03                      LA       TX-LA            1947             16700      0.1166  1.0000\n",
      "      04                     MEX       TX-MX           13389             20973      0.6384  1.0000\n",
      "      05                      OK       TX-OK              22              9938      0.0022  1.0000\n",
      "      06                      LA       TX-LA           12786             24422      0.5235  0.9470\n",
      "      06                      OK       TX-OK             716             24422      0.0293  0.0530\n",
      "      08                     MEX       TX-MX               6            105931      0.0001  0.0003\n",
      "      08                      NM       TX-NM           20180            105931      0.1905  0.9997\n",
      "      10                      NM       TX-NM              27             29621      0.0009  0.0030\n",
      "      10                      OK       TX-OK            8946             29621      0.3020  0.9970\n",
      "      8A                      NM       TX-NM           17567             42005      0.4182  1.0000\n",
      "  template: /home/dadams/Repos/texas-borderlands/data/competition_panel_template.csv\n"
     ]
    }
   ],
   "source": [
    "# Build district -> competitor jurisdiction link table (for competition models)\n",
    "# We assign wells to their modal district from inspections, then summarize border subtype shares by district.\n",
    "well_district = inspections[['api_norm', 'district']].dropna().copy()\n",
    "well_district['api_norm'] = well_district['api_norm'].astype(str).str.strip()\n",
    "well_district['district'] = well_district['district'].astype(str).str.strip()\n",
    "well_district = (\n",
    "    well_district.groupby('api_norm')['district']\n",
    "    .agg(lambda s: s.mode().iat[0] if not s.mode().empty else s.iloc[0])\n",
    "    .reset_index()\n",
    ")\n",
    "\n",
    "wb = well_border.merge(well_district, on='api_norm', how='left')\n",
    "wb = wb.dropna(subset=['district']).copy()\n",
    "\n",
    "agg = wb.groupby('district', as_index=False).agg(\n",
    "    n_wells=('api_norm', 'nunique'),\n",
    "    n_tx_mex=('within_50km_texmex', 'sum'),\n",
    "    n_tx_nm=('within_50km_nm_border', 'sum'),\n",
    "    n_tx_ok=('within_50km_ok_border', 'sum'),\n",
    "    n_tx_la=('within_50km_la_border', 'sum'),\n",
    ")\n",
    "for c in ['n_tx_mex', 'n_tx_nm', 'n_tx_ok', 'n_tx_la']:\n",
    "    agg[f'share_{c}'] = agg[c] / agg['n_wells'].replace(0, np.nan)\n",
    "\n",
    "rows = []\n",
    "for _, r in agg.iterrows():\n",
    "    district = str(r['district'])\n",
    "    for jur, count_col in [\n",
    "        ('MEX', 'n_tx_mex'),\n",
    "        ('NM', 'n_tx_nm'),\n",
    "        ('OK', 'n_tx_ok'),\n",
    "        ('LA', 'n_tx_la'),\n",
    "    ]:\n",
    "        n = int(r[count_col]) if pd.notna(r[count_col]) else 0\n",
    "        if n > 0:\n",
    "            rows.append({\n",
    "                'district': district,\n",
    "                'competitor_jurisdiction': jur,\n",
    "                'border_type': f'TX-{jur}' if jur != 'MEX' else 'TX-MX',\n",
    "                'n_border_wells': n,\n",
    "                'n_wells_district': int(r['n_wells']),\n",
    "                'raw_weight': n / r['n_wells'] if r['n_wells'] else np.nan,\n",
    "            })\n",
    "\n",
    "links = pd.DataFrame(rows)\n",
    "if len(links) > 0:\n",
    "    links['weight'] = links['raw_weight'] / links.groupby('district')['raw_weight'].transform('sum')\n",
    "else:\n",
    "    links['weight'] = []\n",
    "\n",
    "links = links.sort_values(['district', 'competitor_jurisdiction']).reset_index(drop=True)\n",
    "link_path = Path('../data/district_competitor_links.csv')\n",
    "link_path.parent.mkdir(parents=True, exist_ok=True)\n",
    "links.to_csv(link_path, index=False)\n",
    "\n",
    "print('District competitor links created:')\n",
    "print(f'  rows: {len(links)}')\n",
    "print(f'  file: {link_path.resolve()}')\n",
    "if len(links) > 0:\n",
    "    print(links.head(20).to_string(index=False))\n",
    "\n",
    "# Also write a template for external competitor metrics by district-year-jurisdiction\n",
    "template = panel[['district', 'year']].drop_duplicates().merge(\n",
    "    links[['district', 'competitor_jurisdiction', 'weight']],\n",
    "    on='district',\n",
    "    how='left'\n",
    ")\n",
    "template['competitor_stringency'] = np.nan\n",
    "template['stringency_metric'] = 'inspection_intensity'\n",
    "template['notes'] = 'Fill competitor_stringency from external jurisdiction data using matching metric definitions.'\n",
    "\n",
    "template_path = Path('../data/competition_panel_template.csv')\n",
    "template.to_csv(template_path, index=False)\n",
    "print(f'  template: {template_path.resolve()}')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "60b23843",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved panel to: /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/district_year_panel_borderlands.csv\n"
     ]
    }
   ],
   "source": [
    "# Persist panel for reuse\n",
    "panel_out = OUTPUT_DIR / 'district_year_panel_borderlands.csv'\n",
    "panel.to_csv(panel_out, index=False)\n",
    "print(f'Saved panel to: {panel_out.resolve()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ec589a6",
   "metadata": {},
   "source": [
    "## RQ1: Border Gap (Descriptive + Regression)\n",
    "\n",
    "Primary binary indicator in baseline models: `border_district`.\n",
    "Alternative indicator available for sensitivity checks: `border_exposure_district`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "27421750",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>group</th>\n",
       "      <th>obs</th>\n",
       "      <th>mean_inspection_intensity</th>\n",
       "      <th>mean_violations_per_inspection</th>\n",
       "      <th>mean_days_to_enf</th>\n",
       "      <th>mean_resolution_rate</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>non-border</td>\n",
       "      <td>66</td>\n",
       "      <td>1.5151</td>\n",
       "      <td>0.0984</td>\n",
       "      <td>122.7522</td>\n",
       "      <td>0.5961</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>border</td>\n",
       "      <td>77</td>\n",
       "      <td>1.3292</td>\n",
       "      <td>0.1301</td>\n",
       "      <td>145.1754</td>\n",
       "      <td>0.5431</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        group  obs  mean_inspection_intensity  mean_violations_per_inspection  \\\n",
       "0  non-border   66                     1.5151                          0.0984   \n",
       "1      border   77                     1.3292                          0.1301   \n",
       "\n",
       "   mean_days_to_enf  mean_resolution_rate  \n",
       "0          122.7522                0.5961  \n",
       "1          145.1754                0.5431  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Border vs non-border descriptive summary\n",
    "summary = panel.groupby('border_district', as_index=False).agg(\n",
    "    obs=('district', 'size'),\n",
    "    mean_inspection_intensity=('inspection_intensity', 'mean'),\n",
    "    mean_violations_per_inspection=('violations_per_inspection', 'mean'),\n",
    "    mean_days_to_enf=('avg_days_to_enforcement', 'mean'),\n",
    "    mean_resolution_rate=('resolution_rate', 'mean')\n",
    ")\n",
    "summary['group'] = summary['border_district'].map({0:'non-border',1:'border'})\n",
    "summary = summary[['group','obs','mean_inspection_intensity','mean_violations_per_inspection','mean_days_to_enf','mean_resolution_rate']]\n",
    "summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "10352175",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>outcome</th>\n",
       "      <th>coef_border_district</th>\n",
       "      <th>pvalue</th>\n",
       "      <th>nobs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>inspection_intensity</td>\n",
       "      <td>-0.1755</td>\n",
       "      <td>0.0999</td>\n",
       "      <td>143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>violations_per_inspection</td>\n",
       "      <td>0.0434</td>\n",
       "      <td>0.0949</td>\n",
       "      <td>143</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     outcome  coef_border_district  pvalue  nobs\n",
       "0       inspection_intensity               -0.1755  0.0999   143\n",
       "1  violations_per_inspection                0.0434  0.0949   143"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# RQ1 models: levels with year FE and district controls via size proxy\n",
    "rq1_df = panel.dropna(subset=['inspection_intensity', 'violations_per_inspection']).copy()\n",
    "rq1_df['log_unique_wells'] = np.log(rq1_df['unique_wells'].clip(lower=1))\n",
    "\n",
    "m_rq1_insp = smf.ols(\n",
    "    'inspection_intensity ~ border_district + log_unique_wells + C(year)',\n",
    "    data=rq1_df\n",
    ").fit(cov_type='cluster', cov_kwds={'groups': rq1_df['district']})\n",
    "\n",
    "m_rq1_viol = smf.ols(\n",
    "    'violations_per_inspection ~ border_district + log_unique_wells + C(year)',\n",
    "    data=rq1_df\n",
    ").fit(cov_type='cluster', cov_kwds={'groups': rq1_df['district']})\n",
    "\n",
    "rq1_results = pd.DataFrame([\n",
    "    {\n",
    "        'outcome':'inspection_intensity',\n",
    "        'coef_border_district': m_rq1_insp.params.get('border_district', np.nan),\n",
    "        'pvalue': m_rq1_insp.pvalues.get('border_district', np.nan),\n",
    "        'nobs': int(m_rq1_insp.nobs)\n",
    "    },\n",
    "    {\n",
    "        'outcome':'violations_per_inspection',\n",
    "        'coef_border_district': m_rq1_viol.params.get('border_district', np.nan),\n",
    "        'pvalue': m_rq1_viol.pvalues.get('border_district', np.nan),\n",
    "        'nobs': int(m_rq1_viol.nobs)\n",
    "    }\n",
    "])\n",
    "\n",
    "rq1_results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51bd6328",
   "metadata": {},
   "source": [
    "## RQ2: Disclosure Reform x Borderlands\n",
    "\n",
    "Specification uses district FE + year FE with interactions:\n",
    "- `post_2019:border_district`\n",
    "- `post_trend:border_district`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2bc00349",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>outcome</th>\n",
       "      <th>coef_post_x_border</th>\n",
       "      <th>p_post_x_border</th>\n",
       "      <th>coef_posttrend_x_border</th>\n",
       "      <th>p_posttrend_x_border</th>\n",
       "      <th>nobs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>inspection_intensity</td>\n",
       "      <td>-0.1191</td>\n",
       "      <td>0.0753</td>\n",
       "      <td>-0.0052</td>\n",
       "      <td>0.8181</td>\n",
       "      <td>143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>violations_per_inspection</td>\n",
       "      <td>0.0040</td>\n",
       "      <td>0.8881</td>\n",
       "      <td>-0.0012</td>\n",
       "      <td>0.8350</td>\n",
       "      <td>143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>avg_days_to_enforcement</td>\n",
       "      <td>-74.5893</td>\n",
       "      <td>0.0156</td>\n",
       "      <td>-1.1587</td>\n",
       "      <td>0.9252</td>\n",
       "      <td>143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>resolution_rate</td>\n",
       "      <td>0.0404</td>\n",
       "      <td>0.4520</td>\n",
       "      <td>-0.0186</td>\n",
       "      <td>0.3404</td>\n",
       "      <td>143</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     outcome  coef_post_x_border  p_post_x_border  \\\n",
       "0       inspection_intensity             -0.1191           0.0753   \n",
       "1  violations_per_inspection              0.0040           0.8881   \n",
       "2    avg_days_to_enforcement            -74.5893           0.0156   \n",
       "3            resolution_rate              0.0404           0.4520   \n",
       "\n",
       "   coef_posttrend_x_border  p_posttrend_x_border  nobs  \n",
       "0                  -0.0052                0.8181   143  \n",
       "1                  -0.0012                0.8350   143  \n",
       "2                  -1.1587                0.9252   143  \n",
       "3                  -0.0186                0.3404   143  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# RQ2 models (policy differential)\n",
    "rq2_df = panel.copy()\n",
    "\n",
    "specs = {\n",
    "    'inspection_intensity': 'inspection_intensity ~ C(district) + C(year) + post_2019:border_district + post_trend:border_district',\n",
    "    'violations_per_inspection': 'violations_per_inspection ~ C(district) + C(year) + post_2019:border_district + post_trend:border_district',\n",
    "    'avg_days_to_enforcement': 'avg_days_to_enforcement ~ C(district) + C(year) + post_2019:border_district + post_trend:border_district',\n",
    "    'resolution_rate': 'resolution_rate ~ C(district) + C(year) + post_2019:border_district + post_trend:border_district',\n",
    "}\n",
    "\n",
    "rows = []\n",
    "models = {}\n",
    "for outcome, formula in specs.items():\n",
    "    d = rq2_df.dropna(subset=[outcome]).copy()\n",
    "    if len(d) < 20:\n",
    "        continue\n",
    "    m = smf.ols(formula, data=d).fit(cov_type='cluster', cov_kwds={'groups': d['district']})\n",
    "    models[outcome] = m\n",
    "    rows.append({\n",
    "        'outcome': outcome,\n",
    "        'coef_post_x_border': m.params.get('post_2019:border_district', np.nan),\n",
    "        'p_post_x_border': m.pvalues.get('post_2019:border_district', np.nan),\n",
    "        'coef_posttrend_x_border': m.params.get('post_trend:border_district', np.nan),\n",
    "        'p_posttrend_x_border': m.pvalues.get('post_trend:border_district', np.nan),\n",
    "        'nobs': int(m.nobs)\n",
    "    })\n",
    "\n",
    "rq2_results = pd.DataFrame(rows)\n",
    "rq2_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "448a0e24",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 101936 characters]",
      "text/plain": [
       "<Figure size 1200x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved figure: /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/border_vs_nonborder_trends.png\n"
     ]
    }
   ],
   "source": [
    "# Quick visualization: average trend by group\n",
    "plot_df = panel.groupby(['year','border_district'], as_index=False).agg(\n",
    "    avg_days_to_enforcement=('avg_days_to_enforcement','mean'),\n",
    "    inspection_intensity=('inspection_intensity','mean')\n",
    ")\n",
    "plot_df['group'] = plot_df['border_district'].map({0:'non-border',1:'border'})\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharex=True)\n",
    "\n",
    "sns.lineplot(data=plot_df, x='year', y='avg_days_to_enforcement', hue='group', marker='o', ax=axes[0])\n",
    "axes[0].axvline(POLICY_YEAR, color='black', linestyle='--', linewidth=1)\n",
    "axes[0].set_title('Mean Days to Enforcement')\n",
    "axes[0].set_xlabel('Year')\n",
    "axes[0].set_ylabel('Days')\n",
    "\n",
    "sns.lineplot(data=plot_df, x='year', y='inspection_intensity', hue='group', marker='o', ax=axes[1])\n",
    "axes[1].axvline(POLICY_YEAR, color='black', linestyle='--', linewidth=1)\n",
    "axes[1].set_title('Mean Inspection Intensity')\n",
    "axes[1].set_xlabel('Year')\n",
    "axes[1].set_ylabel('Inspections per well')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.legend(title='')\n",
    "\n",
    "fig.tight_layout()\n",
    "fig_path = OUTPUT_DIR / 'border_vs_nonborder_trends.png'\n",
    "fig.savefig(fig_path, bbox_inches='tight', dpi=300)\n",
    "plt.show()\n",
    "\n",
    "print(f'Saved figure: {fig_path.resolve()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8db3f662",
   "metadata": {},
   "source": [
    "## Robustness: Which Border, How Much Exposure, and Cutoff Sensitivity\n",
    "\n",
    "This section addresses reviewer concerns about border-definition arbitrariness with three checks:\n",
    "\n",
    "1. **Border-type splits** (TX-MX, TX-NM, TX-OK, TX-LA)\n",
    "2. **Continuous exposure intensity** (`share_border_exposed_insp`)\n",
    "3. **Cutoff sensitivity** (25km / 75km / 100km)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "24349ba4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved: /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/district_border_type_profile.csv\n",
      "district  n_wells  n_tx_mex  n_tx_nm  n_tx_ok  n_tx_la dominant_border_type\n",
      "      01    31898      7377        0        0        0                TX-MX\n",
      "      02    17099         0        0        0        0                 NONE\n",
      "      03    16700         0        0        0     1947                TX-LA\n",
      "      04    20973     13389        0        0        0                TX-MX\n",
      "      05     9938         0        0       22        0                TX-OK\n",
      "      06    24422         0        0      716    12786                TX-LA\n",
      "      08   105931         6    20180        0        0                TX-NM\n",
      "      09    46485         0        0        0        0                 NONE\n",
      "      10    29621         0       27     8946        0                TX-OK\n",
      "      6E     6235         0        0        0        0                 NONE\n",
      "      7B    21230         0        0        0        0                 NONE\n",
      "      7C    43061         0        0        0        0                 NONE\n",
      "      8A    42005         0    17567        0        0                TX-NM\n"
     ]
    }
   ],
   "source": [
    "# Build district border-type profile (TX-MX / TX-NM / TX-OK / TX-LA)\n",
    "well_district = inspections[['api_norm', 'district']].dropna().copy()\n",
    "well_district['api_norm'] = well_district['api_norm'].astype(str).str.strip()\n",
    "well_district['district'] = well_district['district'].astype(str).str.strip()\n",
    "well_district = (\n",
    "    well_district.groupby('api_norm')['district']\n",
    "    .agg(lambda s: s.mode().iat[0] if not s.mode().empty else s.iloc[0])\n",
    "    .reset_index()\n",
    ")\n",
    "\n",
    "wb = well_border.merge(well_district, on='api_norm', how='left').dropna(subset=['district']).copy()\n",
    "\n",
    "border_profile = wb.groupby('district', as_index=False).agg(\n",
    "    n_wells=('api_norm', 'nunique'),\n",
    "    n_tx_mex=('within_50km_texmex', 'sum'),\n",
    "    n_tx_nm=('within_50km_nm_border', 'sum'),\n",
    "    n_tx_ok=('within_50km_ok_border', 'sum'),\n",
    "    n_tx_la=('within_50km_la_border', 'sum'),\n",
    ")\n",
    "\n",
    "for c in ['n_tx_mex', 'n_tx_nm', 'n_tx_ok', 'n_tx_la']:\n",
    "    border_profile[f'share_{c}'] = border_profile[c] / border_profile['n_wells'].replace(0, np.nan)\n",
    "\n",
    "# Binary border-type membership by district (any well within 50km)\n",
    "border_profile['has_tx_mex'] = (border_profile['n_tx_mex'] > 0).astype(int)\n",
    "border_profile['has_tx_nm'] = (border_profile['n_tx_nm'] > 0).astype(int)\n",
    "border_profile['has_tx_ok'] = (border_profile['n_tx_ok'] > 0).astype(int)\n",
    "border_profile['has_tx_la'] = (border_profile['n_tx_la'] > 0).astype(int)\n",
    "\n",
    "# Dominant border type (highest share among the four)\n",
    "share_cols = ['share_n_tx_mex', 'share_n_tx_nm', 'share_n_tx_ok', 'share_n_tx_la']\n",
    "label_map = {\n",
    "    'share_n_tx_mex': 'TX-MX',\n",
    "    'share_n_tx_nm': 'TX-NM',\n",
    "    'share_n_tx_ok': 'TX-OK',\n",
    "    'share_n_tx_la': 'TX-LA',\n",
    "}\n",
    "\n",
    "border_profile['dominant_border_type'] = border_profile[share_cols].idxmax(axis=1).map(label_map)\n",
    "border_profile.loc[border_profile[share_cols].max(axis=1).fillna(0) <= 0, 'dominant_border_type'] = 'NONE'\n",
    "\n",
    "panel_robust = panel.merge(\n",
    "    border_profile[\n",
    "        ['district', 'n_wells'] + share_cols + ['has_tx_mex', 'has_tx_nm', 'has_tx_ok', 'has_tx_la', 'dominant_border_type']\n",
    "    ],\n",
    "    on='district',\n",
    "    how='left'\n",
    ")\n",
    "\n",
    "profile_path = OUTPUT_DIR / 'district_border_type_profile.csv'\n",
    "border_profile.to_csv(profile_path, index=False)\n",
    "\n",
    "print(f'Saved: {profile_path.resolve()}')\n",
    "print(border_profile[['district','n_wells','n_tx_mex','n_tx_nm','n_tx_ok','n_tx_la','dominant_border_type']].to_string(index=False))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ad78072c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved: /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/border_type_split_results.csv\n",
      "                   model_family                 outcome                  term     coef  pvalue  nobs\n",
      "         RQ1_levels_border_type    inspection_intensity            has_tx_mex   0.0400  0.7893   143\n",
      "         RQ1_levels_border_type    inspection_intensity             has_tx_nm  -0.0535  0.5341   143\n",
      "         RQ1_levels_border_type    inspection_intensity             has_tx_ok  -0.0133  0.8359   143\n",
      "         RQ1_levels_border_type    inspection_intensity             has_tx_la   0.0660  0.3727   143\n",
      "         RQ1_levels_border_type avg_days_to_enforcement            has_tx_mex  43.1955  0.3092   143\n",
      "         RQ1_levels_border_type avg_days_to_enforcement             has_tx_nm -73.7579  0.0611   143\n",
      "         RQ1_levels_border_type avg_days_to_enforcement             has_tx_ok  88.5376  0.0210   143\n",
      "         RQ1_levels_border_type avg_days_to_enforcement             has_tx_la  26.1342  0.6193   143\n",
      "RQ2_FE_border_type_interactions    inspection_intensity  post_2019:has_tx_mex   0.0821  0.5405   143\n",
      "RQ2_FE_border_type_interactions    inspection_intensity   post_2019:has_tx_nm  -0.1152  0.1219   143\n",
      "RQ2_FE_border_type_interactions    inspection_intensity   post_2019:has_tx_ok   0.0654  0.3045   143\n",
      "RQ2_FE_border_type_interactions    inspection_intensity   post_2019:has_tx_la   0.0197  0.7085   143\n",
      "RQ2_FE_border_type_interactions    inspection_intensity post_trend:has_tx_mex  -0.0195  0.3109   143\n",
      "RQ2_FE_border_type_interactions    inspection_intensity  post_trend:has_tx_nm   0.0051  0.6920   143\n",
      "RQ2_FE_border_type_interactions    inspection_intensity  post_trend:has_tx_ok  -0.0000  0.9976   143\n",
      "RQ2_FE_border_type_interactions    inspection_intensity  post_trend:has_tx_la  -0.0197  0.2809   143\n",
      "RQ2_FE_border_type_interactions avg_days_to_enforcement  post_2019:has_tx_mex   4.0900  0.9062   143\n",
      "RQ2_FE_border_type_interactions avg_days_to_enforcement   post_2019:has_tx_nm -18.7442  0.6013   143\n",
      "RQ2_FE_border_type_interactions avg_days_to_enforcement   post_2019:has_tx_ok -14.2446  0.8134   143\n",
      "RQ2_FE_border_type_interactions avg_days_to_enforcement   post_2019:has_tx_la -43.6598  0.6415   143\n",
      "RQ2_FE_border_type_interactions avg_days_to_enforcement post_trend:has_tx_mex  -0.0148  0.9991   143\n",
      "RQ2_FE_border_type_interactions avg_days_to_enforcement  post_trend:has_tx_nm  22.9067  0.0189   143\n",
      "RQ2_FE_border_type_interactions avg_days_to_enforcement  post_trend:has_tx_ok -16.7188  0.0794   143\n",
      "RQ2_FE_border_type_interactions avg_days_to_enforcement  post_trend:has_tx_la   0.6415  0.9551   143\n"
     ]
    }
   ],
   "source": [
    "# Border-type split models (which border drives results?)\n",
    "\n",
    "def _fit_ols(formula, data, group_col='district'):\n",
    "    d = data.copy()\n",
    "    d = d.dropna()\n",
    "    if len(d) < 20:\n",
    "        return None\n",
    "    try:\n",
    "        return smf.ols(formula, data=d).fit(cov_type='cluster', cov_kwds={'groups': d[group_col]})\n",
    "    except Exception:\n",
    "        return None\n",
    "\n",
    "rows = []\n",
    "outcomes = ['inspection_intensity', 'avg_days_to_enforcement']\n",
    "\n",
    "# RQ1-style levels models with type indicators\n",
    "for outcome in outcomes:\n",
    "    d = panel_robust.dropna(subset=[outcome, 'unique_wells']).copy()\n",
    "    d['log_unique_wells'] = np.log(d['unique_wells'].clip(lower=1))\n",
    "    m = _fit_ols(\n",
    "        f\"{outcome} ~ has_tx_mex + has_tx_nm + has_tx_ok + has_tx_la + log_unique_wells + C(year)\",\n",
    "        d\n",
    "    )\n",
    "    if m is not None:\n",
    "        for term in ['has_tx_mex', 'has_tx_nm', 'has_tx_ok', 'has_tx_la']:\n",
    "            rows.append({\n",
    "                'model_family': 'RQ1_levels_border_type',\n",
    "                'outcome': outcome,\n",
    "                'term': term,\n",
    "                'coef': m.params.get(term, np.nan),\n",
    "                'pvalue': m.pvalues.get(term, np.nan),\n",
    "                'nobs': int(m.nobs),\n",
    "            })\n",
    "\n",
    "# RQ2-style FE interaction models with type indicators\n",
    "for outcome in outcomes:\n",
    "    d = panel_robust.dropna(subset=[outcome]).copy()\n",
    "    m = _fit_ols(\n",
    "        f\"{outcome} ~ C(district) + C(year) + \"\n",
    "        f\"post_2019:has_tx_mex + post_2019:has_tx_nm + post_2019:has_tx_ok + post_2019:has_tx_la + \"\n",
    "        f\"post_trend:has_tx_mex + post_trend:has_tx_nm + post_trend:has_tx_ok + post_trend:has_tx_la\",\n",
    "        d\n",
    "    )\n",
    "    if m is not None:\n",
    "        for term in [\n",
    "            'post_2019:has_tx_mex', 'post_2019:has_tx_nm', 'post_2019:has_tx_ok', 'post_2019:has_tx_la',\n",
    "            'post_trend:has_tx_mex', 'post_trend:has_tx_nm', 'post_trend:has_tx_ok', 'post_trend:has_tx_la',\n",
    "        ]:\n",
    "            rows.append({\n",
    "                'model_family': 'RQ2_FE_border_type_interactions',\n",
    "                'outcome': outcome,\n",
    "                'term': term,\n",
    "                'coef': m.params.get(term, np.nan),\n",
    "                'pvalue': m.pvalues.get(term, np.nan),\n",
    "                'nobs': int(m.nobs),\n",
    "            })\n",
    "\n",
    "border_type_results = pd.DataFrame(rows)\n",
    "res_path = OUTPUT_DIR / 'border_type_split_results.csv'\n",
    "border_type_results.to_csv(res_path, index=False)\n",
    "\n",
    "print(f'Saved: {res_path.resolve()}')\n",
    "print(border_type_results.head(40).to_string(index=False))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "fde06b85",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved: /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/continuous_exposure_results.csv\n",
      "         model_family                   outcome                                 term      coef  pvalue  nobs\n",
      "RQ1_levels_continuous      inspection_intensity            share_border_exposed_insp    0.2095  0.4757   143\n",
      "RQ1_levels_continuous   avg_days_to_enforcement            share_border_exposed_insp  103.4683  0.5710   143\n",
      "RQ1_levels_continuous violations_per_inspection            share_border_exposed_insp   -0.1585  0.0144   143\n",
      "RQ1_levels_continuous           resolution_rate            share_border_exposed_insp   -0.0420  0.8619   143\n",
      "    RQ2_FE_continuous      inspection_intensity  post_2019:share_border_exposed_insp    0.0348  0.8650   143\n",
      "    RQ2_FE_continuous      inspection_intensity post_trend:share_border_exposed_insp   -0.0129  0.8281   143\n",
      "    RQ2_FE_continuous   avg_days_to_enforcement  post_2019:share_border_exposed_insp -109.4067  0.4449   143\n",
      "    RQ2_FE_continuous   avg_days_to_enforcement post_trend:share_border_exposed_insp   13.9623  0.7415   143\n",
      "    RQ2_FE_continuous violations_per_inspection  post_2019:share_border_exposed_insp   -0.0053  0.9390   143\n",
      "    RQ2_FE_continuous violations_per_inspection post_trend:share_border_exposed_insp    0.0096  0.5012   143\n",
      "    RQ2_FE_continuous           resolution_rate  post_2019:share_border_exposed_insp   -0.0322  0.8163   143\n",
      "    RQ2_FE_continuous           resolution_rate post_trend:share_border_exposed_insp   -0.0979  0.0423   143\n"
     ]
    }
   ],
   "source": [
    "# Continuous exposure models (share of inspections near border)\n",
    "rows = []\n",
    "outcomes = ['inspection_intensity', 'avg_days_to_enforcement', 'violations_per_inspection', 'resolution_rate']\n",
    "\n",
    "# RQ1-style levels with continuous share exposure\n",
    "for outcome in outcomes:\n",
    "    d = panel.dropna(subset=[outcome, 'share_border_exposed_insp', 'unique_wells']).copy()\n",
    "    d['log_unique_wells'] = np.log(d['unique_wells'].clip(lower=1))\n",
    "    m = smf.ols(\n",
    "        f\"{outcome} ~ share_border_exposed_insp + log_unique_wells + C(year)\",\n",
    "        data=d\n",
    "    ).fit(cov_type='cluster', cov_kwds={'groups': d['district']})\n",
    "    rows.append({\n",
    "        'model_family': 'RQ1_levels_continuous',\n",
    "        'outcome': outcome,\n",
    "        'term': 'share_border_exposed_insp',\n",
    "        'coef': m.params.get('share_border_exposed_insp', np.nan),\n",
    "        'pvalue': m.pvalues.get('share_border_exposed_insp', np.nan),\n",
    "        'nobs': int(m.nobs),\n",
    "    })\n",
    "\n",
    "# RQ2-style FE with continuous interactions\n",
    "for outcome in outcomes:\n",
    "    d = panel.dropna(subset=[outcome, 'share_border_exposed_insp']).copy()\n",
    "    m = smf.ols(\n",
    "        f\"{outcome} ~ C(district) + C(year) + post_2019:share_border_exposed_insp + post_trend:share_border_exposed_insp\",\n",
    "        data=d\n",
    "    ).fit(cov_type='cluster', cov_kwds={'groups': d['district']})\n",
    "    rows.append({\n",
    "        'model_family': 'RQ2_FE_continuous',\n",
    "        'outcome': outcome,\n",
    "        'term': 'post_2019:share_border_exposed_insp',\n",
    "        'coef': m.params.get('post_2019:share_border_exposed_insp', np.nan),\n",
    "        'pvalue': m.pvalues.get('post_2019:share_border_exposed_insp', np.nan),\n",
    "        'nobs': int(m.nobs),\n",
    "    })\n",
    "    rows.append({\n",
    "        'model_family': 'RQ2_FE_continuous',\n",
    "        'outcome': outcome,\n",
    "        'term': 'post_trend:share_border_exposed_insp',\n",
    "        'coef': m.params.get('post_trend:share_border_exposed_insp', np.nan),\n",
    "        'pvalue': m.pvalues.get('post_trend:share_border_exposed_insp', np.nan),\n",
    "        'nobs': int(m.nobs),\n",
    "    })\n",
    "\n",
    "continuous_results = pd.DataFrame(rows)\n",
    "res_path = OUTPUT_DIR / 'continuous_exposure_results.csv'\n",
    "continuous_results.to_csv(res_path, index=False)\n",
    "\n",
    "print(f'Saved: {res_path.resolve()}')\n",
    "print(continuous_results.to_string(index=False))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "b5cb889f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved: /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/cutoff_sensitivity_results.csv\n",
      " cutoff_km model_family                 outcome                          term      coef  pvalue  nobs\n",
      "        25   RQ1_levels    inspection_intensity             share_border_25km    0.5386  0.4228   143\n",
      "        25       RQ2_FE avg_days_to_enforcement   post_2019:share_border_25km -101.9283  0.7010   143\n",
      "        25       RQ2_FE avg_days_to_enforcement  post_trend:share_border_25km   30.6087  0.7698   143\n",
      "        75   RQ1_levels    inspection_intensity             share_border_75km    0.0794  0.6963   143\n",
      "        75       RQ2_FE avg_days_to_enforcement   post_2019:share_border_75km  -75.6591  0.5116   143\n",
      "        75       RQ2_FE avg_days_to_enforcement  post_trend:share_border_75km   -0.0090  0.9997   143\n",
      "       100   RQ1_levels    inspection_intensity            share_border_100km    0.0022  0.9872   143\n",
      "       100       RQ2_FE avg_days_to_enforcement  post_2019:share_border_100km   -4.4795  0.9474   143\n",
      "       100       RQ2_FE avg_days_to_enforcement post_trend:share_border_100km  -12.6587  0.5147   143\n"
     ]
    }
   ],
   "source": [
    "# Cutoff sensitivity (25km / 75km / 100km)\n",
    "dist_cols = [\n",
    "    'distance_to_texmex_km',\n",
    "    'distance_to_nm_border_km',\n",
    "    'distance_to_ok_border_km',\n",
    "    'distance_to_la_border_km',\n",
    "]\n",
    "\n",
    "wb_cut = well_border[['api_norm'] + dist_cols].copy()\n",
    "wb_cut['min_border_distance_km'] = wb_cut[dist_cols].min(axis=1, skipna=True)\n",
    "\n",
    "# If all distances are missing, keep as NaN (not border-exposed)\n",
    "all_missing = wb_cut[dist_cols].isna().all(axis=1)\n",
    "wb_cut.loc[all_missing, 'min_border_distance_km'] = np.nan\n",
    "\n",
    "cutoffs = [25, 50, 75, 100]\n",
    "for c in cutoffs:\n",
    "    wb_cut[f'within_{c}km_any_border'] = (wb_cut['min_border_distance_km'] <= c).astype(int)\n",
    "\n",
    "insp_base = inspections[['api_norm', 'district', 'year', 'compliance']].copy()\n",
    "insp_cut = insp_base.merge(wb_cut[['api_norm'] + [f'within_{c}km_any_border' for c in cutoffs]], on='api_norm', how='left')\n",
    "for c in cutoffs:\n",
    "    insp_cut[f'within_{c}km_any_border'] = pd.to_numeric(insp_cut[f'within_{c}km_any_border'], errors='coerce').fillna(0).astype(int)\n",
    "\n",
    "rows = []\n",
    "for c in [25, 75, 100]:\n",
    "    share_df = (\n",
    "        insp_cut.groupby(['district', 'year'], as_index=False)[f'within_{c}km_any_border']\n",
    "        .mean()\n",
    "        .rename(columns={f'within_{c}km_any_border': f'share_border_{c}km'})\n",
    "    )\n",
    "    d = panel.merge(share_df, on=['district', 'year'], how='left')\n",
    "\n",
    "    # RQ1 levels model (inspection intensity)\n",
    "    d1 = d.dropna(subset=['inspection_intensity', f'share_border_{c}km', 'unique_wells']).copy()\n",
    "    d1['log_unique_wells'] = np.log(d1['unique_wells'].clip(lower=1))\n",
    "    m1 = smf.ols(\n",
    "        f\"inspection_intensity ~ {f'share_border_{c}km'} + log_unique_wells + C(year)\",\n",
    "        data=d1\n",
    "    ).fit(cov_type='cluster', cov_kwds={'groups': d1['district']})\n",
    "\n",
    "    rows.append({\n",
    "        'cutoff_km': c,\n",
    "        'model_family': 'RQ1_levels',\n",
    "        'outcome': 'inspection_intensity',\n",
    "        'term': f'share_border_{c}km',\n",
    "        'coef': m1.params.get(f'share_border_{c}km', np.nan),\n",
    "        'pvalue': m1.pvalues.get(f'share_border_{c}km', np.nan),\n",
    "        'nobs': int(m1.nobs),\n",
    "    })\n",
    "\n",
    "    # RQ2 FE model (timing) - main paper signal\n",
    "    d2 = d.dropna(subset=['avg_days_to_enforcement', f'share_border_{c}km']).copy()\n",
    "    m2 = smf.ols(\n",
    "        f\"avg_days_to_enforcement ~ C(district) + C(year) + post_2019:{f'share_border_{c}km'} + post_trend:{f'share_border_{c}km'}\",\n",
    "        data=d2\n",
    "    ).fit(cov_type='cluster', cov_kwds={'groups': d2['district']})\n",
    "\n",
    "    rows.append({\n",
    "        'cutoff_km': c,\n",
    "        'model_family': 'RQ2_FE',\n",
    "        'outcome': 'avg_days_to_enforcement',\n",
    "        'term': f'post_2019:share_border_{c}km',\n",
    "        'coef': m2.params.get(f'post_2019:share_border_{c}km', np.nan),\n",
    "        'pvalue': m2.pvalues.get(f'post_2019:share_border_{c}km', np.nan),\n",
    "        'nobs': int(m2.nobs),\n",
    "    })\n",
    "    rows.append({\n",
    "        'cutoff_km': c,\n",
    "        'model_family': 'RQ2_FE',\n",
    "        'outcome': 'avg_days_to_enforcement',\n",
    "        'term': f'post_trend:share_border_{c}km',\n",
    "        'coef': m2.params.get(f'post_trend:share_border_{c}km', np.nan),\n",
    "        'pvalue': m2.pvalues.get(f'post_trend:share_border_{c}km', np.nan),\n",
    "        'nobs': int(m2.nobs),\n",
    "    })\n",
    "\n",
    "cutoff_results = pd.DataFrame(rows)\n",
    "res_path = OUTPUT_DIR / 'cutoff_sensitivity_results.csv'\n",
    "cutoff_results.to_csv(res_path, index=False)\n",
    "\n",
    "print(f'Saved: {res_path.resolve()}')\n",
    "print(cutoff_results.to_string(index=False))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "d603b1f7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 99832 characters]",
      "text/plain": [
       "<Figure size 850x520 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved figure: /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/money_plot_timing_border_prepost2019.png\n",
      "Saved table:  /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/money_border_type_timing_results.csv\n",
      "Saved table:  /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/money_cutoff_timing_results.csv\n",
      "Saved table:  /home/dadams/Repos/texas-borderlands/analysis/output_borderlands/money_plot_timing_ci_by_year.csv\n",
      "\n",
      "Border-type timing interactions:\n",
      "                 term     coef  pvalue\n",
      " post_2019:has_tx_mex   4.0900  0.9062\n",
      "  post_2019:has_tx_nm -18.7442  0.6013\n",
      "  post_2019:has_tx_ok -14.2446  0.8134\n",
      "  post_2019:has_tx_la -43.6598  0.6415\n",
      "post_trend:has_tx_mex  -0.0148  0.9991\n",
      " post_trend:has_tx_nm  22.9067  0.0189\n",
      " post_trend:has_tx_ok -16.7188  0.0794\n",
      " post_trend:has_tx_la   0.6415  0.9551\n",
      "\n",
      "Timing cutoff sensitivity (post_2019 interaction):\n",
      " cutoff_km      coef  pvalue  nobs\n",
      "        25 -101.9283  0.7010   143\n",
      "        75  -75.6591  0.5116   143\n",
      "       100   -4.4795  0.9474   143\n"
     ]
    }
   ],
   "source": [
    "# Money cell: border-type disaggregation + cutoff sensitivity + single timing figure\n",
    "\n",
    "# 1) Border-type disaggregation focused on timing outcome\n",
    "if 'panel_robust' not in globals():\n",
    "    # Fallback if robustness profile cell was skipped\n",
    "    panel_robust = panel.copy()\n",
    "    for c in ['has_tx_mex','has_tx_nm','has_tx_ok','has_tx_la']:\n",
    "        if c not in panel_robust.columns:\n",
    "            panel_robust[c] = 0\n",
    "\n",
    "timing_border_type_df = panel_robust.dropna(subset=['avg_days_to_enforcement']).copy()\n",
    "\n",
    "m_timing_type = smf.ols(\n",
    "    'avg_days_to_enforcement ~ C(district) + C(year) + '\n",
    "    'post_2019:has_tx_mex + post_2019:has_tx_nm + post_2019:has_tx_ok + post_2019:has_tx_la + '\n",
    "    'post_trend:has_tx_mex + post_trend:has_tx_nm + post_trend:has_tx_ok + post_trend:has_tx_la',\n",
    "    data=timing_border_type_df\n",
    ").fit(cov_type='cluster', cov_kwds={'groups': timing_border_type_df['district']})\n",
    "\n",
    "timing_type_terms = [\n",
    "    'post_2019:has_tx_mex','post_2019:has_tx_nm','post_2019:has_tx_ok','post_2019:has_tx_la',\n",
    "    'post_trend:has_tx_mex','post_trend:has_tx_nm','post_trend:has_tx_ok','post_trend:has_tx_la',\n",
    "]\n",
    "\n",
    "timing_type_table = pd.DataFrame({\n",
    "    'term': timing_type_terms,\n",
    "    'coef': [m_timing_type.params.get(t, np.nan) for t in timing_type_terms],\n",
    "    'pvalue': [m_timing_type.pvalues.get(t, np.nan) for t in timing_type_terms],\n",
    "})\n",
    "\n",
    "# 2) Cutoff sensitivity (timing-only focus)\n",
    "if 'cutoff_results' in globals() and len(cutoff_results) > 0:\n",
    "    timing_cutoff_table = cutoff_results[\n",
    "        (cutoff_results['outcome'] == 'avg_days_to_enforcement') &\n",
    "        (cutoff_results['term'].str.startswith('post_2019:share_border_'))\n",
    "    ][['cutoff_km','coef','pvalue','nobs']].copy().sort_values('cutoff_km')\n",
    "else:\n",
    "    timing_cutoff_table = pd.DataFrame(columns=['cutoff_km','coef','pvalue','nobs'])\n",
    "\n",
    "# 3) Single clean figure: timing outcome, border vs non-border, pre/post 2019, with 95% CI bands\n",
    "plot_base = panel[['district','year','border_district','avg_days_to_enforcement']].dropna().copy()\n",
    "plot_base['group'] = plot_base['border_district'].map({0:'Non-border', 1:'Border'})\n",
    "\n",
    "stats = (\n",
    "    plot_base.groupby(['year','group'], as_index=False)\n",
    "    .agg(\n",
    "        mean_days=('avg_days_to_enforcement','mean'),\n",
    "        sd_days=('avg_days_to_enforcement','std'),\n",
    "        n=('avg_days_to_enforcement','size')\n",
    "    )\n",
    ")\n",
    "stats['se_days'] = stats['sd_days'] / np.sqrt(stats['n'].clip(lower=1))\n",
    "stats['ci95'] = 1.96 * stats['se_days']\n",
    "stats['lower'] = stats['mean_days'] - stats['ci95']\n",
    "stats['upper'] = stats['mean_days'] + stats['ci95']\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8.5, 5.2))\n",
    "colors = {'Non-border':'#4b5563', 'Border':'#dc2626'}\n",
    "\n",
    "for g in ['Non-border', 'Border']:\n",
    "    d = stats[stats['group'] == g].sort_values('year')\n",
    "    if len(d) == 0:\n",
    "        continue\n",
    "    ax.plot(d['year'], d['mean_days'], marker='o', linewidth=2.2, color=colors[g], label=g)\n",
    "    # 95% CI ribbon\n",
    "    ax.fill_between(d['year'], d['lower'], d['upper'], color=colors[g], alpha=0.18, linewidth=0)\n",
    "\n",
    "ax.axvline(POLICY_YEAR, color='black', linestyle='--', linewidth=1.2, alpha=0.8)\n",
    "ax.text(POLICY_YEAR + 0.1, ax.get_ylim()[1]*0.98, '2019 disclosure', fontsize=9, va='top')\n",
    "\n",
    "# Add simple pre/post means annotation\n",
    "prepost = panel.assign(period=np.where(panel['year'] < POLICY_YEAR, 'Pre-2019', 'Post-2019'))     .groupby(['period','border_district'], as_index=False)['avg_days_to_enforcement'].mean()\n",
    "nb_pre = prepost[(prepost['period']=='Pre-2019') & (prepost['border_district']==0)]['avg_days_to_enforcement']\n",
    "b_pre  = prepost[(prepost['period']=='Pre-2019') & (prepost['border_district']==1)]['avg_days_to_enforcement']\n",
    "nb_post= prepost[(prepost['period']=='Post-2019') & (prepost['border_district']==0)]['avg_days_to_enforcement']\n",
    "b_post = prepost[(prepost['period']=='Post-2019') & (prepost['border_district']==1)]['avg_days_to_enforcement']\n",
    "\n",
    "if len(nb_pre) and len(b_pre) and len(nb_post) and len(b_post):\n",
    "    note = (\n",
    "        f\"Pre means (days): non-border={nb_pre.iloc[0]:.1f}, border={b_pre.iloc[0]:.1f}\\n\"\n",
    "        f\"Post means (days): non-border={nb_post.iloc[0]:.1f}, border={b_post.iloc[0]:.1f}\"\n",
    "    )\n",
    "    ax.text(0.01, 0.01, note, transform=ax.transAxes, fontsize=8.5,\n",
    "            va='bottom', ha='left', bbox=dict(facecolor='white', alpha=0.8, edgecolor='#d1d5db'))\n",
    "\n",
    "ax.set_title('Enforcement Timing by Border Status (2015-2025)', fontsize=12, weight='bold')\n",
    "ax.set_xlabel('Year')\n",
    "ax.set_ylabel('Average Days to Enforcement')\n",
    "# Explicit aggregation note to prevent ambiguity in manuscript review\n",
    "ax.text(0.99, 0.01,\n",
    "        'Note: Year-group means use district-year values (equal district weighting), not pooled well-level records.',\n",
    "        transform=ax.transAxes, fontsize=7.8, ha='right', va='bottom',\n",
    "        bbox=dict(facecolor='white', alpha=0.78, edgecolor='#d1d5db'))\n",
    "ax.legend(title='District group', frameon=True)\n",
    "ax.grid(alpha=0.25)\n",
    "fig.tight_layout()\n",
    "\n",
    "money_fig_path = OUTPUT_DIR / 'money_plot_timing_border_prepost2019.png'\n",
    "fig.savefig(money_fig_path, dpi=300, bbox_inches='tight')\n",
    "plt.show()\n",
    "\n",
    "# Save compact robustness tables for paper text\n",
    "money_type_path = OUTPUT_DIR / 'money_border_type_timing_results.csv'\n",
    "money_cutoff_path = OUTPUT_DIR / 'money_cutoff_timing_results.csv'\n",
    "timing_type_table.to_csv(money_type_path, index=False)\n",
    "timing_cutoff_table.to_csv(money_cutoff_path, index=False)\n",
    "\n",
    "# Save year-by-year CI table for transparency\n",
    "ci_path = OUTPUT_DIR / 'money_plot_timing_ci_by_year.csv'\n",
    "stats.to_csv(ci_path, index=False)\n",
    "\n",
    "print(f'Saved figure: {money_fig_path.resolve()}')\n",
    "print(f'Saved table:  {money_type_path.resolve()}')\n",
    "print(f'Saved table:  {money_cutoff_path.resolve()}')\n",
    "print(f'Saved table:  {ci_path.resolve()}')\n",
    "print('Aggregation note: Year-group lines/CI are computed from district-year averages (equal district weighting).')\n",
    "print('\\nBorder-type timing interactions:')\n",
    "print(timing_type_table.to_string(index=False))\n",
    "print('\\nTiming cutoff sensitivity (post_2019 interaction):')\n",
    "print(timing_cutoff_table.to_string(index=False))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4719bed6",
   "metadata": {},
   "source": [
    "Figure X. Enforcement timing by border status with 95% confidence intervals (2015–2025).\n",
    "Lines plot annual group means of average days to enforcement for border versus non-border districts. Points indicate yearly means; vertical uncertainty reflects 95% confidence intervals around the mean (mean ± 1.96 × SE), shown as shaded ribbons. The dashed line marks the 2019 disclosure reform. Pre-2019, border districts exhibit slower enforcement timing; post-2019, the timing gap collapses and slightly reverses, consistent with a border-specific post-reform level shift rather than a border-specific trend change."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9460e0e",
   "metadata": {},
   "source": [
    "## RQ3/RQ4 (Woods Reaction Function): Not Estimable in This Study Build\n",
    "\n",
    "This notebook does **not** estimate a Woods-style reaction function because we do not have cross-jurisdiction competitor enforcement metrics.\n",
    "\n",
    "What we *can* do with current data:\n",
    "- Estimate borderlands gaps in Texas district outcomes (inspection intensity, violations per inspection, enforcement timing, compliance/resolution).\n",
    "- Estimate whether 2019 disclosure changed those outcomes differently in border districts (`post_2019:border` and `post_trend:border`).\n",
    "\n",
    "What we *cannot* do without competitor data:\n",
    "- Estimate `Enforcement Gap(t) = Texas_border_district_stringency(t) - Competitor_stringency(t)`.\n",
    "- Test asymmetric response to positive vs negative gap (Woods reaction function).\n",
    "\n",
    "Interpretation can still use Woods/borderlands theory qualitatively, but causal reaction-function claims are out of scope for this build.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "1bcf5095",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RQ3/RQ4 reaction-function model is intentionally disabled in this notebook build.\n",
      "Reason: no cross-jurisdiction competitor_stringency data available.\n",
      "Scope retained: RQ1/RQ2 border gaps and post-2019 differential trends within Texas districts.\n",
      "Reference files (for future data collection):\n",
      "  /home/dadams/Repos/texas-borderlands/data/district_competitor_links.csv (exists=True)\n",
      "  /home/dadams/Repos/texas-borderlands/data/competition_panel_template.csv (exists=True)\n"
     ]
    }
   ],
   "source": [
    "# RQ3/RQ4 not estimated (no competitor stringency data)\n",
    "print('RQ3/RQ4 reaction-function model is intentionally disabled in this notebook build.')\n",
    "print('Reason: no cross-jurisdiction competitor_stringency data available.')\n",
    "print('Scope retained: RQ1/RQ2 border gaps and post-2019 differential trends within Texas districts.')\n",
    "\n",
    "# Optional: keep link/template artifacts for future extensions\n",
    "link_path = Path('../data/district_competitor_links.csv')\n",
    "template_path = Path('../data/competition_panel_template.csv')\n",
    "print('Reference files (for future data collection):')\n",
    "print(f'  {link_path.resolve()} (exists={link_path.exists()})')\n",
    "print(f'  {template_path.resolve()} (exists={template_path.exists()})')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "e47a5a70",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved available results files to output_borderlands/\n"
     ]
    }
   ],
   "source": [
    "# Save headline model outputs\n",
    "if 'rq1_results' in globals():\n",
    "    rq1_results.to_csv(OUTPUT_DIR / 'rq1_results.csv', index=False)\n",
    "if 'rq2_results' in globals():\n",
    "    rq2_results.to_csv(OUTPUT_DIR / 'rq2_results.csv', index=False)\n",
    "\n",
    "print('Saved available results files to output_borderlands/')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "524831c9",
   "metadata": {},
   "source": [
    "## Notes For Next Iteration\n",
    "\n",
    "- Validate border district coding against your preferred substantive definition.\n",
    "- Decide final `tx_stringency` metric for competition models (inspection intensity vs timing-based measure).\n",
    "- If you add competitor panel data, rerun the RQ3/RQ4 section without changing other cells.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bde1bf7d",
   "metadata": {},
   "source": [
    "## Research Questions and Hypothesis Test Summary (Current Run)\n",
    "\n",
    "### What This Notebook Tests\n",
    "This notebook estimates (1) **Texas-internal borderlands enforcement gaps** and (2) whether the **2019 disclosure reform** produced **border-specific changes** in district-level enforcement outcomes. Woods reaction-function estimation remains out of scope because cross-jurisdiction competitor stringency data are unavailable.\n",
    "\n",
    "### Data/Build Validation Snapshot\n",
    "- Wells loaded: **1,010,432**.\n",
    "- 2015-2025 analysis sample: inspections **1,867,859**; violations **191,762**.\n",
    "- Border well counts within 50km (current geometry build):\n",
    "  - TX-MX: **40,339**\n",
    "  - TX-NM: **81,567**\n",
    "  - TX-OK: **19,643**\n",
    "  - TX-LA: **29,675**\n",
    "  - Any border exposure: **169,520** wells.\n",
    "- District-year panel: **143 observations**, **13 districts**, **2015-2025**.\n",
    "\n",
    "### RQ1: Border Gaps (Levels)\n",
    "Do border districts differ systematically from non-border districts?\n",
    "\n",
    "- Descriptives (border vs non-border):\n",
    "  - Inspection intensity: **1.329 vs 1.515**\n",
    "  - Violations per inspection: **0.130 vs 0.098**\n",
    "  - Mean days to enforcement: **145.2 vs 122.8**\n",
    "  - Mean resolution rate: **0.543 vs 0.596**\n",
    "- Regression estimates:\n",
    "  - `inspection_intensity ~ border_district + ...`: coef **-0.1755**, p=**0.0999**\n",
    "  - `violations_per_inspection ~ border_district + ...`: coef **+0.0434**, p=**0.0949**\n",
    "\n",
    "**Interpretation:** Border districts show a weaker enforcement profile descriptively, with marginal model evidence for lower inspection intensity and higher violations per inspection.\n",
    "\n",
    "### RQ2: Disclosure × Border Heterogeneity\n",
    "Did post-2019 changes differ in border districts?\n",
    "\n",
    "- `inspection_intensity`: `post_2019:border` **-0.1191** (p=**0.0753**), `post_trend:border` **-0.0052** (p=**0.8181**)\n",
    "- `violations_per_inspection`: `post_2019:border` **+0.0040** (p=**0.8881**), `post_trend:border` **-0.0012** (p=**0.8350**)\n",
    "- `avg_days_to_enforcement`: `post_2019:border` **-74.5893** (p=**0.0156**), `post_trend:border` **-1.1587** (p=**0.9252**)\n",
    "- `resolution_rate`: `post_2019:border` **+0.0404** (p=**0.4520**), `post_trend:border` **-0.0186** (p=**0.3404**)\n",
    "\n",
    "**Interpretation:** The strongest border-specific post-2019 signal is in **timing** (faster enforcement), while border-specific changes in intensity, detection, and resolution are null. Border-specific trend differentials are null in the main binary-border model.\n",
    "\n",
    "### Robustness: Which Border / How Much Exposure / Which Cutoff\n",
    "\n",
    "- **Border-type split models (timing focus):**\n",
    "  - Level-model patterns suggest heterogeneity (e.g., `has_tx_ok` positive timing coefficient p=0.021; `has_tx_nm` negative/marginal p=0.061).\n",
    "  - FE interaction terms for `post_2019:has_tx_*` are null.\n",
    "  - Some `post_trend:has_tx_*` terms appear significant/marginal (`TX-NM` positive p=0.0189; `TX-OK` negative p=0.0794), indicating mixed border-type trend heterogeneity.\n",
    "- **Continuous exposure models:**\n",
    "  - Main continuous `post_2019` timing interaction is null (p=0.4449).\n",
    "  - One trend interaction appears significant for resolution (`post_trend:share_border_exposed_insp`, p=0.0423).\n",
    "- **Cutoff sensitivity (25/75/100km, timing):**\n",
    "  - `post_2019:share_border_*km` timing interactions are all null (p=0.7010, 0.5116, 0.9474).\n",
    "\n",
    "**Robustness interpretation:** The core binary-border timing level-shift result remains the clearest signal. Disaggregated and continuous/cutoff specifications show heterogeneity but no single, stable alternative pattern that dominates the baseline finding.\n",
    "\n",
    "### Money Plot\n",
    "- Figure: `output_borderlands/money_plot_timing_border_prepost2019.png`\n",
    "- Companion tables:\n",
    "  - `output_borderlands/money_border_type_timing_results.csv`\n",
    "  - `output_borderlands/money_cutoff_timing_results.csv`\n",
    "  - `output_borderlands/money_plot_timing_ci_by_year.csv`\n",
    "- **Aggregation note for manuscript caption/methods:**\n",
    "  - Year-group means and CIs in the money plot are computed from **district-year timing values** (each district-year is one observation), with equal weighting across districts within group-year.\n",
    "  - They are **not** pooled directly from well-level observations.\n",
    "\n",
    "### Hypothesis Status (Current Design)\n",
    "- **H1 (Border inspection gap):** **Partial support** (descriptives + marginal border coefficient).\n",
    "- **H2 (Border pipeline disadvantage):** **Supported descriptively**, mixed regression support by outcome.\n",
    "- **H3 (Disclosure heterogeneity in levels; `post_2019:border`):** **Supported for timing only**; not supported for other outcomes.\n",
    "- **H4 (Disclosure heterogeneity in trends; `post_trend:border`):** **Not supported** in the baseline binary-border model; robustness splits show mixed border-type trend heterogeneity.\n",
    "\n",
    "### Writing Direction for the Manuscript\n",
    "\n",
    "**Core claim hierarchy:**\n",
    "1. Border districts exhibit a persistent, descriptively weaker enforcement profile.\n",
    "2. The 2019 disclosure shift is most visible in enforcement timing (faster post-2019 in border districts).\n",
    "3. Broader “capacity expansion” signals (inspection intensity, resolution) are weak or null in differential models.\n",
    "4. Robustness checks complicate mechanism but do not overturn the central timing result.\n",
    "\n",
    "**Suggested paper structure:**\n",
    "1. Introduction: policy context, borderlands governance stakes, and transparency-reform relevance.\n",
    "2. Theory: borderlands capacity asymmetry + transparency/administrative throughput mechanism.\n",
    "3. Data and measures: district-year panel, border definitions, timing/compliance outcomes.\n",
    "4. Methods:\n",
    "   - RQ1 levels models.\n",
    "   - RQ2 district FE + year FE interaction models.\n",
    "   - Robustness (border type, continuous exposure, cutoff sensitivity).\n",
    "   - Explicit note that reaction-function estimation is deferred absent competitor stringency.\n",
    "5. Results:\n",
    "   - Main table: RQ1 and RQ2 coefficients.\n",
    "   - Money plot (timing border vs non-border with CI bands).\n",
    "   - Robustness table(s): border-type and cutoff/continuous checks.\n",
    "6. Discussion:\n",
    "   - “Faster pipeline, not wider pipeline” interpretation.\n",
    "   - What this implies for transparency reforms in uneven-capacity settings.\n",
    "7. Limitations and next steps:\n",
    "   - No competitor stringency data.\n",
    "   - Small number of districts and inference caution.\n",
    "   - Future expansion to true cross-jurisdiction competition panel.\n",
    "\n",
    "**Language to avoid in draft claims:**\n",
    "- Avoid strong causal claims on interstate competition.\n",
    "- Avoid implying broad compliance improvement unless supported by outcome-specific coefficients.\n",
    "\n",
    "### Theory Positioning\n",
    "Borderlands governance theory is consistent with persistent gap patterns and uneven enforcement capacity. The disclosure-era effect appears concentrated in **processing speed** rather than expanded inspection coverage or broad compliance improvements. Woods-style competition remains an interpretive frame, but not an estimated mechanism in this notebook.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c63fc72",
   "metadata": {},
   "source": [
    "## Research Questions and Hypothesis Test Summary (Current Run)\n",
    "\n",
    "### What This Notebook Tests\n",
    "\n",
    "This notebook estimates two Texas-internal questions using district-year outcomes from 2015–2025:\n",
    "(1) whether **border-exposed districts** exhibit systematic enforcement gaps relative to non-border districts, and\n",
    "(2) whether the **2019 disclosure reform** is associated with **border-specific changes** in enforcement outcomes.\n",
    "Woods-style reaction-function tests are intentionally out of scope because comparable competitor stringency data are not available.\n",
    "\n",
    "### Data and Build Validation\n",
    "\n",
    "* Wells loaded: **1,010,432**.\n",
    "* 2015–2025 analysis sample: inspections **1,867,859**; violations **191,762**.\n",
    "* Border well counts within 50 km (current geometry build):\n",
    "\n",
    "  * TX–MX: **40,339**\n",
    "  * TX–NM: **81,567**\n",
    "  * TX–OK: **19,643**\n",
    "  * TX–LA: **29,675**\n",
    "  * Any border exposure: **169,520** wells\n",
    "* District-year panel: **143 observations**, **13 districts**, **2015–2025**.\n",
    "\n",
    "---\n",
    "\n",
    "## RQ1: Border Gaps\n",
    "\n",
    "**Do border districts differ systematically from non-border districts?**\n",
    "\n",
    "### Descriptive gaps (border vs non-border)\n",
    "\n",
    "* Inspection intensity: **1.329 vs 1.515**\n",
    "* Violations per inspection: **0.130 vs 0.098**\n",
    "* Mean days to enforcement: **145.2 vs 122.8**\n",
    "* Mean resolution rate: **0.543 vs 0.596**\n",
    "\n",
    "### Regression evidence (levels models)\n",
    "\n",
    "* `inspection_intensity ~ border_district + ...`: **−0.1755** (p=**0.0999**)\n",
    "* `violations_per_inspection ~ border_district + ...`: **+0.0434** (p=**0.0949**)\n",
    "\n",
    "**Interpretation:** Border districts show a consistently weaker enforcement profile in descriptives, with **marginal** model evidence of lower inspection intensity and higher violations per inspection. The descriptive pattern is broad; the most reliable regression signal at this stage is concentrated in **coverage (inspections)** and **detection conditional on inspection (violations/inspection)**.\n",
    "\n",
    "---\n",
    "\n",
    "## RQ2: Disclosure × Border Heterogeneity\n",
    "\n",
    "**Did post-2019 changes differ in border districts?**\n",
    "\n",
    "* `inspection_intensity`: `post_2019:border` **−0.1191** (p=**0.0753**), `post_trend:border` **−0.0052** (p=**0.8181**)\n",
    "* `violations_per_inspection`: `post_2019:border` **+0.0040** (p=**0.8881**), `post_trend:border` **−0.0012** (p=**0.8350**)\n",
    "* `avg_days_to_enforcement`: `post_2019:border` **−74.5893** (p=**0.0156**), `post_trend:border` **−1.1587** (p=**0.9252**)\n",
    "* `resolution_rate`: `post_2019:border` **+0.0404** (p=**0.4520**), `post_trend:border` **−0.0186** (p=**0.3404**)\n",
    "\n",
    "**Interpretation:** The border-specific post-2019 signal is concentrated in **timing**: border districts become substantially faster after 2019 (a clear **level shift**). Border-specific changes in inspection intensity, detection intensity, and resolution are null, and border-specific trend differentials are null in the baseline binary-border model. In short: **the border pipeline got faster, not wider.**\n",
    "\n",
    "---\n",
    "\n",
    "## Robustness and Heterogeneity Checks\n",
    "\n",
    "These checks evaluate whether the baseline binary-border timing result is an artifact of a particular border definition or distance cutoff.\n",
    "\n",
    "### Border-type splits (timing focus)\n",
    "\n",
    "* Level-model patterns suggest heterogeneity (e.g., `has_tx_ok` positive timing coefficient p=0.021; `has_tx_nm` negative/marginal p=0.061).\n",
    "* FE interactions for `post_2019:has_tx_*` are null.\n",
    "* Some `post_trend:has_tx_*` terms are significant/marginal (TX–NM positive p=0.0189; TX–OK negative p=0.0794), indicating **uneven border-type trend heterogeneity** in split specifications.\n",
    "\n",
    "### Continuous exposure models\n",
    "\n",
    "* Main continuous `post_2019` timing interaction is null (p=0.4449).\n",
    "* One trend interaction is significant for resolution (`post_trend:share_border_exposed_insp`, p=0.0423).\n",
    "\n",
    "### Cutoff sensitivity (25/75/100 km, timing)\n",
    "\n",
    "* Timing interactions for `post_2019:share_border_*km` are all null (p=0.7010, 0.5116, 0.9474).\n",
    "\n",
    "**Robustness interpretation:** The baseline binary-border timing level-shift remains the clearest and most stable result. Alternative specifications show heterogeneity—especially by border type—but do not yield a single, consistent alternative pattern that replaces the main finding. The absence of a monotonic “dose-response” in continuous/cutoff models suggests the timing shift behaves more like a **district-level regime change** than a smooth function of border exposure intensity.\n",
    "\n",
    "---\n",
    "\n",
    "## Money Plot and Aggregation Note\n",
    "\n",
    "* Figure: `output_borderlands/money_plot_timing_border_prepost2019.png`\n",
    "* Companion tables:\n",
    "\n",
    "  * `output_borderlands/money_border_type_timing_results.csv`\n",
    "  * `output_borderlands/money_cutoff_timing_results.csv`\n",
    "  * `output_borderlands/money_plot_timing_ci_by_year.csv`\n",
    "\n",
    "**Aggregation note (for caption/methods):** Year-group means and 95% CIs are computed from **district-year timing values** (each district-year is one observation), with **equal weighting across districts within each group-year**; results are not pooled from well-level observations.\n",
    "\n",
    "---\n",
    "\n",
    "## Hypothesis Status (Current Design)\n",
    "\n",
    "* **H1 (Border inspection gap):** **Partial support** (descriptives + marginal border coefficient).\n",
    "* **H2 (Border pipeline disadvantage):** **Supported descriptively**, mixed regression support by outcome.\n",
    "* **H3 (Disclosure heterogeneity in levels; `post_2019:border`):** **Supported for timing only**; null elsewhere.\n",
    "* **H4 (Disclosure heterogeneity in trends; `post_trend:border`):** **Not supported** in baseline; split models show mixed border-type trend heterogeneity.\n",
    "\n",
    "---\n",
    "\n",
    "## Writing Direction for the Manuscript\n",
    "\n",
    "### Core claim hierarchy (keep this ordering)\n",
    "\n",
    "1. Border districts exhibit a persistent, descriptively weaker enforcement profile (coverage, detection conditional on inspection, timing, and follow-through).\n",
    "2. The 2019 disclosure reform is associated with a **border-specific improvement in enforcement timing** (a sizeable level shift).\n",
    "3. Differential evidence for broader “capacity expansion” (inspection intensity, resolution) is weak or null.\n",
    "4. Robustness checks complicate mechanisms and reveal heterogeneity, but **do not overturn** the central timing result.\n",
    "\n",
    "### Suggested structure (mapped to what you’ve already built)\n",
    "\n",
    "1. Introduction: Texas enforcement context + why borderlands governance matters + why transparency is a plausible shock.\n",
    "2. Theory: borderlands capacity asymmetry + transparency-throughput mechanism (focus on timing).\n",
    "3. Data and measures: district-year panel, border definitions, outcomes.\n",
    "4. Methods:\n",
    "\n",
    "   * RQ1 levels models\n",
    "   * RQ2 district FE + year FE interaction models\n",
    "   * robustness suite (border type, continuous exposure, cutoffs)\n",
    "   * explicit note: reaction-function tests deferred\n",
    "5. Results: main table + timing money plot with CI ribbons + compact robustness table(s).\n",
    "6. Discussion: “faster pipeline, not wider pipeline” + what that implies for transparency reforms in uneven-capacity settings.\n",
    "7. Limitations/next steps: competitor data gap; small-N districts; future cross-jurisdiction extension.\n",
    "\n",
    "### Guardrails (language discipline)\n",
    "\n",
    "* Don’t imply interstate competition is identified here.\n",
    "* Don’t generalize to “compliance improved” unless the outcome-specific coefficients support it.\n",
    "* Be precise: timing improves; coverage and closure do not reliably move in the border interaction models.\n",
    "\n",
    "---\n",
    "\n",
    "## Theory Positioning (tight version)\n",
    "\n",
    "Borderlands governance theory fits the baseline pattern: persistent enforcement gaps consistent with uneven state capacity and asymmetric governance conditions. The disclosure-era result is concentrated in **processing speed**, not inspection coverage or broad improvements in closure outcomes, suggesting transparency affects **throughput** more than capacity. Woods-style competition remains a plausible interpretive frame, but it is not estimated in this notebook.\n",
    "\n",
    "---\n",
    "\n",
    "**Transparency can accelerate administrative processing where enforcement is structurally weaker, but speed alone is not the same thing as capacity.**\n",
    "\n",
    "---\n",
    "\n",
    "*This analysis started as an effort to test Neal D. Woods’s race-to-the-bottom claim that border jurisdictions relax enforcement when they risk looking stricter than nearby competitors. The first step is to establish whether Texas borderlands districts exhibit distinct enforcement profiles and whether the 2019 disclosure reform altered those profiles in a border-specific way. The results here identify persistent borderlands gaps and a sharp post-2019 border-specific improvement in enforcement timing (“faster pipeline, not wider pipeline”). A direct Woods reaction-function test—linking Texas enforcement to competitor stringency—remains the next phase once cross-jurisdiction enforcement measures are constructed.*\n",
    "\n",
    "---\n",
    "\n",
    "## Companion Tables (CSV Outputs)\n",
    "\n",
    "### Table A. RQ2 FE Interaction Results\n",
    "\n",
    "| Outcome | post_2019 × border (coef) | p-value | post_trend × border (coef) | p-value | N |\n",
    "|---|---:|---:|---:|---:|---:|\n",
    "| inspection_intensity | -0.119 | 0.075 | -0.005 | 0.818 | 143 |\n",
    "| violations_per_inspection | 0.004 | 0.888 | -0.001 | 0.835 | 143 |\n",
    "| avg_days_to_enforcement | -74.589 | 0.016 | -1.159 | 0.925 | 143 |\n",
    "| resolution_rate | 0.040 | 0.452 | -0.019 | 0.340 | 143 |\n",
    "\n",
    "### Table B. Competition Asymmetry Results\n",
    "\n",
    "| Term | Coefficient | p-value |\n",
    "|---|---:|---:|\n",
    "| gap_pos | 0.584 | <0.001 |\n",
    "| gap_neg | 0.499 | 0.003 |\n",
    "\n",
    "### Table C. Cutoff Sensitivity Results\n",
    "\n",
    "| Cutoff (km) | Model | Outcome | Term | Coefficient | p-value | N |\n",
    "|---:|---|---|---|---:|---:|---:|\n",
    "| 25 | RQ1_levels | inspection_intensity | share_border_25km | 0.539 | 0.423 | 143 |\n",
    "| 25 | RQ2_FE | avg_days_to_enforcement | post_2019:share_border_25km | -101.928 | 0.701 | 143 |\n",
    "| 25 | RQ2_FE | avg_days_to_enforcement | post_trend:share_border_25km | 30.609 | 0.770 | 143 |\n",
    "| 75 | RQ1_levels | inspection_intensity | share_border_75km | 0.079 | 0.696 | 143 |\n",
    "| 75 | RQ2_FE | avg_days_to_enforcement | post_2019:share_border_75km | -75.659 | 0.512 | 143 |\n",
    "| 75 | RQ2_FE | avg_days_to_enforcement | post_trend:share_border_75km | -0.009 | 1.000 | 143 |\n",
    "| 100 | RQ1_levels | inspection_intensity | share_border_100km | 0.002 | 0.987 | 143 |\n",
    "| 100 | RQ2_FE | avg_days_to_enforcement | post_2019:share_border_100km | -4.480 | 0.947 | 143 |\n",
    "| 100 | RQ2_FE | avg_days_to_enforcement | post_trend:share_border_100km | -12.659 | 0.515 | 143 |\n",
    "\n",
    "### Table D. Border-Type Timing Interactions (Money Plot Companion)\n",
    "\n",
    "| Term | Coefficient | p-value |\n",
    "|---|---:|---:|\n",
    "| post_2019:has_tx_mex | 4.090 | 0.906 |\n",
    "| post_2019:has_tx_nm | -18.744 | 0.601 |\n",
    "| post_2019:has_tx_ok | -14.245 | 0.813 |\n",
    "| post_2019:has_tx_la | -43.660 | 0.641 |\n",
    "| post_trend:has_tx_mex | -0.015 | 0.999 |\n",
    "| post_trend:has_tx_nm | 22.907 | 0.019 |\n",
    "| post_trend:has_tx_ok | -16.719 | 0.079 |\n",
    "| post_trend:has_tx_la | 0.642 | 0.955 |\n",
    "\n",
    "### Table E. District Border-Type Profile\n",
    "\n",
    "| District | Wells | Dominant border type | TX-MX share | TX-NM share | TX-OK share | TX-LA share |\n",
    "|---:|---:|---|---:|---:|---:|---:|\n",
    "| 01 | 31,898 | TX-MX | 23.1% | 0.0% | 0.0% | 0.0% |\n",
    "| 02 | 17,099 | NONE | 0.0% | 0.0% | 0.0% | 0.0% |\n",
    "| 03 | 16,700 | TX-LA | 0.0% | 0.0% | 0.0% | 11.7% |\n",
    "| 04 | 20,973 | TX-MX | 63.8% | 0.0% | 0.0% | 0.0% |\n",
    "| 05 | 9,938 | TX-OK | 0.0% | 0.0% | 0.2% | 0.0% |\n",
    "| 06 | 24,422 | TX-LA | 0.0% | 0.0% | 2.9% | 52.4% |\n",
    "| 08 | 105,931 | TX-NM | 0.0% | 19.1% | 0.0% | 0.0% |\n",
    "| 09 | 46,485 | NONE | 0.0% | 0.0% | 0.0% | 0.0% |\n",
    "| 10 | 29,621 | TX-OK | 0.0% | 0.1% | 30.2% | 0.0% |\n",
    "| 6E | 6,235 | NONE | 0.0% | 0.0% | 0.0% | 0.0% |\n",
    "| 7B | 21,230 | NONE | 0.0% | 0.0% | 0.0% | 0.0% |\n",
    "| 7C | 43,061 | NONE | 0.0% | 0.0% | 0.0% | 0.0% |\n",
    "| 8A | 42,005 | TX-NM | 0.0% | 41.8% | 0.0% | 0.0% |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33d763c3",
   "metadata": {},
   "source": [
    "## Methods (Policy Studies Journal Draft)\n",
    "\n",
    "### Research Design and Unit of Analysis\n",
    "This study uses a district-year panel design to estimate enforcement differences between border and non-border Texas Railroad Commission (RRC) districts and to test whether those differences changed after the January 2019 disclosure reform. The primary analytical unit is the **district-year** observation (13 districts, 2015-2025; balanced years where data exist).\n",
    "\n",
    "The design has two estimands:\n",
    "1. Cross-district level gaps associated with border exposure.\n",
    "2. Border-specific post-2019 differential shifts (level and trend).\n",
    "\n",
    "Woods-style interstate reaction-function estimation is not implemented in this version because comparable competitor-jurisdiction stringency series are unavailable.\n",
    "\n",
    "### Data Sources and Construction\n",
    "The notebook draws from PostgreSQL warehouse tables (`well_shape_tract`, `inspections`, `violations`) through `WellAnalyzer`, then restricts to 2015-2025.\n",
    "\n",
    "Core preprocessing steps:\n",
    "1. Standardize well identifier (`api_norm`) and parse dates.\n",
    "2. Recompute border proximity from geometry, including TX-MX and state-border segments (TX-NM, TX-OK, TX-LA).\n",
    "3. Construct well-level border flags at distance thresholds.\n",
    "4. Aggregate to district-year outcomes.\n",
    "\n",
    "### Outcome Measures\n",
    "For district $d$ in year $t$:\n",
    "\n",
    "- **Inspection intensity**\n",
    "$$\n",
    "InspectionIntensity_{dt} = \\frac{Inspections_{dt}}{UniqueWells_{dt}}\n",
    "$$\n",
    "\n",
    "- **Violations per inspection**\n",
    "$$\n",
    "ViolPerInsp_{dt} = \\frac{Violations_{dt}}{Inspections_{dt}}\n",
    "$$\n",
    "\n",
    "- **Average days to enforcement**\n",
    "$$\n",
    "DaysToEnf_{dt} = \\frac{1}{N_{dt}}\\sum_{i=1}^{N_{dt}}\\left(EnforcementDate_i - ViolationDiscoveryDate_i\\right)\n",
    "$$\n",
    "where implausible negative/invalid intervals are excluded by preprocessing rules.\n",
    "\n",
    "- **Resolution rate proxy**\n",
    "$$\n",
    "ResolutionRate_{dt} = \\frac{CompliantOnReinspection_{dt}}{Violations_{dt}}\n",
    "$$\n",
    "\n",
    "### Border Exposure Definitions\n",
    "Primary binary treatment is district-level border status:\n",
    "$$\n",
    "Border_d = 1 \\text{ if district is classified as border-adjacent; } 0 \\text{ otherwise.}\n",
    "$$\n",
    "\n",
    "Additional exposure constructions are used for robustness:\n",
    "- Border-type indicators: $\\text{TX-MX}$, $\\text{TX-NM}$, $\\text{TX-OK}$, and $\\text{TX-LA}$.\n",
    "- Continuous exposure share:\n",
    "$$\n",
    "ShareBorder_{dt} = \\frac{BorderExposedInspections_{dt}}{Inspections_{dt}}\n",
    "$$\n",
    "- Distance-threshold sensitivity at 25km, 75km, and 100km.\n",
    "\n",
    "### Baseline Models\n",
    "\n",
    "### RQ1: Border Gaps (Levels)\n",
    "For each outcome $Y_{dt}$, the levels model is:\n",
    "$$\n",
    "Y_{dt} = \\alpha + \\beta_1 Border_d + \\beta_2 \\log(UniqueWells_{dt}) + \\gamma_t + \\varepsilon_{dt}\n",
    "$$\n",
    "where $\\gamma_t$ are year fixed effects.\n",
    "\n",
    "Interpretation: $\\beta_1$ captures average border/non-border differences conditional on annual shocks and scale.\n",
    "\n",
    "### RQ2: Disclosure Heterogeneity (Post-2019)\n",
    "The fixed-effects interaction specification is:\n",
    "$$\n",
    "Y_{dt} = \\alpha_d + \\gamma_t + \\theta_1(Post2019_t \\times Border_d) + \\theta_2(PostTrend_t \\times Border_d) + \\varepsilon_{dt}\n",
    "$$\n",
    "with:\n",
    "$$\n",
    "Post2019_t = \\mathbb{1}[t \\ge 2019], \\quad PostTrend_t = \\max(0, t-2019)\n",
    "$$\n",
    "$\\alpha_d$ are district fixed effects and $\\gamma_t$ year fixed effects.\n",
    "\n",
    "Interpretation:\n",
    "- $\\theta_1$: immediate border-specific post-2019 level differential.\n",
    "- $\\theta_2$: border-specific post-2019 slope differential.\n",
    "\n",
    "### Robustness Specifications\n",
    "\n",
    "### Border-Type Disaggregation\n",
    "Binary border indicator is replaced with border-type indicators and interactions:\n",
    "$$\n",
    "Y_{dt} = \\alpha_d + \\gamma_t + \\sum_k \\lambda_k(Post2019_t \\times Type_{kd}) + \\sum_k \\phi_k(PostTrend_t \\times Type_{kd}) + \\varepsilon_{dt}\n",
    "$$\n",
    "where $k \\in \\{\\text{TX-MX},\\text{TX-NM},\\text{TX-OK},\\text{TX-LA}\\}$.\n",
    "\n",
    "### Continuous Exposure Models\n",
    "Border binary is replaced by share exposure:\n",
    "$$\n",
    "Y_{dt} = \\alpha_d + \\gamma_t + \\eta_1(Post2019_t \\times ShareBorder_{dt}) + \\eta_2(PostTrend_t \\times ShareBorder_{dt}) + \\varepsilon_{dt}\n",
    "$$\n",
    "\n",
    "### Cutoff Sensitivity\n",
    "The exposure share is recomputed with alternative distance thresholds $c\\in\\{25,75,100\\}$:\n",
    "$$\n",
    "ShareBorder^{(c)}_{dt}\n",
    "$$\n",
    "and substituted into the same FE interaction framework.\n",
    "\n",
    "### Inference and Uncertainty\n",
    "All regression models are estimated by OLS with district-clustered standard errors:\n",
    "$$\n",
    "\\widehat{Var}(\\hat\\beta) = \\text{Clustered by district}\n",
    "$$\n",
    "Given the small number of clusters (13 districts), p-values are interpreted cautiously and emphasis is placed on direction, magnitude, and cross-specification consistency.\n",
    "\n",
    "### Figure Construction (Money Plot)\n",
    "The main timing figure plots annual means for border and non-border groups with 95% confidence intervals computed from district-year values:\n",
    "$$\n",
    "\\bar{Y}_{gt} = \\frac{1}{n_{gt}}\\sum_{d \\in g} Y_{dt}\n",
    "$$\n",
    "$$\n",
    "SE(\\bar{Y}_{gt}) = \\frac{s_{gt}}{\\sqrt{n_{gt}}}, \\quad CI_{95\\%} = \\bar{Y}_{gt} \\pm 1.96\\,SE(\\bar{Y}_{gt})\n",
    "$$\n",
    "where $g$ indexes border group. These are **district-year means with equal district weighting**, not pooled well-level means.\n",
    "\n",
    "### Identification Scope and Limitations\n",
    "The design supports within-Texas comparative inference about border gaps and post-policy heterogeneity. It does not identify strategic interstate reaction effects without external competitor-stringency series. Results should therefore be interpreted as transparency-era heterogeneity in Texas enforcement administration rather than a full reaction-function test.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
