from pathlib import Path import csv import json from collections import defaultdict import numpy as np 目录 = Path(__file__).resolve().parent def 读表(名称): with (目录/名称).open(encoding="utf-8-sig",newline="") as 文件: return list(csv.DictReader(文件)) def 写表(名称,列名,行): with (目录/名称).open("w",encoding="utf-8-sig",newline="") as 文件: 写入=csv.writer(文件);写入.writerow(列名);写入.writerows(行) # 步骤一：四行是组别与时期的汇总，不是四家门店 原表=读表("两州就业汇总.csv") assert len(原表)==4 assert {(r["州别"],r["调查时期"]) for r in 原表}=={ (州,期) for 州 in ["新泽西州","宾夕法尼亚州东部"] for 期 in ["政策前","政策后"]} 四格={} 汇总=[] for r in 原表: 人数=int(r["有效门店数"]);总量=float(r["全职当量合计"]) assert 人数>0 and np.isfinite(总量) t=int(r["州别"]=="新泽西州");p=int(r["调查时期"]=="政策后") 四格[t,p]=总量/人数 汇总.append([r["州别"],r["调查时期"],人数,总量,四格[t,p]]) # 步骤二：两次相减与反事实 处理变化=四格[1,1]-四格[1,0] 对照变化=四格[0,1]-四格[0,0] 反事实=四格[1,0]+对照变化 DID=处理变化-对照变化 assert abs(DID-2.75360578298058)<1e-10 assert abs(四格[1,1]-反事实-DID)<1e-12 # 步骤三：四格回归仅核对点估计 组合=sorted(四格) X=np.array([[1,t,p,t*p] for t,p in 组合],dtype=float) y=np.array([四格[k] for k in 组合]) 系数=np.linalg.solve(X,y) 自由度=len(y)-np.linalg.matrix_rank(X) assert abs(系数[3]-DID)<1e-12 and 自由度==0 # 步骤四：只改变假设的对照组趋势，原始调查文件保持不变 情景=读表("反事实情景.csv") assert len(情景)==1 假设变化=float(情景[0]["假设对照组变化"]) 情景DID=处理变化-假设变化 情景反事实=四格[1,0]+假设变化 assert np.isfinite(假设变化) and abs(四格[1,1]-情景反事实-情景DID)<1e-12 if 假设变化 == -1: assert abs(情景DID-1.58802136739616)<1e-10 # 步骤五：用另一份明确模拟的门店面板独立迁移 面板=读表("独立练习_模拟门店.csv") assert len(面板)==16 and len({(r["门店编号"],r["时期"]) for r in 面板})==16 门店=defaultdict(set);分格=defaultdict(list) for r in 面板: assert r["组别"] in {"试点组","参照组"} and r["时期"] in {"政策前","政策后"} 门店[r["门店编号"]].add((r["组别"],r["时期"])) 分格[r["组别"],r["时期"]].append(float(r["全职当量"])) assert len(门店)==8 and all(len(v)==2 and len({x[0] for x in v})==1 for v in 门店.values()) assert all(len(v)==4 for v in 分格.values()) 均值={k:sum(v)/len(v) for k,v in 分格.items()} 模拟处理变化=均值["试点组","政策后"]-均值["试点组","政策前"] 模拟对照变化=均值["参照组","政策后"]-均值["参照组","政策前"] 模拟DID=模拟处理变化-模拟对照变化 assert [均值["试点组","政策前"],均值["试点组","政策后"],均值["参照组","政策前"],均值["参照组","政策后"]]==[13,18,11,13] assert 模拟DID==3 and 模拟处理变化==5 and 模拟对照变化==2 # 步骤六：导出核对表和有单位、刻度的趋势图 写表("四格均值结果.csv",["州别","调查时期","有效门店数","全职当量合计","每店平均全职当量"],汇总) 写表("独立练习四格结果.csv",["组别","时期","记录数","每店平均全职当量"],[[*k,len(分格[k]),v] for k,v in sorted(均值.items())]) 核对={"follow":{"treated_change":处理变化,"control_change":对照变化,"did":DID}, "change":{"counterfactual":情景反事实,"did":情景DID}, "transfer":{"observations":len(面板),"treated_change":模拟处理变化,"did":模拟DID}} (目录/"作业数值核对.json").write_text(json.dumps(核对,ensure_ascii=False,indent=2),"utf-8") 写表("DID核对结果.csv",["案例","处理组变化","对照组变化","DID","性质"], [["公开调查四格",处理变化,对照变化,DID,"真实调查汇总"], ["改变反事实趋势",处理变化,假设变化,情景DID,"课堂假设"], ["另一组门店",模拟处理变化,模拟对照变化,模拟DID,"原创模拟"]]) def Y(v):return 330-(v-16)*24 图=['', '', '', '每店平均全职当量'] for v in [16,18,20,22,24,26]: yy=Y(v);图.append(f'{v}') for group,color,dash in [(1,"#222",""),(0,"#777",'stroke-dasharray="7 5"')]: a,b=四格[group,0],四格[group,1] 图.append(f'') for x,v in [(200,a),(550,b)]:图.append(f'') 图 += [f'', f'宾州东部 23.3312', f'新泽西州 20.4394', f'宾州东部 21.1656', f'新泽西州 21.0274', f'假设下的反事实 {反事实:.4f}', '政策前政策后', '点线为平行趋势假设下的反事实；只有一个前期，不能检验平行趋势。'] (目录/"两组趋势与反事实.svg").write_text("".join(图),"utf-8") 结果={"metrics":[{"label":"真实汇总DID","value":round(DID,8)},{"label":"处理组前后变化","value":round(处理变化,8)},{"label":"模拟迁移DID","value":模拟DID}], "columns":["案例","处理组变化","对照组变化","DID"], "rows":[["公开调查四格",处理变化,对照变化,DID],["改变反事实趋势",处理变化,假设变化,情景DID],["原创模拟门店",模拟处理变化,模拟对照变化,模拟DID]], "note":"真实调查按本课已有完整就业分项口径汇总，DID为2.75360578，保留与论文表3的2.76之差。四格回归残差自由度为0，不提供门店层面的标准误或平行趋势检验。情景和迁移数据均非新实证证据。", "checks":[{"name":"四格键与分母有效","ok":True},{"name":"两次差、反事实差和回归系数一致","ok":True},{"name":"四格回归零残差自由度","ok":bool(自由度==0)},{"name":"新模拟面板16行、8家门店且DID为3","ok":True}], "outputs":["四格均值结果.csv","独立练习四格结果.csv","作业数值核对.json","DID核对结果.csv","两组趋势与反事实.svg"]} (目录/"运行结果.json").write_text(json.dumps(结果,ensure_ascii=False,indent=2),"utf-8") print(json.dumps(核对,ensure_ascii=False,indent=2))