非参 / 半参方法：核回归、局部多项式与部分线性模型

*==============================================================*
* 非参 / 半参回归演示：核回归、局部线性、PLM
* 情境：模拟 Y = sin(X) + 0.3*X + eps，看非参能否抓出非线性
*==============================================================*
clear all
set more off
set seed 20260911

* --- 1. 模拟数据 ---
set obs 1000
gen X = runiform()*4*_pi - 2*_pi       // X ~ Uniform[-2pi, 2pi]
gen Z = rnormal()                       // PLM 中的非参变量
gen true_m = sin(X) + 0.3*Z            // 真实 m(X) 与线性 Z 系数
gen Y = true_m + rnormal()*0.5

* --- 2. 非参核回归：lpoly（局部线性） ---
* 语法：lpoly Y X, width(h) gengrid n(100) ci
* 默认 Epanechnikov 核；width 为带宽
lpoly Y X, leg(order(2)) ///
    ci                  ///  画 95% 置信带
    bw cvplot           ///  用 CV 选带宽并画 CV 曲线
    graph(title("局部线性回归：Y 关于 X")) ///
    ytitle("拟合值") xtitle("X")
graph export "lpoly_cv.png", replace

* --- 3. 多档带宽稳健性 ---
local h_default = 0.5
twoway (lpolyci Y X, fcolor(gs12%30) lcolor(navy)), ///
    legend(off) title("带宽 = 0.5")
* 对比 0.25 与 1.0
lpoly Y X, at(X) bwidth(0.25) gen(y025)
lpoly Y X, at(X) bwidth(0.50) gen(y050)
lpoly Y X, at(X) bwidth(1.00) gen(y100)
twoway (line y025 X, sort lcolor(red)) ///
       (line y050 X, sort lcolor(navy)) ///
       (line y100 X, sort lcolor(orange)), ///
       legend(label(1 "h=0.25") label(2 "h=0.50") label(3 "h=1.00")) ///
       title("带宽敏感性：过小/适中/过大")

* --- 4. 部分线性模型 PLM：Robinson (1988) 残差法 ---
* 模型：Y = beta*Z + g(X) + eps
* 步骤：① 非参拟合 E[Y|X]、E[Z|X]；② 残差化；③ OLS
lpoly Y X, at(X) bwidth(0.5) gen(E_Y)
lpoly Z X, at(X) bwidth(0.5) gen(E_Z)
gen tilde_Y = Y - E_Y
gen tilde_Z = Z - E_Z
reg tilde_Y tilde_Z
* tilde_Z 的系数就是 beta 的 Robinson 估计（渐近 sqrt-n 一致）
* 与 OLS 直接对比：
reg Y Z X
* 若 X 与 Z 相关、且 g(X) 非线性，OLS 会偏；PLM 一致

* --- 5. Stata 16+ 官方 npregress（更现代） ---
* npregress kernel Y X, kernel(epanechnikov) bw(cv)
* npregress series  Y X, order(3)  // 序列回归（多项式样条）
* 注意：npregress 是 Stata 16+ 内置，无需 ssc install

# ==============================================================
# 非参 / 半参：statsmodels KernelReg + 手写 Nadaraya-Watson
# ==============================================================
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.nonparametric.kernel_regression import KernelReg
from statsmodels.nonparametric.kde import KDEUnivariate

rng = np.random.default_rng(20260911)

# --- 1. 模拟数据：与 Stata 一致 ---
n = 1000
X = rng.uniform(-2*np.pi, 2*np.pi, n)
Z = rng.normal(size=n)
true_m = np.sin(X) + 0.3*Z
Y = true_m + rng.normal(scale=0.5, size=n)

# --- 2. statsmodels 核回归（局部线性，多维自动带宽） ---
# reg_type='ll' 局部线性；var_type='c' 连续变量
kr = KernelReg(endog=Y, exog=np.column_stack([X, Z]),
               var_type='cc', reg_type='ll', bw='cv_ls')
m_hat, bw = kr.fit()
print("CV 选出的带宽 =", bw)

# --- 3. 手写 Nadaraya-Watson（便于理解） ---
def nw(x0, X, Y, h, kernel='gauss'):
    u = (X - x0)/h
    if kernel == 'gauss':
        w = np.exp(-0.5*u**2)
    else:  # epanechnikov
        w = np.where(np.abs(u) <= 1, 0.75*(1-u**2), 0.0)
    return np.sum(w*Y)/np.sum(w)

grid = np.linspace(X.min(), X.max(), 200)
h = 0.5
nw_fit = np.array([nw(g, X, Y, h) for g in grid])
ll_fit, _ = kr.fit(data_predict=[grid, np.full_like(grid, Z.mean())])

# --- 4. 可视化 ---
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.scatter(X, Y, s=8, alpha=0.3, color='gray', label='观测')
ax.plot(grid, np.sin(grid)+0.3*Z.mean(), 'k--', lw=1.5, label='真实 m(X)')
ax.plot(grid, nw_fit, lw=2, color='#b87333', label=f'NW (h={h})')
ax.plot(grid, ll_fit, lw=2, color='#1c5d6e', label='局部线性 (CV)')
ax.set_xlabel('X'); ax.set_ylabel('Y')
ax.set_title('非参回归：NW vs 局部线性 vs 真实函数')
ax.legend(); plt.tight_layout(); plt.show()

# --- 5. Robinson (1988) PLM 估计 ---
# Y = beta*Z + g(X) + eps；先非参去均值，再 OLS
E_Y = np.array([nw(g, X, Y, h) for g in X])
E_Z = np.array([nw(g, X, Z, h) for g in X])
beta_plm = np.sum((Z-E_Z)*(Y-E_Y)) / np.sum((Z-E_Z)**2)
print(f"Robinson PLM beta_Z = {beta_plm:.3f}  (true = 0.300)")