spb/uqo-research-db

1"""Difference-in-differences event study: COVID-19 and digitalization.23Treatment = firms with above-median *physical-project intensity* in the4pre-period (2015-2019): plant construction, manufacturing/geographic expansion,5infrastructure, capex programs. Outcome = share of a firm-year's project6announcements that are digital (AI / cloud / ERP / data-center / automation /7digital-transformation / technology-deployment). We estimate dynamic8treatment x year effects (2019 omitted) with firm and year fixed effects,9testing whether disrupted incumbents accelerated digital catch-up after 2020.10"""11from __future__ import annotations1213import warnings14import numpy as np15import pandas as pd16import statsmodels.formula.api as smf1718from common import DATA, set_style, savefig, write_table, save_results19import matplotlib.pyplot as plt2021warnings.filterwarnings("ignore")22set_style()2324proj = pd.read_parquet(DATA / "projects.parquet")25proj = proj[proj["ann_year"].between(2013, 2025)].copy()26proj["ann_year"] = proj["ann_year"].astype(int)27PHYS = ("plant_construction", "manufacturing_expansion", "geographic_expansion",28        "infrastructure", "capex_program")2930g = proj.groupby(["cik", "ann_year"])31fy = g.agg(n=("project_id", "count"),32           n_digital=("is_digital", "sum"),33           n_phys=("project_type", lambda s: s.isin(PHYS).sum())).reset_index()34fy = fy.rename(columns={"ann_year": "year"})35fy["digital_share"] = fy["n_digital"] / fy["n"]3637# pre-period physical intensity -> treatment38pre = fy[fy.year.between(2015, 2019)].groupby("cik").apply(39    lambda d: d["n_phys"].sum() / max(d["n"].sum(), 1)).rename("phys_pre")40med = pre.median()41treat = (pre > med).astype(int).rename("treat")42fy = fy.merge(treat, on="cik", how="inner")43fy = fy[fy["year"].between(2015, 2025)].copy()4445# event-study dummies (omit 2019)46for y in sorted(fy["year"].unique()):47    if y == 2019:48        continue49    fy[f"d_{y}"] = ((fy["year"] == y) & (fy["treat"] == 1)).astype(int)5051evt_terms = " + ".join([f"d_{y}" for y in sorted(fy["year"].unique()) if y != 2019])52m = smf.ols(f"digital_share ~ {evt_terms} + C(cik) + C(year)", data=fy)\53       .fit(cov_type="cluster", cov_kwds={"groups": fy["cik"]})5455years = [y for y in sorted(fy["year"].unique())]56coefs, ci = [], []57for y in years:58    if y == 2019:59        coefs.append(0.0); ci.append(0.0)60    else:61        k = f"d_{y}"62        coefs.append(m.params[k]); ci.append(1.96 * m.bse[k])6364fig, ax = plt.subplots(figsize=(7.4, 4.2))65ax.errorbar(years, coefs, yerr=ci, fmt="o-", color="#1b3a5b", capsize=3, ms=4)66ax.axhline(0, color="grey", lw=0.8)67ax.axvline(2019.5, color="#c44e52", ls="--", lw=1.2, label="COVID-19 onset")68ax.set_xlabel("Year"); ax.set_ylabel("Treatment $\\times$ year effect on digital share")69ax.set_title("Event study: digitalization of physically-intensive incumbents")70ax.legend()71savefig(fig, "fig_did_event.pdf")7273# canonical 2x2 DiD74fy["post"] = (fy["year"] >= 2020).astype(int)75did = smf.ols("digital_share ~ treat*post + C(year)", data=fy)\76         .fit(cov_type="cluster", cov_kwds={"groups": fy["cik"]})77b = did.params["treat:post"]; se = did.bse["treat:post"]; pv = did.pvalues["treat:post"]78star = "***" if pv < .01 else "**" if pv < .05 else "*" if pv < .1 else ""7980lt = [r"\begin{tabular}{lc}", r"\toprule", r" & Digital share \\", r"\midrule",81      f"Treat $\\times$ Post & {b:.4f}{star} \\\\", f" & ({se:.4f}) \\\\",82      f"Treat (pre-mean) & {did.params['treat']:.4f} \\\\",83      r"\midrule", r"Year FE & Yes \\",84      f"Observations & {int(did.nobs):,} \\\\",85      f"$R^2$ & {did.rsquared:.3f} \\\\",86      r"\bottomrule", r"\end{tabular}"]87write_table("tab_did", "\n".join(lt))8889# pre-trend joint test (2015-2018 interactions == 0)90pre_terms = [f"d_{y}" for y in years if y < 2019]91ftest = m.f_test(" , ".join([f"{t} = 0" for t in pre_terms]))92save_results("did", {93    "did_coef": float(b), "did_se": float(se), "did_p": float(pv),94    "n_treat_firms": int(treat.sum()), "n_control_firms": int((treat == 0).sum()),95    "pretrend_F": float(np.ravel(ftest.fvalue)[0]), "pretrend_p": float(ftest.pvalue),96    "post_2022_effect": float(np.mean([m.params[f"d_{y}"] for y in (2022, 2023, 2024, 2025)])),97})98print(f"DiD done | treat:post={b:.4f}{star}, pretrend p={float(ftest.pvalue):.3f}")99