assignment4.py

# ============================================================
# BANA 290 — Assignment 4
# Smart Campus Student App Incubator Archive
# Author : Todd Denaro  |  denarot@uci.edu
# ============================================================

import re
import warnings
import numpy as np
import pandas as pd
import requests
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
import statsmodels.formula.api as smf
import statsmodels.api as sm
from bs4 import BeautifulSoup
from scipy import stats
from scipy.stats import gaussian_kde

warnings.filterwarnings("ignore")

BASE_URL = "https://bana290-assignment4.netlify.app"

# ============================================================
# PART 1 — SCRAPE
# ============================================================

# Fetch the index page and collect every archive link stored in
# anchor tags with class "brief-link"; prepend BASE_URL so each
# entry is a fully-qualified URL ready for requests.get()
index_resp = requests.get(BASE_URL)
index_resp.raise_for_status()
index_soup = BeautifulSoup(index_resp.text, "html.parser")

archive_links = [
    BASE_URL + a["href"]
    for a in index_soup.select("a.brief-link")
]

print("Archive pages found:")
for link in archive_links:
    print(" ", link)


# Parse a single archive page: locate the one <table class="archive-table">,
# treat the first <tr> as the header row, and for every data row extract
# the team code from the nested <strong> tag and the app name from
# <span class="app-label"> inside the "team-cell" <td>; return a DataFrame
def scrape_archive_table(url: str) -> pd.DataFrame:
    resp = requests.get(url)
    resp.raise_for_status()
    soup = BeautifulSoup(resp.text, "html.parser")

    table = soup.select_one("table.archive-table")
    rows  = table.find_all("tr")

    headers = [td.get_text(strip=True) for td in rows[0].find_all("td")]

    records = []
    for row in rows[1:]:
        cells  = row.find_all("td")
        record = {}
        for i, (header, cell) in enumerate(zip(headers, cells)):
            if i == 0:
                code     = cell.find("strong").get_text(strip=True)
                name_tag = cell.find("span", class_="app-label")
                name     = name_tag.get_text(strip=True) if name_tag else ""
                record["TEAM_REF"]  = code
                record["TEAM_NAME"] = name
            else:
                record[header] = cell.get_text(strip=True)
        records.append(record)

    df   = pd.DataFrame(records)
    cols = (["TEAM_REF", "TEAM_NAME"]
            + [c for c in df.columns if c not in ("TEAM_REF", "TEAM_NAME")])
    return df[cols]


# Loop over every collected archive URL, call scrape_archive_table to get
# a DataFrame per page, key each result by its URL slug, then inner-merge
# all three tables on TEAM_REF — dropping the duplicate TEAM_NAME columns
# from the second and third tables — to produce one master dataset; save to CSV
raw_dfs = {}
for url in archive_links:
    df_raw = scrape_archive_table(url)
    slug   = url.rstrip("/").split("/")[-1]
    raw_dfs[slug] = df_raw
    print(f"\n--- {slug} ({len(df_raw)} rows) ---")
    print(df_raw.head(3).to_string(index=False))

df_infra   = raw_dfs["fiber-access-bulletin"]
df_metrics = raw_dfs["builder-metrics-ledger"]
df_funding = raw_dfs["anteater-fund-panel"]

drop_name = lambda d: [c for c in d.columns if c != "TEAM_NAME"]

master = (
    df_infra
    .merge(df_metrics[drop_name(df_metrics)], on="TEAM_REF", how="inner")
    .merge(df_funding[drop_name(df_funding)],  on="TEAM_REF", how="inner")
)

master.to_csv("campus_incubator_master.csv", index=False)
print(f"\nMaster dataset: {master.shape[0]} rows x {master.shape[1]} cols")
print("Saved: campus_incubator_master.csv")

# ============================================================
# PART 2 — CLEAN
# ============================================================

df = pd.read_csv("campus_incubator_master.csv")

# Parse DISTANCE_TO_NODE: strip thousands commas, extract the first decimal
# number with a regex, then multiply by 1000 when the string contains "km"
# anywhere so that every value is stored in metres as float64
def parse_distance_m(s: str) -> float:
    s     = s.replace(",", "")
    match = re.search(r"(\d+\.?\d*)", s)
    if not match:
        return np.nan
    value = float(match.group(1))
    if re.search(r"km", s, re.IGNORECASE):
        value = round(value * 1000, 1)
    return value

df["DISTANCE_TO_NODE"] = df["DISTANCE_TO_NODE"].apply(parse_distance_m)

# Extract the first decimal number from ELIGIBILITY_SCORE, AI_INTENSITY,
# and INNOVATION_SCORE to strip every label variant ("Pitch rating = X",
# "panel avg X", "X / 100", "X score", "~X model hrs", etc.) and cast
# each column to float64 so regression models can read them
def parse_first_number(s: str) -> float:
    match = re.search(r"(\d+\.?\d*)", str(s))
    return float(match.group(1)) if match else np.nan

