Bunching / Kink 聚束设计：从拐点处的"扎堆"反推行为弹性

*==============================================================*
* Bunching / Kink 演示：Saez (2010) 弹性估计
* 情境：月收入 z，在 z*=5000 处边际税率从 t0=0.03 跳到 t1=0.10
*==============================================================*
clear all
set more off
set seed 20260911

* --- 首次运行安装 ---
* ssc install bunching, replace       // 社区命令，自动做反事实多项式+bootstrap
* ssc install estout,   replace

* --- 1. 模拟数据：能力 n ~ LogNormal，常数弹性 e=0.3 ---
* 无 kink 时 z0 = n*(1-t0)^e；有 kink 时，z0>5000 的人被"挤"到 5000
set obs 50000
gen n = exp(rnormal(8.5, 0.6))        // 能力分布
local e_true = 0.3
local t0 = 0.03
local t1 = 0.10
local zstar = 5000
gen z0 = n*(1-`t0')^`e_true'          // 无 kink 时的潜在收入
gen z = z0
* 把 z0 落在 [5000, 5000+dz*] 的边际个体挤到拐点
* dz* 由弹性公式反推：dz*/z* = e*(t1-t0)/(1-t0)
local dz = `zstar' * `e_true' * (`t1'-`t0')/(1-`t0')
replace z = `zstar' if z0>=`zstar' & z0<=`zstar'+`dz'

* --- 2. 用社区命令 bunching 跑（推荐主回归）---
* z 是 running variable，a(5000) 指定拐点，nomore= 不截断右侧
bunching z, a(`zstar') ///
    polynomial(4)            ///  4 阶多项式拟合反事实
    required(500)            ///  excluded window 左侧 500
    cf(500)                  ///  excluded window 右侧 500
    binwidth(100)            ///  每 bin 100 元
    start(-2000)             ///  图的左界
    generate(bunch)          ///  输出估计值到 r()
    graph(title("Bunching：个税税率跳档处的聚束"))

* --- 3. 手写版本：histogram + 多项式 + excess mass + bootstrap ---
* 分 bin
local binw = 100
local zmin = 2000
local zmax = 10000
local nb = (`zmax'-`zmin')/`binw'
matrix cnt = J(`nb',1,.)
forvalues j=1/`nb' {
    local lo = `zmin' + (`j'-1)*`binw'
    local hi = `zmin' + `j'*`binw'
    quietly count if z>=`lo' & z<`hi'
    matrix cnt[`j',1] = r(N)
}
clear
set obs `nb'
gen cnt = .
gen zc  = .
forvalues j=1/`nb' {
    quietly replace cnt = cnt[`j',1] in `j'
    quietly replace zc  = `zmin' + (`j'-0.5)*`binw' in `j'
}
gen rel = zc - `zstar'          // 相对拐点的位置

* 多项式拟合反事实：排除拐点左右各 500 元
gen excluded = (abs(rel)<=500)
reg cnt c.rel##c.rel##c.rel##c.rel i.excluded
* 解释：四次多项式 + excluded 虚拟变量；
* excluded 虚拟变量的系数和就是该窗口内被"挤掉"的人数
predict cnt_hat
* excess mass = 实际 - 反事实（在 excluded 窗口内）
gen diff = cnt - cnt_hat if excluded==1
sum diff, detail
scalar B = r(sum)
* 拐点处反事实密度
sum cnt_hat if abs(rel)<=`binw'/2
scalar h0 = r(mean)/`binw'
scalar zstar = `zstar'
* Saez 弹性
scalar e_hat = (B/(h0*zstar)) / ((`t1'-`t0')/(1-`t0'))
display "Excess mass B = " B
display "Saez elasticity e_hat = " e_hat " (true = `e_true')"

* 画直方图 + 反事实线
twoway (bar cnt zc, barwidth(80) color(gs12)) ///
       (line cnt_hat zc, lwidth(2) lcolor(navy)), ///
       xline(`zstar', lpattern(dash) lcolor(maroon)) ///
       legend(label(1 实际人数) label(2 反事实密度)) ///
       xtitle("月应纳税所得额") ytitle("每 bin 人数") ///
       title("Bunching 手写估计：反事实多项式与 excess mass")

*==============================================================*
* Bootstrap 标准误：对个体有放回重抽样，重复全流程
*==============================================================*
* 把上面"模拟数据"段保留为程序，用 bootstrap 调用
cap prog drop bunching_est
prog define bunching_est, rclass
    version 16
    syntax , [zstar(real 5000) binw(real 100) w(real 500) t0(real .03) t1(real .10)]
    tempvar bin rel excl
    gen `bin' = floor((z - (`zstar'-3000))/`binw')
    * 直接调用 bunching 命令的返回值
    bunching z, a(`zstar') polynomial(4) required(`w') cf(`w') binwidth(`binw') ///
        start(-2000) quiet
    return scalar e = r(e)          // 若命令返回弹性，直接取
    return scalar b = r(b)          // 标准化 excess mass
end

* 运行 bootstrap（500 次，演示用；正式文章建议 1000-2000 次）
bootstrap e=r(e) b=r(b), reps(500) seed(20260911) ///
    saving("bunching_boot.dta", replace): bunching_est
estat bootstrap, all
* 输出：e 的偏差、标准误、正态/百分位/BCa 置信区间
* 报告时建议用 BCa 或 percentile 区间，而非正态近似

# ==============================================================
# Bunching 手写估计：模拟数据 + 多项式反事实 + excess mass
# ==============================================================
import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import polynomial as P

rng = np.random.default_rng(20260911)

# --- 1. 模拟数据：与 Stata 版本一致 ---
N = 50000
n = np.exp(rng.normal(8.5, 0.6, N))
e_true, t0, t1, zstar = 0.3, 0.03, 0.10, 5000
z0 = n * (1 - t0) ** e_true
dz = zstar * e_true * (t1 - t0) / (1 - t0)
z = z0.copy()
bunch_mask = (z0 >= zstar) & (z0 <= zstar + dz)
z[bunch_mask] = zstar                      # 把边际个体挤到拐点

# --- 2. 分 bin ---
binw = 100.0
edges = np.arange(2000, 10001, binw)
centers = (edges[:-1] + edges[1:]) / 2
cnt, _ = np.histogram(z, bins=edges)

# --- 3. 多项式拟合反事实（排除拐点 ±500） ---
rel = centers - zstar
excluded = np.abs(rel) <= 500
mask = ~excluded
coef = np.polyfit(rel[mask], cnt[mask], 4)
cnt_hat = np.polyval(coef, rel)

B = (cnt[excluded] - cnt_hat[excluded]).sum()
h0 = cnt_hat[np.abs(rel) <= binw / 2].mean() / binw
e_hat = (B / (h0 * zstar)) / ((t1 - t0) / (1 - t0))
print(f"Excess mass B = {B:.1f}")
print(f"Saez elasticity e_hat = {e_hat:.3f}  (true = {e_true})")

# --- 4. 可视化 ---
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.bar(centers, cnt, width=binw*0.9, color='#cccccc', label='实际人数')
ax.plot(centers, cnt_hat, lw=2, color='#1c5d6e', label='反事实密度（4 阶多项式）')
ax.axvline(zstar, ls='--', color='#a03a2c', label='kink z*=5000')
ax.axvspan(zstar-500, zstar+500, color='#b87333', alpha=0.15, label='excluded window')
ax.set_xlabel('月应纳税所得额'); ax.set_ylabel('每 bin 人数')
ax.set_title(f'Bunching: excess mass B={B:.0f},  e={e_hat:.3f}')
ax.legend(); plt.tight_layout(); plt.show()