第十节：经济类因果分析与双重差分

import pandas as pd
from pathlib import Path

data_path = Path("data/inclusive_finance_did_panel.csv")
df = pd.read_csv(data_path)
print(df.shape)
print(df.dtypes)
print(df.isna().sum())
print(df[["year", "city", "inclusive_finance_index"]].head())

print("城市数：", df["city"].nunique())
print("年份：", sorted(df["year"].unique()))
print("每个城市的年份数：")
print(df.groupby("city")["year"].nunique().value_counts().sort_index())
print("重复的城市—年份：", df.duplicated(["city_code", "year"]).sum())

numeric = ["inclusive_finance_index", "treatment", "post", "did"]
desc = (df[numeric].describe(percentiles=[.25, .5, .75]).T
    .rename(columns={"count":"N", "mean":"Mean", "std":"Std. Dev.", "25%":"P25", "50%":"Median", "75%":"P75"}))
print(desc[["N", "Mean", "Std. Dev.", "Min", "P25", "Median", "P75", "Max"]].round(3))

df["group"] = np.where(df["treatment"].eq(1), "Pilot", "Comparison")
df["period"] = np.where(df["post"].eq(1), "Post", "Pre")
group_stats = (df.groupby(["group", "period"])["inclusive_finance_index"]
    .agg(N="count", Mean="mean", Std="std", Min="min", Max="max").reset_index())
print(group_stats.round(3).to_string(index=False))

corr = df[["inclusive_finance_index", "treatment", "post", "did"]].corr()
print(corr.round(3))
import seaborn as sns
import matplotlib.pyplot as plt
sns.heatmap(corr, annot=True, vmin=-1, vmax=1, cmap="RdBu_r", center=0)
plt.title("变量相关矩阵（描述性）")
plt.tight_layout(); plt.show()

pilot_roots = ["北京市", "天津市", "承德市", "廊坊市", "温州市", "台州市",
               "宁波市", "泉州市", "宜昌市", "佛山市", "东莞市", "重庆市",
               "宜宾市", "玉溪市", "榆林市", "吴忠市"]  # 课堂先展示部分

# 完整名单在附录代码中；这里演示映射的逻辑
df["treatment"] = df["city"].apply(
    lambda city: int(any(root in city or city in root for root in pilot_roots))
)
print(df.groupby("treatment")["city"].nunique())

df["post"] = (df["year"] >= 2020).astype(int)
df["did"] = df["treatment"] * df["post"]
df["analysis_sample"] = df["year"].ne(2019)

print(df[["year", "treatment", "post", "did"]].drop_duplicates().sort_values(
    ["year", "treatment"]
).head(10))
print("主回归行数：", df.loc[df["analysis_sample"]].shape[0])

group_means = (
    df.loc[df["analysis_sample"]]
      .groupby(["year", "treatment"])["inclusive_finance_index"]
      .mean().reset_index()
)
print(group_means.head())

import matplotlib.pyplot as plt
for treated, label in [(0, "比较组"), (1, "试点组")]:
    part = group_means[group_means["treatment"] == treated]
    plt.plot(part["year"], part["inclusive_finance_index"], marker="o", label=label)
plt.axvline(2019.5, linestyle="--", color="darkorange", label="2020 起")
plt.legend(); plt.ylabel("核心指数总分"); plt.show()

import statsmodels.formula.api as smf

sample = df.loc[df["analysis_sample"].eq(1)].copy()
m0 = smf.ols("inclusive_finance_index ~ did", data=sample).fit(cov_type="HC1")
m1 = smf.ols(
    "inclusive_finance_index ~ did + C(city_code) + C(year)", data=sample
).fit(cov_type="cluster", cov_kwds={"groups": sample["city_code"]})

check = df.copy()
check["post_2019"] = (check["year"] >= 2019).astype(int)
check["did_2019"] = check["treatment"] * check["post_2019"]
m2 = smf.ols(
    "inclusive_finance_index ~ did_2019 + C(city_code) + C(year)", data=check
).fit(cov_type="cluster", cov_kwds={"groups": check["city_code"]})

for name, model, term in [("M0", m0, "did"), ("M1", m1, "did"), ("M2", m2, "did_2019")]:
    print(name, round(model.params[term], 4), round(model.bse[term], 4),
          round(model.pvalues[term], 4), int(model.nobs))

# 在拥有 Stata 17+ 授权的环境中运行
%pip install --upgrade stata_setup

import stata_setup
stata_setup.config("STATA_SYSDIR", "mp", splash=False)
from pystata import stata

stata.run('import delimited using "data/inclusive_finance_did_panel.csv", clear')
stata.run('drop if year == 2019')
stata.run('xtset city_code year')
stata.run('xtreg inclusive_finance_index did i.year, fe vce(cluster city_code)')

import statsmodels.formula.api as smf

sample = df.loc[df["analysis_sample"]].copy()
model = smf.ols(
    "inclusive_finance_index ~ did + C(city_code) + C(year)",
    data=sample
).fit(cov_type="cluster", cov_kwds={"groups": sample["city_code"]})

print(model.params["did"])
print(model.bse["did"])
print(model.pvalues["did"])

# 一个简单的敏感性检查：把 2019 也放入样本并定义为政策后
check = df.copy()
check["post_including_2019"] = (check["year"] >= 2019).astype(int)
check["did_including_2019"] = check["treatment"] * check["post_including_2019"]
check_model = smf.ols(
    "inclusive_finance_index ~ did_including_2019 + C(city_code) + C(year)",
    data=check
).fit(cov_type="cluster", cov_kwds={"groups": check["city_code"]})
print(check_model.params["did_including_2019"])