df["ELIGIBILITY_SCORE"] = df["ELIGIBILITY_SCORE"].apply(parse_first_number)
df["AI_INTENSITY"]      = df["AI_INTENSITY"].apply(parse_first_number)
df["INNOVATION_SCORE"]  = df["INNOVATION_SCORE"].apply(parse_first_number)

# Apply IQR capping to AI_INTENSITY and INNOVATION_SCORE: compute Q1, Q3,
# IQR; set lower fence = Q1 - 1.5*IQR and upper fence = Q3 + 1.5*IQR;
# clip both series to those fences so no single extreme team skews the
# causal estimates
def iqr_clip(series: pd.Series) -> pd.Series:
    q1    = series.quantile(0.25)
    q3    = series.quantile(0.75)
    iqr   = q3 - q1
    lower = q1 - 1.5 * iqr
    upper = q3 + 1.5 * iqr
    clipped   = series.clip(lower=lower, upper=upper)
    n_changed = (series != clipped).sum()
    print(f"  {series.name}: fence=[{lower:.2f}, {upper:.2f}]  "
          f"values capped={n_changed}")
    return clipped

print("\nIQR capping:")
df["AI_INTENSITY"]     = iqr_clip(df["AI_INTENSITY"])
df["INNOVATION_SCORE"] = iqr_clip(df["INNOVATION_SCORE"])

df.to_csv("campus_incubator_clean.csv", index=False)
print("\nCleaned dtypes:")
print(df[["DISTANCE_TO_NODE","AI_INTENSITY",
          "INNOVATION_SCORE","ELIGIBILITY_SCORE"]].dtypes.to_string())
print("\nDescriptive statistics:")
print(df[["DISTANCE_TO_NODE","AI_INTENSITY",
          "INNOVATION_SCORE","ELIGIBILITY_SCORE"]].describe().round(2).to_string())
print("\nSaved: campus_incubator_clean.csv")

# ============================================================
# PART 3 — ANALYZE
# ============================================================

df = pd.read_csv("campus_incubator_clean.csv")

# Run a naive OLS of INNOVATION_SCORE on AI_INTENSITY as a baseline
# estimate; this is likely upward-biased because unobserved team quality
# drives both variables simultaneously
ols_model = smf.ols("INNOVATION_SCORE ~ AI_INTENSITY", data=df).fit()
print("\n[OLS — naive estimate]")
print(ols_model.summary2().tables[1][["Coef.", "Std.Err.", "t", "P>|t|"]])
ols_coef = ols_model.params["AI_INTENSITY"]
ols_se   = ols_model.bse["AI_INTENSITY"]

# First stage of 2SLS: regress AI_INTENSITY on DISTANCE_TO_NODE to verify
# instrument relevance; print the F-statistic and flag whether it clears
# the conventional threshold of 10 for a strong instrument
first_stage = smf.ols("AI_INTENSITY ~ DISTANCE_TO_NODE", data=df).fit()
f_stat = first_stage.fvalue
f_pval = first_stage.f_pvalue
r2_fs  = first_stage.rsquared
print("\n[First Stage: AI_INTENSITY ~ DISTANCE_TO_NODE]")
print(first_stage.summary2().tables[1][["Coef.", "Std.Err.", "t", "P>|t|"]])
print(f"  F = {f_stat:.2f}  p = {f_pval:.4f}  R2 = {r2_fs:.4f}  "
      f"{'STRONG' if f_stat > 10 else 'WEAK'}")

df["AI_INTENSITY_hat"] = first_stage.fittedvalues

# Second stage: regress INNOVATION_SCORE on AI_INTENSITY_hat (the fitted
# values from the first stage); the slope coefficient is the 2SLS estimate.
# Correct the standard errors by recomputing residuals with the ORIGINAL
# AI_INTENSITY so the sandwich variance uses the structural equation
stage2_naive = smf.ols("INNOVATION_SCORE ~ AI_INTENSITY_hat", data=df).fit()
beta_0 = stage2_naive.params["Intercept"]
beta_1 = stage2_naive.params["AI_INTENSITY_hat"]

df["resid_structural"] = (df["INNOVATION_SCORE"]
                          - beta_0
                          - beta_1 * df["AI_INTENSITY"])

n   = len(df)
Xh  = sm.add_constant(df["AI_INTENSITY_hat"])
e   = df["resid_structural"].values
XhTXh_inv = np.linalg.inv(Xh.T @ Xh)
meat      = (Xh.T * e**2) @ Xh
V_2sls    = n / (n - 2) * XhTXh_inv @ meat @ XhTXh_inv
se_2sls   = np.sqrt(np.diag(V_2sls))
t_2sls    = beta_1 / se_2sls[1]
p_2sls    = 2 * stats.t.sf(abs(t_2sls), df=n - 2)

print("\n[2SLS — causal estimate]")
print(f"  AI_INTENSITY coef : {beta_1:.4f}  SE={se_2sls[1]:.4f}  "
      f"t={t_2sls:.2f}  p={p_2sls:.4f}")
print(f"\n[OLS vs 2SLS]\n"
      f"  OLS  {ols_coef:.4f} (SE={ols_se:.4f})\n"
      f"  2SLS {beta_1:.4f} (SE={se_2sls[1]:.4f})")

# Sharp RDD: centre ELIGIBILITY_SCORE at the 85-point cutoff, create a
# TREATED dummy equal to 1 for teams that received server credits, then
# fit local-linear regressions on both sides of the threshold; the
# coefficient on TREATED is the Local Average Treatment Effect (LATE)
CUTOFF = 85.0
df["SCORE_C"] = df["ELIGIBILITY_SCORE"] - CUTOFF
df["TREATED"] = (df["ELIGIBILITY_SCORE"] >= CUTOFF).astype(int)

rdd_model = smf.ols(
    "INNOVATION_SCORE ~ TREATED + SCORE_C + TREATED:SCORE_C", data=df
).fit()
rdd_jump = rdd_model.params["TREATED"]
rdd_se   = rdd_model.bse["TREATED"]
rdd_p    = rdd_model.pvalues["TREATED"]

bw     = 10
df_bw  = df[df["SCORE_C"].abs() <= bw]
rdd_bw = smf.ols(
    "INNOVATION_SCORE ~ TREATED + SCORE_C + TREATED:SCORE_C", data=df_bw
).fit()

print("\n[RDD — full sample]")
print(rdd_model.summary2().tables[1][["Coef.", "Std.Err.", "t", "P>|t|"]])
print(f"\n  LATE = {rdd_jump:.3f}  SE={rdd_se:.3f}  p={rdd_p:.4f}")

# Build a four-panel figure: (1) first-stage scatter with regression line
# and F-stat annotation, (2) OLS vs 2SLS fitted lines over the raw scatter,
# (3) RDD scatter with local-linear fits on each side and a labelled jump
# arrow at the cutoff, (4) histogram with KDE overlay for the manipulation
# density test; save the figure to analysis_results.png
sns.set_style("whitegrid")
BLUE   = "#2563EB"
ORANGE = "#F59E0B"
RED    = "#DC2626"
GREY   = "#6B7280"

fig = plt.figure(figsize=(16, 14))
gs  = gridspec.GridSpec(2, 2, figure=fig, hspace=0.40, wspace=0.35)

ax1 = fig.add_subplot(gs[0, 0])
ax1.scatter(df["DISTANCE_TO_NODE"], df["AI_INTENSITY"],
            color=BLUE, alpha=0.65, edgecolors="white", s=60)
x_l = np.linspace(df["DISTANCE_TO_NODE"].min(),
                  df["DISTANCE_TO_NODE"].max(), 200)
ax1.plot(x_l,
         first_stage.params["Intercept"]
         + first_stage.params["DISTANCE_TO_NODE"] * x_l,
         color=RED, linewidth=2)
ax1.set_xlabel("Distance to Fiber Node (metres)", fontsize=11)
ax1.set_ylabel("AI Intensity (compute hrs/wk)", fontsize=11)
ax1.set_title(f"First Stage: AI Intensity ~ Distance\n"
              f"F = {f_stat:.1f}  |  R² = {r2_fs:.3f}",
              fontsize=12, fontweight="bold")
ax1.annotate(f"β = {first_stage.params['DISTANCE_TO_NODE']:.4f}\nF = {f_stat:.1f}",
             xy=(0.05, 0.90), xycoords="axes fraction", fontsize=10,
             color=RED, bbox=dict(boxstyle="round,pad=0.3", fc="white",
                                  ec=RED, alpha=0.8))

ax2 = fig.add_subplot(gs[0, 1])
ax2.scatter(df["AI_INTENSITY"], df["INNOVATION_SCORE"],
            color=GREY, alpha=0.55, edgecolors="white", s=60, label="Observed")
ai_r = np.linspace(df["AI_INTENSITY"].min(), df["AI_INTENSITY"].max(), 200)
ax2.plot(ai_r, ols_model.params["Intercept"] + ols_coef * ai_r,
         color=ORANGE, linewidth=2, label=f"OLS  β={ols_coef:.3f}")
ax2.plot(ai_r, beta_0 + beta_1 * ai_r,
         color=BLUE,   linewidth=2, linestyle="--",
         label=f"2SLS β={beta_1:.3f}")
ax2.set_xlabel("AI Intensity (compute hrs/wk)", fontsize=11)
ax2.set_ylabel("Innovation Score", fontsize=11)
ax2.set_title("OLS vs 2SLS: Innovation ~ AI Intensity",
              fontsize=12, fontweight="bold")
ax2.legend(fontsize=9, framealpha=0.9)

ax3 = fig.add_subplot(gs[1, 0])
lo = df["ELIGIBILITY_SCORE"] < CUTOFF
hi = df["ELIGIBILITY_SCORE"] >= CUTOFF
ax3.scatter(df.loc[lo, "ELIGIBILITY_SCORE"], df.loc[lo, "INNOVATION_SCORE"],
            color=ORANGE, alpha=0.75, edgecolors="white", s=60,
            label="Below cutoff (waitlist)")
ax3.scatter(df.loc[hi, "ELIGIBILITY_SCORE"], df.loc[hi, "INNOVATION_SCORE"],
            color=BLUE,   alpha=0.75, edgecolors="white", s=60,
            label="At/above cutoff (awarded)")

p = rdd_model.params
def rdd_line(x, treated):
    return (p["Intercept"] + p["TREATED"] * treated
            + p["SCORE_C"] * (x - CUTOFF)
            + p["TREATED:SCORE_C"] * treated * (x - CUTOFF))

ax3.plot(np.linspace(df.loc[lo,"ELIGIBILITY_SCORE"].min(), CUTOFF, 100),
         rdd_line(np.linspace(df.loc[lo,"ELIGIBILITY_SCORE"].min(),CUTOFF,100), 0),
         color=ORANGE, linewidth=2.5)
ax3.plot(np.linspace(CUTOFF, df.loc[hi,"ELIGIBILITY_SCORE"].max(), 100),
         rdd_line(np.linspace(CUTOFF,df.loc[hi,"ELIGIBILITY_SCORE"].max(),100), 1),
         color=BLUE,   linewidth=2.5)
ax3.axvline(x=CUTOFF, color=RED, linewidth=1.8, linestyle="--",
            label=f"Cutoff = {CUTOFF}")

jlo = rdd_line(np.array([CUTOFF]), 0)[0]
jhi = rdd_line(np.array([CUTOFF]), 1)[0]
ax3.annotate("", xy=(CUTOFF+0.3, jhi), xytext=(CUTOFF+0.3, jlo),
             arrowprops=dict(arrowstyle="<->", color=RED, lw=1.8))
ax3.text(CUTOFF+0.8, (jlo+jhi)/2,
         f"LATE\n{rdd_jump:+.2f}", color=RED, fontsize=9, va="center")
ax3.set_xlabel("Eligibility Score", fontsize=11)
ax3.set_ylabel("Innovation Score", fontsize=11)
ax3.set_title(f"RDD: Innovation Score around Cutoff = {CUTOFF}\n"
              f"LATE = {rdd_jump:.3f}  (p = {rdd_p:.3f})",
              fontsize=12, fontweight="bold")
ax3.legend(fontsize=9, framealpha=0.9)

ax4 = fig.add_subplot(gs[1, 1])
bins = np.arange(np.floor(df["ELIGIBILITY_SCORE"].min()),
                 np.ceil(df["ELIGIBILITY_SCORE"].max()) + 2, 2)
ax4.hist(df.loc[lo, "ELIGIBILITY_SCORE"], bins=bins,
         color=ORANGE, alpha=0.70, edgecolor="white", label="Below cutoff")
ax4.hist(df.loc[hi, "ELIGIBILITY_SCORE"], bins=bins,
         color=BLUE,   alpha=0.70, edgecolor="white", label="At/above cutoff")
kde   = gaussian_kde(df["ELIGIBILITY_SCORE"].values, bw_method=0.4)
x_kde = np.linspace(df["ELIGIBILITY_SCORE"].min()-2,
                    df["ELIGIBILITY_SCORE"].max()+2, 300)
ax4.plot(x_kde, kde(x_kde) * len(df) * 2, color="black",
         linewidth=2, label="KDE")
ax4.axvline(x=CUTOFF, color=RED, linewidth=1.8, linestyle="--",
            label=f"Cutoff = {CUTOFF}")
ax4.set_xlabel("Eligibility Score", fontsize=11)
ax4.set_ylabel("Count", fontsize=11)
ax4.set_title("Manipulation Test: Density of Eligibility Score\n"
              "(smooth at cutoff = no bunching = continuity holds)",
              fontsize=12, fontweight="bold")
ax4.legend(fontsize=9, framealpha=0.9)

fig.suptitle("Smart Campus Incubator -- IV & RDD Analysis",
             fontsize=15, fontweight="bold", y=1.01)
plt.savefig("analysis_results.png", dpi=150, bbox_inches="tight")
plt.close()
print("\nFigure saved: analysis_results.png")